Merge branch '44-use-full-pivoting-for-computing-site-tensors-for-tci2' into 'main'

A rare event of non-square pivot matrix

Closes #44

See merge request tensors4fields/TensorCrossInterpolation.jl!91

Author: shinaoka

src/matrixlu.jl

@@ -77,11 +77,20 @@ mutable struct rrLU{T}
     npivot::Int
     error::Float64
 
+    function rrLU(rowpermutation::Vector{Int}, colpermutation::Vector{Int}, L::Matrix{T}, U::Matrix{T}, leftorthogonal, npivot, error) where {T}
+        npivot == size(L, 2) || error("L must have the same number of columns as the number of pivots.")
+        npivot == size(U, 1) || error("U must have the same number of rows as the number of pivots.")
+        length(rowpermutation) == size(L, 1) || error("rowpermutation must have length equal to the number of pivots.")
+        length(colpermutation) == size(U, 2) || error("colpermutation must have length equal to the number of pivots.")
+        new{T}(rowpermutation, colpermutation, L, U, leftorthogonal, npivot, error)
+    end
+
     function rrLU{T}(nrows::Int, ncols::Int; leftorthogonal::Bool=true) where {T}
         new{T}(1:nrows, 1:ncols, zeros(nrows, 0), zeros(0, ncols), leftorthogonal, 0, NaN)
     end
 end
 
+
 function rrLU{T}(A::AbstractMatrix{T}; leftorthogonal::Bool=true) where {T}
     rrLU{T}(size(A)...; leftorthogonal=leftorthogonal)
 end
@@ -389,3 +398,70 @@ A special care is taken for a full-rank matrix: the last pivot error is set to z
 function lastpivoterror(lu::rrLU{T})::Float64 where {T}
     return lu.error
 end
+
+
+
+"""
+Solve (LU) x = b
+
+L: lower triangular matrix
+U: upper triangular matrix
+b: right-hand side vector
+
+Return x
+
+Note: Not optimized for performance
+"""
+function solve(L::Matrix{T}, U::Matrix{T}, b::Matrix{T}) where{T}
+    N1, N2, N3 = size(L, 1), size(L, 2), size(U, 2)
+    M = size(b, 2)
+    
+    # Solve Ly = b
+    y = zeros(T, N2, M)
+    for i = 1:N2
+        y[i, :] .= b[i, :]
+        for k in 1:M
+            for j in 1:i-1
+                y[i, k] -= L[i, j] * y[j, k]
+            end
+        end
+        y[i, :] ./= L[i, i]
+    end
+
+    # Solve Ux = y
+    x = zeros(T, N3, M)
+    for i = N3:-1:1
+        x[i, :] .= y[i, :]
+        for k in 1:M
+            for j = i+1:N3
+                x[i, k] -= U[i, j] * x[j, k]
+            end
+        end
+        x[i, :] ./= U[i, i]
+    end
+
+    return x
+end
+
+
+"""
+Override solving Ax = b using LU decomposition
+"""
+function Base.:\(A::rrLU{T}, b::AbstractMatrix{T}) where{T}
+    size(A, 1) == size(A, 2) || error("Matrix must be square.")
+    A.npivot == size(A, 1) || error("rank-deficient matrix is not supportred!")
+    b_perm = b[A.rowpermutation, :]
+    x_perm = solve(A.L, A.U, b_perm)
+    x = similar(x_perm)
+    for i in 1:size(x, 1)
+        x[A.colpermutation[i], :] .= x_perm[i, :]
+    end
+    return x
+end
+
+
+function Base.transpose(A::rrLU{T}) where{T}
+    return rrLU(
+        A.colpermutation, A.rowpermutation,
+        Matrix(transpose(A.U)), Matrix(transpose(A.L)), !A.leftorthogonal, A.npivot, A.error)
+end

src/tensorci2.jl

@@ -290,7 +290,7 @@ function kronecker(
     return MultiIndex[[i, js...] for i in 1:localdim, js in Jset][:]
 end
 
-function setT!(
+function setsitetensor!(
     tci::TensorCI2{ValueType}, b::Int, T::AbstractArray{ValueType,N}
 ) where {ValueType,N}
     tci.sitetensors[b] = reshape(
@@ -328,7 +328,7 @@ end
 # Backward compatibility
 const rmbadpivots! = sweep0site!
 
-function setT!(
+function setsitetensor!(
     tci::TensorCI2{ValueType}, f, b::Int; leftorthogonal=true
 ) where {ValueType}
     leftorthogonal || error("leftorthogonal==false is not supported!")
@@ -342,14 +342,16 @@ function setT!(
 
     if (leftorthogonal && b == length(tci)) ||
         (!leftorthogonal && b == 1)
-        setT!(tci, b, Pi1)
+        setsitetensor!(tci, b, Pi1)
         return tci.sitetensors[b]
     end
 
     P = reshape(
         filltensor(ValueType, f, tci.localdims, tci.Iset[b+1], tci.Jset[b], Val(0)),
         length(tci.Iset[b+1]), length(tci.Jset[b]))
+    length(tci.Iset[b+1]) == length(tci.Jset[b]) || error("Pivot matrix at bond $(b) is not square!")
 
+    #Tmat = transpose(transpose(rrlu(P)) \ transpose(Pi1))
     Tmat = transpose(transpose(P) \ transpose(Pi1))
     tci.sitetensors[b] = reshape(Tmat, length(tci.Iset[b]), tci.localdims[b], length(tci.Iset[b+1]))
     return tci.sitetensors[b]
@@ -396,7 +398,7 @@ function sweep1site!(
         tci.Iset[b+forwardsweep] = Is[rowindices(luci)]
         tci.Jset[b-!forwardsweep] = Js[colindices(luci)]
         if updatetensors
-            setT!(tci, b, forwardsweep ? left(luci) : right(luci))
+            setsitetensor!(tci, b, forwardsweep ? left(luci) : right(luci))
         end
         if any(isnan.(tci.sitetensors[b]))
             error("Error: NaN in tensor T[$b]")
@@ -417,7 +419,7 @@ function sweep1site!(
         end
         localtensor = reshape(filltensor(
             ValueType, f, tci.localdims, tci.Iset[lastupdateindex], tci.Jset[lastupdateindex], Val(1)), shape)
-        setT!(tci, lastupdateindex, localtensor)
+        setsitetensor!(tci, lastupdateindex, localtensor)
     end
     nothing
 end
@@ -560,8 +562,8 @@ function updatepivots!(
     tci.Iset[b+1] = Icombined[rowindices(luci)]
     tci.Jset[b] = Jcombined[colindices(luci)]
     if length(extraIset) == 0 && length(extraJset) == 0
-        setT!(tci, b, left(luci))
-        setT!(tci, b + 1, right(luci))
+        setsitetensor!(tci, b, left(luci))
+        setsitetensor!(tci, b + 1, right(luci))
     end
     updateerrors!(tci, b, pivoterrors(luci))
     nothing
@@ -744,24 +746,22 @@ function optimize!(
         end
     end
 
-    if rank(tci) > maxbonddim
-        errornormalization = normalizeerror ? tci.maxsamplevalue : 1.0
-        abstol = pivottolerance * errornormalization;
-        sweep2site!(
-            tci, f, 1;
-            abstol=abstol,
-            maxbonddim=maxbonddim,
-            pivotsearch=pivotsearch,
-            strictlynested=strictlynested,
-            verbosity=verbosity
-            )
-    end
+    # Extra one sweep by the 1-site update to
+    #  (1) Remove unnecessary pivots added by global pivots
+    #      Note: a pivot matrix can be non-square after adding global pivots,
+    #            or the bond dimension exceeds maxbonddim
+    #  (2) Compute site tensors
+    errornormalization = normalizeerror ? tci.maxsamplevalue : 1.0
+    abstol = pivottolerance * errornormalization;
+    sweep1site!(
+        tci,
+        f,
+        abstol=abstol,
+        maxbonddim=maxbonddim,
+    )
 
-    if !issitetensorsavailable(tci)
-        fillsitetensors!(tci, f)
-    end
+    _sanitycheck(tci)
 
-    errornormalization = normalizeerror ? tci.maxsamplevalue : 1.0
     return ranks, errors ./ errornormalization
 end
 
@@ -993,13 +993,17 @@ end
 
 function fillsitetensors!(
     tci::TensorCI2{ValueType}, f) where {ValueType}
-    #==
-    for b in 1:length(tci)-1
-       rmbadpivots!(tci, f, b)
-    end
-    ==#
     for b in 1:length(tci)
-       setT!(tci, f, b)
+       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

test/runtests.jl

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

test/test_matrixlu.jl

@@ -1,6 +1,7 @@
 using Test
 import TensorCrossInterpolation as TCI
 using LinearAlgebra
+using Random
 
 @testset "LU decomposition" begin
     @testset "Argmax finder" begin
@@ -151,9 +152,9 @@ using LinearAlgebra
             0.490269 0.810266 0.7946
         ]
         q = [
-            0.239552   0.306094  0.299063  0.0382492  0.185462  0.0334971  0.697561   0.389596  0.105665  0.0912763
-            0.0570609  0.56623   0.97183   0.994184   0.371695  0.284437   0.993251   0.902347  0.572944  0.0531369
-            0.45002    0.461168  0.6086    0.613702   0.543997  0.759954   0.0959818  0.638499  0.407382  0.482592
+            0.239552 0.306094 0.299063 0.0382492 0.185462 0.0334971 0.697561 0.389596 0.105665 0.0912763
+            0.0570609 0.56623 0.97183 0.994184 0.371695 0.284437 0.993251 0.902347 0.572944 0.0531369
+            0.45002 0.461168 0.6086 0.613702 0.543997 0.759954 0.0959818 0.638499 0.407382 0.482592
         ]
 
         A = p * q
@@ -175,11 +176,11 @@ using LinearAlgebra
 
     @testset "lastpivoterror for limited maxrank or tolerance" begin
         A = [
-            0.433088   0.956638   0.0907974  0.0447859  0.0196053
-            0.855517   0.782503   0.291197   0.540828   0.358579
-            0.37455    0.536457   0.205479   0.75896    0.701206
-            0.47272    0.0172539  0.518177   0.242864   0.461635
-            0.0676373  0.450878   0.672335   0.77726    0.540691
+            0.433088 0.956638 0.0907974 0.0447859 0.0196053
+            0.855517 0.782503 0.291197 0.540828 0.358579
+            0.37455 0.536457 0.205479 0.75896 0.701206
+            0.47272 0.0172539 0.518177 0.242864 0.461635
+            0.0676373 0.450878 0.672335 0.77726 0.540691
         ]
 
         lu = TCI.rrlu(A, maxrank=2)
@@ -208,4 +209,46 @@ using LinearAlgebra
         @test size(lu) == size(A)
         @test maximum(abs.(TCI.left(lu) * TCI.right(lu) .- A)) < 1e-3
     end
+
+
+    @testset "transpose" begin
+        Random.seed!(1234)
+        N1, N2, N3 = 5, 10, 3
+        A = rand(N1, N2)
+
+        tlu = transpose(TCI.rrlu(A))
+
+        @test TCI.left(tlu) * TCI.right(tlu) ≈ A'
+    end
+
+    @testset "solve by rrLU" begin
+        Random.seed!(1234)
+        N1, N2, N3 = 5, 5, 5
+        M = 2
+        L = tril(rand(N1, N2))
+        U = triu(rand(N2, N3))
+        b = rand(N1, M)
+
+        A = L * U
+        lua = TCI.rrlu(A)
+        @test TCI.left(lua) * TCI.right(lua) ≈ A
+
+        @test A * (lua \ b) ≈ b
+    end
+
+    #==
+    @testset "solve by rrLU (large matrix)" begin
+        M = 1000
+        N = 1000
+
+        L = qr(rand(N, N)).Q # Well-behaved matrix
+        U = qr(rand(N, N)).Q
+        b = rand(N, M)
+
+        A = L * U
+        lua = TCI.rrlu(A)
+        breconst = A * (lua \ b)
+        @test breconst ≈ b
+    end
+    ==#
 end