[WIP] Boundary-MPS approximate contraction of an iPEPO window with CTMRG environment

[Yue-Zhengyuan]

This PR uses the finite-MPO/finite-MPS zip-up compression introduced in QuantumKitHub/MPSKit.jl#470.

Summary

This PR adds approximate expectation-value contractions for a single-layer InfinitePEPO using finite boundary MPSs and a CTMRG environment:

• MPOObservable represents a physical open-boundary MPO acting along a path in the 2D network.
• expectation_value_approx contracts an MPOObservable in its enclosing rectangular window.
• correlator_approx measures a dense two-site operator between one fixed first site and many second sites, reusing partial contractions with adaptive window depth and fixed width.

Each PEPO row becomes a finite MPO, which is applied to a boundary MPS with zip-up compression and optional one-site DMRG refinement.
The implementation includes virtual-space flips and twists, with U(1) and FermionParity regression coverage.

Design

Path-based MPO observables

MPOObservable stores:

• sites: lattice sites on which physical MPO tensors act;
• mpo: the open-boundary MPO tensors, with mpo[k] acting on sites[k];
• path: a non-self-intersecting nearest-neighbor path containing the operator sites in MPO order and any intermediate sites carrying the virtual string.

Constructors accept explicit MPO tensors, an explicitly routed path, or a dense AbstractTensorMap.
The dense constructor orders the sites and operator legs, decomposes the operator, and connects consecutive sites by horizontal-first paths.
Automatic routing is limited: it rejects paths that revisit a site or pass through another operator site out of order.
Explicit paths allow other valid routes, including paths that cross between two rows more than once at different columns.

The functions mpo_path_first, mpo_path_middle, mpo_path_last, and mpo_path_string insert the operator or route its string, trace the PEPO physical legs, and fuse the MPO string with the appropriate PEPO virtual legs.
This avoids materializing a separate BraidingTensor at each routing site.
Fusing these legs is not necessarily the most efficient contraction strategy, particularly at routing sites.

Finite-window expectation values

expectation_value_approx(rho, observable, env; trunc, maxiter = 1, direction = :auto)

For north-to-south contraction, CTMRG corners and north edges form the initial FiniteMPS.
Each row includes the west/east CTMRG edges and the traced or operator-modified PEPO tensors.
After zip-up compression and optional DMRG refinement, the final MPS is contracted against the south CTMRG boundary.
The expectation value is normalized by a separate observable-free contraction of the same window.

Column contraction reuses this implementation by rotating the state, environment, and observable.
Automatic direction selection uses rows for wide or square windows and columns for tall windows; callers can override it with :rows or :columns.
The south boundary is stored in the adjointed representation expected by the MPS contractions.

WindowApprox is an internal wrapper for the zip-up and optional DMRG algorithms.
For both public measurement functions, the default trunc limits the rank to the largest CTMRG boundary dimension, and maxiter = 0 disables DMRG refinement.

Two-site correlators with adaptive depth

correlator_approx(rho, op, i, js, env; trunc, maxiter = 1)

The interface follows correlator: i is a fixed CartesianIndex{2}, and js is a vector of Cartesian indices, a CartesianIndices collection, or a single Cartesian index.
A collection returns a vector in vec(js) order; a single second site returns a scalar.
Targets must be nonempty, unique, and distinct from i.
Operator leg 1 always acts at i, and leg 2 at each j.

The contraction direction is automatic, with no direction keyword:

• Row contraction is valid if all targets satisfy j[1] ≥ i[1].
• Column contraction is valid if all targets satisfy j[2] ≤ i[2].
• If both are valid, the rectangle enclosing all sites determines the choice: rows for wide or square rectangles, columns for tall rectangles.
• If neither is valid, the call raises an ArgumentError.

For row contraction, every window has the same column range, covering i and all targets, but ends at its own target row.
The operator string travels south along the source column and then horizontally along the target row.
Column contraction uses the rotated equivalent: a fixed row range and a window ending at each target column.

A single forward loop maintains two north states: one without observables for normalization, and one carrying the source insertion and open MPO string.
At each target row, the algorithm constructs the CTMRG south boundary directly below it, closes the targets using shared horizontal environments, and normalizes them with the observable-free contraction using the same row closure.
The target row is closed without further truncation in both numerator and denominator.
This avoids mixing a target-row numerator with a differently truncated full-sweep denominator, which can otherwise produce identity expectations far from one.

The implementation reuses the operator decomposition, forward propagation, and horizontal contractions within each target row.
It does not store previous row MPOs or boundary histories and requires no reverse sweep or adjoint row MPOs.
Adding deeper targets preserves earlier results as long as the fixed width and automatically selected direction remain unchanged.

Validation

Tests cover independent adaptive-window references, target ordering, column rotation on rectangular unit cells, propagation through rows without targets, and preservation of earlier values when the depth is extended.
Truncated identity checks cover U(1) and FermionParity, both contraction directions, and optional DMRG refinement.
A separate physical-state regression compares against expectation_value.
MPO routing and finite-window tests also cover turns, reversed ordering, and a U-shaped path.

The two correlator test files passed locally with 27 assertions, run serially with --jobs=1.

Remaining work

• Generalize the contraction backend to iPEPS and purified iPEPO networks.
• Improve automatic path routing and explore alternatives to fusing MPO strings into PEPO virtual legs.

[codecov]

Codecov Report

❌ Patch coverage is 12.92035% with 492 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
...rc/algorithms/contractions/mpo_path/pepo_1layer.jl	0.00%	115 Missing ⚠️
...orithms/contractions/window/twosite/pepo_1layer.jl	0.00%	106 Missing ⚠️
src/algorithms/contractions/mpo_path/routing.jl	0.00%	105 Missing ⚠️
src/algorithms/contractions/window/pepo_1layer.jl	33.33%	54 Missing ⚠️
src/algorithms/contractions/transfer.jl	0.00%	41 Missing ⚠️
src/algorithms/expval_approx.jl	0.00%	30 Missing ⚠️
src/algorithms/contractions/window/tools.jl	0.00%	23 Missing ⚠️
src/algorithms/correlator_approx.jl	0.00%	18 Missing ⚠️

Files with missing lines	Coverage Δ
src/PEPSKit.jl	100.00% <ø> (ø)
src/operators/localoperator.jl	89.11% <100.00%> (+15.72%)	⬆️
src/utility/indexing.jl	92.50% <ø> (+2.50%)	⬆️
src/algorithms/correlator_approx.jl	0.00% <0.00%> (ø)
src/algorithms/contractions/window/tools.jl	0.00% <0.00%> (ø)
src/algorithms/expval_approx.jl	0.00% <0.00%> (ø)
src/algorithms/contractions/transfer.jl	49.43% <0.00%> (-42.23%)	⬇️
src/algorithms/contractions/window/pepo_1layer.jl	33.33% <33.33%> (ø)
src/algorithms/contractions/mpo_path/routing.jl	0.00% <0.00%> (ø)
...orithms/contractions/window/twosite/pepo_1layer.jl	0.00% <0.00%> (ø)
... and 1 more

... and 13 files 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]

@lkdvos @leburgel Does the MPSKit support for fermions (via planar contraction) still work if the MPS/MPO arrows do not follow the standard directions?

[leburgel]

@lkdvos @leburgel Does the MPSKit support for fermions (via planar contraction) still work if the MPS/MPO arrows do not follow the standard directions?

Yes, I think it should just work regardless of which spaces exactly are dual.

[Yue-Zhengyuan]

Fermion support is now here, but the way to achieve it is not ideal. I first put the virtual arrow directions in the iPEPO in the standard direction. But the on the right boundary, the "physical space" is dual, as required by the CTMRGEnv constructor (see C₂ below):

    [1; 2]      [1 2; 3]        [1; 2]
    C₁-←-2      1-←-E₁-←-3      1-←-C₂
    ↓               ↓               ↑
    1               2               2

Thus I still need to use many planar operations (some of which I still don't fully understand) when constructing the left and right ends of each row, and the south boundary.

(Edit: now I also flip the east envspace to standardize all arrows. Things then make more sense.)

Otherwise, if I dispatch away from planar operations, it seems that I need to overload lots of things in MPSKit (not only in ZipUp but also DMRG for further refinement), which may not be worthwhile. (But anyway we still need to do it when generalizing to approximate contraction of multi-layer networks?)

[leburgel]

I was going to take some time to go over this, but it seems this PR is rather large and I completely misjudged the time I allocated. I'll try to go over the details as soon as I can, but maybe I can start with a few general remarks already before I submit a full review:

---

As far as I can tell, the newly introduced MPOObservable could in spirit serve as a term in a LocalOperator, where instead of associating a dense TensorMap to a group of sites we associate a finite MPO (encoded here as a vector of MPO tensors). The only difference is the path field that is additionally included in the MPOObservable struct. This leads me to the question: what is the benefit of being able to manually specify a path, over just determining this automatically? Also, what exactly is the difficulty in this automatic path finding (as indicated in the note above)? If it's not more expensive to always use a default path, then we can directly use MPOs as LocalOperator terms, whose expectation value we could then evaluate exactly or approximately. This would be very nice to have, hence my question.

---

For the finite window contraction, would it be possible to directly pass an algorithm for this? I can imagine wanting to just use DMRG2 (or even just DMRG) to do this, rather than starting from ZipUp and then optionally refining. Probably the way it is now is fine for practical purposes, at least the ones you had in mind, but it feels like it would be useful to have direct control over the algorithm.

---

Regarding fusing the MPO path, instead of fusing everything into the virtual spaces, we could (at some point) try to keep the virtual MPO legs around explicitly as auxiliary legs for MPS tensors when they are encountered. The problem is then to contract two "matching" auxiliary legs corresponding to the same virtual MPO bond when they encounter each other in two tensors being contracted. This is a lot more tedious to do in our framework, since we don't have "labeled" tensor indices that recognize each other automatically. Here we can just pass around explicit labels for all of the auxilary legs that come from a virtual MPO bond, and keep track of them. Not sure how feasible this is, and definitely not sure if this is something we want to do here, but just wanted to bring it up already.

---

Regarding fermion support, I'm not sure if I'm entirely following the discussion. I think the MPSKit routines should just be able to handle physical spaces with a "non-standard" duality. Using planar operations where possible should ensure this just works, in the way that is consistent with MPSKit conventions. Therefore, I don't really see the need to consider dispatching away from planar operations. The only reason for this as far as I'm concerned would be to ensure differentiability, which is not a goal at all for now I think. So even flipping this east env space shouldn't be necessary, I think.

The only reason to not use planar operations should be if we simply can't, e.g. for multi-layer networks. There, we have to overload the core contractions used in MPSKit anyway, but we additionally have to put additional twists to ensure the result is what is expected by MPSKit conventions. For the MPSKit contractions already overloaded in PEPSKit, this should have been done in a way that they can also handle arbitrary arrow directions.

---

For multiplying in the bottom boundary, I would really like to avoid (doubly) conjugating edge tensors, this always leads to issues sooner or later in my experience. The question then becomes if we need to twist anything on the resulting top MPS (after all of the MPO applications) to ensure we can just contract in south edges in the usual way and actually get what we want. I'm not entire sure about that.

[Yue-Zhengyuan]

In principle I should split this PR in two: first focusing only on expectation_value_approx, with a follow up with correlator_approx. But these two share much infrastructure, so I put them together...

• Path routing for MPO terms

I agree that needing an additional path to specify the MPO is not ideal. The difficulty is in writing a clever algorithm that creates a path that go through the MPO sites in order, while also being non-self-intersecting. You can consider the sites (1, 1) - (6, 1) - (2, 0) - (3, 2) - (4, 0) - (5, 2) (which is unlikely to appear in practice, though). It is not obvious to me which should be the "canonical" path joining them without self-intersection.

• Algorithm struct for window contraction

It is already there, called WindowApprox. The problem with skipping ZipUp is how you initialize a ψ0 for DMRG/DMRG2 without randomly guessing.

Currently this struct is not exported. Instead, we provide kwargs trunc and maxiter in the API to build the struct internally:

WindowApprox(Zipup(; trunc), _approx_dmrg(maxiter))

• Fusing MPO virtual space with the PEPO

I want to postpone this. Currently ZipUp (and DMRG?) only support(s) MPS with one physical leg per site.

• Fermion support

Changing arrow directions in MPSKit is like planar contraction with flippers without twists. But this is not the case in PEPSKit, where we should use flip that preserves @tensor calls. First flipping all arrows to the standard direction appears to automatically put the required twists back in the network; the change of arrows is more like a by-product that is unimportant to MPSKit. I suggest we keep this pre-processing for now, and only revisit it at the end.

• Conjugation of south boundary

In the cache WindowRowCache used by correlator_approx, we not only need the south CTM boundary, but also need the contraction of the south boundary with some rows above it. So avoiding conjugating south boundary twice will require more than a custom version of dot. We will also need to write ψ * O, where the MPO O acts on the bra-state ψ. This is a bit redundant to me.

[leburgel]

Sorry it took me so long to get through this, I couldn't seem to find the time to sit down and get through the whole thing in one go.

I have a first round of comments, most of which are related to the general remarks I had a while back. An overview:

• Path routing for MPO terms: I think we can make this problem much easier by giving up on ordering the sites according to the matrix-like linear ordering we use for dense local operators. There the choice doesn't matter much so imposing a canonical ordering makes sense. But for MPOs, choosing the site ordering in a way that makes the routing easy makes much more sense to me.

• Algorithm struct for window contraction: I left a more in-detail comment below, but I think that if we use DMRG with an expansion we can merge things quite cleanly actually.

• Fusing MPO virtual space with the PEPO: definitely fair, fusing is the quickest way to actually get something up and running.

• Fermion support: as always I am still confused by this. If it's verified to work properly, there's not many objections I can raise, so that would be fine for me. We can always revisit how we do things later, but I don't have any conceptual arguments right now.

• Conjugation of the south boundary: again there's a more detailed comment below, but I think we can avoid conjugations by rotating instead. This comes with its own issues, but I think it's conceptually cleaner than taking adjoints. To implement ψ * O, we would just rotate O, compute the normal thing, and take care to reverse ψ back and forth in the appropriate way. This is the same way I would do infinite boundary MPS contractions, so I think that could be a viable alternative.

[leburgel]

I would personally implement this through rotation without conjugation, which reverses the left-right ordering of the south edge tensors and collects them into a finite MPS. This respects the natural virtual orientation of the south edge tensors, and avoids having to take the "adjoint" of individual MPS tensors.

[Yue-Zhengyuan]

I think we still need to handle the physical legs if we use the rotation approach? South boundary will have dual physical spaces.

[leburgel]

I somehow thought dual physical spaces weren't really an issue for MPSKit, I though we only made hard assumptions on the virtual space duality. I might be wrong or be misremembering, I just thought I tried this out at some point a long time ago.

[leburgel]

The main reason why I originally asked if we could get rid of this path field is that this would then just become a specific term in a LocalOperator. Conceptually, I quite like just being able to call expectation_value_approx on any LocalOperator, where we then evaluate every term individually in the appropriate way (using some kind of local_expectation_value_aproximate hook).

For dense tensors, we could just convert them to an MPO when calling the approximate version. The path routing could then happen inside this local expectation value computation, concerning only the relevant path.

The non-intersecting path routing problem is indeed hard in general, if we're handed an MPO acting on specific sites. But usually, we would want to evaluate some operator which we can represent as a dense tensor map (I assume this is the most common case?), which we then decompose into an MPO and evaluate approximately. If we choose the site order in this decomposition in a way that takes into account the routing that we'll have to do, this makes the problem much easier. So I think we could get rid of the routing difficulties by simply allowing sites that are not ordered according to the lattice ordering. This anyway makes a lot of sense for MPO observables I think, where intuitively we would order the sites to end up with the virtual string of smallest length, which will generally be quite easy to route.

Does this make sense, and sound like a reasonable thing to try?

In the end I would like to not do any path routing and virtual fusing at all, and just pass around the MPO indices appropriately automatically. This is not at all necessary now, but it would be good to already get ahead of the routing issue.

[Yue-Zhengyuan]

For expectation_value, we may add an additional algorithm struct like ExactContract for the existing exact contraction, and dispatch to boundary MPS approximation using the contraction algorithm:

expectation_value(state, O, env, [alg = ExactContract()])
expectation_value(state, O, env, alg = WindowApprox(zipup, dmrg))

So extra names like expectation_value_approx will be unnecessary.

correlator_approx may still be kept, since it accepts a list of 2-site bonds instead a fixed first site + a list of second sites. (But I may also make the transition?)

I think even in common cases, the MPO is not always first constructed as a dense map? (e.g. some are just a TensorProductTerm) And it is not easy either to come up with a algorithm that reorder the sites acted on by the operator so that the string length is minimized.

[leburgel]

If we let go of the operator site ordering being the same as for dense tensor local terms, I think we could simplify things quite a bit?

[Yue-Zhengyuan]

Major changes made today (mainly on correlator_approx):

• Changed correlator_approx to accept one fixed site i and second sites js, matching correlator. Direction is selected automatically from the target positions; square windows prefer row-by-row contraction.
• Introduced adaptive window depth with fixed width: each contraction ends at its target row or column.
• Removed WindowRowCache, the reverse sweep that needs _adjoint_mpo, and operator swapping.

[Yue-Zhengyuan]

Now the only big issue is designing the MPO terms. Ideally this should also include PBC-MPOs that form a closed loop. I'll take a closer look at #425 and see if I can get some inspirations.