Zip-up compression for FiniteMPO-FiniteMPS product

docs/src/assets/mpskit.bib

@@ -498,6 +498,22 @@ @misc{shen2025
   keywords = {Condensed Matter - Strongly Correlated Electrons,High Energy Physics - Theory}
 }
 
+@article{stoudenmire2010,
+  title = {Minimally Entangled Typical Thermal State Algorithms},
+  author = {Stoudenmire, E. M. and White, Steven R.},
+  year = {2010},
+  month = may,
+  journal = {New Journal of Physics},
+  volume = {12},
+  number = {5},
+  pages = {055026},
+  doi = {10.1088/1367-2630/12/5/055026},
+  url = {https://doi.org/10.1088/1367-2630/12/5/055026},
+  eprint = {1002.1305},
+  archiveprefix = {arXiv},
+  primaryclass = {cond-mat.str-el}
+}
+
 % NOTE: lecture notes; MPSKit citation not independently confirmed during review.
 @misc{sinha2025,
   title = {Lectures on Quantum Field Theory on a Quantum Computer},

docs/src/changelog.md

@@ -23,6 +23,11 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 
 - Addition of `FiniteMPS`/`FiniteMPO` with different scalar types, through a new
   `Base.similar(ψ, ::Type{S})` for `S <: Number` on `FiniteMPS`. ([#484](https://github.com/QuantumKitHub/MPSKit.jl/pull/484))
+- `Zipup`, an algorithm for `approximate`/`approximate!` that compresses a finite MPO-MPS product in
+  a single sweep, optionally followed by a sweep in the opposite direction that imposes the final
+  truncation. The sweep direction is selected by the `left_to_right` keyword. Both
+  `approximate((O, ϕ), alg)` and `approximate!(ψ, (O, ϕ), alg)` are supported, where the destination
+  `ψ` is a write target rather than an initial guess and may alias `ϕ`; they return `(ψ, ϵ)`.
 
 ### Changed
 

src/MPSKit.jl

@@ -39,6 +39,7 @@ export TDVP, TDVP2, WI, WII, TaylorCluster
 export changebonds, changebonds!
 export VUMPSSvdCut, OptimalExpand, SvdCut, RandExpand, SketchedExpand
 export NoiseSchedule, FunctionalSchedule, ExponentialDecay, Warmup, DMRG3S
+export Zipup
 export propagator
 export DynamicalDMRG, NaiveInvert, Jeckelmann
 export exact_diagonalization, fidelity_susceptibility
@@ -192,6 +193,7 @@ include("algorithms/statmech/idmrg.jl")
 include("algorithms/fidelity_susceptibility.jl")
 
 include("algorithms/approximate/approximate.jl")
+include("algorithms/approximate/zipup.jl")
 include("algorithms/approximate/vomps.jl")
 include("algorithms/approximate/fvomps.jl")
 include("algorithms/approximate/idmrg.jl")

src/algorithms/approximate/approximate.jl

@@ -4,10 +4,12 @@
     approximate!(ψ₀, (O, ψ), algorithm, [environments]) -> (ψ, environments, ϵ)
     approximate(ψ₀, ψ, algorithm, [environments]) -> (ψ, environments, ϵ)
     approximate!(ψ₀, ψ, algorithm, [environments]) -> (ψ, environments, ϵ)
+    approximate((O, ψ), algorithm) -> (ψ′, ϵ)
+    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.
@@ -35,12 +37,15 @@ struct itself instead (e.g. `DMRG(; tol, maxiter, verbosity)`).
 Each algorithm below only supports a subset of the general interface. Check this table before
 picking one — in particular, note that **only `DMRG`/`DMRG2` accept a bare state `ψ`**; the
 infinite algorithms always require an explicit `(O, ψ)` tuple, and **`VOMPS` has no in-place
-`approximate!`** at all.
+`approximate!`** at all. `Zipup` is a single sweep rather than an iterative optimization, so it uses
+no environments and returns `(ψ, ϵ)`; its `ψ₀` is a write destination, not an initial guess, and it
+may be omitted.
 
 | Algorithm | Scheme                        | State `ψ₀`                        | bare `ψ` allowed? | `approximate!` |
 |:--------- |:----------------------------- |:---------------------------------- |:------------------:|:--------------:|
 | `DMRG`    | single-site, fixes bond dim    | `AbstractFiniteMPS`                | ✅                  | ✅              |
 | `DMRG2`   | two-site, truncates via `trunc` | `AbstractFiniteMPS`            | ✅                  | ✅              |
+| `Zipup`   | streaming MPO-MPS compression | `FiniteMPS` destination, optional | ❌ (tuple only)    | ✅              |
 | `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/zipup.jl

@@ -0,0 +1,182 @@
+"""
+$(TYPEDEF)
+
+Algorithm that approximates an open-boundary finite MPO-MPS product using a zip-up sweep, optionally
+followed by a zip-down sweep in the opposite direction. The MPO and MPS are contracted one site at a
+time, and the enlarged virtual bond is truncated immediately. The sweep direction is selected by
+`left_to_right`.
+
+    approximate((O, ϕ), alg::Zipup) -> ψ, ϵ
+    approximate!(ψ, (O, ϕ), alg::Zipup) -> ψ, ϵ
+
+Contrary to the variational algorithms, this algorithm requires no initial guess:
+the in-place version simply uses `ψ` as the destination of the sweep, overwriting its contents, and may alias `ϕ`.
+The out-of-place version allocates a destination with the promoted scalar type of `O` and `ϕ`.
+Both return the truncation error `ϵ` alongside the approximated state.
+
+## Fields
+
+$(TYPEDFIELDS)
+
+## Constructors
+
+    Zipup(; trunc, alg_svd=Defaults.alg_svd(), left_to_right=true)
+    Zipup(alg_zipup, [alg_zipdown]; left_to_right=true)
+
+Create a `Zipup` algorithm with the given truncated gauge algorithm, or by passing a truncation scheme and singular value decomposition algorithm.
+The keyword `trunc` can be either one truncation strategy for a single zip-up sweep, or a tuple `(zipup_trunc, zipdown_trunc)` for a zip-up sweep followed by a zip-down sweep.
+Equivalently, one can pass the corresponding truncated gauge algorithms directly as `alg_zipup` and `alg_zipdown`.
+The keyword `left_to_right` selects the direction of the zip-up sweep, the zip-down sweep always running in the opposite direction.
+
+Following Paeckel et al., if the desired final bond dimension is `D`, one can use a more permissive zip-up truncation, e.g. rank `2D` with stricter tolerances, and use `alg_zipdown` to impose the final truncation.
+
+## References
+
+- [Stoudenmire and White New J. Phys. 12 (2010)](@cite stoudenmire2010)
+- [Paeckel et al. Ann. of Phys. 411 (2019)](@cite paeckel2019)
+"""
+struct Zipup{
+        U <: MatrixAlgebraKit.TruncatedAlgorithm,
+        D <: Union{Nothing, MatrixAlgebraKit.TruncatedAlgorithm},
+    } <: Algorithm
+    "algorithm used for gauging and truncating the local tensors during the zip-up sweep"
+    alg_zipup::U
+    "algorithm used for the final locally gauged truncation pass; `nothing` skips this pass"
+    alg_zipdown::D
+    "if `true`, zip up from left to right and truncate from right to left, and vice versa"
+    left_to_right::Bool
+end
+
+function Zipup(alg_zipup, alg_zipdown = nothing; left_to_right::Bool = true)
+    return Zipup(alg_zipup, alg_zipdown, left_to_right)
+end
+
+function Zipup(; trunc, alg_svd = Defaults.alg_svd(), left_to_right::Bool = true)
+    if trunc isa TruncationStrategy
+        return Zipup(MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trunc); left_to_right)
+    elseif trunc isa Tuple{<:TruncationStrategy, <:TruncationStrategy}
+        alg_zipup = MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trunc[1])
+        alg_zipdown = MatrixAlgebraKit.TruncatedAlgorithm(alg_svd, trunc[2])
+        return Zipup(alg_zipup, alg_zipdown; left_to_right)
+    else
+        throw(ArgumentError("`trunc` should be a truncation strategy or a tuple of two truncation strategies"))
+    end
+end
+
+function approximate!(ψ::FiniteMPS, (O, ϕ)::Tuple{Any, <:FiniteMPS}, alg::Zipup)
+    N = check_length(ψ, O, ϕ)
+    T = TensorOperations.promote_contract(scalartype(O), scalartype(ϕ))
+    promote_type(T, scalartype(ψ)) === scalartype(ψ) ||
+        throw(ArgumentError("destination state with scalartype $(scalartype(ψ)) cannot hold the result with scalartype $T"))
+    for i in 1:N
+        physicalspace(ϕ, i) == _input_physicalspace(O[i]) ||
+            throw(SpaceMismatch("MPO input physical space does not match MPS physical space at site $i"))
+    end
+
+    return if alg.left_to_right
+        zip_left_right!(ψ, O, ϕ, alg.alg_zipup, alg.alg_zipdown)
+    else
+        zip_right_left!(ψ, O, ϕ, alg.alg_zipup, alg.alg_zipdown)
+    end
+end
+
+function approximate(Oϕ::Tuple{Any, <:FiniteMPS}, alg::Zipup)
+    O, ϕ = Oϕ
+    T = TensorOperations.promote_contract(scalartype(O), scalartype(ϕ))
+    return approximate!(similar(ϕ, T), Oϕ, alg)
+end
+
+@doc """
+    zip_left_right!(ψ, O, ϕ, alg_zipup, [alg_zipdown]) -> ψ, ϵ
+    zip_right_left!(ψ, O, ϕ, alg_zipup, [alg_zipdown]) -> ψ, ϵ
+
+Contract the MPO `O` with the MPS `ϕ` in a single sweep, truncating the enlarged virtual bond at every
+site with `alg_zipup`, and write the result into `ψ`. `zip_left_right!` zips up from left to right,
+`zip_right_left!` from right to left. Unless `alg_zipdown` is `nothing`, a second sweep in the
+opposite direction imposes a final truncation with `alg_zipdown` in a locally gauged basis, leaving
+the gauge center of `ψ` at the far end. The destination may alias `ϕ`.
+
+Also returns the truncation error `ϵ`, the largest 2-norm of the discarded singular values over all
+bonds and both sweeps.
+"""
+zip_left_right!
+@doc (@doc zip_left_right!) zip_right_left!
+
+function zip_left_right!(ψ::FiniteMPS, O, ϕ::FiniteMPS, alg_zipup, alg_zipdown = nothing)
+    N = length(ψ)
+
+    # obtain all input tensors before overwriting the destination, such that `ψ === ϕ` is allowed:
+    # from here on, the input is only queried through `Aϕs`, never through `ϕ` itself
+    Aϕs = map(i -> i == 1 ? ϕ.AC[1] : ϕ.AR[i], 1:N)
+
+    # the sweep re-derives the entire state: discard all cached tensors, as their spaces are stale
+    # TODO: "reallocate" tensors?"
+    foreach(f -> fill!(f, missing), (ψ.ALs, ψ.ARs, ψ.ACs, ψ.Cs))
+
+    A = storagetype(eltype(ψ))
+    Fₗ = fuser(A, left_virtualspace(Aϕs[1]), left_virtualspace(O, 1))
+    ϵ = zero(real(scalartype(ψ)))
+
+    # zip up from left to right, leaving the gauge center on the last site
+    for i in 1:(N - 1)
+        Aᶻ = _fuse_mpo_mps_left(O[i], Aϕs[i], Fₗ)
+        AL, Fₗ, ϵᵢ = left_gauge(Aᶻ, alg_zipup) # right factor doubles as the next left fuser
+        ψ.ALs[i] = AL
+        ϵ = max(ϵ, ϵᵢ)
+    end
+    Fᵣ = fuser(A, right_virtualspace(Aϕs[N]), right_virtualspace(O, N))
+    ψ.ACs[N] = _fuse_mpo_mps(O[N], Aϕs[N], Fₗ, Fᵣ)
+
+    # zip down from right to left, truncating in a locally gauged basis
+    if !isnothing(alg_zipdown)
+        for i in N:-1:2
+            ψ, ϵᵢ = right_gauge!(ψ, i, ψ.AC[i], alg_zipdown)
+            ϵ = max(ϵ, ϵᵢ)
+        end
+    end
+
+    return ψ, ϵ
+end
+
+function zip_right_left!(ψ::FiniteMPS, O, ϕ::FiniteMPS, alg_zipup, alg_zipdown = nothing)
+    N = length(ψ)
+
+    Aϕs = map(i -> i == N ? ϕ.AC[N] : ϕ.AL[i], 1:N)
+    foreach(f -> fill!(f, missing), (ψ.ALs, ψ.ARs, ψ.ACs, ψ.Cs))
+
+    A = storagetype(eltype(ψ))
+    # the right-hand fusers are oriented as `(Vmps ⊗ Vmpo) ← Vfused`, matching the factor that
+    # replaces them on the next site
+    Vᵣ = right_virtualspace(Aϕs[N]) ⊗ right_virtualspace(O, N)
+    Fᵣ = isomorphism(A, Vᵣ, fuse(Vᵣ))
+    ϵ = zero(real(scalartype(ψ)))
+
+    # zip up from right to left, leaving the gauge center on the first site
+    for i in N:-1:2
+        Aᶻ = _fuse_mpo_mps_right(O[i], Aϕs[i], Fᵣ)
+        Fᵣ, AR, ϵᵢ = _right_gauge_zip(Aᶻ, alg_zipup) # left factor doubles as the next right fuser
+        ψ.ARs[i] = AR
+        ϵ = max(ϵ, ϵᵢ)
+    end
+    # the carry is oriented such that it can simply be composed with the last local tensor
+    Fₗ = fuser(A, left_virtualspace(Aϕs[1]), left_virtualspace(O, 1))
+    ψ.ACs[1] = _fuse_mpo_mps_left(O[1], Aϕs[1], Fₗ) * Fᵣ
+
+    # zip down from left to right, truncating in a locally gauged basis
+    if !isnothing(alg_zipdown)
+        for i in 1:(N - 1)
+            ψ, ϵᵢ = left_gauge!(ψ, i, ψ.AC[i], alg_zipdown)
+            ϵ = max(ϵ, ϵᵢ)
+        end
+    end
+
+    return ψ, ϵ
+end
+
+# `right_gauge` for the tensors of a right-to-left zip-up sweep: these are already partitioned across
+# the new bond, so the leg permutation that `right_gauge` applies to MPS tensors has to be skipped
+function _right_gauge_zip(Aᶻ, alg::MatrixAlgebraKit.TruncatedAlgorithm)
+    U, S, Vᴴ, ϵ = svd_trunc(Aᶻ, alg)
+    C = LinearAlgebra.rmul!(U, S) # C = U * S, matching `RightOrthViaSVD`
+    return C, _transpose_front(Vᴴ), ϵ
+end

src/operators/mpo.jl

@@ -278,10 +278,50 @@ function Base.:*(mpo::InfiniteMPO, mps::InfiniteMPS)
     return changebonds(InfiniteMPS(As), SvdCut(; trunc = notrunc()))
 end
 
+"""
+Fuse the left and right virtual legs of the product of MPO-MPS tensors
+```
+        ┌---A---┐
+    1 --Fl  |  Fr-- 3   =>  A′[1 2; 3]
+        └---O---┘
+            |
+            2
+```
+"""
 function _fuse_mpo_mps(O::MPOTensor, A::MPSTensor, Fₗ, Fᵣ)
     @plansor A′[-1 -2; -3] := Fₗ[-1; 1 3] * A[1 2; 4] * O[3 -2; 2 5] * conj(Fᵣ[-3; 4 5])
     return A′ isa AbstractBlockTensorMap ? TensorMap(A′) : A′
 end
+"""
+Fuse the left virtual legs of the product of MPO-MPS tensors
+```
+        ┌---A--- 3
+    1 --Fl  |       =>  A′[1 2; 3 4]
+        └---O--- 4
+            |
+            2
+```
+"""
+function _fuse_mpo_mps_left(O::MPOTensor, A::MPSTensor, Fₗ)
+    @plansor A′[-1 -2; -3 -4] := Fₗ[-1; 1 3] * A[1 2; -3] * O[3 -2; 2 -4]
+    return A′ isa AbstractBlockTensorMap ? TensorMap(A′) : A′
+end
+"""
+Fuse the right virtual legs of the product of MPO-MPS tensors
+```
+    1 --A---┐
+        |  Fr-- 3   =>  A′[1 2; 3 4]
+    2 --O---┘
+        |
+        4
+```
+"""
+function _fuse_mpo_mps_right(O::MPOTensor, A::MPSTensor, Fᵣ)
+    # the resulting tensor is partitioned across the new bond
+    # `_transpose_front` of the right factor is again an MPS tensor
+    @plansor A′[-1 -2; -3 -4] := A[-1 1; 2] * O[-2 -4; 1 3] * Fᵣ[2 3; -3]
+    return A′ isa AbstractBlockTensorMap ? TensorMap(A′) : A′
+end
 
 function Base.:*(mpo::FiniteMPO{<:MPOTensor}, x::AbstractTensorMap)
     @assert length(mpo) > 1

src/states/abstractmps.jl

@@ -169,6 +169,11 @@ physicalspace(O::MPOTensor) = space(O, 2)
 physicalspace(O::AbstractBlockTensorMap{<:Any, <:Any, 2, 2}) = only(space(O, 2))
 physicalspace(ψ::AbstractMPS) = map(Base.Fix1(physicalspace, ψ), eachsite(ψ))
 
+# the input physical space of an MPO tensor, i.e. the space of the state it can be applied to.
+# This need not equal its output physical space `physicalspace`.
+_input_physicalspace(O::MPOTensor) = dual(space(O, 3))
+_input_physicalspace(O::AbstractBlockTensorMap{<:Any, <:Any, 2, 2}) = dual(only(space(O, 3)))
+
 """
     eachsite(state::AbstractMPS)