[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?)