Revert TaylorCluster Optimizations of #199

Project.toml

@@ -22,7 +22,6 @@ RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
 TensorKitManifolds = "11fa318c-39cb-4a83-b1ed-cdc7ba1e3684"
 Transducers = "28d57a85-8fef-5791-bfe6-a80928e7c999"
-TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
 VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 
 [compat]
@@ -50,7 +49,6 @@ TensorOperations = "5"
 Test = "1"
 TestExtras = "0.3"
 Transducers = "0.4"
-TupleTools = "1.6.0"
 VectorInterface = "0.2, 0.3, 0.4, 0.5"
 julia = "1.10"
 

src/MPSKit.jl

@@ -12,7 +12,6 @@ using RecipesBase
 using VectorInterface
 using Accessors
 using HalfIntegers
-import TupleTools as TT
 using DocStringExtensions
 
 using LinearAlgebra: diag, Diagonal

src/algorithms/timestep/timeevmpo.jl

@@ -38,12 +38,13 @@
     "include higher-order corrections"
     extension::Bool = true
     "approximate compression of corrections, accurate up to order `N`"
-    compression::Bool = true
+    compression::Bool = false
 end
 
 const WI = TaylorCluster(; N=1, extension=false, compression=false)
 
-function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster)
+function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster;
+                       tol=eps(real(scalartype(H)))^(3 / 4))
     N = alg.N
     τ = -1im * dt
 
@@ -80,73 +81,65 @@
     if alg.extension
         H_next = H_n * H′
         linds_next = LinearIndices(ntuple(i -> V, N + 1))
-        cinds_next = CartesianIndices(linds_next)
         for (i, slice) in enumerate(parent(H_n))
-            for Inext in nonzero_keys(H_next[i])
-                aₑ = cinds_next[Inext[1]]
-                bₑ = cinds_next[Inext[4]]
-                for c in findall(==(1), aₑ.I), d in findall(==(V), bₑ.I)
-                    a = TT.deleteat(Tuple(aₑ), c)
-                    b = TT.deleteat(Tuple(bₑ), d)
-                    if all(>(1), b) && !(all(in((1, V), a)) && any(==(V), a))
-                        n1 = count(==(1), a) + 1
-                        n3 = count(==(V), b) + 1
-                        factor = τ * factorial(N) / (factorial(N + 1) * n1 * n3)
-                        slice[linds[a...], 1, 1, linds[b...]] += factor * H_next[i][Inext]
-                    end
+            for a in cinds, b in cinds
+                all(>(1), b.I) || continue
+                all(in((1, V)), a.I) && any(==(V), a.I) && continue
+
+                n1 = count(==(1), a.I) + 1
+                n3 = count(==(V), b.I) + 1
+                factor = τ * factorial(N) / (factorial(N + 1) * n1 * n3)
+
+                for c in 1:(N + 1), d in 1:(N + 1)
+                    aₑ = insert!([a.I...], c, 1)
+                    bₑ = insert!([b.I...], d, V)
+
+                    # TODO: use VectorInterface for memory efficiency
+                    slice[linds[a], 1, 1, linds[b]] += factor *
+                                                       H_next[i][linds_next[aₑ...], 1, 1,
+                                                                 linds_next[bₑ...]]
                 end
             end
         end
     end
 
     # loopback step: Algorithm 1
     # constructing the Nth order time evolution MPO
-    if H isa FiniteMPOHamiltonian
-        mpo = FiniteMPO(parent(H_n))
-    else
-        mpo = InfiniteMPO(parent(H_n))
-    end
+    mpo = MPO(parent(H_n))
     for slice in parent(mpo)
-        for (I, v) in nonzero_pairs(slice)
-            a = cinds[I[1]]
-            b = cinds[I[4]]
-            if all(in((1, V)), a.I) && any(!=(1), a.I)
-                delete!(slice, I)
-            elseif any(!=(1), b.I) && all(in((1, V)), b.I)
-                exponent = count(==(V), b.I)
-                factor = τ^exponent * factorial(N - exponent) / factorial(N)
-                Idst = CartesianIndex(Base.setindex(I.I, 1, 4))
-                slice[Idst] += factor * v
-                delete!(slice, I)
+        for b in cinds[2:end]
+            all(in((1, V)), b.I) || continue
+
+            b_lin = linds[b]
+            a = count(==(V), b.I)
+            factor = τ^a * factorial(N - a) / factorial(N)
+            slice[:, 1, 1, 1] = slice[:, 1, 1, 1] + factor * slice[:, 1, 1, b_lin]
+            for I in nonzero_keys(slice)
+                (I[1] == b_lin || I[4] == b_lin) && delete!(slice, I)
             end
         end
     end
 
     # Remove equivalent rows and columns: Algorithm 2
     for slice in parent(mpo)
-        for I in nonzero_keys(slice)
-            a = cinds[I[1]]
-            a_sort = CartesianIndex(sort(collect(a.I); by=(!=(1)))...)
-            n1_a = count(==(1), a.I)
-            n3_a = count(==(V), a.I)
-            if n3_a <= n1_a && a_sort != a
-                Idst = CartesianIndex(Base.setindex(I.I, linds[a_sort], 1))
-                slice[Idst] += slice[I]
-                delete!(slice, I)
-            elseif a != a_sort
-                delete!(slice, I)
-            end
+        for c in cinds
+            c_lin = linds[c]
+            s_c = CartesianIndex(sort(collect(c.I); by=(!=(1)))...)
+            s_r = CartesianIndex(sort(collect(c.I); by=(!=(V)))...)
+
+            n1 = count(==(1), c.I)
+            n3 = count(==(V), c.I)
 
-            b = cinds[I[4]]
-            b_sort = CartesianIndex(sort(collect(b.I); by=(!=(V)))...)
-            n1_b = count(==(1), b.I)
-            n3_b = count(==(V), b.I)
-            if n3_b > n1_b && b_sort != b
-                Idst = CartesianIndex(Base.setindex(I.I, linds[b_sort], 4))
-                slice[Idst] += slice[I]
-                delete!(slice, I)
-            elseif b != b_sort
-                delete!(slice, I)
+            if n3 <= n1 && s_c != c
+                slice[linds[s_c], 1, 1, :] += slice[c_lin, 1, 1, :]
+                for I in nonzero_keys(slice)
+                    (I[1] == c_lin || I[4] == c_lin) && delete!(slice, I)
+                end
+            elseif n3 > n1 && s_r != c
+                slice[:, 1, 1, linds[s_r]] += slice[:, 1, 1, c_lin]
+                for I in nonzero_keys(slice)
+                    (I[1] == c_lin || I[4] == c_lin) && delete!(slice, I)
+                end
             end
         end
     end
@@ -161,14 +154,10 @@
                 b = CartesianIndex(replace(a.I, V => 1))
                 b_lin = linds[b]
                 factor = τ^n1 * factorial(N - n1) / factorial(N)
+                slice[:, 1, 1, b_lin] += factor * slice[:, 1, 1, a_lin]
+
                 for I in nonzero_keys(slice)
-                    if I[4] == a_lin
-                        Idst = CartesianIndex(Base.setindex(I.I, b_lin, 4))
-                        slice[Idst] += factor * slice[I]
-                        delete!(slice, I)
-                    elseif I[1] == a_lin
-                        delete!(slice, I)
-                    end
+                    (I[1] == a_lin || I[4] == a_lin) && delete!(slice, I)
                 end
             end
         end
@@ -180,7 +169,7 @@
         mpo[end] = mpo[end][:, :, :, 1]
     end
 
-    return remove_orphans!(mpo)
+    return remove_orphans!(mpo; tol=tol)
 end
 
 has_prod_elem(slice, t1, t2) = all(map(x -> contains(slice, x...), zip(t1, t2)))
@@ -223,65 +212,6 @@
     end
 end
 
-# function make_time_mpo(H::MPOHamiltonian{S,T}, dt, alg::WII) where {S,T}
-#     WA = PeriodicArray{T,3}(undef, H.period, H.odim - 2, H.odim - 2)
-#     WB = PeriodicArray{T,2}(undef, H.period, H.odim - 2)
-#     WC = PeriodicArray{T,2}(undef, H.period, H.odim - 2)
-#     WD = PeriodicArray{T,1}(undef, H.period)
-#
-#     δ = dt * (-1im)
-#
-#     for i in 1:(H.period), j in 2:(H.odim - 1), k in 2:(H.odim - 1)
-#         init_1 = isometry(storagetype(H[i][1, H.odim]), codomain(H[i][1, H.odim]),
-#                           domain(H[i][1, H.odim]))
-#         init = [init_1, zero(H[i][1, k]), zero(H[i][j, H.odim]), zero(H[i][j, k])]
-#
-#         (y, convhist) = exponentiate(1.0, RecursiveVec(init),
-#                                      Arnoldi(; tol=alg.tol, maxiter=alg.maxiter)) do x
-#             out = similar(x.vecs)
-#
-#             @plansor out[1][-1 -2; -3 -4] := δ * x[1][-1 1; -3 -4] *
-#                                              H[i][1, H.odim][2 3; 1 4] * τ[-2 4; 2 3]
-#
-#             @plansor out[2][-1 -2; -3 -4] := δ * x[2][-1 1; -3 -4] *
-#                                              H[i][1, H.odim][2 3; 1 4] * τ[-2 4; 2 3]
-#             @plansor out[2][-1 -2; -3 -4] += sqrt(δ) * x[1][1 2; -3 4] *
-#                                              H[i][1, k][-1 -2; 3 -4] * τ[3 4; 1 2]
-#
-#             @plansor out[3][-1 -2; -3 -4] := δ * x[3][-1 1; -3 -4] *
-#                                              H[i][1, H.odim][2 3; 1 4] * τ[-2 4; 2 3]
-#             @plansor out[3][-1 -2; -3 -4] += sqrt(δ) * x[1][1 2; -3 4] *
-#                                              H[i][j, H.odim][-1 -2; 3 -4] * τ[3 4; 1 2]
-#
-#             @plansor out[4][-1 -2; -3 -4] := δ * x[4][-1 1; -3 -4] *
-#                                              H[i][1, H.odim][2 3; 1 4] * τ[-2 4; 2 3]
-#             @plansor out[4][-1 -2; -3 -4] += x[1][1 2; -3 4] * H[i][j, k][-1 -2; 3 -4] *
-#                                              τ[3 4; 1 2]
-#             @plansor out[4][-1 -2; -3 -4] += sqrt(δ) * x[2][1 2; -3 -4] *
-#                                              H[i][j, H.odim][-1 -2; 3 4] * τ[3 4; 1 2]
-#             @plansor out[4][-1 -2; -3 -4] += sqrt(δ) * x[3][-1 4; -3 3] *
-#                                              H[i][1, k][2 -2; 1 -4] * τ[3 4; 1 2]
-#
-#             return RecursiveVec(out)
-#         end
-#         convhist.converged == 0 &&
-#             @warn "exponentiate failed to converge: normres = $(convhist.normres)"
-#
-#         WA[i, j - 1, k - 1] = y[4]
-#         WB[i, j - 1] = y[3]
-#         WC[i, k - 1] = y[2]
-#         WD[i] = y[1]
-#     end
-#
-#     W2 = PeriodicArray{Union{T,Missing},3}(missing, H.period, H.odim - 1, H.odim - 1)
-#     W2[:, 2:end, 2:end] = WA
-#     W2[:, 2:end, 1] = WB
-#     W2[:, 1, 2:end] = WC
-#     W2[:, 1, 1] = WD
-#
-#     return SparseMPO(W2)
-# end
-
 function make_time_mpo(H::InfiniteMPOHamiltonian{T}, dt, alg::WII) where {T}
     WA = H.A
     WB = H.B

test/algorithms.jl

@@ -769,4 +769,41 @@ end
     end
 end
 
+@testset "TaylorCluster time evolution" begin
+    L = 4
+    Hs = [transverse_field_ising(; L=L), heisenberg_XXX(; L=L)]
+
+    Ns = [1, 2, 3]
+    dts = [1e-2, 1e-3]
+    for H in Hs
+        ψ = FiniteMPS(L, physicalspace(H, 1), ℂ^16)
+        for N in Ns
+            εs = zeros(ComplexF64, 2)
+            for (i, dt) in enumerate(dts)
+                ψ₀, _ = find_groundstate(ψ, H, DMRG(; verbosity=0))
+                E₀ = expectation_value(ψ₀, H)
+
+                O = make_time_mpo(H, dt, TaylorCluster(; N=N))
+
+                ψ₁, _ = approximate(ψ₀, (O, ψ₀), DMRG(; verbosity=0))
+                εs[i] = norm(dot(ψ₀, ψ₁) - exp(-im * E₀ * dt))
+            end
+            @test (log(εs[2]) - log(εs[1])) / (log(dts[2]) - log(dts[1])) ≈ N + 1 atol = 0.1
+        end
+
+        for N in Ns
+            εs = zeros(ComplexF64, 2)
+            for (i, dt) in enumerate(dts)
+                ψ₀, _ = find_groundstate(ψ, H, DMRG(; verbosity=0))
+                E₀ = expectation_value(ψ₀, H)
+
+                O = make_time_mpo(H, dt, TaylorCluster(; N=N, compression=true))
+
+                ψ₁, _ = approximate(ψ₀, (O, ψ₀), DMRG(; verbosity=0))
+                εs[i] = norm(dot(ψ₀, ψ₁) - exp(-im * E₀ * dt))
+            end
+            @test (log(εs[2]) - log(εs[1])) / (log(dts[2]) - log(dts[1])) ≈ N atol = 0.1
+        end
+    end
+end
 end