Merge branch '37-add-missing-features-described-in-the-code-release-paper' into 'main'

add GK nodes & weights

See merge request tensors4fields/TensorCrossInterpolation.jl!68

Author: shinaoka

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
@@ -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/tensorci2.jl

@@ -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
 

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)

test/test_integration.jl

@@ -0,0 +1,36 @@
+using Test
+import TensorCrossInterpolation as TCI
+import Random
+
+@testset "Integrate polynomials" begin
+    Random.seed!(1234)
+
+    coefficients = [
+        0.23637074801483304, 0.20661524945577847, 0.1850826417895819, 0.8433788714289417,
+        0.5801482873508491, 0.20339438932656262, 0.21593267492457668, 0.8052490409622802,
+        0.7189346124875339, 0.9400806688257749, 0.355210845205325, 0.5251561513473092,
+        0.6819965273401778, 0.9221987248861162, 0.04166444723413998
+    ]
+
+    polynomial(x::Float64) = sum(c * x^(i-1) for (i, c) in enumerate(coefficients))
+    polynomialintegral(x::Float64) = sum(c * x^i / i for (i, c) in enumerate(coefficients))
+    f(x::Vector{Float64}) = prod(polynomial.(x))
+
+    N = 5
+    exactval = polynomialintegral(1.0)^N
+    @test TCI.integrate(Float64, f, fill(0.0, N), fill(1.0, N)) ≈ exactval
+
+    b = rand(N)
+    a = rand(N)
+    exactval = prod(polynomialintegral.(b) .- polynomialintegral.(a))
+    @test TCI.integrate(Float64, f, a, b) ≈ exactval
+end
+
+@testset "Integrate 10d function" begin
+    function f(x)
+        return 1000 * cos(10 * sum(x .^ 2)) * exp(-sum(x)^4 / 1000)
+    end
+    I15 = TCI.integrate(Float64, f, fill(-1.0, 10), fill(+1.0, 10); GKorder=15, tolerance=1e-8)
+    Iref = -5.4960415218049
+    @test abs(I15 - Iref) < 1e-3
+end