Merge branch 'main' into 37-recompression-interface

Author: rittermarc

Project.toml

@@ -6,15 +6,16 @@ version = "0.9.0"
 [deps]
 EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 
 [compat]
 EllipsisNotation = "1"
 julia = "1.6"
 
 [extras]
-JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
+JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 QuanticsGrids = "634c7f73-3e90-4749-a1bd-001b8efc642d"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

docs/src/documentation.md

@@ -2,9 +2,6 @@
 
 Documentation of all types and methods in module [TensorCrossInterpolation](https://gitlab.com/tensors4fields/TensorCrossInterpolation.jl).
 
-```@index
-```
-
 ## Matrix approximation
 
 
@@ -47,6 +44,12 @@ Modules = [TensorCrossInterpolation]
 Pages = ["tensorci2.jl"]
 ```
 
+### Integration
+```@autodocs
+Modules = [TensorCrossInterpolation]
+Pages = ["integration.jl"]
+```
+
 ## Helpers and utility methods
 ```@autodocs
 Modules = [TensorCrossInterpolation]

docs/src/index.md

@@ -8,10 +8,6 @@ This is the documentation for [TensorCrossInterpolation](https://gitlab.com/tens
 
 With the user manual and usage examples below, users should be able to use this library as a "black box" in most cases. Detailed documentation of (almost) all methods can be found in the [Documentation](@ref) section, and [Implementation](@ref) contains a detailed explanation of this implementation of TCI.
 
-```@contents
-Pages = ["index.md", "documentation.md", "implementation.md"]
-```
-
 ## Interpolating functions
 
 The most convenient way to create a TCI is [`crossinterpolate2`](@ref). For example, consider the lorentzian in 5 dimensions, i.e. $f(\mathbf v) = 1/(1 + \mathbf v^2)$ on a mesh $\mathbf{v} \in \{1, 2, ..., 10\}^5$.
@@ -34,7 +30,7 @@ println("TCI approximation: $(tci([1, 2, 3, 4, 5]))")
 For easy integration into tensor network algorithms, the tensor train can be converted to ITensors MPS format. If you're using julia version 1.9 or later, an extension is automatically loaded if both `TensorCrossInterpolation.jl` and `ITensors.jl` are present.
 For older versions of julia, use the package using [TCIITensorConversion.jl](https://gitlab.com/tensors4fields/tciitensorconversion.jl).
 
-## Sums
+## Sums and Integrals
 
 Tensor trains are a way to efficiently obtain sums over all lattice sites, since this sum can be factorized:
 ```@example simple
@@ -47,6 +43,22 @@ println("Sum of tensor train: $sumtt")
 ```
 For further information, see [`sum`](@ref).
 
+This factorized sum can be used for efficient evaluation of high-dimensional integrals. This is implemented with Gauss-Kronrod quadrature rules in [`integrate`](@ref). For example, the integral
+```math
+I = 10^3 \int\limits_{[-1, +1]^{10}} d^{10} \vec{x} \,
+     \cos\!\left(10 \textstyle\sum_{n=1}^{10} x_n^2 \right)
+     \exp\!\left[-10^{-3}\left(\textstyle\sum_{n=1}^{10} x_n\right)^4\right]
+```
+is evaluated by the following code:
+```@example simple
+function f(x)
+    return 1e3 * cos(10 * sum(x .^ 2)) * exp(-sum(x)^4 / 1e3)
+end
+I = TCI.integrate(Float64, f, fill(-1.0, 10), fill(+1.0, 10); GKorder=15, tolerance=1e-8)
+println("GK15 integral value: $I")
+```
+The argument `GKorder` controls the Gauss-Kronrod quadrature rule used for the integration, and `tolerance` controls the tolerance in the TCI approximation, which is distinct from the tolerance in the integral. For complicated functions, it is recommended to integrate using two different GK rules and to compare the results to get a good estimate of the discretization error.
+
 ## Properties of the TCI object
 
 After running the code above, `tci` is a [`TensorCI2`](@ref) object that can be interrogated for various properties. The most important ones are the rank (i.e. maximum bond dimension) and the link dimensions:

src/TensorCrossInterpolation.jl

@@ -2,6 +2,7 @@ module TensorCrossInterpolation
 
 using LinearAlgebra
 using EllipsisNotation
+import QuadGK
 
 # To add a method for rank(tci)
 import LinearAlgebra: rank, diag
@@ -12,7 +13,7 @@ import Base: ==, +
 import Base: isempty, iterate, getindex, lastindex, broadcastable
 import Base: length, size, sum
 
-export crossinterpolate, crossinterpolate2, optfirstpivot
+export crossinterpolate1, crossinterpolate2, optfirstpivot
 export tensortrain
 
 include("util.jl")
@@ -31,5 +32,6 @@ include("tensorci1.jl")
 include("tensorci2.jl")
 include("tensortrain.jl")
 include("conversion.jl")
+include("integration.jl")
 
 end

src/integration.jl

@@ -0,0 +1,57 @@
+"""
+    function integrate(
+        ::Type{ValueType},
+        f,
+        a::Vector{ValueType},
+        b::Vector{ValueType};
+        tolerance=1e-8,
+        GKorder::Int=15
+    ) where {ValueType}
+
+Integrate the function `f` using TCI and Gauss--Kronrod quadrature rules.
+
+Arguments:
+- ValueType: return type of `f`.
+- a: Vector of lower bounds in each dimension. Effectively, the lower corner of the hypercube that is being integrated over.
+- b: Vector of upper bounds in each dimension.
+- tolerance: tolerance of the TCI approximation for the values of f.
+- GKorder: Order of the Gauss--Kronrod rule, e.g. 15.
+"""
+function integrate(
+    ::Type{ValueType},
+    f,
+    a::Vector{ValueType},
+    b::Vector{ValueType};
+    tolerance=1e-8,
+    GKorder::Int=15
+) where {ValueType}
+    if iseven(GKorder)
+        error("Gauss--Kronrod order must be odd, e.g. 15 or 61.")
+    end
+
+    if length(a) != length(b)
+        error("Integral bounds must have the same dimensionality, but got $(length(a)) lower bounds and $(length(b)) upper bounds.")
+    end
+
+    nodes1d, weights1d, _ = QuadGK.kronrod(GKorder ÷ 2, -1, +1)
+    nodes = @. (b - a) * (nodes1d' + 1) / 2 + a
+    weights = @. (b - a) * weights1d' / 2
+    normalization = GKorder^length(a)
+
+    localdims = fill(length(nodes1d), length(a))
+
+    function F(indices)
+        x = [nodes[n, i] for (n, i) in enumerate(indices)]
+        w = prod(weights[n, i] for (n, i) in enumerate(indices))
+        return w * f(x) * normalization
+    end
+
+    tci2, ranks, errors = crossinterpolate2(
+        ValueType,
+        F,
+        localdims;
+        tolerance
+    )
+
+    return sum(tci2) / normalization
+end

src/tensorci1.jl

@@ -1,7 +1,7 @@
 """
     mutable struct TensorCI1{ValueType} <: AbstractTensorTrain{ValueType}
 
-Type that represents tensor cross interpolations created using the TCI1 algorithm. Users may want to create these using [`crossinterpolate`](@ref) rather than calling a constructor directly.
+Type that represents tensor cross interpolations created using the TCI1 algorithm. Users may want to create these using [`crossinterpolate1`](@ref) rather than calling a constructor directly.
 """
 mutable struct TensorCI1{ValueType} <: AbstractTensorTrain{ValueType}
     Iset::Vector{IndexSet{MultiIndex}}
@@ -469,7 +469,7 @@ end
 
 
 @doc raw"""
-    function crossinterpolate(
+    function crossinterpolate1(
         ::Type{ValueType},
         f,
         localdims::Union{Vector{Int},NTuple{N,Int}},
@@ -504,7 +504,7 @@ Notes:
 
 See also: [`optfirstpivot`](@ref), [`CachedFunction`](@ref), [`crossinterpolate2`](@ref)
 """
-function crossinterpolate(
+function crossinterpolate1(
     ::Type{ValueType},
     f,
     localdims::Union{Vector{Int},NTuple{N,Int}},
@@ -551,3 +551,30 @@ function crossinterpolate(
     errornormalization = normalizeerror ? tci.maxsamplevalue : 1.0
     return tci, ranks, errors ./ errornormalization
 end
+
+@doc raw"""
+    function crossinterpolate(
+        ::Type{ValueType},
+        f,
+        localdims::Union{Vector{Int},NTuple{N,Int}},
+        firstpivot::MultiIndex=ones(Int, length(localdims));
+        tolerance::Float64=1e-8,
+        maxiter::Int=200,
+        sweepstrategy::Symbol=:backandforth,
+        pivottolerance::Float64=1e-12,
+        verbosity::Int=0,
+        additionalpivots::Vector{MultiIndex}=MultiIndex[],
+        normalizeerror::Bool=true
+    ) where {ValueType, N}
+
+Deprecated, and only included for backward compatibility. Please use [`crossinterpolate1`](@ref) instead.
+"""
+function crossinterpolate(
+    ::Type{ValueType},
+    f,
+    localdims::Union{Vector{Int},NTuple{N,Int}},
+    firstpivot::MultiIndex=ones(Int, length(localdims));
+    kwargs...
+) where {ValueType, N}
+    return crossinterpolate1(ValueType, f, localdims, firstpivot; kwargs...)
+end

src/tensorci2.jl

@@ -632,7 +632,7 @@ Notes:
 - By default, no caching takes place. Use the [`CachedFunction`](@ref) wrapper if your function is expensive to evaluate.
 
 
-See also: [`crossinterpolate2`](@ref), [`optfirstpivot`](@ref), [`CachedFunction`](@ref), [`crossinterpolate`](@ref)
+See also: [`crossinterpolate2`](@ref), [`optfirstpivot`](@ref), [`CachedFunction`](@ref), [`crossinterpolate1`](@ref)
 """
 function optimize!(
     tci::TensorCI2{ValueType},
@@ -778,8 +778,8 @@ function sweep2site!(
         extraIset = tci.Iset
         extraJset = tci.Jset
         if length(tci.Iset_history) > 0
-            extraIset = union(extraIset, tci.Iset_history[end])
-            extraJset = union(extraJset, tci.Jset_history[end])
+            extraIset = union.(extraIset, tci.Iset_history[end])
+            extraJset = union.(extraJset, tci.Jset_history[end])
         end
     end
 
@@ -871,7 +871,7 @@ Notes:
 - By default, no caching takes place. Use the [`CachedFunction`](@ref) wrapper if your function is expensive to evaluate.
 
 
-See also: [`optimize!`](@ref), [`optfirstpivot`](@ref), [`CachedFunction`](@ref), [`crossinterpolate`](@ref)
+See also: [`optimize!`](@ref), [`optfirstpivot`](@ref), [`CachedFunction`](@ref), [`crossinterpolate1`](@ref)
 """
 function crossinterpolate2(
     ::Type{ValueType},

src/tensortrain.jl

@@ -146,7 +146,9 @@ function to_tensors(obj::TensorTrainFit{ValueType}, x::Vector{ValueType}) where
     ]
 end
 
-_evaluate(tt, indexset) = only(prod(T[:, i, :] for (T, i) in zip(tt, indexset)))
+function _evaluate(tt::Vector{Array{V, 3}}, indexset) where {V}
+    only(prod(T[:, i, :] for (T, i) in zip(tt, indexset)))
+end
 
 function (obj::TensorTrainFit{ValueType})(x::Vector{ValueType}) where {ValueType}
     tensors = to_tensors(obj, x)

src/util.jl

@@ -69,11 +69,11 @@ Optimize the first pivot for a tensor cross interpolation.
 
 Arguments:
 - `f` is function to be interpolated.
-- `localdims::Union{Vector{Int},NTuple{N,Int}}` determines the local dimensions of the function parameters (see [`crossinterpolate`](@ref)).
+- `localdims::Union{Vector{Int},NTuple{N,Int}}` determines the local dimensions of the function parameters (see [`crossinterpolate1`](@ref)).
 - `fistpivot::MultiIndex=ones(Int, length(localdims))` is the starting point for the optimization. It is advantageous to choose it close to a global maximum of the function.
 - `maxsweep` is the maximum number of optimization sweeps. Default: `1000`.
 
-See also: [`crossinterpolate`](@ref)
+See also: [`crossinterpolate1`](@ref)
 """
 function optfirstpivot(
     f,

test/runtests.jl

@@ -1,7 +1,6 @@
 import TensorCrossInterpolation as TCI
 using Test
 using LinearAlgebra
-import ITensors
 
 include("test_with_aqua.jl")
 include("test_with_jet.jl")
@@ -17,6 +16,7 @@ include("test_batcheval.jl")
 include("test_tensorci1.jl")
 include("test_tensorci2.jl")
 include("test_tensortrain.jl")
+include("test_integration.jl")
 
 #==
 if VERSION.major >= 2 || (VERSION.major == 1 && VERSION.minor >= 9)