Skip to content

Commit

Permalink
Revert TaylorCluster Optimizations of #199 (#220)
Browse files Browse the repository at this point in the history 
* revert broken TaylorCluster optimizations

* make mpo using `MPO(parent(H))` instead of checking for finite or infinite case

* remove commented-out code, copy master branch changes

* remove TupleTools

* add scaling tests for finite time evolution with TaylorClusters
  • Loading branch information
@VictorVanthilt
VictorVanthilt authored Jan 17, 2025
1 parent 925edc4 commit e42e898
Show file tree
Hide file tree
Showing 4 changed files with 87 additions and 123 deletions.
2 changes: 0 additions & 2 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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"

Expand Down
1 change: 0 additions & 1 deletion src/MPSKit.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@ using RecipesBase
using VectorInterface
using Accessors
using HalfIntegers
import TupleTools as TT
using DocStringExtensions

using LinearAlgebra: diag, Diagonal
Expand Down
170 changes: 50 additions & 120 deletions src/algorithms/timestep/timeevmpo.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,12 +38,13 @@ $(TYPEDFIELDS)
"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

Expand Down Expand Up @@ -80,73 +81,65 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster)
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
Expand All @@ -161,14 +154,10 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster)
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
Expand All @@ -180,7 +169,7 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster)
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)))
Expand Down Expand Up @@ -223,65 +212,6 @@ function interweave(a, b)
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
Expand Down
37 changes: 37 additions & 0 deletions test/algorithms.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

0 comments on commit e42e898

Please sign in to comment.