Quantics TCI of univariate function

Quantics TCI of univariate function#

using Plots
gr() # Use GR backend for plotting

import QuanticsGrids as QG
using QuanticsTCI: quanticscrossinterpolate, integral

# defines mutable struct `SemiLogy` and sets shorthands `semilogy` and `semilogy!`
@userplot SemiLogy
@recipe function f(t::SemiLogy)
    x = t.args[begin]
    y = t.args[end]
    ε = nextfloat(0.0)

    yscale := :log10
    # Warning: Invalid negative or zero value 0.0 found at series index 16 for log10 based yscale
    # prevent log10(0) from being -Inf
    (x, ε .+ y)
end

Example 1#

The first example is taken from Fig. 1 in Ritter2024.

We are going to compute the integral $I [f] = \int_{0}^{\ln 20} d x f (x)$ of the function

f (x) = \cos (\frac{x}{B}) \cos (\frac{x}{4 \sqrt{5} B}) e^{- x^{2}} + 2 e^{- x},

where $B = 2^{- 30}$ . The integral evaluates to $I [f] = 19 / 10 + O (e^{- 1 / (4 B^{2})})$ .

We first construct a QTT representation of the function $f (x)$ as follows:

B = 2^(-30) # global variable
function f(x)
    return cos(x / B) * cos(x / (4 * sqrt(5) * B)) * exp(-x^2) + 2 * exp(-x)
end

println(f(0.2))

1.7990930849034745

Let’s examine the behaviour of $f (x)$ . This function involves structure on widely different scales: rapid, incommensurate oscillations and a slowly decaying envelope. We’ll use PythonPlot.jl visualisation library which uses Python library matplotlib behind the scenes.

For small $x$ we have:

xs = LinRange(0, 2.0^(-23), 1000)

plt = plot(title="$(nameof(f))")
plot!(plt, xs, f.(xs), label="$(nameof(f))", legend=true)
plt

For $x \in (0, 3]$ we will get:

xs2 = LinRange(2.0^(-23), 3, 100000)
plt = plot(title="$(nameof(f))")
plot!(plt, xs2, f.(xs2), label="$(nameof(f))", legend=true)
plt

QTT representation#

We construct a QTT representation of this function on the domain $[0, \ln 20]$ , discretized on a quantics grid of size $2^{R}$ with $R = 40$ bits:

R = 40 # number of bits
xmin = 0.0
xmax = log(20.0)
N = 2^R # size of the grid
# * Uniform grid (includeendpoint=false, default):
#   -xmin, -xmin+dx, ...., -xmin + (2^R-1)*dx
#     where dx = (xmax - xmin)/2^R.
#   Note that the grid does not include the end point xmin.
#
# * Uniform grid (includeendpoint=true):
#   -xmin, -xmin+dx, ...., xmin-dx, xmin,
#     where dx = (xmax - xmin)/(2^R-1).
qgrid = QG.DiscretizedGrid{1}(R, xmin, xmax; includeendpoint=true)
ci, ranks, errors = quanticscrossinterpolate(Float64, f, qgrid; maxbonddim=15)

(QuanticsTCI.QuanticsTensorCI2{Float64}(TensorCrossInterpolation.TensorCI2{Float64} with rank 15, QuanticsGrids.DiscretizedGrid{1}(40, (0.0,), (2.995732273553991,), 2, :fused, true), TensorCrossInterpolation.CachedFunction{Float64, UInt128} with 45662 entries), [17, 15, 15], [3.003930472475819e-8, 2.4814437584894508e-8, 3.086415959710106e-8])

length(ci.tci)

Here, we’ve created the object ci of type QuanticsTensorCI2{Float64}. This can be evaluated at an linear index $i$ ( $1 \leq i \leq 2^{R}$ ) as follows:

for i in [1, 2, 3, 2^R] # Linear indices
    # restore original coordinate `x` from linear index `i`
    x = QG.grididx_to_origcoord(qgrid, i)
    println("x: $(x), i: $(i), tci: $(ci(i)), ref: $(f(x))")
end

x: 0.0, i: 1, tci: 3.0, ref: 3.0
x: 2.724602630730001e-12, i: 2, tci: 2.999995667173124, ref: 2.999995667173124
x: 5.449205261460002e-12, i: 3, tci: 2.9999826687427715, ref: 2.9999826687427715
x: 2.995732273553991, i: 1099511627776, tci: 0.10000652643362676, ref: 0.1000065283433056

We see that ci(i) approximates the original f at x = QG.grididx_to_origcoord(qgrid, i). Let’s plot them together.

maxindex = QG.origcoord_to_grididx(qgrid, 2.0^(-23))
testindices = Int.(round.(LinRange(1, maxindex, 1000)))

xs = [QG.grididx_to_origcoord(qgrid, i) for i in testindices]
ys = f.(xs)
yci = ci.(testindices)

plt = plot(title="$(nameof(f)) and TCI", xlabel="x", ylabel="y")
plot!(plt, xs, ys, label="$(nameof(f))", legend=true)
plot!(plt, xs, yci, label="tci", linestyle=:dash, alpha=0.7, legend=true)
plt

Above, one can see that the original function is interpolated very accurately.

Let’s plot of $x$ vs interpolation error $| f (x) - ci (x) |$ for small $x$

ys = f.(xs)
yci = ci.(testindices)
plt = plot(title="x vs interpolation error: $(nameof(f))", xlabel="x", ylabel="interpolation error")
semilogy!(xs, abs.(ys .- yci), label="log(|f(x) - ci(x)|)", yscale=:log10, legend=:bottomright, ylim=(1e-16, 1e-7), yticks=10.0 .^ collect(-16:1:-7))
plt

… and for all $x$ :

plt = plot(title="x vs interpolation error: $(nameof(f))", xlabel="x", ylabel="interpolation error")

testindices = Int.(round.(LinRange(1, 2^R, 1000)))
xs = [QG.grididx_to_origcoord(qgrid, i) for i in testindices]
ys = f.(xs)
yci = ci.(testindices)
semilogy!(xs, abs.(ys .- yci), label="log(|f(x) - ci(x)|)", legend=true, ylim=(1e-16, 1e-6), yticks=10.0 .^ collect(-16:1:-6))
plt

The function is approximated with an accuracy $\approx 10^{- 7}$ over the entire domain.

We are now ready to compute the integral $I [f] = \int_{0}^{\ln 20} d x f (x) ≃ 19 / 10$ using the QTT representation of $f (x)$ .

integral(ci), 19 / 10

(1.8999999945916517, 1.9)

integral(ci) is equivalent to calling QuanticsTCI.sum(ci) and multiplying the result by the interval length divided by $2^{R}$ .

sum(ci) * (log(20) - 0) / 2^R, 19 / 10

(1.8999999945899237, 1.9)

About `ci::QuanticsTensorCI2{Float64}`#

Let’s dive into the ci object:

println(typeof(ci))

QuanticsTCI.QuanticsTensorCI2{Float64}

As we’ve seen before, ci is an object of QuanticsTensorCI2{Float64} in QuanticsTCI.jl, which is a thin wrapper of TensorCI2{Float64} in TensorCrossInterpolation.jl. The undering object of TensorCI2{Float64} type can be accessed as ci.tci. This will be useful for obtaining more detailed information on the TCI results.

For instance, ci.tci.maxsamplevalue is an estimate of the abosolute maximum value of the function, and ci.tci.pivoterrors stores the error as function of the bond dimension computed by prrLU. In the following figure, we plot the normalized error vs. bond dimension, showing an exponential decay.

# Plot error vs bond dimension obtained by prrLU
plt = plot(title="normalized error vs. bond dimension: $(nameof(f))", xlabel="Bond dimension", ylabel="Normalization error")
semilogy!(1:length(ci.tci.pivoterrors), ci.tci.pivoterrors ./ ci.tci.maxsamplevalue, marker=:x, ylim=(1e-8, 10), yticks=(10.0 .^ (-10:1:0)), legend=false)
plt

Function evaluations#

Our TCI algorithm does not call elements of the entire tensor, but constructs the TT (Tensor Train) from some elements chosen adaptively. On which points $x \in [0, 3]$ was the function evaluated to construct a QTT representation of the function $f (x)$ ? Let’s find out. One can retrieve the information on the function evaluations as follows.

import QuanticsTCI
# Dict{Float64,Float64}
# key: `x`
# value: function value at `x`
evaluated = QuanticsTCI.cachedata(ci)

Dict{Float64, Float64} with 45662 entries:
  0.000878537 => 2.31392
  0.714978    => 0.379167
  1.46693e-8  => 2.1889
  0.738286    => 0.925032
  0.0726419   => 2.47049
  0.74619     => 0.877457
  0.714913    => 0.759495
  0.714654    => 1.07975
  0.714678    => 1.04103
  1.50007     => 0.447397
  0.714913    => 1.01016
  0.67287     => 1.39331
  0.0471769   => 1.13029
  0.715207    => 1.1412
  2.80889     => 0.120594
  0.000808344 => 1.155
  0.714681    => 0.454686
  0.714913    => 1.35309
  0.000879206 => 2.04983
  1.53746e-6  => 2.05443
  0.00123981  => 2.87518
  0.583997    => 0.754877
  0.000868523 => 2.18544
  0.738656    => 0.499682
  0.00100509  => 1.08149
  ⋮           => ⋮

Let’s plot f and the evaluated points together.

f̂(x) = ci(QG.origcoord_to_quantics(qgrid, x))
xs = LinRange(0, 2.0^(-23), 1000)

xs_evaluated = collect(keys(evaluated))
fs_evaluated = [evaluated[x] for x in xs_evaluated]

plt = plot(title="$(nameof(f)) and TCI", xlabel="x", ylabel="y", xlim=(0, maximum(xs)))
plot!(plt, xs, f.(xs), label="$(nameof(f))")
scatter!(plt, xs_evaluated, fs_evaluated, marker=:x, label="evaluated points")
plt

Example 2#

We now consider the function:

\begin{array}{r} \begin{aligned} f (x) & = sinc (x) + 3 e^{- 0.3 (x - 4)^{2}} sinc (x - 4) \\ - \cos (4 x)^{2} - 2 sinc (x + 10) e^{- 0.6 (x + 9)} + 4 \cos (2 x) e^{- | x + 5 |} \\ + \frac{6}{x - 11} + \sqrt{(| x |)} \arctan (x / 15) . \end{aligned} \end{array}

One can construct a QTT representation of this function on the domain $[- 10, 10)$ using a quantics grid of size $2^{R}$ ( $R = 20$ ):

import QuanticsGrids as QG
using QuanticsTCI

R = 20 # number of bits
N = 2^R  # size of the grid

qgrid = QG.DiscretizedGrid{1}(R, -10, 10; includeendpoint=false)

# Function of interest
function oscillation_fn(x)
    return (
        sinc(x) + 3 * exp(-0.3 * (x - 4)^2) * sinc(x - 4) - cos(4 * x)^2 -
        2 * sinc(x + 10) * exp(-0.6 * (x + 9)) + 4 * cos(2 * x) * exp(-abs(x + 5)) +
        6 * 1 / (x - 11) + sqrt(abs(x)) * atan(x / 15))
end

# Convert to quantics format and sweep
ci, ranks, errors = quanticscrossinterpolate(Float64, oscillation_fn, qgrid; maxbonddim=15)

(QuanticsTCI.QuanticsTensorCI2{Float64}(TensorCrossInterpolation.TensorCI2{Float64} with rank 14, QuanticsGrids.DiscretizedGrid{1}(20, (-10.0,), (10.0,), 2, :fused, false), TensorCrossInterpolation.CachedFunction{Float64, UInt128} with 3414 entries), [12, 14, 14, 14], [4.56720264014334e-9, 7.150033486208613e-9, 7.150033486208613e-9, 7.150033486208613e-9])

for i in [1, 2, 2^R] # Linear indices
    x = QG.grididx_to_origcoord(qgrid, i)
    println("x: $(x), tci: $(ci(i)), ref: $(oscillation_fn(x))")
end

x: -10.0, tci: -6.223187220679993, ref: -6.223187220679993
x: -9.999980926513672, tci: -6.223066148054395, ref: -6.223066148054395
x: 9.999980926513672, tci: -4.585194580210248, ref: -4.585194580210247

Above, one can see that the original function is interpolated very accurately. The function grididx_to_origcoord transforms a linear index to a coordinate point $x$ in the original domain ( $- 10 \leq x < 10$ ).

In the following figure, we plot the normalized error vs. bond dimension, showing an exponential decay.

# Plot error vs bond dimension obtained by prrLU
plt = plot(xlabel="Bond dimension", ylabel="Normalization error", title="normalized error vs. bond dimension")
semilogy!(1:length(ci.tci.pivoterrors), ci.tci.pivoterrors ./ ci.tci.maxsamplevalue, marker=:x, ylim=(1e-8, 10), yticks=(10.0 .^ (-10:1:0)), legend=false)
plt

Example 3#

Control the error of the TCI by a tolerance#

We interpolate the same function as in Example 2, but this time we use a tolerance to control the error of the TCI. The tolerance is a positive number that determines the maximum error of the TCI, which is scaled by an estimate of the abosolute maximum of the function. The TCI algorithm will adaptively increase the bond dimension until the error is below the tolerance.

tol = 1e-8 # Tolerance for the error

# Convert to quantics format and sweep
ci_tol, ranks_tol, errors_tol = quanticscrossinterpolate(
    Float64, oscillation_fn, qgrid;
    tolerance=tol,
    normalizeerror=true, # Normalize the error by the maximum sample value,
    verbosity=1, loginterval=1, # Log the error every `loginterval` iterations
)

iteration = 1, rank = 12, error= 4.380931188449161e-8, maxsamplevalue= 6.223187220679993, nglobalpivot=3
Rejected 3 global pivots added in the previous iteration, errors are [3.270628212703741e-10, 1.4377037338419996e-9, 4.3811990146913615e-9]
iteration = 2, rank = 14, error= 4.4476743568461786e-8, maxsamplevalue= 6.223187220679993, nglobalpivot=0
iteration = 3, rank = 14, error= 4.4476743568461786e-8, maxsamplevalue= 6.223187220679993, nglobalpivot=0
iteration = 4, rank = 14, error= 4.4476743568461786e-8, maxsamplevalue= 6.223187220679993, nglobalpivot=0

(QuanticsTCI.QuanticsTensorCI2{Float64}(TensorCrossInterpolation.TensorCI2{Float64} with rank 14, QuanticsGrids.DiscretizedGrid{1}(20, (-10.0,), (10.0,), 2, :fused, false), TensorCrossInterpolation.CachedFunction{Float64, UInt128} with 3581 entries), [12, 14, 14, 14], [7.039690488325799e-9, 7.146939661507069e-9, 7.146939661507069e-9, 7.146939661507069e-9])

println("Max abs sampled value is $(ci_tol.tci.maxsamplevalue)")

Max abs sampled value is 6.223187220679993

errors_tol ./ ci_tol.tci.maxsamplevalue

4-element Vector{Float64}:
 1.131203391235366e-9
 1.14843719272295e-9
 1.14843719272295e-9
 1.14843719272295e-9

Estimate the error of the TCI#

Wait! Since we did not sample the function over the entire domain, we do not know the true error of the TCI. In theory, we can estimate the error of the TCI by comparing the function values at the sampled points with the TCI values at the same points. But, it is not practical to compare the function values with the TCI values at all points in the domain. The function estimatetrueerror in TensorCrossInterpolation.jl provides a good estimate of the error of the TCI. The algorithm finds indices (points) where the error is large by a randomized global search algorithm starting with a set of random initial points.

import TensorCrossInterpolation as TCI
pivoterror_global = TCI.estimatetrueerror(TCI.TensorTrain(ci.tci), ci.quanticsfunction; nsearch=100) # Results are sorted in descending order of the error

81-element Vector{Tuple{Vector{Int64}, Float64}}:
 ([1, 2, 1, 1, 2, 1, 1, 2, 2, 1, 1, 2, 1, 2, 1, 1, 2, 1, 2, 1], 7.228325804575775e-8)
 ([1, 2, 1, 1, 2, 1, 1, 2, 2, 1, 1, 2, 1, 2, 1, 1, 2, 1, 1, 2], 7.218350006610308e-8)
 ([1, 2, 1, 1, 2, 1, 1, 2, 2, 1, 1, 2, 1, 2, 1, 2, 1, 1, 1, 2], 7.108977806424832e-8)
 ([1, 1, 1, 1, 2, 2, 1, 1, 2, 1, 1, 2, 1, 2, 1, 1, 2, 1, 2, 1], 7.040376348577126e-8)
 ([1, 1, 1, 1, 2, 2, 1, 1, 2, 1, 1, 2, 1, 2, 1, 1, 2, 1, 1, 2], 7.033973470349508e-8)
 ([1, 2, 1, 1, 2, 1, 2, 1, 1, 1, 1, 2, 1, 2, 1, 1, 2, 1, 1, 2], 7.01437330441479e-8)
 ([1, 1, 1, 1, 2, 2, 1, 2, 1, 1, 1, 2, 1, 2, 1, 1, 2, 1, 2, 1], 6.850878553343875e-8)
 ([1, 1, 1, 1, 2, 2, 1, 2, 1, 1, 1, 2, 1, 2, 1, 1, 2, 1, 1, 2], 6.844331190691832e-8)
 ([1, 1, 1, 1, 2, 2, 1, 2, 1, 1, 1, 2, 1, 2, 1, 2, 1, 1, 1, 2], 6.733362822863853e-8)
 ([1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 2, 1, 2, 1, 1, 2, 1, 2, 1], 6.38076578240998e-8)
 ([1, 2, 1, 1, 2, 1, 2, 1, 1, 1, 1, 2, 1, 1, 2, 2, 2, 2, 2, 1], 6.373489380706587e-8)
 ([1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 2, 1, 2, 1, 1, 2, 2, 1, 2], 6.370827598800588e-8)
 ([1, 1, 1, 1, 2, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 2, 2, 2, 1, 2], 6.356614790092863e-8)
 ⋮
 ([1, 2, 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 2], 2.323339959309223e-8)
 ([2, 1, 1, 2, 1, 1, 2, 2, 2, 1, 2, 1, 1, 2, 1, 2, 1, 1, 1, 2], 2.2546057243388873e-8)
 ([2, 1, 2, 2, 1, 2, 1, 2, 1, 1, 1, 2, 2, 1, 2, 2, 1, 1, 2, 1], 1.9336143797232808e-8)
 ([2, 1, 2, 2, 1, 1, 1, 1, 1, 2, 2, 1, 2, 1, 2, 2, 1, 1, 2, 1], 1.8184483030481147e-8)
 ([2, 2, 1, 1, 2, 2, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 2, 2, 1], 1.6522432977339463e-8)
 ([2, 2, 1, 2, 1, 2, 1, 2, 2, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 1], 1.5915256862397698e-8)
 ([1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 2, 2, 2, 2], 1.5276709008915645e-8)
 ([2, 2, 2, 1, 1, 2, 2, 1, 1, 2, 2, 1, 2, 1, 2, 2, 1, 1, 1, 1], 1.445382435960596e-8)
 ([2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2], 1.3113687913346439e-8)
 ([2, 2, 2, 2, 2, 1, 1, 1, 2, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1], 1.3005156063172763e-8)
 ([2, 1, 1, 2, 1, 1, 2, 1, 2, 1, 2, 1, 2, 1, 1, 2, 2, 2, 1, 2], 1.2944330496367229e-8)
 ([1, 1, 2, 2, 1, 2, 2, 1, 1, 2, 2, 1, 1, 2, 2, 2, 2, 2, 1, 1], 5.7397369079836835e-9)

Now, you can see the error estimate of the TCI is below the tolerance of $10^{- 8}$ (or close to it).

println("The largest error found is $(pivoterror_global[1][2]) and the corresponding pivot is $(pivoterror_global[1][1]).")
println("The tolerance used is $(tol * ci_tol.tci.maxsamplevalue).")

The largest error found is 7.228325804575775e-8 and the corresponding pivot is [1, 2, 1, 1, 2, 1, 1, 2, 2, 1, 1, 2, 1, 2, 1, 1, 2, 1, 2, 1].
The tolerance used is 6.223187220679993e-8.