Merge branch '47-better-global-error-estimate' into 'main'

Implement estimatetrueerror

Closes #47

See merge request tensors4fields/TensorCrossInterpolation.jl!90

Author: shinaoka

docs/src/documentation.md

@@ -53,6 +53,6 @@ Pages = ["integration.jl"]
 ## Helpers and utility methods
 ```@autodocs
 Modules = [TensorCrossInterpolation]
-Pages = ["cachedfunction.jl", "batcheval.jl", "util.jl"]
+Pages = ["cachedfunction.jl", "batcheval.jl", "util.jl", "globalsearch.jl"]
 ```
 

docs/src/index.md

@@ -117,6 +117,16 @@ tci, ranks, errors = TCI.crossinterpolate2(
 ```
 This algorithm optimizes the given index set (in this case `[1, 2, 3, 4, 5]`) by searching for a maximum absolute value, alternating through the dimensions. If no starting point is given, `[1, 1, ...]` is used.
 
+## Estiamte true interpolation error by random global search
+Since the TCI update algorithms are local, the true interpolation error is not known. However, the error can be estimated by global searches. This is implemented in the function `estimatetrueerror`:
+
+```julia
+pivoterrors = TCI.estimatetrueerror(TCI.TensorTrain(tci), f)
+```
+
+This function approximately estimates the error that would be reached by repeating a greedy search from a random initial point.
+The result is a vector of a found indexset and the corresponding error, sorted by error. The error is the maximum absolute difference between the function and the TT approximation.
+
 ## Caching
 During constructing a TCI, the function to be interpolated can be evaluated for the same index set multiple times.
 If an evaluation of the function to be interpolated is costly, i.e., takes more than 100 ns,

src/TensorCrossInterpolation.jl

@@ -36,5 +36,6 @@ include("tensortrain.jl")
 include("conversion.jl")
 include("integration.jl")
 include("contraction.jl")
+include("globalsearch.jl")
 
 end

src/globalsearch.jl

@@ -0,0 +1,81 @@
+"""
+    function estimatetrueerror(
+        tt::TensorTrain{ValueType,3}, f;
+        nsearch::Int = 100,
+        initialpoints::Union{Nothing,AbstractVector{MultiIndex}} = nothing,
+    )::Vector{Tuple{MultiIndex,Float64}} where {ValueType}
+
+Estimate the global error by comparing the exact function value and
+the TT approximation using a greedy algorithm and returns a unique list of pairs of the pivot and the error (error is the absolute difference between the exact function value and the TT approximation).
+On return, the list is sorted in descending order of the error.
+
+Arguments:
+- `tt::TensorTrain{ValueType,3}`: The tensor train to be compared with `f`
+- `f`: The function to be compared with `tt`.
+- `nsearch::Int`: The number of initial points to be used in the search (defaults to 100).
+- `initialpoints::Union{Nothing,AbstractVector{MultiIndex}}`: The initial points to be used in the search (defaults to `nothing`). If `initialpoints` is not `nothing`, `nsearch` is ignored.
+"""
+function estimatetrueerror(
+    tt::TensorTrain{ValueType,3}, f;
+    nsearch::Int = 100,
+    initialpoints::Union{Nothing,AbstractVector{MultiIndex}} = nothing,
+)::Vector{Tuple{MultiIndex,Float64}} where {ValueType}
+    if nsearch <= 0 && initialpoints === nothing
+        error("No search is performed")
+    end
+    nsearch >= 0 || error("nsearch must be non-negative")
+
+    if nsearch > 0 && initialpoints === nothing
+        # Use random initial points
+        initialpoints = [[rand(1:first(d)) for d in sitedims(tt)] for _ in 1:nsearch]
+    end
+
+    ttcache = TTCache(tt)
+
+    pivoterror = [_floatingzone(ttcache, f, initp) for initp in initialpoints]
+
+    p = sortperm([e for (_, e) in pivoterror], rev=true)
+
+    return unique(pivoterror[p])
+end
+
+
+function _floatingzone(
+    ttcache::TTCache{ValueType}, f, pivot
+)::Tuple{MultiIndex,Float64} where {ValueType}
+    n = length(ttcache)
+
+    maxerror = abs(f(pivot) - ttcache(pivot))
+    localdims = first.(sitedims(ttcache))
+
+    while true
+        prev_maxerror = maxerror
+        for ipos in 1:n
+            exactdata = filltensor(
+                ValueType,
+                f,
+                localdims,
+                [pivot[1:ipos-1]],
+                [pivot[ipos+1:end]],
+                Val(1)
+            )
+            prediction = filltensor(
+                ValueType,
+                ttcache,
+                localdims,
+                [pivot[1:ipos-1]],
+                [pivot[ipos+1:end]],
+                Val(1)
+            )
+            err = vec(abs.(exactdata .- prediction))
+            pivot[ipos] = argmax(err)
+            maxerror = maximum(err)
+        end
+
+        if maxerror == prev_maxerror
+            break
+        end
+    end
+
+    return pivot, maxerror
+end

test/runtests.jl

@@ -2,6 +2,7 @@ import TensorCrossInterpolation as TCI
 using Test
 using LinearAlgebra
 
+#==
 include("test_with_aqua.jl")
 include("test_with_jet.jl")
 include("test_util.jl")
@@ -19,3 +20,5 @@ include("test_tensortrain.jl")
 include("test_conversion.jl")
 include("test_contraction.jl")
 include("test_integration.jl")
+==#
+include("test_globalsearch.jl")

test/test_globalsearch.jl

@@ -0,0 +1,37 @@
+using Test
+import TensorCrossInterpolation as TCI
+import TensorCrossInterpolation: crossinterpolate2, MultiIndex
+import Random
+import QuanticsGrids as QD
+
+@testset "globalsearch" begin
+    Random.seed!(1240)
+    R = 20
+    abstol = 1e-10
+
+    grid = QD.DiscretizedGrid{1}(R, (0.0,), (1.0,))
+
+    fx(x) = exp(-x) + 1e-3 * sin(1000 * x)
+    f(bitlist::MultiIndex) = fx(QD.quantics_to_origcoord(grid, bitlist)[1])
+
+    abstol = 1e-4
+    localdims = fill(2, R)
+    firstpivots = [ones(Int, R), vcat(1, fill(2, R - 1))]
+    tci, ranks, errors = crossinterpolate2(
+        Float64,
+        f,
+        localdims,
+        firstpivots;
+        tolerance=abstol,
+        maxbonddim=1,
+        verbosity=1,
+        normalizeerror=false,
+    )
+
+
+    pivoterrors = TCI.estimatetrueerror(TCI.TensorTrain(tci), f)
+
+    errors = [e for (_, e) in pivoterrors] 
+    @test all([abs(f(p) - tci(p)) for (p, _) in pivoterrors] .== errors)
+    @test all(errors[1:end-1] .>= errors[2:end]) # check if errors are sorted in descending order
+end