[Feature] FiniteMPOHamiltonian compression

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"

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.},

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")

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"
@@ -105,3 +112,22 @@ end
 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
+        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

src/operators/jordanmpotensor.jl

@@ -67,26 +67,37 @@ 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)])
-    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
+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

src/operators/ortho.jl

@@ -0,0 +1,190 @@
+function left_canonicalize!(
+        H::FiniteMPOHamiltonian, i::Int;
+        alg = QRpos(), trscheme::TruncationScheme = notrunc()
+    )
+    1 ≤ i < length(H) || throw(ArgumentError("Bounds error in canonicalize"))
+
+    W = H[i]
+    S = spacetype(W)
+    d = sqrt(dim(physicalspace(W)))
+
+    # 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)
+
+    # 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()
+            Q, Σ, Vᴴ = tsvd!(tmp; alg, trunc = trscheme)
+            R = Σ * Vᴴ
+        end
+
+        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()
+            Q, Σ, Vᴴ = tsvd!(tmp; alg, trunc = trscheme)
+            R = Σ * Vᴴ
+        end
+        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
+
+    # absorb into next site
+    W′ = H[i + 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
+        C′ = W′.C # empty
+    end
+
+    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(
+        right_virtualspace(H[i]) ⊗ physicalspace(W′) ← domain(W′),
+        A′, B′, C′, D′
+    )
+    return H
+end
+
+function right_canonicalize!(
+        H::FiniteMPOHamiltonian, i::Int;
+        alg = LQpos(), trscheme::TruncationScheme = notrunc()
+    )
+    1 < i ≤ length(H) || throw(ArgumentError("Bounds error in canonicalize"))
+
+    W = H[i]
+    S = spacetype(W)
+    d = sqrt(dim(physicalspace(W)))
+
+    # 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)
+
+    # LQ of second row
+    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()
+            U, Σ, Q = tsvd!(tmp; alg, trunc = trscheme)
+            R = U * Σ
+        end
+
+        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()
+            U, Σ, Q = tsvd!(tmp; alg, trunc = trscheme)
+            R = U * Σ
+        end
+        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
+
+    # absorb into previous site
+    W′ = H[i - 1]
+
+    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
+        B′ = W′.B
+    end
+
+    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

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 = 1.0e-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