Zip-up compression for FiniteMPO-FiniteMPS product

[Yue-Zhengyuan]

This PR adds approximate((O, ψ), ::Zipup) for open-boundary FiniteMPO acting on FiniteMPS, without normalization, first introduced in 1002.1305:

approximate((O, ψ)::Tuple{Any, <:FiniteMPS}, alg::Zipup)

This provides a streaming MPO-MPS application plus SVD compression path. The algorithm struct Zipup controls the SVD algorithm and truncation strategy.

Motivation

This is useful for PEPS boundary-MPS style contractions (QuantumKitHub/PEPSKit.jl#396), where applying a double-layer transfer MPO to a boundary MPS can temporarily create very large virtual bonds. The existing initialization pattern,

changebonds(O * ψ, SvdCut(; trscheme); normalize=false)

first materializes the full enlarged O * ψ across all sites, then compresses it. That is unnecessarily memory-heavy when the product is only needed as an initialization that will be further optimized with approximate.

The new Zipup algorithm contracts one MPO/MPS site pair at a time from right to left, absorbs the "carry" (residue tensor produced by right-orth gauge conversion) from the already-compressed right block, and immediately truncates before moving left. This avoids storing the fully enlarged product on every site.

Notes

The Zipup approximation matches changebonds(O * ψ, SvdCut(; trscheme); normalize=false) when no meaningful singular values are discarded, but can differ under real truncation because Zipup truncates while building the product rather than after materializing the full product.

[Yue-Zhengyuan]

Please let me know if there are actually better functions in MPSKit that does more or less the same thing.

[Yue-Zhengyuan]

@lkdvos I cannot run the format check on GitHub either. So there's something to be fixed.

[lkdvos]

Thanks for getting this started! I'll look into the formatter, let's not worry about that here.

For consistency, it would probably make the most sense to match this into approximate!(out, (O, in), alg, envs = nothing), but I guess that stumbles back on the issue we had before where it is not obvious how to construct the out beforehand. It might be reasonable to simply allow approximate((O, in), alg), and construct a FiniteMPS without any tensors initialized, so just the exterior vectors get reused?

I'm a bit confused about the necessity for GenericMPO.
Usually when we have two legs, the main point is that the MPO is structured in such a way that it acts only on one of the MPS legs, by keeping the layers separate. This then leads to memory and computational benefits, as opposed to just fusing everything together.
Here, there isn't really a lot of benefit to support the GenericMPO, since it requires us to maintain more contractions, and is generally expected to be slightly slower because the fusiontree overhead grows with the number of legs.

As an alternative, I would suggest writing this in terms of primitive pieces that can be overloaded for the PEPS case, which keep the PEPS layers separate to keep the computational benefits in. Here that would probably just be fuse_mpo_mps?

[codecov]

Codecov Report

❌ Patch coverage is 91.66667% with 2 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/algorithms/approximate/zipup.jl	91.66%	2 Missing ⚠️

Files with missing lines	Coverage Δ
src/MPSKit.jl	100.00% <ø> (ø)
src/algorithms/approximate/approximate.jl	56.00% <ø> (ø)
src/algorithms/approximate/zipup.jl	91.66% <91.66%> (ø)

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[Yue-Zhengyuan]

It might be reasonable to simply allow approximate((O, in), alg)

Ah, the same interface debate as compress in QuantumKitHub/PEPSKit.jl#311 ... For alg, should I introduce a new Zipper (something like that) instead of reusing SvdCut? I need to distinguish it from changebonds(O * in, alg).

construct a FiniteMPS without any tensors initialized, so just the exterior vectors get reused?

What is meant by "exterior vectors"?

I'm a bit confused about the necessity for GenericMPO.

Indeed I realized this later. We just overload fuse_mpo_mps for PEPS/PEPO applications. Here in MPSKit I should just write for the simplest one-physical-leg case.

[lkdvos]

What is meant by "exterior vectors"?

Here I mean that a FiniteMPS is at its core just a set of vectors (ALs, ARs, ACs, Cs), which are each filled with tensors. I was thinking whether it could be worth it to just initialize the vectors as undef (or perhaps just missing), so no tensors are present but we have something to instantiate into. This bypasses the need to allocate output tensors in advance, but still keeps the same interface?
Let me know if that sounds too far-fetched though, I don't want to overcomplicate things.

Should I introduce a new Zipper (something like that) instead of reusing SvdCut?

I think this is probably sensible, it might additionally store something like a direction or a backend and allocator in the future.

[Yue-Zhengyuan]

Initializing with undef just for matching interface is a bit unnecessary for me. What about also calling the function compress? Might be more descriptive than the current zipper.

[lkdvos]

It's just really annoying to have both approximate and compress for the same thing though, just based on the distinction of not needing an initial guess, in that case I would much rather make the interface of approximate more adapted to that

[Yue-Zhengyuan]

I feel that the default _fuse_mpo_mps, which uses @plansor, will cause troubles for iPEPO measurements. For column-by-column contraction from right to left, the default iPEPO arrow convention (top to bottom) will make the MPO virtual spaces dual.

[lkdvos]

I'm not sure if that really is a problem, since you will presumably have to overload that anyways?

[Yue-Zhengyuan]

The case of one-layer iPEPO tr(ρO) actually matches the current signature, so it is a bit difficult to overload.

[lkdvos]

I'm not sure if that really is a problem, since you will presumably have to overload that anyways?

[lkdvos]

Is there a reason we need this copy here? And if so, could you restructure it to be approximate((O, psi), alg) = approximate(copy(psi), (O, psi), alg)?

[Yue-Zhengyuan]

I was being over-cautious here. I'll remove the copy.

[VinceNeede]

Thanks for putting this together — the zip-up contraction is a nice addition. A few questions/comments while reading through it, cross-checking against the original references:

1. Naming

The algorithm is implemented here as Zipper, but the original method (Stoudenmire & White, New J. Phys. 12, 055026 (2010), DOI: 10.1088/1367-2630/12/5/055026 calls it "zip-up." Was the rename intentional (e.g. to avoid a name clash, or because the implementation deviates from the original enough to warrant a different name), or would it be worth aligning the naming with the literature so people searching for "zip-up" can find it?

2. Missing truncation-only sweep

As I understand the original algorithm, it consists of two sweeps: one sweep that does the contraction with only a loose/minimal truncation (to keep the intermediate bond dimension under control without discarding meaningful weight), followed by a second sweep that performs the actual truncation down to the target maxdim/cutoff. In this PR I only see the first sweep (compression with minimal truncation) — is the second, final-truncation sweep intentionally left out for now (e.g. deferred to a follow-up PR), or should it be included here?

Related to this: Paeckel et al. (Annals of Physics 411, 167998 (2019), DOI: 10.1016/j.aop.2019.167998 give a concrete prescription for the intermediate/minimal truncation stage — using a looser bond dimension (2× the target maxdim) and a looser cutoff (cutoff/10) relative to the final truncation values, I think this would be worth implementing as default behavior if the second sweep is added.

3. Orthogonality center of the MPO

Stoudenmire & White also move the orthogonality center of the Hamiltonian/MPO to the last site (given this is done right-to-left here). I don't see that happening in this implementation — is the MPO's orthogonality center supposed to be moved before by the user, or is this something that still needs to be added?

Happy to be corrected on any of this if I'm misreading the code or the references — just want to make sure the implementation matches the papers it's based on before it lands.

[Yue-Zhengyuan]

@VinceNeede Hope the following answers your questions.

Naming

I'm referring to the algorithm introduced in https://arxiv.org/abs/2310.08533 (Section IV), where it is indeed called "zipper". I have now renamed to Zipup. Thanks!

Missing truncation-only sweep

From you description, I think this is done by changebonds(O * psi, ::SvdCut), where O * psi does a slight truncation, and changebonds do the big truncation to desired dimension. Here zipper directly do the big truncation: as put by the original authors, "we abstain from applying the whole row of t’s at once".

Orthogonality center of the MPO

By construction the zipper algorithm directly produce an MPS with the orthogonality center at the first/last site (depending on whether you sweep from right to left, or left to right), so no need to move it again manually. Also, I think the orthogonality center is automatically handled by the FiniteMPS machinery, so users usually don't need to explicitly manipulate it.

[Yue-Zhengyuan]

Well, after some reading, the "zip-up" algorithm and the zipper seems to be the same thing... maybe I should cite the 2010 work instead?

[VinceNeede]

Thanks for the updates, glad the naming got sorted!

On point 2: You're right that changebonds!(::SvdCut) sweeps right-to-left, so with Zipup now going left-to-right, chaining the two does reproduce the paper's two-sweep structure directionally. My concern is less about whether it's possible to chain them, and more that a user calling Zipup alone gets only the first sweep, with no indication that a second call is needed to match the published algorithm. Since this only works correctly if the two calls use matched, deliberately different truncation parameters (Paeckel et al., loose 2×maxdim/cutoff/10 on the first pass, target values on the second), I think this is worth surfacing rather than leaving implicit — either by adding the second sweep as an opt-in inside Zipup itself, or at minimum documenting the exact chain (approximate with loosened trscheme, then changebonds(..., SvdCut) with the final one) directly in the docstring, so the two-sweep behavior isn't something a user has to reconstruct from the original paper.

On point 3 (orthogonality center): To clarify, I meant the gauge of the MPO, not the resulting MPS. In the original algorithm, this gauge-fixing step isn't primarily about canonical form for its own sake — it's there to control conditioning. Quoting Stoudenmire & White directly: the arrangement guarantees that the basis to the right of the current site, from the product of the MPO and the MPS, is "not drastically ill-conditioned," meaning no basis state ends up with a norm much bigger than one.

I don't think the Jordan-block structure of a Hamiltonian MPO gets you this for free. The triangular structure only constrains the flow of information between channels (no cycling back to an earlier block) — it says nothing about the magnitude of the entries. For short-range Hamiltonians with O(1) couplings, the channels are probably all reasonably scaled, so this may not show up in practice. But for long-range or multi-scale Hamiltonians (e.g. power-law couplings, or terms spanning very different energy scales), the "already-closed" channels can accumulate sums of very different magnitude as the sweep progresses — exactly the kind of basis ill-conditioning the gauge-fixing step is meant to prevent. So I think this is a real, if perhaps rare in typical test cases, gap rather than something the Jordan form makes moot.

That said, I don't know whether MPSKit has any existing notion of a gauge/orthogonality center for FiniteMPO at all (independent of the Jordan-form Hamiltonian case) — is there such machinery, or would this be genuinely new territory? If it's out of scope for this PR, a docstring note that the MPO's own conditioning isn't handled by Zipup would at least flag it for anyone applying this to long-range or multi-scale Hamiltonians.

[Yue-Zhengyuan]

I would refrain from incorporating SvdCut changebonds into the zip-up algorithm to keep things focused. Even in the original paper the second step is referred to as another "fitting algorithm".

For gauging MPO, I guess you can always convert it to an MPS with 2 physical legs per site. FiniteMPO appears to be just a thin wrapper over the vector of MPO tensors, so it does not have the complete machinery to handle the gauges like FiniteMPS.

[VinceNeede]

One more clarification on point 2, since I think we might be talking about two different things: the "fitting algorithm" the paper refers to (e.g. a DMRG-style variational sweep) is a separate, optional refinement that operates once you're already at the target bond dimension — it improves the overlap with the exact result but doesn't change the bond dimension itself.

The two-sweep structure I meant is different: it's how zip-up itself reaches the target bond dimension without losing information along the way. The intermediate sweep has to use a loose truncation precisely because it isn't reaching the target bond dimension yet — cutting straight to the target on a single pass risks discarding weight that the second sweep would have recovered. The second sweep is the one actually responsible for reaching the desired final bond dimension, using the better-conditioned information available after the first pass. So this two-sweep step is intrinsic to the zip-up compression itself, and only after it's done would a separate fitting algorithm optionally be applied. Just wanted to make that distinction explicit, since I think that's where the disagreement comes from.

On the MPO point: given FiniteMPO doesn't currently have the gauge machinery to support this, would it be reasonable to at least add a note in the docstring that the MPO's own conditioning isn't handled by Zipup? That way anyone applying this to long-range or multi-scale Hamiltonians, where it's more likely to matter, knows it's a known limitation rather than an oversight.

[Yue-Zhengyuan]

What would you suggest about the interface if we do two successive truncations? Do we pass two alg_gauges to Zipup?

[VinceNeede]

For the interface, I'd suggest two separate keyword arguments, something like alg_gauge for the intermediate sweep and alg_gauge_final for the final truncation, (or, construct the alg internally and expose trscheme and trscheme_final as keywords) with the following fallback behavior:

• If only alg_gauge_final is supplied, alg_gauge is constructed automatically from it using the Paeckel prescription (2 * maxdim, cutoff / 10).
• If only alg_gauge is supplied and alg_gauge_final is not, no final truncation sweep is performed (preserving the current single-sweep behavior for anyone who wants it).

This also matches what ITensor (maintained by Stoudenmire and White themselves) does in practice — their contract docs for alg="zipup" describe exactly this two-stage structure: an initial sweep with loose truncation (not locally optimal), followed by a call to truncate! with the actual target parameters, explicitly citing the same Paeckel et al. prescription for the loosened intermediate cutoff/maxdim. So this interface would bring Zipup in line with both the reference paper and the existing reference implementation.

[Yue-Zhengyuan]

ITensor reference is very helpful! I'll soon update to match the behavior there.

Not all truncation strategies have a dimension or error cutoff. So for automatic determination of alg_gauge for the first sweep, I will make an external constructor like Zipper(; D, atol, rtol, two_sweep::Bool).

[lkdvos]

Let me try and weigh in on this, although I might have missed some pieces:

I have no preference about the name, feel free to pick whichever, we can refer to both in the documentation.

Considering the gauge, and how MPSKit handles this, and what it means: In the state-truncation case, we want to truncate in the center gauge, because that way the truncation of the local tensor is provably the optimal truncation that best approximates the original state. This optimality is slightly lost when you consider doing a truncation on all sites, as these things don't commute and earlier truncations might affect later ones. This can be mediated by afterwards doing a sweeping algorithm that improves on this, although typically the SVD sweep is already quite good, because at least the local truncations are optimal.

Going to MPOMPS compression schemes, things get a little bit more difficult. In principle, in order to do an optimal truncation, you have to do this in the center gauge, but to do this you effectively have to do the full mpomps contraction without truncating, then bring this in canonical gauge, and then truncate. (+ sweeps afterwards if you want to improve the accuracy). Of course, this is quite expensive, so people have a lot of different ways and algorithms that try to improve on this, with various benefits and drawbacks.

The zipper algorithm variants are like you implemented here and like the references, namely that simply assume the local gauge of the MPO*MPS contraction is not too far wrong (for example when the MPO is close to the identity), and if we don't have to truncate too heavily this is a much cheaper approach. This can then be slightly improved by indeed doing a first sweep that only lightly truncates (in the wrong gauge, but with minor truncation), followed by a reverse sweep that can truncate more heavily (since we can now bring this in the correct gauge). While in principle we could do this by just chaining this with a changebonds(::SvdCut) call, I think here it is quite natural and it would indeed be convenient to bundle these two options (although we should probably look into sharing some of the code).

Considering the gauge of the MPO, this is a completely different story. While there is functionality to bring this into a gauged form, I am not sure you can prove that this is helpful, since that only tells us the optimal way to truncate the MPO by itself (in the frobenius norm, not the operator norm), which doesn't really say anything about the gauge of MPO*MPS. There is probably something to be said about doing at least some kind of gauging procedure on the MPO to avoid pathological cases, but I would argue that this requires quite a lot of information about what kind of MPO it is, so I'd rather not have that hardcoded into this algorithm.

A final note here is that (at least according to me) the main spirit of the Zipper algorithm is that it is fast, and meant for cases where the optimal accuracy might not be the most important. In that sense, I'd argue that doing too much work trying to make everything optimal (such as gauging the MPO) is somewhat against that spirit, because if this is really required, you're probably better off with following up with a sweeping algorithm instead.

---

A small sidenote, @Yue-Zhengyuan is right that MPSKit already handles gauging of the states: calling mps.AC[i] will bring the ith site to the center gauge and give you that tensor.

[VinceNeede]

Thanks for the detailed explanation, the distinction between optimality in the Frobenius norm (for the MPO alone) versus the operator norm (which is what actually matters for MPO*MPS) makes sense, and I agree that's a good reason not to hardcode MPO gauging into this algorithm without knowing more about the specific MPO structure. Glad the two-sweep truncation will be bundled in, happy to leave the MPO conditioning question as a scope call for now, with a sweeping algorithm as the natural fallback if it turns out to matter in practice.