From 4046fadb8f1a4d827a544f7f63f1e1381df2e782 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 3 Aug 2025 19:43:52 -0400 Subject: [PATCH 01/17] Start setting up ortho for mpohamiltonian --- src/operators/ortho.jl | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) create mode 100644 src/operators/ortho.jl diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl new file mode 100644 index 000000000..8fed7be18 --- /dev/null +++ b/src/operators/ortho.jl @@ -0,0 +1,29 @@ +function TensorKit.leftorth!(W::JordanMPOTensor; alg = QRpos()) + # orthogonalize second column against first + WI = removeunit(W[1, 1, 1, 1], 1) + @tensor t[l; r] := conj(WI[p; p' l]) * W.C[p; p' r] + I = sectortype(W) + S = spacetype(W) + + @tensor t[l; r] := conj(removeunit(W[1, 1, 1, 1], 1)[p; p' l]) * W.C[p; p' r] + @tensor C′[p; p' r] := W.C[p; p' r] - WI[p; p' l] * t[l; r] + + # QR of second column + CA = cat(insertleftunit(C′, 1), W.A; dims = 1) + Q, r = leftorth(CA, ((3, 1, 2), (4,)); alg) + Q′ = transpose(Q, ((2, 3), (1, 4))) + V = codomain(W) ← physicalspace(W) ⊗ BlockTensorKit.oplus(oneunit(S), right_virtualspace(Q′), oneunit(S)) + Q1 = SparseBlockTensorMap(Q′[2:end, 1, 1, 1]) + Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1) + + # Assemble output + W′ = JordanMPOTensor(V, Q1, W.B, Q2, W.D) + + R = similar(W, right_virtualspace(W′) ← right_virtualspace(W′)) + R[1, 1] = id!(R[1, 1]) + R[1, 2] = t + R[2, 2] = r + R[3, 3] = id!(R[3, 3]) + + return W′, R +end From b5b32c866b0b80eae4693a099a5d6dfacda0253b Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 3 Aug 2025 19:44:05 -0400 Subject: [PATCH 02/17] Add left_canonicalize! --- src/operators/ortho.jl | 46 +++++++++++++++++++++++++++++++++++++----- 1 file changed, 41 insertions(+), 5 deletions(-) diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl index 8fed7be18..a11b5082b 100644 --- a/src/operators/ortho.jl +++ b/src/operators/ortho.jl @@ -1,16 +1,18 @@ -function TensorKit.leftorth!(W::JordanMPOTensor; alg = QRpos()) +function TensorKit.leftorth(W::JordanMPOTensor; alg = QRpos()) # orthogonalize second column against first WI = removeunit(W[1, 1, 1, 1], 1) @tensor t[l; r] := conj(WI[p; p' l]) * W.C[p; p' r] I = sectortype(W) S = spacetype(W) - @tensor t[l; r] := conj(removeunit(W[1, 1, 1, 1], 1)[p; p' l]) * W.C[p; p' r] - @tensor C′[p; p' r] := W.C[p; p' r] - WI[p; p' l] * t[l; r] + @plansor t[l; r] := conj(removeunit(W[1, 1, 1, 1], 1)[p; p' l]) * W.C[p; p' r] + # @plansor C′[p; p' r] := W.C[p; p' r] - WI[p; p' l] * t[l; r] + @plansor C′[p; p' r] := -WI[p; p' l] * t[l; r] + add!(C′, W.C) # QR of second column - CA = cat(insertleftunit(C′, 1), W.A; dims = 1) - Q, r = leftorth(CA, ((3, 1, 2), (4,)); alg) + CA = transpose(cat(insertleftunit(C′, 1), W.A; dims = 1), ((3, 1, 2), (4,))) + Q, r = leftorth!(CA; alg) Q′ = transpose(Q, ((2, 3), (1, 4))) V = codomain(W) ← physicalspace(W) ⊗ BlockTensorKit.oplus(oneunit(S), right_virtualspace(Q′), oneunit(S)) Q1 = SparseBlockTensorMap(Q′[2:end, 1, 1, 1]) @@ -27,3 +29,37 @@ function TensorKit.leftorth!(W::JordanMPOTensor; alg = QRpos()) return W′, R end + +function left_canonicalize!(H::MPOHamilonian, i::Int; alg = QRPos()) + @assert i != 1 "TBA" + + W = H[i] + + # orthogonalize second column against first + WI = removeunit(W[1, 1, 1, 1], 1) + @tensor t[l; r] := conj(WI[p; p' l]) * W.C[p; p' r] + # @plansor C′[p; p' r] := W.C[p; p' r] - WI[p; p' l] * t[l; r] + @plansor C′[p; p' r] := -WI[p; p' l] * t[l; r] + add!(C′, W.C) + + # QR of second column + CA = transpose(cat(insertleftunit(C′, 1), W.A; dims = 1), ((3, 1, 2), (4,))) + Q, R = leftorth!(CA; alg) + Q′ = transpose(Q, ((2, 3), (1, 4))) + Q1 = SparseBlockTensorMap(Q′[2:end, 1, 1, 1]) + Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1) + V = BlockTensorKit.oplus(oneunit(spacetype(W)), right_virtualspace(Q′), oneunit(spacetype(W))) + H[i] = JordanMPOTensor(codomain(W) ← physicalspace(W) ⊗ V, Q1, W.B, Q2, W.D) + + # absorb into next site + W′ = H[i + 1] + @plansor A′[l p; p' r] := R[l; r'] * W′.A[r' p; p' r] + @plansor B′[l p; p'] := R[l; r] * W′.B[r p; p'] + @plansor C′[l p; p' r] := t[l; r'] * W′.A[r' p; p' r] + C′ = add!(removeunit(C, 1), W′.C) + @plansor D′[l p; p'] := t[l; r] * W′.B[r p; p'] + D′ = add!(removeunit(D, 1), W′.D) + + H[i + 1] = JordanMPOTensor(V ⊗ physicalspace(W′) ← domain(W′), A′, B′, C′, D′) + return H +end From 0df7c4df96485ac7c15627c579208456cfba9959 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 3 Aug 2025 20:21:32 -0400 Subject: [PATCH 03/17] Add right_canonicalize! --- src/operators/ortho.jl | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl index a11b5082b..0eccededd 100644 --- a/src/operators/ortho.jl +++ b/src/operators/ortho.jl @@ -63,3 +63,36 @@ function left_canonicalize!(H::MPOHamilonian, i::Int; alg = QRPos()) H[i + 1] = JordanMPOTensor(V ⊗ physicalspace(W′) ← domain(W′), A′, B′, C′, D′) return H end + +function right_canonicalize!(H::MPOHamilonian, i::Int; alg = RQpos()) + @assert i != length(H) "TBA" + + W = H[i] + + # orthogonalize second row against last + WI = removeunit(W[end, 1, 1, end], 4) + @tensor t[l; r] := conj(WI[r p; p']) * W.B[l p; p'] + @plansor B′[l p; p'] := -WI[r p; p'] * t[l; r] + add!(B′, W.B) + + # LQ of second row + AB = transpose(cat(insertleftunit(B′, 4), W.A; dims = 4), ((1,), (3, 4, 2))) + R, Q = rightorth!(AB; alg) + Q′ = transpose(Q, ((1, 4), (2, 3))) + Q1 = SparseBlockTensorMap(Q′[1, 1, 1, 2:end]) + Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 4) + V = BlockTensorKit.oplus(oneunit(spacetype(W)), left_virtualspace(Q′), oneunit(spacetype(W))) + H[i] = JordanMPOTensor(V ⊗ physicalspace(W) ← domain(W), Q1, Q2, W.C, W.D) + + # absorb into previous site + W′ = H[i - 1] + @plansor A′[l p; p' r] := W′.A[l p; p' r] * R[r; r] + @plansor B′[l p; p' r'] := W′.A[l p; p' r'] * t[r'; r] R[l; r] * W′.B[r p; p'] + B′ = add!(removeunit(B′, 4), W′.B) + @plansor C′[p; p' r] := W′.C[p; p' r'] * R[r'; r] + @plansor D′[p; p' r] := W′.C[p; p' l] *t[l; r] + D′ = add!(removeunit(D, 3), W′.D) + + H[i - 1] = JordanMPOTensor(V ⊗ physicalspace(W′) ← domain(W′), A′, B′, C′, D′) + return H +end From 622f4cb6ed7484cd7ccd919cb543a8cb6619100b Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 4 Aug 2025 11:49:22 -0400 Subject: [PATCH 04/17] Add JordanMPOTensor constructor that promotes types --- src/operators/jordanmpotensor.jl | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl index e46356f35..939b47d12 100644 --- a/src/operators/jordanmpotensor.jl +++ b/src/operators/jordanmpotensor.jl @@ -67,10 +67,8 @@ end function JordanMPOTensor( V::TensorMapSumSpace{S, 2, 2}, - A::SparseBlockTensorMap{TA, E, S, 2, 2}, - B::SparseBlockTensorMap{TB, E, S, 2, 1}, - C::SparseBlockTensorMap{TC, E, S, 1, 2}, - D::SparseBlockTensorMap{TD, E, S, 1, 1} + A::SparseBlockTensorMap{TA, E, S, 2, 2}, B::SparseBlockTensorMap{TB, E, S, 2, 1}, + C::SparseBlockTensorMap{TC, E, S, 1, 2}, D::SparseBlockTensorMap{TD, E, S, 1, 1} ) where {E, S, TA, TB, TC, TD} allVs = eachspace(V) VA = space(allVs[2:(end - 1), 1, 1, 2:(end - 1)]) @@ -87,6 +85,18 @@ function JordanMPOTensor( return JordanMPOTensor{E, S, TA, TB, TC, TD}(V, A, B, C, D) end +function JordanMPOTensor( + V::TensorMapSumSpace{S, 2, 2}, + A::AbstractTensorMap{E, S, 2, 2}, B::AbstractTensorMap{E, S, 2, 1}, + C::AbstractTensorMap{E, S, 1, 2}, D::AbstractTensorMap{E, S, 1, 1} + ) where {E, S} + return JordanMPOTensor(V, + A isa SparseBlockTensorMap ? A : SparseBlockTensorMap(A), + B isa SparseBlockTensorMap ? B : SparseBlockTensorMap(B), + C isa SparseBlockTensorMap ? C : SparseBlockTensorMap(C), + D isa SparseBlockTensorMap ? D : SparseBlockTensorMap(D) + ) +end function JordanMPOTensor(W::SparseBlockTensorMap{TT, E, S, 2, 2}) where {TT, E, S} @assert W[1, 1, 1, 1] isa BraidingTensor && W[end, 1, 1, end] isa BraidingTensor From ad84949f535436bdbe3aa8f46fac28dc84228e9b Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 4 Aug 2025 11:49:32 -0400 Subject: [PATCH 05/17] Add error messages --- src/operators/jordanmpotensor.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl index 939b47d12..ba011ff3e 100644 --- a/src/operators/jordanmpotensor.jl +++ b/src/operators/jordanmpotensor.jl @@ -72,16 +72,16 @@ function JordanMPOTensor( ) where {E, S, TA, TB, TC, TD} allVs = eachspace(V) VA = space(allVs[2:(end - 1), 1, 1, 2:(end - 1)]) - VA == space(A) || throw(SpaceMismatch("A-block has incompatible spaces")) + VA == space(A) || throw(SpaceMismatch(lazy"A-block has incompatible spaces:\n$VA\n$(space(A))")) VB = removeunit(space(allVs[2:(end - 1), 1, 1, end]), 4) - VB == space(B) || throw(SpaceMismatch("B-block has incompatible spaces")) + VB == space(B) || throw(SpaceMismatch(lazy"B-block has incompatible spaces:\n$VB\n$(space(B))")) VC = removeunit(space(allVs[1, 1, 1, 2:(end - 1)]), 1) - VC == space(C) || throw(SpaceMismatch("C-block has incompatible spaces")) + VC == space(C) || throw(SpaceMismatch(lazy"C-block has incompatible spaces:\n$VC\n$(space(C))")) VD = removeunit(removeunit(space(allVs[1, 1, 1, end:end]), 4), 1) - VD == space(D) || throw(SpaceMismatch("D-block has incompatible spaces")) + VD == space(D) || throw(SpaceMismatch(lazy"D-block has incompatible spaces:\n$VD\n$(space(D))")) return JordanMPOTensor{E, S, TA, TB, TC, TD}(V, A, B, C, D) end From fa7a45056fb77affe729c875a4958c9b785ad014 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 4 Aug 2025 11:50:03 -0400 Subject: [PATCH 06/17] Add SvdCut --- src/MPSKit.jl | 1 + src/algorithms/changebonds/svdcut.jl | 18 ++++ src/operators/ortho.jl | 120 +++++++++++++++++++++------ 3 files changed, 114 insertions(+), 25 deletions(-) diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 4bab6bf82..51ee82e32 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -115,6 +115,7 @@ include("operators/abstractmpo.jl") include("operators/mpo.jl") include("operators/jordanmpotensor.jl") include("operators/mpohamiltonian.jl") # the mpohamiltonian objects +include("operators/ortho.jl") include("operators/multilinempo.jl") include("operators/projection.jl") include("operators/timedependence.jl") diff --git a/src/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl index cd60bf1b8..a73133344 100644 --- a/src/algorithms/changebonds/svdcut.jl +++ b/src/algorithms/changebonds/svdcut.jl @@ -105,3 +105,21 @@ end function changebonds(ψ, H, alg::SvdCut, envs = environments(ψ, H)) return changebonds(ψ, alg), envs end + +function changebonds!(H::FiniteMPOHamiltonian, alg::SvdCut) + # orthogonality center to the left + for i in length(H):-1:2 + H = right_canonicalize!(H, i) + end + + # swipe right + for i in 1:(length(H) - 1) + H = left_canonicalize!(H, i; alg = alg.alg_svd, alg.trscheme) + end + # swipe left -- TODO: do we really need this double sweep? + for i in length(H):-1:2 + H = right_canonicalize!(H, i; alg = alg.alg_svd, alg.trscheme) + end + + return H +end diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl index 0eccededd..d34b0a6c5 100644 --- a/src/operators/ortho.jl +++ b/src/operators/ortho.jl @@ -30,10 +30,17 @@ function TensorKit.leftorth(W::JordanMPOTensor; alg = QRpos()) return W′, R end -function left_canonicalize!(H::MPOHamilonian, i::Int; alg = QRPos()) - @assert i != 1 "TBA" +function left_canonicalize!( + H::MPOHamiltonian, i::Int; + alg = QRpos(), trscheme::TruncationScheme = notrunc() + ) + if H isa FiniteMPOHamiltonian + 1 ≤ i < length(H) || throw(ArgumentError("Bounds error in canonicalize")) + end W = H[i] + S = spacetype(W) + d = sqrt(dim(physicalspace(W))) # orthogonalize second column against first WI = removeunit(W[1, 1, 1, 1], 1) @@ -42,13 +49,41 @@ function left_canonicalize!(H::MPOHamilonian, i::Int; alg = QRPos()) @plansor C′[p; p' r] := -WI[p; p' l] * t[l; r] add!(C′, W.C) - # QR of second column - CA = transpose(cat(insertleftunit(C′, 1), W.A; dims = 1), ((3, 1, 2), (4,))) - Q, R = leftorth!(CA; alg) - Q′ = transpose(Q, ((2, 3), (1, 4))) - Q1 = SparseBlockTensorMap(Q′[2:end, 1, 1, 1]) - Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1) - V = BlockTensorKit.oplus(oneunit(spacetype(W)), right_virtualspace(Q′), oneunit(spacetype(W))) + # QR of second colum + if i == 1 + tmp = transpose(C′, ((2, 1), (3,))) + if trscheme == notrunc() + Q, R = leftorth!(tmp; alg) + else + @assert alg == SVD() || alg == SDD() + # TODO: truncation might not be using normalized tensors? + Q, Σ, Vᴴ = tsvd!(tmp; alg, trunc = trscheme) + R = Σ * Vᴴ + end + scale!(Q, d) + scale!(R, inv(d)) + + V = BlockTensorKit.oplus(oneunit(S), space(R, 1), oneunit(S)) + + Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ space(R, 1)) + Q2 = transpose(Q, ((2,), (1, 3))) + else + tmp = transpose(cat(insertleftunit(C′, 1), W.A; dims = 1), ((3, 1, 2), (4,))) + if trscheme == notrunc() + Q, R = leftorth!(tmp; alg) + else + @assert alg == SVD() || alg == SDD() + # TODO: truncation might not be using normalized tensors? + Q, Σ, Vᴴ = tsvd!(tmp; alg, trunc = trscheme) + R = Σ * Vᴴ + end + scale!(Q, d) + scale!(R, inv(d)) + Q′ = transpose(Q, ((2, 3), (1, 4))) + Q1 = Q′[2:end, 1, 1, 1] + Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1) + V = BlockTensorKit.oplus(oneunit(S), right_virtualspace(Q′), oneunit(S)) + end H[i] = JordanMPOTensor(codomain(W) ← physicalspace(W) ⊗ V, Q1, W.B, Q2, W.D) # absorb into next site @@ -56,43 +91,78 @@ function left_canonicalize!(H::MPOHamilonian, i::Int; alg = QRPos()) @plansor A′[l p; p' r] := R[l; r'] * W′.A[r' p; p' r] @plansor B′[l p; p'] := R[l; r] * W′.B[r p; p'] @plansor C′[l p; p' r] := t[l; r'] * W′.A[r' p; p' r] - C′ = add!(removeunit(C, 1), W′.C) + C′ = add!(removeunit(C′, 1), W′.C) @plansor D′[l p; p'] := t[l; r] * W′.B[r p; p'] - D′ = add!(removeunit(D, 1), W′.D) + D′ = add!(removeunit(D′, 1), W′.D) H[i + 1] = JordanMPOTensor(V ⊗ physicalspace(W′) ← domain(W′), A′, B′, C′, D′) return H end -function right_canonicalize!(H::MPOHamilonian, i::Int; alg = RQpos()) - @assert i != length(H) "TBA" +function right_canonicalize!( + H::MPOHamiltonian, i::Int; + alg = RQpos(), trscheme::TruncationScheme = notrunc() + ) + if H isa FiniteMPOHamiltonian + 1 < i ≤ length(H) || throw(ArgumentError("Bounds error in canonicalize")) + end W = H[i] + S = spacetype(W) + d = sqrt(dim(physicalspace(W))) # orthogonalize second row against last WI = removeunit(W[end, 1, 1, end], 4) - @tensor t[l; r] := conj(WI[r p; p']) * W.B[l p; p'] + @plansor t[l; r] := conj(WI[r p; p']) * W.B[l p; p'] @plansor B′[l p; p'] := -WI[r p; p'] * t[l; r] add!(B′, W.B) # LQ of second row - AB = transpose(cat(insertleftunit(B′, 4), W.A; dims = 4), ((1,), (3, 4, 2))) - R, Q = rightorth!(AB; alg) - Q′ = transpose(Q, ((1, 4), (2, 3))) - Q1 = SparseBlockTensorMap(Q′[1, 1, 1, 2:end]) - Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 4) - V = BlockTensorKit.oplus(oneunit(spacetype(W)), left_virtualspace(Q′), oneunit(spacetype(W))) + if i == length(H) + tmp = transpose(B′, ((1,), (3, 2))) + if trscheme == notrunc() + R, Q = rightorth!(tmp; alg) + else + @assert alg == SVD() || alg == SDD() + # TODO: truncation might not be using normalized tensors? + U, Σ, Q = tsvd!(tmp; alg, trunc = trscheme) + R = U * Σ + end + scale!(Q, d) + scale!(R, inv(d)) + + R, Q = rightorth!(transpose(B′, ((1,), (3, 2))); alg) + V = BlockTensorKit.oplus(oneunit(S), space(Q, 1), oneunit(S)) + Q1 = typeof(W.A)(undef, space(Q, 1) ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}()) + Q2 = transpose(Q, ((1, 3), (2,))) + else + tmp = transpose(cat(insertleftunit(B′, 4), W.A; dims = 4), ((1,), (3, 4, 2))) + if trscheme == notrunc() + R, Q = rightorth!(tmp; alg) + else + @assert alg == SVD() || alg == SDD() + # TODO: truncation might not be using normalized tensors? + U, Σ, Q = tsvd!(tmp; alg, trunc = trscheme) + R = U * Σ + end + scale!(Q, d) + scale!(R, inv(d)) + Q′ = transpose(Q, ((1, 4), (2, 3))) + Q1 = SparseBlockTensorMap(Q′[1, 1, 1, 2:end]) + Q2 = removeunit(SparseBlockTensorMap(Q′[1, 1, 1, 1:1]), 4) + V = BlockTensorKit.oplus(oneunit(S), left_virtualspace(Q′), oneunit(S)) + end H[i] = JordanMPOTensor(V ⊗ physicalspace(W) ← domain(W), Q1, Q2, W.C, W.D) # absorb into previous site W′ = H[i - 1] - @plansor A′[l p; p' r] := W′.A[l p; p' r] * R[r; r] - @plansor B′[l p; p' r'] := W′.A[l p; p' r'] * t[r'; r] R[l; r] * W′.B[r p; p'] + @plansor A′[l p; p' r] := W′.A[l p; p' r'] * R[r'; r] + @plansor B′[l p; p' r] := W′.A[l p; p' r'] * t[r'; r] B′ = add!(removeunit(B′, 4), W′.B) @plansor C′[p; p' r] := W′.C[p; p' r'] * R[r'; r] - @plansor D′[p; p' r] := W′.C[p; p' l] *t[l; r] - D′ = add!(removeunit(D, 3), W′.D) + @plansor D′[p; p' r] := W′.C[p; p' r'] * t[r'; r] + D′ = add!(removeunit(D′, 3), W′.D) - H[i - 1] = JordanMPOTensor(V ⊗ physicalspace(W′) ← domain(W′), A′, B′, C′, D′) + H[i - 1] = JordanMPOTensor(codomain(W′) ← physicalspace(W′) ⊗ V, A′, B′, C′, D′) return H end From e89a521c0f9ae4cd952bcdf12f14740b9701d3b7 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 4 Aug 2025 11:50:03 -0400 Subject: [PATCH 07/17] Add SvdCut --- src/MPSKit.jl | 1 + src/algorithms/changebonds/svdcut.jl | 19 +++++ src/operators/ortho.jl | 120 +++++++++++++++++++++------ 3 files changed, 115 insertions(+), 25 deletions(-) diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 4bab6bf82..51ee82e32 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -115,6 +115,7 @@ include("operators/abstractmpo.jl") include("operators/mpo.jl") include("operators/jordanmpotensor.jl") include("operators/mpohamiltonian.jl") # the mpohamiltonian objects +include("operators/ortho.jl") include("operators/multilinempo.jl") include("operators/projection.jl") include("operators/timedependence.jl") diff --git a/src/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl index cd60bf1b8..c87c50d19 100644 --- a/src/algorithms/changebonds/svdcut.jl +++ b/src/algorithms/changebonds/svdcut.jl @@ -105,3 +105,22 @@ end function changebonds(ψ, H, alg::SvdCut, envs = environments(ψ, H)) return changebonds(ψ, alg), envs end + +function changebonds!(H::FiniteMPOHamiltonian, alg::SvdCut) + # orthogonality center to the left + for i in length(H):-1:2 + H = right_canonicalize!(H, i) + end + + # swipe right + for i in 1:(length(H) - 1) + @info i + H = left_canonicalize!(H, i; alg = alg.alg_svd, alg.trscheme) + end + # swipe left -- TODO: do we really need this double sweep? + for i in length(H):-1:2 + H = right_canonicalize!(H, i; alg = alg.alg_svd, alg.trscheme) + end + + return H +end diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl index 0eccededd..d34b0a6c5 100644 --- a/src/operators/ortho.jl +++ b/src/operators/ortho.jl @@ -30,10 +30,17 @@ function TensorKit.leftorth(W::JordanMPOTensor; alg = QRpos()) return W′, R end -function left_canonicalize!(H::MPOHamilonian, i::Int; alg = QRPos()) - @assert i != 1 "TBA" +function left_canonicalize!( + H::MPOHamiltonian, i::Int; + alg = QRpos(), trscheme::TruncationScheme = notrunc() + ) + if H isa FiniteMPOHamiltonian + 1 ≤ i < length(H) || throw(ArgumentError("Bounds error in canonicalize")) + end W = H[i] + S = spacetype(W) + d = sqrt(dim(physicalspace(W))) # orthogonalize second column against first WI = removeunit(W[1, 1, 1, 1], 1) @@ -42,13 +49,41 @@ function left_canonicalize!(H::MPOHamilonian, i::Int; alg = QRPos()) @plansor C′[p; p' r] := -WI[p; p' l] * t[l; r] add!(C′, W.C) - # QR of second column - CA = transpose(cat(insertleftunit(C′, 1), W.A; dims = 1), ((3, 1, 2), (4,))) - Q, R = leftorth!(CA; alg) - Q′ = transpose(Q, ((2, 3), (1, 4))) - Q1 = SparseBlockTensorMap(Q′[2:end, 1, 1, 1]) - Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1) - V = BlockTensorKit.oplus(oneunit(spacetype(W)), right_virtualspace(Q′), oneunit(spacetype(W))) + # QR of second colum + if i == 1 + tmp = transpose(C′, ((2, 1), (3,))) + if trscheme == notrunc() + Q, R = leftorth!(tmp; alg) + else + @assert alg == SVD() || alg == SDD() + # TODO: truncation might not be using normalized tensors? + Q, Σ, Vᴴ = tsvd!(tmp; alg, trunc = trscheme) + R = Σ * Vᴴ + end + scale!(Q, d) + scale!(R, inv(d)) + + V = BlockTensorKit.oplus(oneunit(S), space(R, 1), oneunit(S)) + + Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ space(R, 1)) + Q2 = transpose(Q, ((2,), (1, 3))) + else + tmp = transpose(cat(insertleftunit(C′, 1), W.A; dims = 1), ((3, 1, 2), (4,))) + if trscheme == notrunc() + Q, R = leftorth!(tmp; alg) + else + @assert alg == SVD() || alg == SDD() + # TODO: truncation might not be using normalized tensors? + Q, Σ, Vᴴ = tsvd!(tmp; alg, trunc = trscheme) + R = Σ * Vᴴ + end + scale!(Q, d) + scale!(R, inv(d)) + Q′ = transpose(Q, ((2, 3), (1, 4))) + Q1 = Q′[2:end, 1, 1, 1] + Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1) + V = BlockTensorKit.oplus(oneunit(S), right_virtualspace(Q′), oneunit(S)) + end H[i] = JordanMPOTensor(codomain(W) ← physicalspace(W) ⊗ V, Q1, W.B, Q2, W.D) # absorb into next site @@ -56,43 +91,78 @@ function left_canonicalize!(H::MPOHamilonian, i::Int; alg = QRPos()) @plansor A′[l p; p' r] := R[l; r'] * W′.A[r' p; p' r] @plansor B′[l p; p'] := R[l; r] * W′.B[r p; p'] @plansor C′[l p; p' r] := t[l; r'] * W′.A[r' p; p' r] - C′ = add!(removeunit(C, 1), W′.C) + C′ = add!(removeunit(C′, 1), W′.C) @plansor D′[l p; p'] := t[l; r] * W′.B[r p; p'] - D′ = add!(removeunit(D, 1), W′.D) + D′ = add!(removeunit(D′, 1), W′.D) H[i + 1] = JordanMPOTensor(V ⊗ physicalspace(W′) ← domain(W′), A′, B′, C′, D′) return H end -function right_canonicalize!(H::MPOHamilonian, i::Int; alg = RQpos()) - @assert i != length(H) "TBA" +function right_canonicalize!( + H::MPOHamiltonian, i::Int; + alg = RQpos(), trscheme::TruncationScheme = notrunc() + ) + if H isa FiniteMPOHamiltonian + 1 < i ≤ length(H) || throw(ArgumentError("Bounds error in canonicalize")) + end W = H[i] + S = spacetype(W) + d = sqrt(dim(physicalspace(W))) # orthogonalize second row against last WI = removeunit(W[end, 1, 1, end], 4) - @tensor t[l; r] := conj(WI[r p; p']) * W.B[l p; p'] + @plansor t[l; r] := conj(WI[r p; p']) * W.B[l p; p'] @plansor B′[l p; p'] := -WI[r p; p'] * t[l; r] add!(B′, W.B) # LQ of second row - AB = transpose(cat(insertleftunit(B′, 4), W.A; dims = 4), ((1,), (3, 4, 2))) - R, Q = rightorth!(AB; alg) - Q′ = transpose(Q, ((1, 4), (2, 3))) - Q1 = SparseBlockTensorMap(Q′[1, 1, 1, 2:end]) - Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 4) - V = BlockTensorKit.oplus(oneunit(spacetype(W)), left_virtualspace(Q′), oneunit(spacetype(W))) + if i == length(H) + tmp = transpose(B′, ((1,), (3, 2))) + if trscheme == notrunc() + R, Q = rightorth!(tmp; alg) + else + @assert alg == SVD() || alg == SDD() + # TODO: truncation might not be using normalized tensors? + U, Σ, Q = tsvd!(tmp; alg, trunc = trscheme) + R = U * Σ + end + scale!(Q, d) + scale!(R, inv(d)) + + R, Q = rightorth!(transpose(B′, ((1,), (3, 2))); alg) + V = BlockTensorKit.oplus(oneunit(S), space(Q, 1), oneunit(S)) + Q1 = typeof(W.A)(undef, space(Q, 1) ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}()) + Q2 = transpose(Q, ((1, 3), (2,))) + else + tmp = transpose(cat(insertleftunit(B′, 4), W.A; dims = 4), ((1,), (3, 4, 2))) + if trscheme == notrunc() + R, Q = rightorth!(tmp; alg) + else + @assert alg == SVD() || alg == SDD() + # TODO: truncation might not be using normalized tensors? + U, Σ, Q = tsvd!(tmp; alg, trunc = trscheme) + R = U * Σ + end + scale!(Q, d) + scale!(R, inv(d)) + Q′ = transpose(Q, ((1, 4), (2, 3))) + Q1 = SparseBlockTensorMap(Q′[1, 1, 1, 2:end]) + Q2 = removeunit(SparseBlockTensorMap(Q′[1, 1, 1, 1:1]), 4) + V = BlockTensorKit.oplus(oneunit(S), left_virtualspace(Q′), oneunit(S)) + end H[i] = JordanMPOTensor(V ⊗ physicalspace(W) ← domain(W), Q1, Q2, W.C, W.D) # absorb into previous site W′ = H[i - 1] - @plansor A′[l p; p' r] := W′.A[l p; p' r] * R[r; r] - @plansor B′[l p; p' r'] := W′.A[l p; p' r'] * t[r'; r] R[l; r] * W′.B[r p; p'] + @plansor A′[l p; p' r] := W′.A[l p; p' r'] * R[r'; r] + @plansor B′[l p; p' r] := W′.A[l p; p' r'] * t[r'; r] B′ = add!(removeunit(B′, 4), W′.B) @plansor C′[p; p' r] := W′.C[p; p' r'] * R[r'; r] - @plansor D′[p; p' r] := W′.C[p; p' l] *t[l; r] - D′ = add!(removeunit(D, 3), W′.D) + @plansor D′[p; p' r] := W′.C[p; p' r'] * t[r'; r] + D′ = add!(removeunit(D′, 3), W′.D) - H[i - 1] = JordanMPOTensor(V ⊗ physicalspace(W′) ← domain(W′), A′, B′, C′, D′) + H[i - 1] = JordanMPOTensor(codomain(W′) ← physicalspace(W′) ⊗ V, A′, B′, C′, D′) return H end From 602b398fdd6ef64f9493d63493f74a94f3c55e5b Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 4 Aug 2025 13:41:34 -0400 Subject: [PATCH 08/17] fix more edge cases --- src/operators/ortho.jl | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl index d34b0a6c5..4f076cdb2 100644 --- a/src/operators/ortho.jl +++ b/src/operators/ortho.jl @@ -88,10 +88,15 @@ function left_canonicalize!( # absorb into next site W′ = H[i + 1] - @plansor A′[l p; p' r] := R[l; r'] * W′.A[r' p; p' r] + if i != length(H) - 1 + @plansor A′[l p; p' r] := R[l; r'] * W′.A[r' p; p' r] + @plansor C′[l p; p' r] := t[l; r'] * W′.A[r' p; p' r] + C′ = add!(removeunit(C′, 1), W′.C) + else + A′ = similar(W′.A, space(R, 1) ⊗ physicalspace(W′) ← domain(W′.A)) + C′ = W′.C + end @plansor B′[l p; p'] := R[l; r] * W′.B[r p; p'] - @plansor C′[l p; p' r] := t[l; r'] * W′.A[r' p; p' r] - C′ = add!(removeunit(C′, 1), W′.C) @plansor D′[l p; p'] := t[l; r] * W′.B[r p; p'] D′ = add!(removeunit(D′, 1), W′.D) @@ -101,7 +106,7 @@ end function right_canonicalize!( H::MPOHamiltonian, i::Int; - alg = RQpos(), trscheme::TruncationScheme = notrunc() + alg = LQpos(), trscheme::TruncationScheme = notrunc() ) if H isa FiniteMPOHamiltonian 1 < i ≤ length(H) || throw(ArgumentError("Bounds error in canonicalize")) @@ -156,9 +161,14 @@ function right_canonicalize!( # absorb into previous site W′ = H[i - 1] - @plansor A′[l p; p' r] := W′.A[l p; p' r'] * R[r'; r] - @plansor B′[l p; p' r] := W′.A[l p; p' r'] * t[r'; r] - B′ = add!(removeunit(B′, 4), W′.B) + if i != 2 + @plansor A′[l p; p' r] := W′.A[l p; p' r'] * R[r'; r] + @plansor B′[l p; p' r] := W′.A[l p; p' r'] * t[r'; r] + B′ = add!(removeunit(B′, 4), W′.B) + else + A′ = similar(W′.A, codomain(W′.A) ← physicalspace(W′.A) ⊗ space(R, 2)') + B′ = W′.B + end @plansor C′[p; p' r] := W′.C[p; p' r'] * R[r'; r] @plansor D′[p; p' r] := W′.C[p; p' r'] * t[r'; r] D′ = add!(removeunit(D′, 3), W′.D) From 5565f287d87b772568e661b523741dd22dd43b51 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 4 Aug 2025 12:22:09 -0400 Subject: [PATCH 09/17] Add copying version --- src/algorithms/changebonds/svdcut.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl index c87c50d19..feb2f9936 100644 --- a/src/algorithms/changebonds/svdcut.jl +++ b/src/algorithms/changebonds/svdcut.jl @@ -106,6 +106,7 @@ function changebonds(ψ, H, alg::SvdCut, envs = environments(ψ, H)) return changebonds(ψ, alg), envs end +changebonds(mpo::FiniteMPOHamiltonian, alg::SvdCut) = changebonds!(copy(mpo), alg) function changebonds!(H::FiniteMPOHamiltonian, alg::SvdCut) # orthogonality center to the left for i in length(H):-1:2 From 3bbf163c0cf48777e63c123164326ba9a2c8a0bb Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 20 Aug 2025 14:49:50 +0200 Subject: [PATCH 10/17] formatter --- src/operators/jordanmpotensor.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl index ba011ff3e..7fcc7c6e8 100644 --- a/src/operators/jordanmpotensor.jl +++ b/src/operators/jordanmpotensor.jl @@ -90,7 +90,8 @@ function JordanMPOTensor( A::AbstractTensorMap{E, S, 2, 2}, B::AbstractTensorMap{E, S, 2, 1}, C::AbstractTensorMap{E, S, 1, 2}, D::AbstractTensorMap{E, S, 1, 1} ) where {E, S} - return JordanMPOTensor(V, + return JordanMPOTensor( + V, A isa SparseBlockTensorMap ? A : SparseBlockTensorMap(A), B isa SparseBlockTensorMap ? B : SparseBlockTensorMap(B), C isa SparseBlockTensorMap ? C : SparseBlockTensorMap(C), From 84273bf542f5a710923c175d6133c9a0b9aeb251 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 21 Aug 2025 11:25:55 +0200 Subject: [PATCH 11/17] little bit of cleanup --- src/operators/ortho.jl | 44 ++++-------------------------------------- 1 file changed, 4 insertions(+), 40 deletions(-) diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl index 4f076cdb2..93befc6ce 100644 --- a/src/operators/ortho.jl +++ b/src/operators/ortho.jl @@ -1,42 +1,8 @@ -function TensorKit.leftorth(W::JordanMPOTensor; alg = QRpos()) - # orthogonalize second column against first - WI = removeunit(W[1, 1, 1, 1], 1) - @tensor t[l; r] := conj(WI[p; p' l]) * W.C[p; p' r] - I = sectortype(W) - S = spacetype(W) - - @plansor t[l; r] := conj(removeunit(W[1, 1, 1, 1], 1)[p; p' l]) * W.C[p; p' r] - # @plansor C′[p; p' r] := W.C[p; p' r] - WI[p; p' l] * t[l; r] - @plansor C′[p; p' r] := -WI[p; p' l] * t[l; r] - add!(C′, W.C) - - # QR of second column - CA = transpose(cat(insertleftunit(C′, 1), W.A; dims = 1), ((3, 1, 2), (4,))) - Q, r = leftorth!(CA; alg) - Q′ = transpose(Q, ((2, 3), (1, 4))) - V = codomain(W) ← physicalspace(W) ⊗ BlockTensorKit.oplus(oneunit(S), right_virtualspace(Q′), oneunit(S)) - Q1 = SparseBlockTensorMap(Q′[2:end, 1, 1, 1]) - Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1) - - # Assemble output - W′ = JordanMPOTensor(V, Q1, W.B, Q2, W.D) - - R = similar(W, right_virtualspace(W′) ← right_virtualspace(W′)) - R[1, 1] = id!(R[1, 1]) - R[1, 2] = t - R[2, 2] = r - R[3, 3] = id!(R[3, 3]) - - return W′, R -end - function left_canonicalize!( - H::MPOHamiltonian, i::Int; + H::FiniteMPOHamiltonian, i::Int; alg = QRpos(), trscheme::TruncationScheme = notrunc() ) - if H isa FiniteMPOHamiltonian - 1 ≤ i < length(H) || throw(ArgumentError("Bounds error in canonicalize")) - end + 1 ≤ i < length(H) || throw(ArgumentError("Bounds error in canonicalize")) W = H[i] S = spacetype(W) @@ -105,12 +71,10 @@ function left_canonicalize!( end function right_canonicalize!( - H::MPOHamiltonian, i::Int; + H::FiniteMPOHamiltonian, i::Int; alg = LQpos(), trscheme::TruncationScheme = notrunc() ) - if H isa FiniteMPOHamiltonian - 1 < i ≤ length(H) || throw(ArgumentError("Bounds error in canonicalize")) - end + 1 < i ≤ length(H) || throw(ArgumentError("Bounds error in canonicalize")) W = H[i] S = spacetype(W) From b0586268cbfd59f798b4b2017673c237c8c95a1b Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 21 Aug 2025 12:09:03 +0200 Subject: [PATCH 12/17] Add some tests --- test/operators.jl | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/test/operators.jl b/test/operators.jl index 12f3a9706..69a84d5db 100644 --- a/test/operators.jl +++ b/test/operators.jl @@ -124,11 +124,13 @@ module TestOperators @testset "Finite MPOHamiltonian" begin L = 3 T = ComplexF64 - for V in (ℂ^2, U1Space(0 => 1, 1 => 1)) + for V in (ℂ^2, U1Space(-1 => 1, 0 => 1, 1 => 1)) lattice = fill(V, L) - O₁ = rand(T, V, V) + O₁ = randn(T, V, V) + O₁ += O₁' E = id(storagetype(O₁), domain(O₁)) - O₂ = rand(T, V^2 ← V^2) + O₂ = randn(T, V^2 ← V^2) + O₂ += O₂' H1 = FiniteMPOHamiltonian(lattice, i => O₁ for i in 1:L) H2 = FiniteMPOHamiltonian(lattice, (i, i + 1) => O₂ for i in 1:(L - 1)) @@ -164,6 +166,9 @@ module TestOperators FiniteMPOHamiltonian(lattice, 3 => O₁) @test 0.8 * H1 + 0.2 * H1 ≈ H1 atol = 1.0e-6 @test convert(TensorMap, H1 + H2) ≈ convert(TensorMap, H1) + convert(TensorMap, H2) atol = 1.0e-6 + H1_trunc = changebonds(H1, SvdCut(; trscheme=truncdim(0))) + @test H1_trunc ≈ H1 + @test all(left_virtualspace(H1_trunc) .== left_virtualspace(H1)) # test dot and application state = rand(T, prod(lattice)) @@ -181,12 +186,21 @@ module TestOperators vertical_operators = Dict( (I, I + I_vertical) => O₂ for I in eachindex(IndexCartesian(), square) if I[1] < size(square, 1) ) - operators = merge(local_operators, vertical_operators) + I_horizontal = CartesianIndex(0, 1) + horizontal_operators = Dict( + (I, I + I_horizontal) => O₂ for I in eachindex(IndexCartesian(), square) if I[2] < size(square, 1) + ) + operators = merge(local_operators, vertical_operators, horizontal_operators) H4 = FiniteMPOHamiltonian(grid, operators) @test H4 ≈ FiniteMPOHamiltonian(grid, local_operators) + - FiniteMPOHamiltonian(grid, vertical_operators) + FiniteMPOHamiltonian(grid, vertical_operators) + + FiniteMPOHamiltonian(grid, horizontal_operators) atol=1e-4 + + H5 = changebonds(H4 / 3 + 2H4 / 3, SvdCut(; trscheme = truncbelow(1.0e-12))) + psi = FiniteMPS(physicalspace(H5), V ⊕ oneunit(V)) + @test expectation_value(psi, H4) ≈ expectation_value(psi, H5) end end From dcb2587e208e95da610e439a889bd1cac2c4e631 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 21 Aug 2025 15:16:55 +0200 Subject: [PATCH 13/17] Handle more edge cases --- src/operators/ortho.jl | 129 +++++++++++++++++++++++++++-------------- 1 file changed, 87 insertions(+), 42 deletions(-) diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl index 93befc6ce..c3ef8cdd3 100644 --- a/src/operators/ortho.jl +++ b/src/operators/ortho.jl @@ -15,58 +15,83 @@ function left_canonicalize!( @plansor C′[p; p' r] := -WI[p; p' l] * t[l; r] add!(C′, W.C) - # QR of second colum - if i == 1 + # QR of second column + if size(W, 1) == 1 tmp = transpose(C′, ((2, 1), (3,))) + if trscheme == notrunc() Q, R = leftorth!(tmp; alg) else @assert alg == SVD() || alg == SDD() - # TODO: truncation might not be using normalized tensors? Q, Σ, Vᴴ = tsvd!(tmp; alg, trunc = trscheme) R = Σ * Vᴴ end - scale!(Q, d) - scale!(R, inv(d)) - - V = BlockTensorKit.oplus(oneunit(S), space(R, 1), oneunit(S)) - Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ space(R, 1)) - Q2 = transpose(Q, ((2,), (1, 3))) + if dim(R) == 0 # fully truncated + V = BlockTensorKit.oplus(oneunit(S), oneunit(S)) + Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}()) + Q2 = typeof(W.C)(undef, physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}()) + else + V = BlockTensorKit.oplus(oneunit(S), space(R, 1), oneunit(S)) + scale!(Q, d) + scale!(R, inv(d)) + Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ space(R, 1)) + Q2 = transpose(Q, ((2,), (1, 3))) + end + H[i] = JordanMPOTensor(codomain(W) ← physicalspace(W) ⊗ V, Q1, W.B, Q2, W.D) else tmp = transpose(cat(insertleftunit(C′, 1), W.A; dims = 1), ((3, 1, 2), (4,))) if trscheme == notrunc() Q, R = leftorth!(tmp; alg) else @assert alg == SVD() || alg == SDD() - # TODO: truncation might not be using normalized tensors? Q, Σ, Vᴴ = tsvd!(tmp; alg, trunc = trscheme) R = Σ * Vᴴ end - scale!(Q, d) - scale!(R, inv(d)) - Q′ = transpose(Q, ((2, 3), (1, 4))) - Q1 = Q′[2:end, 1, 1, 1] - Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1) - V = BlockTensorKit.oplus(oneunit(S), right_virtualspace(Q′), oneunit(S)) + if dim(R) == 0 # fully truncated + V = BlockTensorKit.oplus(oneunit(S), oneunit(S)) + Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}()) + Q2 = typeof(W.C)(undef, physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}()) + else + scale!(Q, d) + scale!(R, inv(d)) + Q′ = transpose(Q, ((2, 3), (1, 4))) + Q1 = Q′[2:end, 1, 1, 1] + Q2 = removeunit(SparseBlockTensorMap(Q′[1:1, 1, 1, 1]), 1) + V = BlockTensorKit.oplus(oneunit(S), right_virtualspace(Q′), oneunit(S)) + end + H[i] = JordanMPOTensor(codomain(W) ← physicalspace(W) ⊗ V, Q1, W.B, Q2, W.D) end - H[i] = JordanMPOTensor(codomain(W) ← physicalspace(W) ⊗ V, Q1, W.B, Q2, W.D) # absorb into next site W′ = H[i + 1] - if i != length(H) - 1 + + if size(W′, 4) > 1 && dim(R) != 0 @plansor A′[l p; p' r] := R[l; r'] * W′.A[r' p; p' r] + else + A′ = similar(W′.A, right_virtualspace(H[i].A) ⊗ physicalspace(W′) ← domain(W′.A)) + end + + if size(W′, 4) > 1 @plansor C′[l p; p' r] := t[l; r'] * W′.A[r' p; p' r] C′ = add!(removeunit(C′, 1), W′.C) else - A′ = similar(W′.A, space(R, 1) ⊗ physicalspace(W′) ← domain(W′.A)) - C′ = W′.C + C′ = W′.C # empty end - @plansor B′[l p; p'] := R[l; r] * W′.B[r p; p'] + + if dim(R) != 0 + @plansor B′[l p; p'] := R[l; r] * W′.B[r p; p'] + else + B′ = similar(W′.B, right_virtualspace(H[i].A) ⊗ physicalspace(W′) ← domain(W′.B)) + end + @plansor D′[l p; p'] := t[l; r] * W′.B[r p; p'] D′ = add!(removeunit(D′, 1), W′.D) - H[i + 1] = JordanMPOTensor(V ⊗ physicalspace(W′) ← domain(W′), A′, B′, C′, D′) + H[i + 1] = JordanMPOTensor( + right_virtualspace(H[i]) ⊗ physicalspace(W′) ← domain(W′), + A′, B′, C′, D′ + ) return H end @@ -87,56 +112,76 @@ function right_canonicalize!( add!(B′, W.B) # LQ of second row - if i == length(H) + if size(W, 4) == 1 tmp = transpose(B′, ((1,), (3, 2))) if trscheme == notrunc() R, Q = rightorth!(tmp; alg) else @assert alg == SVD() || alg == SDD() - # TODO: truncation might not be using normalized tensors? U, Σ, Q = tsvd!(tmp; alg, trunc = trscheme) R = U * Σ end - scale!(Q, d) - scale!(R, inv(d)) - R, Q = rightorth!(transpose(B′, ((1,), (3, 2))); alg) - V = BlockTensorKit.oplus(oneunit(S), space(Q, 1), oneunit(S)) - Q1 = typeof(W.A)(undef, space(Q, 1) ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}()) - Q2 = transpose(Q, ((1, 3), (2,))) + if dim(R) == 0 + V = BlockTensorKit.oplus(oneunit(S), oneunit(S)) + Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}()) + Q2 = typeof(W.B)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W)) + else + V = BlockTensorKit.oplus(oneunit(S), space(Q, 1), oneunit(S)) + scale!(Q, d) + scale!(R, inv(d)) + Q1 = typeof(W.A)(undef, space(Q, 1) ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}()) + Q2 = transpose(Q, ((1, 3), (2,))) + end + H[i] = JordanMPOTensor(V ⊗ physicalspace(W) ← domain(W), Q1, Q2, W.C, W.D) else tmp = transpose(cat(insertleftunit(B′, 4), W.A; dims = 4), ((1,), (3, 4, 2))) if trscheme == notrunc() R, Q = rightorth!(tmp; alg) else @assert alg == SVD() || alg == SDD() - # TODO: truncation might not be using normalized tensors? U, Σ, Q = tsvd!(tmp; alg, trunc = trscheme) R = U * Σ end - scale!(Q, d) - scale!(R, inv(d)) - Q′ = transpose(Q, ((1, 4), (2, 3))) - Q1 = SparseBlockTensorMap(Q′[1, 1, 1, 2:end]) - Q2 = removeunit(SparseBlockTensorMap(Q′[1, 1, 1, 1:1]), 4) - V = BlockTensorKit.oplus(oneunit(S), left_virtualspace(Q′), oneunit(S)) + if dim(R) == 0 + V = BlockTensorKit.oplus(oneunit(S), oneunit(S)) + Q1 = typeof(W.A)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W) ⊗ SumSpace{S}()) + Q2 = typeof(W.B)(undef, SumSpace{S}() ⊗ physicalspace(W) ← physicalspace(W)) + else + scale!(Q, d) + scale!(R, inv(d)) + Q′ = transpose(Q, ((1, 4), (2, 3))) + Q1 = SparseBlockTensorMap(Q′[1, 1, 1, 2:end]) + Q2 = removeunit(SparseBlockTensorMap(Q′[1, 1, 1, 1:1]), 4) + V = BlockTensorKit.oplus(oneunit(S), left_virtualspace(Q′), oneunit(S)) + end + H[i] = JordanMPOTensor(V ⊗ physicalspace(W) ← domain(W), Q1, Q2, W.C, W.D) end - H[i] = JordanMPOTensor(V ⊗ physicalspace(W) ← domain(W), Q1, Q2, W.C, W.D) # absorb into previous site W′ = H[i - 1] - if i != 2 + + if size(W′, 1) > 1 && dim(R) != 0 @plansor A′[l p; p' r] := W′.A[l p; p' r'] * R[r'; r] + else + A′ = similar(W′.A, codomain(W′.A) ← physicalspace(W′.A) ⊗ left_virtualspace(H[i].A)) + end + + if size(W′, 1) > 1 @plansor B′[l p; p' r] := W′.A[l p; p' r'] * t[r'; r] B′ = add!(removeunit(B′, 4), W′.B) else - A′ = similar(W′.A, codomain(W′.A) ← physicalspace(W′.A) ⊗ space(R, 2)') B′ = W′.B end - @plansor C′[p; p' r] := W′.C[p; p' r'] * R[r'; r] + + if dim(R) != 0 + @plansor C′[p; p' r] := W′.C[p; p' r'] * R[r'; r] + else + C′ = similar(W′.C, codomain(W′.C) ← physicalspace(W′) ⊗ left_virtualspace(H[i].A)) + end + @plansor D′[p; p' r] := W′.C[p; p' r'] * t[r'; r] D′ = add!(removeunit(D′, 3), W′.D) - H[i - 1] = JordanMPOTensor(codomain(W′) ← physicalspace(W′) ⊗ V, A′, B′, C′, D′) return H end From a3bc86614577fddb279147fa958ad9989c3f3262 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 22 Aug 2025 16:51:28 +0200 Subject: [PATCH 14/17] Bump BlockTensorKit compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f5a00ec63..27b180f8b 100644 --- a/Project.toml +++ b/Project.toml @@ -25,7 +25,7 @@ VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8" [compat] Accessors = "0.1" Aqua = "0.8.9" -BlockTensorKit = "0.1.6" +BlockTensorKit = "0.2" Combinatorics = "1" Compat = "3.47, 4.10" DocStringExtensions = "0.9.3" From efde24ea705a241281a0e97da1e5a2be6fb0b99d Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 28 Aug 2025 21:41:07 +0200 Subject: [PATCH 15/17] Formatter --- test/operators.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/operators.jl b/test/operators.jl index 69a84d5db..3ea1ae003 100644 --- a/test/operators.jl +++ b/test/operators.jl @@ -166,7 +166,7 @@ module TestOperators FiniteMPOHamiltonian(lattice, 3 => O₁) @test 0.8 * H1 + 0.2 * H1 ≈ H1 atol = 1.0e-6 @test convert(TensorMap, H1 + H2) ≈ convert(TensorMap, H1) + convert(TensorMap, H2) atol = 1.0e-6 - H1_trunc = changebonds(H1, SvdCut(; trscheme=truncdim(0))) + H1_trunc = changebonds(H1, SvdCut(; trscheme = truncdim(0))) @test H1_trunc ≈ H1 @test all(left_virtualspace(H1_trunc) .== left_virtualspace(H1)) @@ -196,7 +196,7 @@ module TestOperators @test H4 ≈ FiniteMPOHamiltonian(grid, local_operators) + FiniteMPOHamiltonian(grid, vertical_operators) + - FiniteMPOHamiltonian(grid, horizontal_operators) atol=1e-4 + FiniteMPOHamiltonian(grid, horizontal_operators) atol = 1.0e-4 H5 = changebonds(H4 / 3 + 2H4 / 3, SvdCut(; trscheme = truncbelow(1.0e-12))) psi = FiniteMPS(physicalspace(H5), V ⊕ oneunit(V)) From 5728f5c032e2ddd8aded5f724ae0ad5162cc01e3 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 31 Aug 2025 13:25:49 +0200 Subject: [PATCH 16/17] update docstring and reference --- docs/src/assets/mpskit.bib | 16 ++++++++++++++++ src/algorithms/changebonds/svdcut.jl | 9 ++++++++- 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/docs/src/assets/mpskit.bib b/docs/src/assets/mpskit.bib index 11b22006a..41320c7ed 100644 --- a/docs/src/assets/mpskit.bib +++ b/docs/src/assets/mpskit.bib @@ -196,6 +196,22 @@ @article{paeckel2019 keywords = {Density matrix renormalization group (DMRG),Matrix-product states (MPS),Strongly-correlated systems,Time-dependent variational principle (TDVP),Time-evolution methods,Time-evolving block decimation (TEBD)} } +@article{parker2020, + title = {Local Matrix Product Operators: {{Canonical}} Form, Compression, and Control Theory}, + shorttitle = {Local Matrix Product Operators}, + author = {Parker, Daniel E. and Cao, Xiangyu and Zaletel, Michael P.}, + year = {2020}, + month = jul, + journal = {Physical Review B}, + volume = {102}, + number = {3}, + pages = {035147}, + publisher = {American Physical Society}, + doi = {10.1103/PhysRevB.102.035147}, + url = {https://link.aps.org/doi/10.1103/PhysRevB.102.035147}, + abstract = {We present a new method for compressing matrix product operators (MPOs) which represent sums of local terms, such as Hamiltonians. Just as with area law states, such local operators may be fully specified with a small amount of information per site. Standard matrix product state (MPS) tools are ill-suited to this case, due to extensive Schmidt values that coexist with intensive ones, and Jordan blocks in the transfer matrix. We ameliorate these issues by introducing an ``almost Schmidt decomposition'' that respects locality. Our method is ``{$\varepsilon$}-close'' to the accuracy of MPS-based methods for finite MPOs, and extends seamlessly to the thermodynamic limit, where MPS techniques are inapplicable. In the framework of control theory, our method generalizes Kung's algorithm for model order reduction. Several examples are provided, including an all-MPO version of the operator recursion method (Lanczos algorithm) directly in the thermodynamic limit. All results are accompanied by practical algorithms, well-suited for the large MPOs that arise in DMRG for long-range or quasi-2D models.} +} + @article{pollmann2009, title = {Theory of {{Finite-Entanglement Scaling}} at {{One-Dimensional Quantum Critical Points}}}, author = {Pollmann, Frank and Mukerjee, Subroto and Turner, Ari M. and Moore, Joel E.}, diff --git a/src/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl index d9340cf42..7ad26e706 100644 --- a/src/algorithms/changebonds/svdcut.jl +++ b/src/algorithms/changebonds/svdcut.jl @@ -1,11 +1,18 @@ """ $(TYPEDEF) -An algorithm that uses truncated SVD to change the bond dimension of a ψ. +An algorithm that uses truncated SVD to change the bond dimension of a state or operator. +This is achieved by a sweeping algorithm that locally performs (optimal) truncations in a gauged basis. + +See also [`changebonds(!)`](@ref changebonds) ## Fields $(TYPEDFIELDS) + +## References + +* [Parker et al. Phys. Rev. B 102 (2020)](@cite parker2020) """ @kwdef struct SvdCut{S} <: Algorithm "algorithm used for the singular value decomposition" From 22d1a3574b65c2e44c2db4ef51321fe379562fcf Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 5 Sep 2025 12:48:55 -0400 Subject: [PATCH 17/17] Add comment explaining todo --- src/operators/ortho.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/operators/ortho.jl b/src/operators/ortho.jl index c3ef8cdd3..1b7ecb924 100644 --- a/src/operators/ortho.jl +++ b/src/operators/ortho.jl @@ -11,6 +11,7 @@ function left_canonicalize!( # orthogonalize second column against first WI = removeunit(W[1, 1, 1, 1], 1) @tensor t[l; r] := conj(WI[p; p' l]) * W.C[p; p' r] + # TODO: the following is currently broken due to a TensorKit bug # @plansor C′[p; p' r] := W.C[p; p' r] - WI[p; p' l] * t[l; r] @plansor C′[p; p' r] := -WI[p; p' l] * t[l; r] add!(C′, W.C) @@ -108,6 +109,8 @@ function right_canonicalize!( # orthogonalize second row against last WI = removeunit(W[end, 1, 1, end], 4) @plansor t[l; r] := conj(WI[r p; p']) * W.B[l p; p'] + # TODO: the following is currently broken due to a TensorKit bug + # @plansor B′[l p; p'] := W.B[l p; p'] - WI[r p; p'] * t[l; r] @plansor B′[l p; p'] := -WI[r p; p'] * t[l; r] add!(B′, W.B)