Skip to content

Commit cb9ffcc

Browse files
mtfishmanmtfishman
authored
Normalize exponentiate to fix TDVP, add more tests (#43)
* Fix imaginary TDVP by normalizing exponentiate output * Add tests to loop over unit cell size, QN conservation, time step
1 parent 76c32a2 commit cb9ffcc

File tree

4 files changed

+261
-123
lines changed

4 files changed

+261
-123
lines changed

‎examples/vumps/transfer_matrix_spectrum.jl‎

Lines changed: 197 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,197 @@
1+
using ITensors
2+
using ITensorInfiniteMPS
3+
4+
##############################################################################
5+
# VUMPS parameters
6+
#
7+
8+
maxdim = 10 # Maximum bond dimension
9+
cutoff = 1e-6 # Singular value cutoff when increasing the bond dimension
10+
max_vumps_iters = 20 # Maximum number of iterations of the VUMPS algorithm at a fixed bond dimension
11+
tol = 1e-5 # Precision error tolerance for outer loop of VUMPS or TDVP
12+
outer_iters = 5 # Number of times to increase the bond dimension
13+
time_step = -Inf # -Inf corresponds to VUMPS
14+
solver_tol = (x -> x / 100) # Tolerance for the local solver (eigsolve in VUMPS and exponentiate in TDVP)
15+
multisite_update_alg = "parallel" # Choose between ["sequential", "parallel"]
16+
conserve_qns = true
17+
N = 2 # Number of sites in the unit cell (1 site unit cell is currently broken)
18+
19+
# Parameters of the transverse field Ising model
20+
model_params = (J=1.0, h=0.9)
21+
22+
##############################################################################
23+
# CODE BELOW HERE DOES NOT NEED TO BE MODIFIED
24+
#
25+
26+
model = Model"ising"()
27+
28+
function space_shifted(::Model"ising", q̃sz; conserve_qns=true)
29+
if conserve_qns
30+
return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1]
31+
else
32+
return [QN() => 2]
33+
end
34+
end
35+
36+
space_ = fill(space_shifted(model, 0; conserve_qns=conserve_qns), N)
37+
s = infsiteinds("S=1/2", N; space=space_)
38+
initstate(n) = "↑"
39+
ψ = InfMPS(s, initstate)
40+
41+
# Form the Hamiltonian
42+
H = InfiniteITensorSum(model, s; model_params...)
43+
44+
# Check translational invariance
45+
@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))
46+
47+
vumps_kwargs = (
48+
tol=tol,
49+
maxiter=max_vumps_iters,
50+
solver_tol=solver_tol,
51+
multisite_update_alg=multisite_update_alg,
52+
)
53+
subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
54+
#ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...)
55+
56+
# Alternate steps of running VUMPS and increasing the bond dimension
57+
@time for _ in 1:outer_iters
58+
println("\nIncrease bond dimension")
59+
global ψ = subspace_expansion(ψ, H; subspace_expansion_kwargs...)
60+
println("Run VUMPS with new bond dimension")
61+
global ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...)
62+
end
63+
64+
# Check translational invariance
65+
@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))
66+
67+
#
68+
# Compare to DMRG
69+
#
70+
71+
Nfinite = 100
72+
sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
73+
Hfinite = MPO(model, sfinite; model_params...)
74+
ψfinite = randomMPS(sfinite, initstate)
75+
@show flux(ψfinite)
76+
sweeps = Sweeps(10)
77+
setmaxdim!(sweeps, maxdim)
78+
setcutoff!(sweeps, cutoff)
79+
energy_finite_total, ψfinite = @time dmrg(Hfinite, ψfinite, sweeps)
80+
@show energy_finite_total / Nfinite
81+
82+
function energy_local(ψ1, ψ2, h)
83+
ϕ = ψ1 * ψ2
84+
return (noprime(ϕ * h) * dag(ϕ))[]
85+
end
86+
87+
function ITensors.expect(ψ, o)
88+
return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[]
89+
end
90+
91+
# Exact energy at criticality: 4/pi = 1.2732395447351628
92+
93+
nfinite = Nfinite ÷ 2
94+
orthogonalize!(ψfinite, nfinite)
95+
hnfinite = ITensor(model, sfinite[nfinite], sfinite[nfinite + 1]; model_params...)
96+
energy_finite = energy_local(ψfinite[nfinite], ψfinite[nfinite + 1], hnfinite)
97+
energy_infinite = energy_local(ψ.AL[1], ψ.AL[2] * ψ.C[2], H[(1, 2)])
98+
@show energy_finite, energy_infinite
99+
@show abs(energy_finite - energy_infinite)
100+
101+
energy_exact = reference(model, Observable("energy"); model_params...)
102+
@show energy_exact
103+
104+
Sz1_finite = expect(ψfinite[nfinite], "Sz")
105+
orthogonalize!(ψfinite, nfinite + 1)
106+
Sz2_finite = expect(ψfinite[nfinite + 1], "Sz")
107+
Sz1_infinite = expect(ψ.AL[1] * ψ.C[1], "Sz")
108+
Sz2_infinite = expect(ψ.AL[2] * ψ.C[2], "Sz")
109+
110+
@show Sz1_finite, Sz2_finite
111+
@show Sz1_infinite, Sz2_infinite
112+
113+
##############################################################################
114+
# Compute eigenspace of the transfer matrix
115+
#
116+
117+
using Arpack
118+
using KrylovKit
119+
using LinearAlgebra
120+
121+
T = TransferMatrix(ψ.AL)
122+
Tᵀ = transpose(T)
123+
vⁱᴿ = randomITensor(dag(input_inds(T)))
124+
vⁱᴸ = randomITensor(dag(input_inds(Tᵀ)))
125+
126+
neigs = 10
127+
tol = 1e-10
128+
λ⃗ᴿ, v⃗ᴿ, right_info = eigsolve(T, vⁱᴿ, neigs, :LM; tol=tol)
129+
λ⃗ᴸ, v⃗ᴸ, left_info = eigsolve(Tᵀ, vⁱᴸ, neigs, :LM; tol=tol)
130+
131+
@show norm(T(v⃗ᴿ[1]) - λ⃗ᴿ[1] * v⃗ᴿ[1])
132+
@show norm(Tᵀ(v⃗ᴸ[1]) - λ⃗ᴸ[1] * v⃗ᴸ[1])
133+
134+
@show λ⃗ᴿ
135+
@show λ⃗ᴸ
136+
@show flux.(v⃗ᴿ)
137+
138+
neigs = length(v⃗ᴿ)
139+
140+
# Normalize the vectors
141+
N⃗ = [(translatecell(v⃗ᴸ[n], 1) * v⃗ᴿ[n])[] for n in 1:neigs]
142+
143+
v⃗ᴿ ./= sqrt.(N⃗)
144+
v⃗ᴸ ./= sqrt.(N⃗)
145+
146+
# Form a second starting vector orthogonal to v⃗ᴿ[1]
147+
# This doesn't work. TODO: project out v⃗ᴿ[1], v⃗ᴸ[1] from T
148+
#λ⃗ᴿ², v⃗ᴿ², right_info_2 = eigsolve(T, vⁱᴿ², neigs, :LM; tol=tol)
149+
150+
# Projector onto the n-th eigenstate
151+
function proj(v⃗ᴸ, v⃗ᴿ, n)
152+
Lⁿ = v⃗ᴸ[n]
153+
Rⁿ = v⃗ᴿ[n]
154+
return ITensorMap(
155+
[translatecell(Lⁿ, 1), translatecell(Rⁿ, -1)]; input_inds=inds(Rⁿ), output_inds=inds(Lⁿ)
156+
)
157+
end
158+
159+
P⃗ = [proj(v⃗ᴸ, v⃗ᴿ, n) for n in 1:neigs]
160+
T⁻P = T - sum(P⃗)
161+
162+
#vⁱᴿ² = vⁱᴿ - (translatecell(v⃗ᴸ[1], 1) * vⁱᴿ)[] / norm(v⃗ᴿ[1]) * v⃗ᴿ[1]
163+
#@show norm(dag(vⁱᴿ²) * v⃗ᴿ[1])
164+
165+
λ⃗ᴾᴿ, v⃗ᴾᴿ, right_info = eigsolve(T⁻P, vⁱᴿ, neigs, :LM; tol=tol)
166+
@show λ⃗ᴾᴿ
167+
168+
vⁱᴿ⁻ᵈᵃᵗᵃ = vec(array(vⁱᴿ))
169+
λ⃗ᴿᴬ, v⃗ᴿ⁻ᵈᵃᵗᵃ = Arpack.eigs(T; v0=vⁱᴿ⁻ᵈᵃᵗᵃ, nev=neigs)
170+
171+
## XXX: this is giving an error about trying to set the element of the wrong QN block for:
172+
## maxdim = 5
173+
## cutoff = 1e-12
174+
## max_vumps_iters = 10
175+
## outer_iters = 10
176+
## model_params = (J=1.0, h=0.8)
177+
##
178+
## v⃗ᴿᴬ = [itensor(v⃗ᴿ⁻ᵈᵃᵗᵃ[:, n], input_inds(T); tol=1e-4) for n in 1:length(λ⃗ᴿᴬ)]
179+
## @show flux.(v⃗ᴿᴬ)
180+
181+
@show λ⃗ᴿᴬ
182+
183+
# Full eigendecomposition
184+
185+
Tfull = prod(T)
186+
DV = eigen(Tfull, input_inds(T), output_inds(T))
187+
188+
@show norm(Tfull * DV.V - DV.Vt * DV.D)
189+
190+
d = diag(array(DV.D))
191+
192+
p = sortperm(d; by=abs, rev=true)
193+
@show p[1:neigs]
194+
@show d[p[1:neigs]]
195+
196+
println("Error if ED with Arpack")
197+
@show d[p[1:neigs]] - λ⃗ᴿᴬ

‎examples/vumps/vumps_ising.jl‎

Lines changed: 26 additions & 100 deletions
Original file line numberDiff line numberDiff line change
@@ -1,15 +1,20 @@
11
using ITensorInfiniteMPS
2-
using ITensorInfiniteMPS.ITensors
2+
using ITensors
33

44
##############################################################################
5-
# VUMPS parameters
5+
# VUMPS/TDVP parameters
66
#
77

8-
maxdim = 5 # Maximum bond dimension
9-
cutoff = 1e-12 # Singular value cutoff when increasing the bond dimension
10-
max_vumps_iters = 10 # Maximum number of iterations of the VUMPS algorithm at a fixed bond dimension
11-
outer_iters = 10 # Number of times to increase the bond dimension
12-
time_step = -Inf # -Inf corresponds to VUMPS
8+
maxdim = 20 # Maximum bond dimension
9+
cutoff = 1e-6 # Singular value cutoff when increasing the bond dimension
10+
max_vumps_iters = 10 # Maximum number of iterations of the VUMPS/TDVP algorithm at a fixed bond dimension
11+
tol = 1e-5 # Precision error tolerance for outer loop of VUMPS or TDVP
12+
outer_iters = 5 # Number of times to increase the bond dimension
13+
time_step = -Inf # -Inf corresponds to VUMPS, finite time_step corresponds to TDVP
14+
solver_tol = (x -> x / 100) # Tolerance for the local solver (eigsolve in VUMPS and exponentiate in TDVP)
15+
multisite_update_alg = "parallel" # Choose between ["sequential", "parallel"]. Only parallel works with TDVP.
16+
conserve_qns = false
17+
N = 2 # Number of sites in the unit cell (1-site unit cell is currently broken)
1318

1419
# Parameters of the transverse field Ising model
1520
model_params = (J=1.0, h=0.9)
@@ -18,15 +23,17 @@ model_params = (J=1.0, h=0.9)
1823
# CODE BELOW HERE DOES NOT NEED TO BE MODIFIED
1924
#
2025

21-
N = 2 # Number of sites in the unit cell
22-
2326
model = Model"ising"()
2427

25-
function space_shifted(::Model"ising", q̃sz)
26-
return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1]
28+
function space_shifted(::Model"ising", q̃sz; conserve_qns=true)
29+
if conserve_qns
30+
return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1]
31+
else
32+
return [QN() => 2]
33+
end
2734
end
2835

29-
space_ = fill(space_shifted(model, 0), N)
36+
space_ = fill(space_shifted(model, 0; conserve_qns=conserve_qns), N)
3037
s = infsiteinds("S=1/2", N; space=space_)
3138
initstate(n) = "↑"
3239
ψ = InfMPS(s, initstate)
@@ -37,9 +44,14 @@ H = InfiniteITensorSum(model, s; model_params...)
3744
# Check translational invariance
3845
@show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))
3946

40-
vumps_kwargs = (tol=1e-5, maxiter=max_vumps_iters)
47+
vumps_kwargs = (
48+
tol=tol,
49+
maxiter=max_vumps_iters,
50+
solver_tol=solver_tol,
51+
multisite_update_alg=multisite_update_alg,
52+
)
4153
subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
42-
ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...)
54+
#ψ = tdvp(H, ψ; time_step=time_step, vumps_kwargs...)
4355

4456
# Alternate steps of running VUMPS and increasing the bond dimension
4557
@time for _ in 1:outer_iters
@@ -97,89 +109,3 @@ Sz2_infinite = expect(ψ.AL[2] * ψ.C[2], "Sz")
97109

98110
@show Sz1_finite, Sz2_finite
99111
@show Sz1_infinite, Sz2_infinite
100-
101-
##############################################################################
102-
# Compute eigenspace of the transfer matrix
103-
#
104-
105-
using Arpack
106-
using KrylovKit
107-
using LinearAlgebra
108-
109-
T = TransferMatrix(ψ.AL)
110-
Tᵀ = transpose(T)
111-
vⁱᴿ = randomITensor(dag(input_inds(T)))
112-
vⁱᴸ = randomITensor(dag(input_inds(Tᵀ)))
113-
114-
neigs = 10
115-
tol = 1e-10
116-
λ⃗ᴿ, v⃗ᴿ, right_info = eigsolve(T, vⁱᴿ, neigs, :LM; tol=tol)
117-
λ⃗ᴸ, v⃗ᴸ, left_info = eigsolve(Tᵀ, vⁱᴸ, neigs, :LM; tol=tol)
118-
119-
@show norm(T(v⃗ᴿ[1]) - λ⃗ᴿ[1] * v⃗ᴿ[1])
120-
@show norm(Tᵀ(v⃗ᴸ[1]) - λ⃗ᴸ[1] * v⃗ᴸ[1])
121-
122-
@show λ⃗ᴿ
123-
@show λ⃗ᴸ
124-
@show flux.(v⃗ᴿ)
125-
126-
neigs = length(v⃗ᴿ)
127-
128-
# Normalize the vectors
129-
N⃗ = [(translatecell(v⃗ᴸ[n], 1) * v⃗ᴿ[n])[] for n in 1:neigs]
130-
131-
v⃗ᴿ ./= sqrt.(N⃗)
132-
v⃗ᴸ ./= sqrt.(N⃗)
133-
134-
# Form a second starting vector orthogonal to v⃗ᴿ[1]
135-
# This doesn't work. TODO: project out v⃗ᴿ[1], v⃗ᴸ[1] from T
136-
#λ⃗ᴿ², v⃗ᴿ², right_info_2 = eigsolve(T, vⁱᴿ², neigs, :LM; tol=tol)
137-
138-
# Projector onto the n-th eigenstate
139-
function proj(v⃗ᴸ, v⃗ᴿ, n)
140-
Lⁿ = v⃗ᴸ[n]
141-
Rⁿ = v⃗ᴿ[n]
142-
return ITensorMap(
143-
[translatecell(Lⁿ, 1), translatecell(Rⁿ, -1)]; input_inds=inds(Rⁿ), output_inds=inds(Lⁿ)
144-
)
145-
end
146-
147-
P⃗ = [proj(v⃗ᴸ, v⃗ᴿ, n) for n in 1:neigs]
148-
T⁻P = T - sum(P⃗)
149-
150-
#vⁱᴿ² = vⁱᴿ - (translatecell(v⃗ᴸ[1], 1) * vⁱᴿ)[] / norm(v⃗ᴿ[1]) * v⃗ᴿ[1]
151-
#@show norm(dag(vⁱᴿ²) * v⃗ᴿ[1])
152-
153-
λ⃗ᴾᴿ, v⃗ᴾᴿ, right_info = eigsolve(T⁻P, vⁱᴿ, neigs, :LM; tol=tol)
154-
@show λ⃗ᴾᴿ
155-
156-
vⁱᴿ⁻ᵈᵃᵗᵃ = vec(array(vⁱᴿ))
157-
λ⃗ᴿᴬ, v⃗ᴿ⁻ᵈᵃᵗᵃ = Arpack.eigs(T; v0=vⁱᴿ⁻ᵈᵃᵗᵃ, nev=neigs)
158-
159-
## XXX: this is giving an error about trying to set the element of the wrong QN block for:
160-
## maxdim = 5
161-
## cutoff = 1e-12
162-
## max_vumps_iters = 10
163-
## outer_iters = 10
164-
## model_params = (J=1.0, h=0.8)
165-
##
166-
## v⃗ᴿᴬ = [itensor(v⃗ᴿ⁻ᵈᵃᵗᵃ[:, n], input_inds(T); tol=1e-4) for n in 1:length(λ⃗ᴿᴬ)]
167-
## @show flux.(v⃗ᴿᴬ)
168-
169-
@show λ⃗ᴿᴬ
170-
171-
# Full eigendecomposition
172-
173-
Tfull = prod(T)
174-
DV = eigen(Tfull, input_inds(T), output_inds(T))
175-
176-
@show norm(Tfull * DV.V - DV.Vt * DV.D)
177-
178-
d = diag(array(DV.D))
179-
180-
p = sortperm(d; by=abs, rev=true)
181-
@show p[1:neigs]
182-
@show d[p[1:neigs]]
183-
184-
println("Error if ED with Arpack")
185-
@show d[p[1:neigs]] - λ⃗ᴿᴬ

‎src/vumps_localham.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -676,6 +676,7 @@ end
676676

677677
return function tdvp_solver(M, time_step, v₀, solver_tol)
678678
v, info = exponentiate(M, time_step, v₀; ishermitian=true, tol=solver_tol)
679+
v = v / norm(v)
679680
return nothing, v, info
680681
end
681682

0 commit comments

Comments
 (0)