Merge branch '48-document-global-pivot-search-and-activate-it-by-default' into 'main'

Resolve "Document global pivot search and activate it by default"

See merge request tensors4fields/TensorCrossInterpolation.jl!92

Author: shinaoka

docs/src/index.md

@@ -117,8 +117,24 @@ 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.
 
+## Combing TCI2 and global pivot search
+The main algorithm for adding new pivots in TCI2 is the 2-site algorithm, which is local.
+The 2-site algorithm alone may miss some regions with high interpolation error.
+
+The current TCI2 implementation provides the combination of the 2-site algorithm and a global search algorithm to find such regions.
+This functionality is activated by default.
+In the function [`crossinterpolate2`](@ref), we alternate between a 2-site-update sweep and a global pivot insertion.
+After a 2-site-update sweep, we search for index sets with high interpolation errors (> the given tolerance multiplied by the parameter `tolmarginglobalsearch`) and add them to the TCI2 object, and then we continue with a 2-site-update sweep.
+
+The number of initial points used in one global search is controlled by the parameter `nsearchglobalpivot`.
+You may consider increasing this number if the global search is not effective (check the number of pivots found and timings of the global search by setting `verbosity` to a higher value!).
+The maximum number of global pivots inserted at once is controlled by the parameter `maxnglobalpivot`.
+
+A rare failure case is that the global search find the index sets with high interpolation errors, but the 2-site algorithm fails to add these pivots into the TCI2 object.
+This will end up adding the same index sets in the next global search, leading to an endless loop.
+
 ## 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`:
+Since most of 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`](@ref):
 
 ```julia
 pivoterrors = TCI.estimatetrueerror(TCI.TensorTrain(tci), f)

src/globalsearch.jl

@@ -32,7 +32,7 @@ function estimatetrueerror(
 
     ttcache = TTCache(tt)
 
-    pivoterror = [_floatingzone(ttcache, f, initp) for initp in initialpoints]
+    pivoterror = [_floatingzone(ttcache, f; initp=initp) for initp in initialpoints]
 
     p = sortperm([e for (_, e) in pivoterror], rev=true)
 
@@ -41,14 +41,25 @@ end
 
 
 function _floatingzone(
-    ttcache::TTCache{ValueType}, f, pivot
+    ttcache::TTCache{ValueType}, f;
+    earlystoptol::Float64 = typemax(Float64),
+    nsweeps=typemax(Int), initp::Union{Nothing,MultiIndex} = nothing
 )::Tuple{MultiIndex,Float64} where {ValueType}
+    nsweeps > 0 || error("nsweeps should be positive!")
+
+    localdims = first.(sitedims(ttcache))
+
     n = length(ttcache)
 
+    if initp === nothing
+        pivot = [rand(1:d) for d in localdims]
+    else
+        pivot = initp
+    end
+
     maxerror = abs(f(pivot) - ttcache(pivot))
-    localdims = first.(sitedims(ttcache))
 
-    while true
+    for isweep in 1:nsweeps
         prev_maxerror = maxerror
         for ipos in 1:n
             exactdata = filltensor(
@@ -72,10 +83,28 @@ function _floatingzone(
             maxerror = maximum(err)
         end
 
-        if maxerror == prev_maxerror
+        if maxerror == prev_maxerror || maxerror > earlystoptol # early stop
             break
         end
     end
 
     return pivot, maxerror
+end
+
+
+function fillsitetensors!(
+    tci::TensorCI2{ValueType}, f) where {ValueType}
+    for b in 1:length(tci)
+       setsitetensor!(tci, f, b)
+    end
+    nothing
+end
+
+
+function _sanitycheck(tci::TensorCI2{ValueType})::Bool where {ValueType}
+    for b in 1:length(tci)-1
+        length(tci.Iset[b+1]) == length(tci.Jset[b]) || error("Pivot matrix at bond $(b) is not square!")
+    end
+
+    return true
 end

src/integration.jl

@@ -23,7 +23,6 @@ function integrate(
     a::Vector{ValueType},
     b::Vector{ValueType};
     GKorder::Int=15,
-    normalizeerror=false,
     kwargs...
 ) where {ValueType}
     if iseven(GKorder)
@@ -52,7 +51,6 @@ function integrate(
         F,
         localdims;
         nsearchglobalpivot=10,
-        normalizeerror,
         kwargs...
     )
 

src/tensorci2.jl

@@ -214,7 +214,7 @@ function addglobalpivots2sitesweep!(
     pivotsearch::Symbol=:full,
     verbosity::Int=0,
     ntry::Int=10,
-    strictlynested::Bool=true
+    strictlynested::Bool=false
 )::Int where {F,ValueType}
     if any(length(tci) .!= length.(pivots))
         throw(DimensionMismatch("Please specify a pivot as one index per leg of the MPS."))
@@ -606,9 +606,9 @@ end
         normalizeerror::Bool=true,
         ncheckhistory=3,
         maxnglobalpivot::Int=5,
-        nsearchglobalpivot::Int=0,
+        nsearchglobalpivot::Int=5,
         tolmarginglobalsearch::Float64=10.0,
-        strictlynested::Bool=true
+        strictlynested::Bool=false
     ) where {ValueType}
 
 Perform optimization sweeps using the TCI2 algorithm. This will sucessively improve the TCI approximation of a function until it fits `f` with an error smaller than `tolerance`, or until the maximum bond dimension (`maxbonddim`) is reached.
@@ -629,7 +629,7 @@ Arguments:
 - `normalizeerror::Bool` determines whether to scale the error by the maximum absolute value of `f` found during sampling. If set to `false`, the algorithm continues until the *absolute* error is below `tolerance`. If set to `true`, the algorithm uses the absolute error divided by the maximum sample instead. This is helpful if the magnitude of the function is not known in advance. Default: `true`.
 - `ncheckhistory::Int` is the number of history points to use for convergence checks. Default: `3`.
 - `maxnglobalpivot::Int` can be set to `>= 0`. Default: `5`.
-- `nsearchglobalpivot::Int` can be set to `>= 0`. Default: `0`.
+- `nsearchglobalpivot::Int` can be set to `>= 0`. Default: `5`.
 - `tolmarginglobalsearch` can be set to `>= 1.0`. Seach global pivots where the interpolation error is larger than the tolerance by `tolmarginglobalsearch`.  Default: `10.0`.
 - `strictlynested::Bool` determines whether to preserve partial nesting in the TCI algorithm. Default: `false`.
 - `checkbatchevaluatable::Bool` Check if the function `f` is batch evaluatable. Default: `false`.
@@ -655,7 +655,7 @@ function optimize!(
     normalizeerror::Bool=true,
     ncheckhistory::Int=3,
     maxnglobalpivot::Int=5,
-    nsearchglobalpivot::Int=0,
+    nsearchglobalpivot::Int=5,
     tolmarginglobalsearch::Float64=10.0,
     strictlynested::Bool=false,
     checkbatchevaluatable::Bool=false
@@ -704,10 +704,11 @@ function optimize!(
             sweepstrategy=sweepstrategy,
             fillsitetensors=true
             )
-        if verbosity > 0 && length(globalpivots) > 0
-            nrejections = length([p for p in globalpivots if abs(evaluate(tci, p) - f(p)) > abstol])
+        if verbosity > 0 && length(globalpivots) > 0 && mod(iter, loginterval) == 0
+            abserr = [abs(evaluate(tci, p) - f(p)) for p in globalpivots]
+            nrejections = length(abserr .> abstol)
             if nrejections > 0
-                println("  Rejected $(nrejections) global pivots added in the previous iteration")
+                println("  Rejected $(nrejections) global pivots added in the previous iteration, errors are $(abserr)")
                 flush(stdout)
             end
         end
@@ -776,7 +777,7 @@ function sweep2site!(
     sweepstrategy::Symbol=:backandforth,
     pivotsearch::Symbol=:full,
     verbosity::Int=0,
-    strictlynested::Bool=true,
+    strictlynested::Bool=false,
     fillsitetensors::Bool=true
 ) where {ValueType}
     invalidatesitetensors!(tci)
@@ -850,9 +851,9 @@ end
         normalizeerror::Bool=true,
         ncheckhistory=3,
         maxnglobalpivot::Int=5,
-        nsearchglobalpivot::Int=0,
+        nsearchglobalpivot::Int=5,
         tolmarginglobalsearch::Float64=10.0,
-        strictlynested::Bool=true
+        strictlynested::Bool=false
     ) where {ValueType,N}
 
 Cross interpolate a function ``f(\mathbf{u})`` using the TCI2 algorithm. Here, the domain of ``f`` is ``\mathbf{u} \in [1, \ldots, d_1] \times [1, \ldots, d_2] \times \ldots \times [1, \ldots, d_{\mathscr{L}}]`` and ``d_1 \ldots d_{\mathscr{L}}`` are the local dimensions.
@@ -873,9 +874,9 @@ Arguments:
 - `normalizeerror::Bool` determines whether to scale the error by the maximum absolute value of `f` found during sampling. If set to `false`, the algorithm continues until the *absolute* error is below `tolerance`. If set to `true`, the algorithm uses the absolute error divided by the maximum sample instead. This is helpful if the magnitude of the function is not known in advance. Default: `true`.
 - `ncheckhistory::Int` is the number of history points to use for convergence checks. Default: `3`.
 - `maxnglobalpivot::Int` can be set to `>= 0`. Default: `5`.
-- `nsearchglobalpivot::Int` can be set to `>= 0`. Default: `0`.
+- `nsearchglobalpivot::Int` can be set to `>= 0`. Default: `5`.
 - `tolmarginglobalsearch` can be set to `>= 1.0`. Seach global pivots where the interpolation error is larger than the tolerance by `tolmarginglobalsearch`.  Default: `10.0`.
-- `strictlynested::Bool=true` determines whether to preserve partial nesting in the TCI algorithm. Default: `true`.
+- `strictlynested::Bool=false` determines whether to preserve partial nesting in the TCI algorithm. Default: `true`.
 - `checkbatchevaluatable::Bool` Check if the function `f` is batch evaluatable. Default: `false`.
 
 Notes:
@@ -916,8 +917,9 @@ function searchglobalpivots(
     end
 
     pivots = Dict{Float64,MultiIndex}()
+    ttcache = TTCache(tci)
     for _ in 1:nsearch
-        pivot, error = _floatingzone(tci, f, 10 * abstol)
+        pivot, error = _floatingzone(ttcache, f; earlystoptol = 10 * abstol, nsweeps=100)
         if error > abstol
             pivots[error] = pivot
         end
@@ -941,69 +943,3 @@ function searchglobalpivots(
     return [p for (_,p) in pivots]
 end
 
-
-function _floatingzone(
-    tci::TensorCI2{ValueType}, f, abstol;
-    nsweeps=100
-)::Tuple{MultiIndex,Float64} where {ValueType}
-    nsweeps > 0 || error("nsweeps should be positive!")
-
-    localdims = tci.localdims
-
-    n = length(tci)
-
-    ttcache = TTCache(tci)
-
-    pivot = [rand(1:d) for d in localdims]
-
-    maxerror = abs(f(pivot) - ttcache(pivot))
-
-    for isweep in 1:nsweeps
-        prev_maxerror = maxerror
-        for ipos in 1:n
-            exactdata = filltensor(
-                ValueType,
-                f,
-                tci.localdims,
-                [pivot[1:ipos-1]],
-                [pivot[ipos+1:end]],
-                Val(1)
-            )
-            prediction = filltensor(
-                ValueType,
-                ttcache,
-                tci.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 || maxerror > abstol # early stop
-            break
-        end
-    end
-
-    return pivot, maxerror
-end
-
-
-function fillsitetensors!(
-    tci::TensorCI2{ValueType}, f) where {ValueType}
-    for b in 1:length(tci)
-       setsitetensor!(tci, f, b)
-    end
-    nothing
-end
-
-
-function _sanitycheck(tci::TensorCI2{ValueType})::Bool where {ValueType}
-    for b in 1:length(tci)-1
-        length(tci.Iset[b+1]) == length(tci.Jset[b]) || error("Pivot matrix at bond $(b) is not square!")
-    end
-
-    return true
-end