Quantics TCI of univariate funciton (advanced topics)

Click here to download the notebook locally.

Quantics TCI of univariate funciton (advanced topics)#

import QuanticsGrids as QG
import TensorCrossInterpolation as TCI

using Plots

Example 1 (continuation)#

Performance tips#

Let’s recall again the function \(f(x)\) from the previous page.

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

where \(B = 2^{-30}\).

In the previous page, we implemented \(f(x)\) as

B = 2^(-30) # global variable
# This is simple, but not recommended
function f(x)
    return cos(x/B) * cos(x/(4*sqrt(5)*B)) * exp(-x^2) + 2 * exp(-x)
end

However, since the implementation treats B as a untyped global variables, Julia compiler won’t generate efficient code. See Performance tips at the official Julia documentatoin to learn more. In Julia, this can be written as below:

# Define callable struct
Base.@kwdef struct Ritter2024 <: Function
    B::Float64 = 2^(-30)
end

# Make Ritter2024 be "callable" object.
function (obj::Ritter2024)(x)
    B = obj.B
    return cos(x / B) * cos(x / (4 * sqrt(5) * B)) * exp(-x^2) + 2 * exp(-x)
end

f = Ritter2024()
nothing # hide

Such “callable” objects are sometimes called “functors.” See Function like objects at the official Julia documentation to learn more.

The inside work of quanticscrossinterpolate#

Let us see the inside work of the quanticscrossinterpolate function in QuanticsTCI.jl by implementing the previous example using QuanticsGrids.jl and TensorCrossInterpolation.jl.

We first define a grid as in the previous example:

import QuanticsGrids as QG

R = 40 # number of bits
xmin = 0.0
xmax = log(20)
qgrid = QG.DiscretizedGrid{1}(R, xmin, xmax; includeendpoint=false)

QuanticsGrids.DiscretizedGrid{1}(40, (0.0,), (2.995732273553991,), 2, :fused, false)

We now define a function that takes a quantics index, which is named qf.

localdims = fill(2, R) # Sizes of local indices
qf(q) = f(QG.quantics_to_origcoord(qgrid, q))
cf = TCI.CachedFunction{Float64}(qf, localdims)

(::TensorCrossInterpolation.CachedFunction{Float64, UInt128}) (generic function with 3 methods)

Here, we’ve defined a function qf that accepts quantics q as an argument. Then we wrap qf using TCI.CachedFunction. The function QuanticsGrids.quantics_to_origcoord converts a quantics index to a point (in this example, of type Float64) in the original coordinate system. TCI.CachedFunction caches function evaluations; cf caches the pair of input and output that are used during constructing a QTT representation. This is useful if a function evaluation takes more than 100 ns or if you want to record function evaluations.

Choosing good initial pivots is critical for numerical stability. In the following code, we generate initial pivots by finding local maxima of \(|f(x)|\). The function optfirstpivot maximizes the given function qf using single-index updates (zero-temperature Monte Carlo).

nrandominitpivot = 5

# random initial pivot
initialpivots = [
        TCI.optfirstpivot(qf, localdims, [rand(1:d) for d in localdims]) for _ in 1:nrandominitpivot
]

qf.(initialpivots) # Function values at initial pivots

5-element Vector{Float64}:
 2.859709450326715
 2.964309079411982
 2.8752388302345873
 2.7620042373840032
 2.8668685832700014

We are ready to construct a QTT representation of cf via TCI.crossinterpolate2:

ci, ranks, errors = TCI.crossinterpolate2(Float64, cf, localdims, initialpivots; maxbonddim=15)

(TensorCrossInterpolation.TensorCI2{Float64} with rank 15, [20, 15, 15], [2.4824872103285508e-8, 4.388253825686919e-8, 3.405657211732417e-8])

The integral of \(f(x)\) can be computed by summing all the elements in the QTT representation and multiplying by the interval length divided by \(2^\mathcal{R}\).

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

(1.8999999960268668, 1.9)

You can retrieve the results of the function evaluations during the TCI construction as follows.

cache = TCI.cachedata(cf) # Dict{Vector{Int},Float64}

_qs = collect(keys(cache)) # quantics indices
_xs = QG.quantics_to_origcoord.(Ref(qgrid), _qs) # original coordinates
xs_evaluated = sort(_xs);

The object xs_evaluated stores points what we wanted to know.

We expect ci::TensorCrossInterpolation.TensorCI2{Float64} approximates cf(and f): Just in case, we will check qf approximates f and qf and cf functions output the same result:

x = 0.2
# convert `x` in the original coordinate system to the corresponding `q` of quantics
q = QG.origcoord_to_quantics(qgrid, x)
println("f(x) = $(f(x)), qf(q) = $(qf(q)), cf(q) = $(cf(q)), ci(q) = $(ci(q))")

f(x) = 1.7990930849034745, qf(q) = 1.798322002608539, cf(q) = 1.798322002608539, ci(q) = 1.7983219680357108

@show length(TCI.cachedata(cf))

length(TCI.cachedata(cf)) = 59220

59220
Contents