Sparsity

.gitignore

@@ -1 +1 @@
-/Manifest.toml
+/Manifest.toml

examples/highly_oscillatory.jl

@@ -1,41 +1,43 @@
 using PolynomialQTT
 using CairoMakie
 import QuanticsGrids as QG
-import TensorCrossInterpolation as TCI 
+import TensorCrossInterpolation as TCI
+using LinearAlgebra
 
 f(x) = cos(x^2) + sin(π * x)
-a,b = -2.0, sqrt(2)
+a, b = -2.0, sqrt(2)
 K = 10
 R = 8
 
 tt = PolynomialQTT.interpolatesinglescale(f, a, b, R, K)
+tt_adaptive = interpolateadaptive(f, a, b, R, K)
 
 grid = QG.DiscretizedGrid(R, a, b)
 
-plotquantics = QG.grididx_to_quantics.(Ref(grid), 1:2^R)
-plotx = QG.grididx_to_origcoord.(Ref(grid), 1:2^R)
+plotquantics = QG.grididx_to_quantics.(Ref(grid), 1:(2^R))
+plotx = QG.grididx_to_origcoord.(Ref(grid), 1:(2^R))
 origdata = f.(plotx)
 ttdata = tt.(plotquantics)
 
-let 
+let
     fig = Figure()
-    ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"f(x)", title="Single-scale interpolation")
-    lines!(ax, plotx, f.(plotx), label="Original function", linewidth=2.0, linestyle=:solid)
-    lines!(ax, plotx, ttdata, label=L"K=%$K", linewidth=2.0,linestyle=:dash)
-    axislegend(ax, pos=:best)
+    ax = Axis(fig[1, 1], xlabel = L"x", ylabel = L"f(x)", title = "Single-scale interpolation")
+    lines!(ax, plotx, f.(plotx), label = "Original function", linewidth = 2.0, linestyle = :solid)
+    lines!(ax, plotx, ttdata, label = L"K=%$K", linewidth = 2.0, linestyle = :dash)
+    axislegend(ax, pos = :best)
     fig
 end
 
-let 
+let
     fig = Figure()
-    ax = Axis(fig[1, 1], xlabel=L"x", ylabel="Error", title="Error single-scale interpolation")
+    ax = Axis(fig[1, 1], xlabel = L"x", ylabel = "Error", title = "Error single-scale interpolation")
     lines!(ax, plotx, abs.(ttdata .- origdata))
     fig
 end
 
-let 
+let
     fig = Figure()
-    ax = Axis(fig[1, 1], xlabel=L"\ell", ylabel=L"\chi_\ell", title="Bond dimension")
-    scatterlines!(ax,  1:R-1, TCI.linkdims(tt))
+    ax = Axis(fig[1, 1], xlabel = L"\ell", ylabel = L"\chi_\ell", title = "Bond dimension")
+    scatterlines!(ax, 1:(R - 1), TCI.linkdims(tt))
     fig
-end
+end

examples/multiscale.jl

@@ -1,41 +1,42 @@
 using PolynomialQTT
 using CairoMakie
 import QuanticsGrids as QG
-import TensorCrossInterpolation as TCI 
+import TensorCrossInterpolation as TCI
 
 α = 0.01
-f(x) = exp(-0.5*(x/α)^2)
-a,b = 0.0, 1.0
+f(x) = exp(-0.5 * (x / α)^2)
+a, b = 0.0, 1.0
 K = 25
 R = 8
 
 tt = PolynomialQTT.interpolatesinglescale(f, a, b, R, K)
+tt_adaptive = interpolateadaptive(f, a, b, R, K)
 
 grid = QG.DiscretizedGrid(R, a, b)
 
-plotquantics = QG.grididx_to_quantics.(Ref(grid), 1:2^R)
-plotx = QG.grididx_to_origcoord.(Ref(grid), 1:2^R)
+plotquantics = QG.grididx_to_quantics.(Ref(grid), 1:(2^R))
+plotx = QG.grididx_to_origcoord.(Ref(grid), 1:(2^R))
 origdata = f.(plotx)
 ttdata = tt.(plotquantics)
 
-let 
+let
     fig = Figure()
-    ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"f(x)", title="Single-scale interpolation")
-    lines!(ax, plotx, f.(plotx), label=L"Origignal function", linewidth=2.0, linestyle=:solid)
-    lines!(ax, plotx, ttdata, label=L"K=%$K", linewidth=2.0,linestyle=:dash)
+    ax = Axis(fig[1, 1], xlabel = L"x", ylabel = L"f(x)", title = "Single-scale interpolation")
+    lines!(ax, plotx, f.(plotx), label = L"Origignal function", linewidth = 2.0, linestyle = :solid)
+    lines!(ax, plotx, ttdata, label = L"K=%$K", linewidth = 2.0, linestyle = :dash)
     fig
 end
 
-let 
+let
     fig = Figure()
-    ax = Axis(fig[1, 1], xlabel=L"x", ylabel="Error", title="Error single-scale interpolation")
+    ax = Axis(fig[1, 1], xlabel = L"x", ylabel = "Error", title = "Error single-scale interpolation")
     lines!(ax, plotx, abs.(ttdata .- origdata))
     fig
 end
 
-let 
+let
     fig = Figure()
-    ax = Axis(fig[1, 1], xlabel=L"\ell", ylabel=L"\chi_\ell", title="Bond dimension")
-    scatterlines!(ax,  1:R-1, TCI.linkdims(tt))
+    ax = Axis(fig[1, 1], xlabel = L"\ell", ylabel = L"\chi_\ell", title = "Bond dimension")
+    scatterlines!(ax, 1:(R - 1), TCI.linkdims(tt))
     fig
-end
+end

examples/sparse_interpolation.jl

@@ -0,0 +1,71 @@
+using PolynomialQTT
+using CairoMakie
+import QuanticsGrids as QG
+import TensorCrossInterpolation as TCI
+
+α   = 0.1
+f(x) = α / sqrt(α^2 + (x - 0.5)^2)
+a, b = 0.0, 1.0
+R    = 10          # depth: 2^R = 1024 grid points
+N_dense  = 100     # polynomial degree for dense construction
+N_sparse = 200     # polynomial degree for sparse construction (larger N is affordable)
+
+# Evaluation grid
+grid        = QG.DiscretizedGrid{1}(R, a, b)
+quanticsinds = QG.grididx_to_quantics.(Ref(grid), 1:(2^R))
+plotx       = QG.grididx_to_origcoord.(Ref(grid), 1:(2^R))
+origdata    = f.(plotx)
+
+# Dense construction at N_dense
+tt_dense = PolynomialQTT.interpolatesinglescale(f, a, b, R, N_dense)
+err_dense = maximum(abs.(tt_dense.(quanticsinds) .- origdata))
+
+# Sparse constructions at N_sparse for varying M
+M_values = [5, 10, 15, 20, 30]
+tt_sparse = [PolynomialQTT.interpolatesinglescale_sparse(f, a, b, R, N_sparse, M) for M in M_values]
+errs_sparse = [maximum(abs.(tt.(quanticsinds) .- origdata)) for tt in tt_sparse]
+
+# ── Plot 1: the function ──────────────────────────────────────────────────────
+let
+    fig = Figure()
+    ax = Axis(fig[1, 1], xlabel = L"x", ylabel = L"f(x)",
+        title = L"Lorentzian peak: $f(x) = \alpha / \sqrt{\alpha^2 + (x-\frac{1}{2})^2}$, $\alpha = %$α$")
+    lines!(ax, plotx, origdata, linewidth = 2)
+    fig
+end
+
+# ── Plot 2: interpolation error vs bandwidth M ────────────────────────────────
+let
+    fig = Figure()
+    ax = Axis(fig[1, 1], xlabel = L"M \text{ (bandwidth)}", ylabel = "Max absolute error",
+        title = "Sparse construction error vs bandwidth  (N=$N_sparse)",
+        yscale = log10)
+    scatterlines!(ax, M_values, errs_sparse, label = "sparse  N=$N_sparse", linewidth = 2)
+    hlines!(ax, [err_dense], color = :black, linestyle = :dash,
+        label = "dense  N=$N_dense")
+    axislegend(ax, pos = :topright)
+    fig
+end
+
+# ── Plot 3: bond dimensions ───────────────────────────────────────────────────
+# Pick the best-performing sparse result (M = 10)
+M_best = 10
+tt_best = tt_sparse[findfirst(==(M_best), M_values)]
+
+let
+    fig = Figure()
+    ax = Axis(fig[1, 1], xlabel = L"\ell", ylabel = L"\chi_\ell",
+        title = "Bond dimensions  (α=$α)")
+    scatterlines!(ax, 1:(R - 1), TCI.linkdims(tt_dense),
+        label = "dense  N=$N_dense", linewidth = 2, marker = :circle)
+    scatterlines!(ax, 1:(R - 1), TCI.linkdims(tt_best),
+        label = "sparse  N=$N_sparse, M=$M_best", linewidth = 2, marker = :diamond)
+    axislegend(ax, pos = :topright)
+    fig
+end
+
+# ── Print summary ─────────────────────────────────────────────────────────────
+println("Dense  N=$N_dense:  max err = $err_dense  rank = $(TCI.rank(tt_dense))")
+for (M, tt, err) in zip(M_values, tt_sparse, errs_sparse)
+    println("Sparse N=$N_sparse M=$M:  max err = $err  rank = $(TCI.rank(tt))")
+end

src/PolynomialQTT.jl

@@ -7,5 +7,6 @@ import Base: in
 
 include("interval.jl")
 include("interpolation.jl")
+include("angular_lagrange.jl")
 
 end

src/angular_lagrange.jl

@@ -0,0 +1,123 @@
+
+"""
+For sparse interpolation
+"""
+function angular_local_lagrange(P::LagrangePolynomials{Float64}, M::Int)
+    N = length(P.grid) - 1
+    @assert N >= 2M "Need N ≥ 2M (got N=$N, M=$M)"
+
+    Ā = zeros(N + 1, 2, N + 1)   # indices: [α+1, σ+1, β+1]
+
+    for σ in 0:1, β in 0:N
+        x   = (Float64(σ) + P.grid[β + 1]) / 2
+        θ_x = acos(clamp(1 - 2x, -1.0, 1.0))
+
+        # Nearest angular grid index (grid is uniform: θ^γ = πγ/N)
+        ι = clamp(round(Int, θ_x * N / π), 0, N)
+
+        # Shift window to stay within [0, N] while keeping size 2M+1
+        ι_lo = clamp(ι - M, 0, N - 2M)
+        window = ι_lo:(ι_lo + 2M)          # exactly 2M+1 indices
+
+        # Angular positions of the window nodes
+        θ_win = [π * γ / N for γ in window]
+
+        # Evaluate the local Lagrange basis L^α(θ_x) for each α in the window
+        for (k, α) in enumerate(window)
+            L = one(Float64)
+            for j in eachindex(θ_win)
+                j == k && continue
+                L *= (θ_x - θ_win[j]) / (θ_win[k] - θ_win[j])
+            end
+            Ā[α + 1, σ + 1, β + 1] = L
+        end
+    end
+
+    return Ā
+end
+
+"""
+    interpolatesinglescale_sparse(f, a, b, numbits, polynomialdegree, M; ...)
+
+Sparse variant of `interpolatesinglescale` (Section 4.3).
+"""
+function interpolatesinglescale_sparse(
+        f,
+        a::Float64, b::Float64,
+        numbits::Int,
+        polynomialdegree::Int,
+        M::Int;
+        tolerance  = 1.0e-12,
+        maxbonddim = typemax(Int)
+    )
+    _scale(x) = a + (b - a) * x
+    P = getChebyshevGrid(polynomialdegree)
+
+    Aleft = [
+        f(_scale((sigma + P.grid[beta + 1]) / 2))
+            for sigma in [0, 1], beta in 0:polynomialdegree
+    ]
+    Acenter = angular_local_lagrange(P, M)
+    Aright  = [
+        P(alpha, sigma / 2)
+            for alpha in 0:polynomialdegree, sigma in [0, 1]
+    ]
+
+    train_ = [
+        reshape(Aleft, 1, 2, size(Aleft, 2)),
+        fill(Acenter, numbits - 2)...,
+        reshape(Aright, size(Aright, 1), 2, 1),
+    ]
+
+    return _compress_train(train_, tolerance, maxbonddim)
+end
+
+
+"""
+    interpolatesinglescale_sparse(f, a, b, numbits, polynomialdegree, M; ...)
+
+N-dimensional version of the sparse single-scale construction.  The sparse
+1-D core replaces `interpolationtensor(P)` inside the Kronecker product built
+by `_direct_product_coretensors`.
+"""
+function interpolatesinglescale_sparse(
+        f,
+        a::NTuple{D, Float64}, b::NTuple{D, Float64},
+        numbits::Int,
+        polynomialdegree::Int,
+        M::Int;
+        tolerance  = 1.0e-12,
+        maxbonddim = typemax(Int)
+    ) where {D}
+    _scale(x::NTuple{D, Float64})::NTuple{D, Float64} = tuple((a .+ (b .- a) .* x)...)
+
+    P = getChebyshevGrid(polynomialdegree)
+
+    Aleft = Array{Float64, 2D}(
+        undef, ntuple(_ -> 2, D)..., ntuple(_ -> polynomialdegree + 1, D)...
+    )
+    for sigmas in Iterators.product(ntuple(_ -> 0:1, D)...)
+        for betas in Iterators.product(ntuple(_ -> 0:polynomialdegree, D)...)
+            point = tuple(((sigmas[n] + P.grid[betas[n] + 1]) / 2 for n in 1:D)...)
+            Aleft[(sigmas .+ 1)..., (betas .+ 1)...] = f(_scale(point)...)
+        end
+    end
+
+    # Only this line differs from interpolatesinglescale: sparse core instead of dense
+    Acenter_ = angular_local_lagrange(P, M)
+    Acenter  = _direct_product_coretensors(fill(Acenter_, D))
+
+    Aright_ = [
+        P(alpha, sigma / 2)
+            for alpha in 0:polynomialdegree, sigma in [0, 1]
+    ]
+    Aright = _direct_product_coretensors(fill(reshape(Aright_, size(Aright_)..., 1), D))
+
+    train_ = [
+        reshape(Aleft, 1, 2^D, :),
+        fill(Acenter, numbits - 2)...,
+        Aright,
+    ]
+
+    return _compress_train(train_, tolerance, maxbonddim)
+end