Zipper compression for FiniteMPO-FiniteMPS product

src/MPSKit.jl

@@ -38,6 +38,7 @@ export time_evolve, timestep, timestep!, make_time_mpo
 export TDVP, TDVP2, WI, WII, TaylorCluster
 export changebonds, changebonds!
 export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand
+export Zipper
 export propagator
 export DynamicalDMRG, NaiveInvert, Jeckelmann
 export exact_diagonalization, fidelity_susceptibility
@@ -188,6 +189,7 @@ include("algorithms/statmech/idmrg.jl")
 include("algorithms/fidelity_susceptibility.jl")
 
 include("algorithms/approximate/approximate.jl")
+include("algorithms/approximate/zipper.jl")
 include("algorithms/approximate/vomps.jl")
 include("algorithms/approximate/fvomps.jl")
 include("algorithms/approximate/idmrg.jl")

src/algorithms/approximate/approximate.jl

@@ -4,10 +4,11 @@
     approximate!(ψ₀, (O, ψ), algorithm, [environments]) -> (ψ, environments, ϵ)
     approximate(ψ₀, ψ, algorithm, [environments]) -> (ψ, environments, ϵ)
     approximate!(ψ₀, ψ, algorithm, [environments]) -> (ψ, environments, ϵ)
+    approximate((O, ψ), algorithm) -> ψ′
 
 Compute an approximation to the application of an operator `O` to the state `ψ` in the form
-of an MPS `ψ₀`. If only a state `ψ` is supplied instead of the `(O, ψ)` pair, `ψ₀` is
-approximated directly to `ψ` (i.e. `O` is taken to be the identity).
+of an MPS, using initial guess `ψ₀`. If only a state `ψ` is supplied instead of the `(O, ψ)` pair,
+`ψ₀` is approximated directly to `ψ` (i.e. `O` is taken to be the identity).
 
 **Not every algorithm supports every combination of arguments below** — see the per-algorithm
 notes at the end of this docstring before picking one.
@@ -41,6 +42,7 @@ infinite algorithms always require an explicit `(O, ψ)` tuple, and **`VOMPS` ha
 |:--------- |:----------------------------- |:---------------------------------- |:------------------:|:--------------:|
 | `DMRG`    | single-site, fixes bond dim    | `AbstractFiniteMPS`                | ✅                  | ✅              |
 | `DMRG2`   | two-site, truncates via `trscheme` | `AbstractFiniteMPS`            | ✅                  | ✅              |
+| `Zipper`  | right-to-left streaming MPO-MPS compression | none, uses `(O, ψ)` directly | ❌                  | ❌              |
 | `IDMRG`   | single-site, thermodynamic limit | `InfiniteMPS` / `MultilineMPS`  | ❌ (tuple only)     | ✅              |
 | `IDMRG2`  | two-site, thermodynamic limit, needs unit cell ≥ 2 | `InfiniteMPS` / `MultilineMPS` | ❌ (tuple only) | ✅ |
 | `VOMPS`   | tangent-space truncation       | `InfiniteMPS` / `MultilineMPS`     | ❌ (tuple only)     | ❌ (out-of-place only) |

src/algorithms/approximate/zipper.jl

@@ -0,0 +1,57 @@
+"""
+    Zipper(; alg_svd=Defaults.alg_svd(), trscheme)
+
+Algorithm that approximates an open-boundary finite MPO-MPS product using a right-to-left
+zipper sweep. The MPO and MPS are contracted one site at a time, and the enlarged virtual
+bond is truncated immediately using `trscheme`.
+
+Use as:
+
+    approximate((O, ψ), Zipper(; trscheme))
+
+This returns an unnormalized compression of `O * ψ`, comparable to
+`changebonds(O * ψ, SvdCut(; trscheme); normalize=false)`, but without storing
+the fully enlarged product on every site.
+"""
+@kwdef struct Zipper{S} <: Algorithm
+    "algorithm used for the singular value decomposition"
+    alg_svd::S = Defaults.alg_svd()
+
+    "algorithm used for truncation of the local gauge tensors"
+    trscheme::TruncationStrategy
+end
+
+function approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipper)
+    N = check_length(O, ψ)
+    if !isunitspace(left_virtualspace(O, 1)) || !isunitspace(right_virtualspace(O, N))
+        throw(ArgumentError("Zipper is only implemented for open-boundary MPOs"))
+    end
+
+    ψ′ = copy(ψ)
+    T = TensorOperations.promote_contract(scalartype(O), scalartype(ψ))
+    A = TensorKit.similarstoragetype(eltype(ψ), T)
+
+    Fᵣ = fuser(A, right_virtualspace(ψ′, N), right_virtualspace(O, N))
+    local carry
+    alg_gauge = MatrixAlgebraKit.TruncatedAlgorithm(alg.alg_svd, alg.trscheme)
+
+    As = map(N:-1:1) do i
+        Aψ = i == 1 ? ψ′.AC[1] : ψ′.AR[i]
+        physicalspace(Aψ) == physicalspace(O[i]) ||
+            throw(SpaceMismatch("MPO input physical space does not match MPS physical space at site $i"))
+        Fₗ = fuser(A, left_virtualspace(ψ′, i), left_virtualspace(O, i))
+        Aᶻ = _fuse_mpo_mps(O[i], Aψ, Fₗ, Fᵣ)
+        i < N && (Aᶻ = Aᶻ * carry)
+
+        if i == 1
+            return Aᶻ
+        else
+            C, AR = right_gauge(Aᶻ, alg_gauge)
+            carry = C
+            Fᵣ = Fₗ
+            return AR
+        end
+    end
+
+    return FiniteMPS(reverse(As); normalize = false, overwrite = true)
+end

test/algorithms/dynamical_dmrg.jl

@@ -60,4 +60,3 @@ end
         @test data ≈ predicted atol = 1.0e-8
     end
 end
-

test/algorithms/zipper.jl

@@ -0,0 +1,39 @@
+println("
+-----------------------------
+|      Zipper tests          |
+-----------------------------
+")
+
+using .TestSetup
+using Test
+using MPSKit
+using TensorKit
+using TensorKit: ℙ
+using Random
+
+spacelist = [(ℙ^4, ℙ^3), (Rep[SU₂](1 => 1), Rep[SU₂](0 => 2, 1 => 2, 2 => 1))]
+
+@testset "Finite MPO-MPS zipper $(spacetype(pspace))" for (pspace, Dspace) in spacelist
+    Random.seed!(1357)
+    L = 6
+    Wspace = Dspace
+    Vspaces = [oneunit(Wspace); fill(Wspace, L - 1); oneunit(Wspace)]
+    O = FiniteMPO(
+        [rand(ComplexF64, Vspaces[i] ⊗ pspace ← pspace ⊗ Vspaces[i + 1]) for i in 1:L]
+    )
+    ψ = FiniteMPS(rand, ComplexF64, L, pspace, Dspace)
+    O_copy = copy(O)
+    ψ_copy = copy(ψ)
+
+    trscheme = trunctol(; atol = 1.0e-10)
+    ref = changebonds(O * ψ, SvdCut(; trscheme); normalize = false)
+    got = approximate((O, ψ), Zipper(; trscheme))
+
+    @test norm(ref - got) / norm(ref) < 1.0e-10
+    @test norm(ψ - ψ_copy) < 1.0e-12
+    @test all(i -> norm(O[i] - O_copy[i]) < 1.0e-12, 1:length(O))
+
+    Dcut = 4
+    got_tr = approximate((O, ψ), Zipper(; trscheme = truncrank(Dcut)))
+    @test maximum(i -> dim(left_virtualspace(got_tr, i)), 2:length(got_tr)) <= Dcut
+end