Skip to content

Commit e89a521

Browse files
lkdvoslkdvos
committed
Add SvdCut
1 parent ad84949 commit e89a521

File tree

3 files changed

+115
-25
lines changed

3 files changed

+115
-25
lines changed

‎src/MPSKit.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -115,6 +115,7 @@ include("operators/abstractmpo.jl")
115115
include("operators/mpo.jl")
116116
include("operators/jordanmpotensor.jl")
117117
include("operators/mpohamiltonian.jl") # the mpohamiltonian objects
118+
include("operators/ortho.jl")
118119
include("operators/multilinempo.jl")
119120
include("operators/projection.jl")
120121
include("operators/timedependence.jl")

‎src/algorithms/changebonds/svdcut.jl‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -105,3 +105,22 @@ end
105105
function changebonds(ψ, H, alg::SvdCut, envs = environments(ψ, H))
106106
return changebonds(ψ, alg), envs
107107
end
108+
109+
function changebonds!(H::FiniteMPOHamiltonian, alg::SvdCut)
110+
# orthogonality center to the left
111+
for i in length(H):-1:2
112+
H = right_canonicalize!(H, i)
113+
end
114+
115+
# swipe right
116+
for i in 1:(length(H) - 1)
117+
@info i
118+
H = left_canonicalize!(H, i; alg = alg.alg_svd, alg.trscheme)
119+
end
120+
# swipe left -- TODO: do we really need this double sweep?
121+
for i in length(H):-1:2
122+
H = right_canonicalize!(H, i; alg = alg.alg_svd, alg.trscheme)
123+
end
124+
125+
return H
126+
end

‎src/operators/ortho.jl‎

Lines changed: 95 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -30,10 +30,17 @@ function TensorKit.leftorth(W::JordanMPOTensor; alg = QRpos())
3030
return W′, R
3131
end
3232

33-
function left_canonicalize!(H::MPOHamilonian, i::Int; alg = QRPos())
34-
@assert i != 1 "TBA"
33+
function left_canonicalize!(
34+
H::MPOHamiltonian, i::Int;
35+
alg = QRpos(), trscheme::TruncationScheme = notrunc()
36+
)
37+
if H isa FiniteMPOHamiltonian
38+
1 ≤ i < length(H) || throw(ArgumentError("Bounds error in canonicalize"))
39+
end
3540

3641
W = H[i]
42+
S = spacetype(W)
43+
d = sqrt(dim(physicalspace(W)))
3744

3845
# orthogonalize second column against first
3946
WI = removeunit(W[1, 1, 1, 1], 1)
@@ -42,57 +49,120 @@ function left_canonicalize!(H::MPOHamilonian, i::Int; alg = QRPos())
4249
@plansor C′[p; p' r] := -WI[p; p' l] * t[l; r]
4350
add!(C′, W.C)
4451

45-
# QR of second column
46-
CA = transpose(cat(insertleftunit(C′, 1), W.A; dims = 1), ((3, 1, 2), (4,)))
47-
Q, R = leftorth!(CA; alg)
48-
Q′ = transpose(Q, ((2, 3), (1, 4)))
49-
Q1 = SparseBlockTensorMap(Q′[2:end, 1, 1, 1])
50-
Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1)
51-
V = BlockTensorKit.oplus(oneunit(spacetype(W)), right_virtualspace(Q′), oneunit(spacetype(W)))
52+
# QR of second colum
53+
if i == 1
54+
tmp = transpose(C′, ((2, 1), (3,)))
55+
if trscheme == notrunc()
56+
Q, R = leftorth!(tmp; alg)
57+
else
58+
@assert alg == SVD() || alg == SDD()
59+
# TODO: truncation might not be using normalized tensors?
60+
Q, Σ, Vᴴ = tsvd!(tmp; alg, trunc = trscheme)
61+
R = Σ * Vᴴ
62+
end
63+
scale!(Q, d)
64+
scale!(R, inv(d))
65+
66+
V = BlockTensorKit.oplus(oneunit(S), space(R, 1), oneunit(S))
67+
68+
Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ space(R, 1))
69+
Q2 = transpose(Q, ((2,), (1, 3)))
70+
else
71+
tmp = transpose(cat(insertleftunit(C′, 1), W.A; dims = 1), ((3, 1, 2), (4,)))
72+
if trscheme == notrunc()
73+
Q, R = leftorth!(tmp; alg)
74+
else
75+
@assert alg == SVD() || alg == SDD()
76+
# TODO: truncation might not be using normalized tensors?
77+
Q, Σ, Vᴴ = tsvd!(tmp; alg, trunc = trscheme)
78+
R = Σ * Vᴴ
79+
end
80+
scale!(Q, d)
81+
scale!(R, inv(d))
82+
Q′ = transpose(Q, ((2, 3), (1, 4)))
83+
Q1 = Q′[2:end, 1, 1, 1]
84+
Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1)
85+
V = BlockTensorKit.oplus(oneunit(S), right_virtualspace(Q′), oneunit(S))
86+
end
5287
H[i] = JordanMPOTensor(codomain(W) ← physicalspace(W) ⊗ V, Q1, W.B, Q2, W.D)
5388

5489
# absorb into next site
5590
W′ = H[i + 1]
5691
@plansor A′[l p; p' r] := R[l; r'] * W′.A[r' p; p' r]
5792
@plansor B′[l p; p'] := R[l; r] * W′.B[r p; p']
5893
@plansor C′[l p; p' r] := t[l; r'] * W′.A[r' p; p' r]
59-
C′ = add!(removeunit(C, 1), W′.C)
94+
C′ = add!(removeunit(C′, 1), W′.C)
6095
@plansor D′[l p; p'] := t[l; r] * W′.B[r p; p']
61-
D′ = add!(removeunit(D, 1), W′.D)
96+
D′ = add!(removeunit(D′, 1), W′.D)
6297

6398
H[i + 1] = JordanMPOTensor(V ⊗ physicalspace(W′) ← domain(W′), A′, B′, C′, D′)
6499
return H
65100
end
66101

67-
function right_canonicalize!(H::MPOHamilonian, i::Int; alg = RQpos())
68-
@assert i != length(H) "TBA"
102+
function right_canonicalize!(
103+
H::MPOHamiltonian, i::Int;
104+
alg = RQpos(), trscheme::TruncationScheme = notrunc()
105+
)
106+
if H isa FiniteMPOHamiltonian
107+
1 < i ≤ length(H) || throw(ArgumentError("Bounds error in canonicalize"))
108+
end
69109

70110
W = H[i]
111+
S = spacetype(W)
112+
d = sqrt(dim(physicalspace(W)))
71113

72114
# orthogonalize second row against last
73115
WI = removeunit(W[end, 1, 1, end], 4)
74-
@tensor t[l; r] := conj(WI[r p; p']) * W.B[l p; p']
116+
@plansor t[l; r] := conj(WI[r p; p']) * W.B[l p; p']
75117
@plansor B′[l p; p'] := -WI[r p; p'] * t[l; r]
76118
add!(B′, W.B)
77119

78120
# LQ of second row
79-
AB = transpose(cat(insertleftunit(B′, 4), W.A; dims = 4), ((1,), (3, 4, 2)))
80-
R, Q = rightorth!(AB; alg)
81-
Q′ = transpose(Q, ((1, 4), (2, 3)))
82-
Q1 = SparseBlockTensorMap(Q′[1, 1, 1, 2:end])
83-
Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 4)
84-
V = BlockTensorKit.oplus(oneunit(spacetype(W)), left_virtualspace(Q′), oneunit(spacetype(W)))
121+
if i == length(H)
122+
tmp = transpose(B′, ((1,), (3, 2)))
123+
if trscheme == notrunc()
124+
R, Q = rightorth!(tmp; alg)
125+
else
126+
@assert alg == SVD() || alg == SDD()
127+
# TODO: truncation might not be using normalized tensors?
128+
U, Σ, Q = tsvd!(tmp; alg, trunc = trscheme)
129+
R = U * Σ
130+
end
131+
scale!(Q, d)
132+
scale!(R, inv(d))
133+
134+
R, Q = rightorth!(transpose(B′, ((1,), (3, 2))); alg)
135+
V = BlockTensorKit.oplus(oneunit(S), space(Q, 1), oneunit(S))
136+
Q1 = typeof(W.A)(undef, space(Q, 1) ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}())
137+
Q2 = transpose(Q, ((1, 3), (2,)))
138+
else
139+
tmp = transpose(cat(insertleftunit(B′, 4), W.A; dims = 4), ((1,), (3, 4, 2)))
140+
if trscheme == notrunc()
141+
R, Q = rightorth!(tmp; alg)
142+
else
143+
@assert alg == SVD() || alg == SDD()
144+
# TODO: truncation might not be using normalized tensors?
145+
U, Σ, Q = tsvd!(tmp; alg, trunc = trscheme)
146+
R = U * Σ
147+
end
148+
scale!(Q, d)
149+
scale!(R, inv(d))
150+
Q′ = transpose(Q, ((1, 4), (2, 3)))
151+
Q1 = SparseBlockTensorMap(Q′[1, 1, 1, 2:end])
152+
Q2 = removeunit(SparseBlockTensorMap(Q′[1, 1, 1, 1:1]), 4)
153+
V = BlockTensorKit.oplus(oneunit(S), left_virtualspace(Q′), oneunit(S))
154+
end
85155
H[i] = JordanMPOTensor(V ⊗ physicalspace(W) ← domain(W), Q1, Q2, W.C, W.D)
86156

87157
# absorb into previous site
88158
W′ = H[i - 1]
89-
@plansor A′[l p; p' r] := W′.A[l p; p' r] * R[r; r]
90-
@plansor B′[l p; p' r'] := W′.A[l p; p' r'] * t[r'; r] R[l; r] * W′.B[r p; p']
159+
@plansor A′[l p; p' r] := W′.A[l p; p' r'] * R[r'; r]
160+
@plansor B′[l p; p' r] := W′.A[l p; p' r'] * t[r'; r]
91161
B′ = add!(removeunit(B′, 4), W′.B)
92162
@plansor C′[p; p' r] := W′.C[p; p' r'] * R[r'; r]
93-
@plansor D′[p; p' r] := W′.C[p; p' l] *t[l; r]
94-
D′ = add!(removeunit(D, 3), W′.D)
163+
@plansor D′[p; p' r] := W′.C[p; p' r'] * t[r'; r]
164+
D′ = add!(removeunit(D′, 3), W′.D)
95165

96-
H[i - 1] = JordanMPOTensor(V ⊗ physicalspace(W′) ← domain(W′), A′, B′, C′, D′)
166+
H[i - 1] = JordanMPOTensor(codomain(W′) ← physicalspace(W′) ⊗ V, A′, B′, C′, D′)
97167
return H
98168
end

0 commit comments

Comments
 (0)