Compress 2-layer InfinitePEPO with local truncation

[Yue-Zhengyuan]

This PR introduces the most naive algorithm to compress two 1-layer iPEPOs into a new 1-layer one with smaller virtual bond dimension:

compress((ρ1, ρ2), alg::LocalTruncation)

It is a polished version of @sanderdemeyer's work in ClusterExpansions.jl, with the extra support for arbitrary unit cell sizes (ClusterExpansions.jl focused on the trivial unit cell case). At least it can serve as initialization for other advanced approximation methods.

The algorithm minimizes the local cost function at each virtual bond, which is the norm of

        ↓ ╱      ↓ ╱            ↓ ╱                 ↓ ╱
    ----A2---←---B2---      ----A2-←-|╲       ╱|--←-B2---
      ╱ |      ╱ |            ╱ |    | ╲     ╱ |  ╱ |
        ↓        ↓       -      ↓    |P1├-←-┤P2|    ↓
        | ╱      | ╱            | ╱  | ╱     ╲ |    | ╱
    ----A1---←---B1---      ----A1-←-|╱       ╲|--←-B1---
      ╱ ↓      ╱ ↓            ╱ ↓                 ╱ ↓

The projectors are found following steps described in Appendix A, 1905.02351 (Boundary Tensor Renormalization Group). To handle fermions, the virtual arrows are required to follow the leftwards/downwards directions.

In addition, as suggested by @leburgel, I introduced lazy multiplication of two iPEPOs ρ1 * ρ2 by stacking layers of ρ1 below ρ2 without doing any contraction. Since we restrict all codomain physical spaces to be non-dual, we don't need to worry about fermionic twists here.

An immediate generalization is truncating the product of an InfinitePEPO and an InfinitePEPS, but currently I don't have concrete test cases.

[codecov]

Codecov Report

❌ Patch coverage is 98.50746% with 1 line in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/operators/infinitepepo.jl	90.00%	1 Missing ⚠️

Files with missing lines	Coverage Δ
src/algorithms/compress/local.jl	100.00% <100.00%> (ø)
src/states/infinitepartitionfunction.jl	75.00% <100.00%> (+0.42%)	⬆️
src/states/infinitepeps.jl	66.95% <100.00%> (+0.28%)	⬆️
src/utility/util.jl	68.81% <ø> (ø)
src/operators/infinitepepo.jl	77.96% <90.00%> (+1.11%)	⬆️

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

[sanderdemeyer]

Cool!

What you could do in the PEPS-PEPO case is to construct the PEPO as $$ exp(-d\tau H) $$, and see whether many applications of this PEPO on a random initial PEPS make it converge to the ground state. That is what I have as one of the tests in my ClusterExpansions code, which you can find here. This PEPO could be made using the Trotter decomposition, or via Simple Update, like you did in the previous tests. This should also allow for a fermionic test case.

[Yue-Zhengyuan]

@sanderdemeyer I don't know if you knew it before, but recently someone mentioned to me that we can optimize the HOTRG projectors (in your terminology, the "oblique projectors" in "no-environment truncation") with automatic differentiation, which is equivalent to second RG (see this paper). If you haven't tried it before, maybe I can help test this idea (also as an exercise to learn AD for me), which is also a great extension to this PR to actually make it appealing.

[Yue-Zhengyuan]

From the experience of 2508.05406 and 2602.22113, this simple truncation method can already yield some decent results. So we may just merge this after settling on some technical details like the interface.

[Yue-Zhengyuan]

@lkdvos Do you have time to take a look? I have changed from overloading approximate to changebonds. If this still looks weird, I can give up overloading MPSKit functions and just call it compress.

[lkdvos]

I'm still not entirely convinced by trying to fit it in changebonds, but I don't think it matters all that much to be honest. I think I would lean in the direction of just giving it a separate name, but if you prefer to keep it like this that definitely also works for me!

I left some other comments through the code that are mostly minor, otherwise I think this is probably fine!

[leburgel]

I like the rename to compress, this takes away all the annoyance that comes with trying to fit things into MPSKit signatures. My only real comment is on using the two PEPOs in the product as separate arguments, the simplest issue with this being that it's not very future-proof against any changes in the number of layers being approximated.

If there's no need for a more-layer version in the foreseeable future I would be fine with it as is, just wanted to bring it up.

[Yue-Zhengyuan]

I now think that ρ1 * ρ2 is not good. It is mainly because ρ1 * ρ2 will make a copy of the input in cat.

Some immediate generalizations are (ψ_bra, ψ_ket) (the bra-ket network for two iPEPS) and (ρ, ψ). They don't fit with * for iPEPO only.

We never compress more than 2 layers at once. Only the cost function may have more layers. So I guess even tuples are not necessary. @Adwait-Naravane if we ever want to write this, how should we pass the additional layers to the function? Maybe inside the truncation algorithm struct?

I remember testing the layer-wise QR approach myself, and found that it wasn't minimizing the cost function. But anyway, if it can reduce cost, I can also provide an option to use it.

BTW, QuantumKitHub/TNRKit.jl#172, which is doing essentially the same thing, took only one month to get merged. We spend too much time arguing about the interface.

[leburgel]

I now think that ρ1 * ρ2 is not good. It is mainly because ρ1 * ρ2 will make a copy of the input in cat.
Some immediate generalizations are (ψ_bra, ψ_ket) (the bra-ket network for two iPEPS) and (ρ, ψ). They don't fit with * for iPEPO only.

Sounds fair, I like the tuple option anyway.

We never compress more than 2 layers at once. Only the cost function may have more layers. So I guess even tuples are not necessary. @Adwait-Naravane if we ever want to write this, how should we pass the additional layers to the function? Maybe inside the truncation algorithm struct?

Even if you only compress two layers at a time, you can still pass all the layers in 1 tuple and then just apply the 2-layer compression method recursively adding one extra layer in every time. So maybe I'm misunderstanding, but I don't think there's a need to pass the additional layers elsewhere or add them to the truncation algorithm struct. If we ever want to do multi-layer compression we just add 1 additional method with the same interface, which dispatches to the ::Tuple{<:InfinitePEPO, <:InfinitePEPO} recursively internally. So this seems fine to me.

BTW, QuantumKitHub/TNRKit.jl#172, which is doing essentially the same thing, took only one month to get merged. We spend too much time arguing about the interface.

We already have the tendency to break things quite often, to the dismay of quite a few users, so I wouldn't say we necessarily spend too much time discussing interfaces. However, it is true that we're just too slow. Having to wait days or even weeks for a review is just not nice. Being annoying about interface choices while also not making time to quickly follow up on reviews is a bad combination. I will try to improve on the speed at least and try to make some more time, that I can definitely get behind.

[Yue-Zhengyuan]

I'm currently testing the layerwise-QR approach and see whether it really minimizes the cost. Let's aim for merging in this week.

[leburgel]

Final comments, but otherwise looks good to me now.

[Yue-Zhengyuan]

Just checked more carefully on the layerwise-QR approach currently used in TNRKit. It indeed minimizes the local cost function, so I was wrong in this aspect. However, the computational cost is higher for layerwise-QR (if we trust the analysis done by Codex, attached below). So let's only keep the current HOTRG-like implementation.

@Adwait-Naravane @VictorVanthilt You may want to revisit the two-layer QR in TTNR and switch to 3D-HOTRG to compress the temperature (imaginary time) direction.

Details

layerwise_qr=false path (current PEPSKit implementation):

This forms the two-layer Gram matrix MdagM / MMdag, whose matrix size is

D^2 x D^2

With a good contraction order, each Gram construction costs roughly

O(d^3 D^5)

then eigh_full! on the D^2 x D^2 matrix costs

O((D^2)^3) = O(D^6)

So one left/right orthogonalization costs

O(d^3 D^5 + D^6)

and the full virtual_projector has the same scaling, plus the common final R1 * R2 and SVD cost O(D^6).

layerwise_qr=true path (TNRKit implementation):

Each single-layer QR treats a PEPO tensor as a matrix of size approximately

(d D^3) x (d D)

so each QR costs

O(d D^3 * (d D)^2) = O(d^3 D^5)

After two such QRs, the intermediate M has matrix size roughly

(d^2 D^2) x (D^2)

and QR of that costs

O(d^2 D^2 * (D^2)^2) = O(d^2 D^6)

So one left/right orthogonalization costs

O(d^3 D^5 + d^2 D^6)

and the whole virtual_projector costs

O(d^3 D^5 + d^2 D^6)

up to the common final O(D^6) SVD.

Bottom line

Cost-wise, layerwise_qr=false is cheaper asymptotically in the physical dimension:

false: O(d^3 D^5 + D^6)
true:  O(d^3 D^5 + d^2 D^6)

For fixed small d, both are effectively O(D^6) in virtual bond dimension, but layerwise_qr=true carries an extra d^2 prefactor in the final QR stage. So the layerwise path is not computationally cheaper in the dense uniform-dimension estimate; its advantage would have to come from numerical behavior, block/symmetry structure, or practical contraction/QR implementation details.