To be able to edit code and run cells, you need to run the notebook yourself. Where would you like to run the notebook?

This notebook takes about 50 seconds to run.

In the cloud (experimental)

Binder is a free, open source service that runs scientific notebooks in the cloud! It will take a while, usually 2-7 minutes to get a session.

On your computer

(Recommended if you want to store your changes.)

  1. Copy the notebook URL:
  2. Run Pluto

    (Also see: How to install Julia and Pluto)

  3. Paste URL in the Open box

Frontmatter

If you are publishing this notebook on the web, you can set the parameters below to provide HTML metadata. This is useful for search engines and social media.

Author 1

To run this notebook locally, copy and paste the following into your Julia REPL:

using Pkg; Pkg.activate(temp=true); Pkg.add("Pluto")
BASE_URL = "https://raw.githubusercontent.com/tensor4all/T4APlutoExamples/refs/heads/main/pluto_notebooks/"
notebook = "quantics1d.jl"
url = joinpath(BASE_URL, notebook)
using Pluto; Pluto.run(notebook=download(url))
👀 Reading hidden code
193 μs

Quantics TCI of univariate function

👀 Reading hidden code
179 μs
begin
using Plots
gr() # Use GR backend for plotting

import QuanticsGrids as QG
import TensorCrossInterpolation as TCI

using QuanticsTCI: QuanticsTCI, quanticscrossinterpolate, integral
end
👀 Reading hidden code
1.5 s
begin
# 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
end
👀 Reading hidden code
186 ms

Example 1

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

We are going to compute the integral \ (generic function with 175 methods)mathrm{I}[f] = \int_0^{\ln 20} \mathrm{d}x f(x) $ of the function

$$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}$. The integral evaluates to $\mathrm{I}[f] = 19/10 + O(e^{-1/(4B^2)})$.

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

👀 Reading hidden code
70.5 ms
begin
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))
end
👀 Reading hidden code
❔
70.7 ms

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:

👀 Reading hidden code
337 μs
let
xs = LinRange(0, 2.0^(-23), 1000)

plt = plot(title = "$(nameof(f))")
plot!(plt, xs, f.(xs), label = "$(nameof(f))", legend = true)
plt
end
👀 Reading hidden code
409 ms

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

👀 Reading hidden code
188 μs
let
xs2 = LinRange(2.0^(-23), 3, 100000)
plt = plot(title = "$(nameof(f))")
plot!(plt, xs2, f.(xs2), label = "$(nameof(f))", legend = true)
plt
end
👀 Reading hidden code
36.1 ms

QTT representation

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

👀 Reading hidden code
276 μs
begin
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)
end
👀 Reading hidden code
11.7 s
40
length(ci.tci)
👀 Reading hidden code
16.3 μs

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

👀 Reading hidden code
211 μs
let
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
end
👀 Reading hidden code
❔
92.1 ms

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

👀 Reading hidden code
221 μs
let
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
end
👀 Reading hidden code
263 ms

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

Let's plot of $x$ vs interpolation error $|f(x) - \mathrm{ci}(x)|$ for small $x$

👀 Reading hidden code
251 μs
let
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 = "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
end
👀 Reading hidden code
290 ms

... and for all $x$:

👀 Reading hidden code
187 μs
let
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
end
👀 Reading hidden code
289 ms
Loading more cells...