Infinite MPOHamiltonian compression via SvdCut

docs/src/changelog.md

@@ -24,9 +24,16 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 - `ClusterExpansion` time-evolution algorithm for `make_time_mpo`, constructing evolution MPOs
   for finite and infinite nearest-neighbor Hamiltonians through a nonperturbative cluster expansion
   that matches exact exponentials on clusters of up to `N` sites. ([#252](https://github.com/QuantumKitHub/MPSKit.jl/pull/252))
+- `changebonds(::InfiniteMPOHamiltonian, ::SvdCut)` compresses infinite Hamiltonians with the
+  almost-Schmidt decomposition of [Parker et al.](https://arxiv.org/abs/1909.06341), and
+  `left_canonicalize!`/`right_canonicalize!` bring an `InfiniteMPOHamiltonian` into canonical form. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539))
 
 ### Changed
 
+- `left_canonicalize!`/`right_canonicalize!(::FiniteMPOHamiltonian)` now sweep over the whole chain
+  instead of a single site, and return the gauge transformations alongside the operator, matching
+  the infinite versions. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539))
+
 ### Deprecated
 
 ### Removed
@@ -35,6 +42,8 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 
 - `convert(TensorMap, H)` for single-site `FiniteMPOHamiltonian`s no longer contracts the site tensor with itself. ([#252](https://github.com/QuantumKitHub/MPSKit.jl/pull/252))
 - `show` methods for finite MPS/MPOs no longer check unit space condition through `oneunit`, which is not defined for multifusion spaces. ([#538](https://github.com/QuantumKitHub/MPSKit.jl/pull/538))
+- The canonical forms of `MPOHamiltonian`s, and therefore `SvdCut` compression, are now correct when
+  interactions contain identity components or redundant channels. ([#539](https://github.com/QuantumKitHub/MPSKit.jl/pull/539))
 
 ### Performance
 

src/algorithms/changebonds/svdcut.jl

@@ -153,22 +153,15 @@ 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
+    right_canonicalize!(H)
 
     # a Jordan bond has a start and a finish level on top of its interaction channels
     channels_before = [dim(left_virtualspace(H, i)) - 2 for i in 2:length(H)]
 
-    # swipe right
     alg_trunc = MatrixAlgebraKit.TruncatedAlgorithm(alg.alg_svd, alg.trunc)
-    for i in 1:(length(H) - 1)
-        H = left_canonicalize!(H, i; alg = MatrixAlgebraKit.LeftOrthViaSVD(alg_trunc))
-    end
-    # swipe left -- TODO: do we really need this double sweep?
-    for i in length(H):-1:2
-        H = right_canonicalize!(H, i; alg = MatrixAlgebraKit.RightOrthViaSVD(alg_trunc))
-    end
+    left_canonicalize!(H; alg = MatrixAlgebraKit.LeftOrthViaSVD(alg_trunc))
+    # TODO: do we really need this double sweep?
+    right_canonicalize!(H; alg = MatrixAlgebraKit.RightOrthViaSVD(alg_trunc))
 
     emptied = findall(
         i -> channels_before[i] > 0 && dim(left_virtualspace(H, i + 1)) - 2 == 0,
@@ -178,3 +171,55 @@ function changebonds!(H::FiniteMPOHamiltonian, alg::SvdCut)
 
     return H
 end
+
+changebonds(H::InfiniteMPOHamiltonian, alg::SvdCut) = changebonds!(copy(H), alg)
+function changebonds!(H::InfiniteMPOHamiltonian, alg::SvdCut)
+    channels_before = [dim(right_virtualspace(H[i])) - 2 for i in 1:length(H)]
+
+    # the interaction channels of the gauge between the right- and left-canonical forms
+    right_canonicalize!(H)
+    _, Cs = left_canonicalize!(H)
+
+    Us = map(C -> first(svd_trunc(C; trunc = alg.trunc, alg = alg.alg_svd)), Cs)
+    for i in 1:length(H)
+        H[i] = _project_channels(Us[i - 1], H[i], Us[i])
+    end
+
+    emptied = findall(
+        i -> channels_before[i] > 0 && dim(right_virtualspace(H[i])) - 2 == 0,
+        eachindex(channels_before)
+    )
+    isempty(emptied) || @warn "`SvdCut` discarded every interaction crossing bond(s) $emptied; those terms are gone from the compressed operator. Loosen `trunc` or check the scale of the input Hamiltonian."
+
+    return H
+end
+
+# Compute `diag(1, Uₗ', 1) * W * diag(1, Uᵣ, 1)`.
+function _project_channels(Uₗ, W::JordanMPOTensor, Uᵣ)
+    P = physicalspace(W)
+    Vₗ, Vᵣ = _channelspace(Uₗ), _channelspace(Uᵣ)
+    WA, WB, WC = W.A, W.B, W.C
+    Uₗᴴ = adjoint(Uₗ)
+
+    if dim(Vₗ) == 0 || dim(Vᵣ) == 0
+        A = similar(WA, Vₗ ⊗ P ← P ⊗ Vᵣ)
+    else
+        @plansor A[l p; p' r] := Uₗᴴ[l; l'] * WA[l' p; p' r'] * Uᵣ[r'; r]
+    end
+    if dim(Vₗ) == 0
+        B = similar(WB, Vₗ ⊗ P ← P)
+    else
+        @plansor B[l p; p'] := Uₗᴴ[l; l'] * WB[l' p; p']
+    end
+    if dim(Vᵣ) == 0
+        C = similar(WC, P ← P ⊗ Vᵣ)
+    else
+        @plansor C[p; p' r] := WC[p; p' r'] * Uᵣ[r'; r]
+    end
+
+    unit = rightunitspace(P)
+    V = ⊞(unit, Vₗ, unit) ⊗ P ← P ⊗ ⊞(unit, Vᵣ, unit)
+    return JordanMPOTensor(V, A, B, C, W.D)
+end
+
+_channelspace(U) = dim(U) == 0 ? SumSpace{spacetype(U)}() : domain(U)[1]