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

[Yue-Zhengyuan]

(This PR requires QuantumKitHub/MPSKit.jl#470 which has not been released yet. @lkdvos)

Summary

This PR adds approximate expectation-value contractions using boundary MPS method for finite windows of a single-layer InfinitePEPO:

• MPOObservable represents a physical OBC-MPO acting on a path of the 2D network.
• expectation_value_approx measures the expectation value of MPOObservable.
• correlator_approx specializes on 2-site correlators on many bonds, with a caching mechanism to reuse intermediate contraction results.

The implementation converts CTMRG boundaries to finite MPSs, applies PEPO rows as finite MPOs with MPSKit's zip-up contraction, optionally refines each MPO-MPS product with one-site DMRG, and closes the window against the opposite CTMRG boundary.

Design

Note

This is definitely not finalized, and up to further improvements.

Path-based MPO observables

MPOObservable stores three aligned pieces of information:

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

The purpose of path is to avoid materializing a BraidingTensor at each of the routing sites.

An MPOObservable can be built from explicit MPO tensors, from an manually routed path, or from a dense AbstractTensorMap. The dense constructor will first order the sites to be acted, decompose the operator to an MPO, and automatically choose a path to connect the sites.

Note

Currently the automatic path routing does not work for all sites.

Fusing an MPO path into a PEPO

The functions mpo_path_first, mpo_path_middle, mpo_path_last, and mpo_path_string applies the physical MPO on the iPEPO, and then trace out the physical legs. The virtual strings of the MPO are fused with the iPEPO tensor virtual legs. It is assumed that each site can only be passed by the MPO once, hence the requirement that the MPO path is not self-intersecting.

Note

Fusing the MPO strings with the iPEPO virtual legs is not the optimal way to contract, especially for the routing sites.

Finite-window contraction

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

For row-by-row contraction of a rectangular window in the iPEPO from north to south, expectation_value_approx proceeds in the following steps.

1. builds a north FiniteMPS from CTMRG corners and north edges;
2. builds each FiniteMPO row from the west/east CTMRG edges and traced or observable-modified PEPO tensors;
3. applies the row with approximate((W, psi), Zipup(...));
4. optionally refines that result with one-site DMRG;
5. contracts the final state with a prepared south-boundary MPS using dot.

The south tensors are conjugated and permuted when the boundary MPS is built, in order to cancel the conjugation applied to the first argument of dot.

Column sweeps reuse the same backend by rotating the PEPO, CTMRG environment, coordinates, and observable. With direction = :auto, wide windows use row sweeps and tall or square windows use rotated column sweeps.

WindowApprox is an internal, unexported wrapper for the zip-up and optional DMRG algorithms. Public callers only choose trunc, maxiter, and direction. By default, trunc is truncrank(chi), where chi is the largest CTMRG corner-space dimension; maxiter = 0 disables DMRG refinement.

Two-site correlators with caching mechanism

correlator_approx(
    rho::InfinitePEPO, op::AbstractTensorMap, bonds::AbstractVector, env::CTMRGEnv;
    trunc, maxiter = 1, direction = :auto,
)

correlator_approx is specialized to measure the same two-site operator op on many bonds (using one window that covers all bonds, even if some bonds can be covered with a smaller window) with a caching mechanism to reduce repeated contractions. The cache stores:

• every ordinary row MPO in the shared window;
• "north prefixes" obtained by absorbing rows from north to south;
• "south suffixes" obtained by absorbing adjoint row MPOs from south to north;
• the common window norm.

bonds should not contain duplicates. First, each bond is ordered geometrically. After ordering, the first and the second site is referred to as the "source" and the "target", respectively.
A swapped flag records if the bond is reversed by the ordering.

Bonds with the same source and swapped flag are put in one evaluation group to share common contractions. For row-by-row contraction, correlator_approx proceeds as follows:

• closes all targets in the same row as the source using shared horizontal environments;
• routes the MPO string along the source column towards later rows;
• for later rows containing targets, closes all targets in them incrementally, similar to the 1-row correlator.

Thus rows before the source and after the target are reused through north/south caches, the vertical string is propagated once per source, and horizontal segments are shared between ordered targets.

TODO

The most important task is adding the support for fermions.

• Use EdgeTransferMatrix to avoid double conjugation of the south boundary and easier generalization to the case of iPEPS and purified iPEPO.
• Carefully check the twists that should be added to the fusers when fusing the MPO strings with the iPEPO.
• Dispatch away from the planar contractions used by MPSKit to properly add the twists for fermions.

[codecov]

Codecov Report

❌ Patch coverage is 23.12139% with 399 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
...orithms/contractions/window/twosite/pepo_1layer.jl	0.00%	125 Missing ⚠️
...rc/algorithms/contractions/mpo_path/pepo_1layer.jl	0.00%	110 Missing ⚠️
src/algorithms/contractions/window/pepo_1layer.jl	30.68%	61 Missing ⚠️
.../algorithms/contractions/window/twosite/caching.jl	0.00%	34 Missing ⚠️
src/algorithms/contractions/window/tools.jl	0.00%	27 Missing ⚠️
src/algorithms/correlator_approx.jl	0.00%	20 Missing ⚠️
src/algorithms/expval_approx.jl	0.00%	18 Missing ⚠️
src/operators/mpo_observable.jl	95.55%	4 Missing ⚠️

Files with missing lines	Coverage Δ
src/PEPSKit.jl	100.00% <ø> (ø)
src/operators/localoperator.jl	75.59% <100.00%> (+0.59%)	⬆️
src/operators/mpo_observable.jl	95.55% <95.55%> (ø)
src/algorithms/expval_approx.jl	0.00% <0.00%> (ø)
src/algorithms/correlator_approx.jl	0.00% <0.00%> (ø)
src/algorithms/contractions/window/tools.jl	0.00% <0.00%> (ø)
.../algorithms/contractions/window/twosite/caching.jl	0.00% <0.00%> (ø)
src/algorithms/contractions/window/pepo_1layer.jl	30.68% <30.68%> (ø)
...rc/algorithms/contractions/mpo_path/pepo_1layer.jl	0.00% <0.00%> (ø)
...orithms/contractions/window/twosite/pepo_1layer.jl	0.00% <0.00%> (ø)

... and 6 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.