Local truncation for two layers of InfinitePEPO

[Yue-Zhengyuan]

This PR introduces the most naive algorithm to truncate the virtual space of the product of two layers of InfinitePEPOs. 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.

The function doing the truncation overloads MPSKit.changebonds, and dispatches on the algorithm struct LocalTruncation <: BondChangeAlgorithm.

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.30508% with 1 line in your changes missing coverage. Please review.

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

Files with missing lines	Coverage Δ
src/algorithms/compress/local.jl	100.00% <100.00%> (ø)
src/utility/util.jl	68.81% <ø> (ø)
src/operators/infinitepepo.jl	77.19% <83.33%> (+0.34%)	⬆️

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

[leburgel]

Can we split off isa(state, InfinitePEPS) ? 1 : 2 into a dedicated num_physical_ind method (or something similar)? That way it's clearer that that's what's being probed here, and also there's probably other places where this an be reused.

[leburgel]

Probably adding the trunctol here is unavoidable, but I think it would be better to control this with some global value in Defaults rather than hardcoding it here. Also, would an rtol not be safer, for tensors with non-standard normalizations?

[Yue-Zhengyuan]

If we are confident about the implementation (even for fermions about the twists), we can actually manually set all negative eigenvalues in D.data to 0 similar to positive_approx. In this way a truncation tolerance can be avoided.

[leburgel]

Do we have an appropriate test to check this, so we can then just manually discard all negative eigenvalues? If this is too much effort that's fine for now, I would just like to keep track of hardcoded tolerances as much as possible by collecting them all in the same place.

[Yue-Zhengyuan]

In TNRKit (in the context of 3D HOTRG) I have a test that uses random fermionic PEPOTensors to ensure that no negative eigenvalues are produced due to wrong twists. I can port that here.

[leburgel]

This signature is a bit tricky I think. Here, it makes total sense to do it like this. However, for other uses of approximate (which are not implemented yet, but we probably should soon), things would probably look like

approximate(n0, to_approximate, alg::ApproximateAlgorithm, optional_args...)

where n0 is some initial guess specifying the desired structure (e.g. PEPS with the virtual spaces you want), to_approximate is the thing you want to approximate (e.g. a PEPO and a PEPS, representing the application of the PEPO to the PEPS), and then your algorithm and any other things you might potentially need. Here, all of the structure of n0 is in the trunc field of the LocalApprox algorithm, so there's no need for that one. Then to_approximate = (ρ1, ρ2) can just be added as two separate arguments.

I think this is fine as it is, but just to say that we might have to add quite different signatures for different algorithms. That being said, since there's really no n0 here, it would be very artificial to match the example signature I gave here (it would become approximate(trunc, (ρ1, ρ2), alg), which is a bit ridiculous).

So I don't really know how to do this better for now, just wanted to make note of this.

[Yue-Zhengyuan]

Also adding to this list: one may also want to approximate the inner product of two PEPSs (or PEPOs regarded as purified PEPSs), obtaining a partition function with a moderate bond dimension. This is encountered when one wants to feed it into TNR.

[leburgel]

Good point. If there's no good general solution we can just continue with what you suggested, but maybe we can discuss how we should approach this in more detail soon.

Thinking about it a bit more, the case implemented here and the the ones we both brought up feel like some kind of mixture of MPSKit.approximate (implying some variational approximation) and MPSKit.changebonds (implying a pure bond dimension truncation or expansion). In that sense, this use case could be written as

MPSKit.changebonds(ρ1 * ρ2, alg)

where we can perform the multiplication lazily as a stacking of the PEPO unit cells. That might actually be the most fitting I think. However, if we have too much trouble fitting the use cases here into the MPSKit algorithms and their usual signatures, we can always consider moving away from those altogether.

[Yue-Zhengyuan]

I took a look at MPSKit. There is also the syntax

changebonds(ψ, H, alg, ...)

which I think is optimizing the MPO-MPS product H * ψ. I guess it is also OK to directly use changebonds(ρ1, ρ2, alg) without even performing the lazy stacking first?

[Yue-Zhengyuan]

I'm curious why the gradients test needs so much more time (and allocates much more memory) when running with min (Julia v1.10; around 110 min) than with v1.12 (around 40 min). There's not such a huge gap for other test groups.

[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!

[lkdvos]

This might be a bit nitpicky, but can you maybe change the name of these functions, since you aren't actually performing a QR, but rather some form of square root?

Just out of curiosity, do you have any idea how this compares to a BP-style or SU-style approach? I can't help but feel like this is some rank-1 approximation of the double-layer object, but somehow in a symmetrized/hermitianized manner. (This is completely irrelevant for this PR, so feel free to ignore, I'm just curious)

[Yue-Zhengyuan]

change the name of these functions

This is a bit like QuantumKitHub/MatrixAlgebraKit.jl#109, where I wanted to do left/right_orth but only need L/R, and don't care which algorithm is used to obtain Q. Maybe just call it left_orth_twolayer, or do you have any suggestions in the context of MAK?

how this compares to a BP-style or SU-style approach

I have to say I don't have a clear answer... One thing this can do but BP/SU cannot is 1x1 unit cell.

rank-1 approximation

Naively I can't see where a truncrank(1) is involved. What do you mean exactly?

[lkdvos]

I think the idea I had is that this feels like the surroundings of the local bond are approximated by some kind of local tensor, i.e. there is no edge between the qr_twolayer and lq_twolayer part, but it might also just be a bit of a stretch so this probably doesn't buy us much.

How about reduced_twolayer_pepo_left/right? Maybe twolayer_left_residual, to indicate it's only the R factor?

It almost looks like a square root too, so I guess this is a cholesky decomposition? Honestly naming things is hard, but calling it qr or lq seems a bit misleading.
For the MatrixAlgebraKit issue, I have to admit that this slipped my attention a bit, for which I apologize. There as well, it is kind of hard to come up with a name, and there is not too much to gain from not forming the Q I think.

[Yue-Zhengyuan]

Below is copied from Gemini. A below is the two-layer A1 A2, and M is the MdagM tensor.

---

1. The Structural Difference (Full vs. Triangular)

When you compute $M = A^\dagger A$, $M$ is a Hermitian positive semi-definite matrix.

• Cholesky Decomposition factors $M$ explicitly as $M = R_{chol}^\dagger R_{chol}$, where $R_{chol}$ is strictly an upper triangular matrix. (This $R_{chol}$ is identical to the $R$ you would get if you did a standard QR decomposition directly on $A$).
• Your Eigendecomposition Method factors $M = U^\dagger D U$ and defines $R_{eig} = \sqrt{D} U$. However, $R_{eig}$ is generally a full matrix, not an upper triangular one.

Mathematically, because both matrices satisfy $R^\dagger R = M$, they are related by a unitary transformation. There exists some unitary matrix $V$ such that $R_{chol} = V R_{eig}$.

2. Connection to Polar Decomposition

While standard QR yields the Cholesky factor, your eigendecomposition method is actually computing a variant of the Polar Decomposition.

In polar decomposition, $A = W P$, where $W$ is an isometry ($W^\dagger W = I$) and $P$ is a unique Hermitian positive semi-definite matrix given by $P = \sqrt{A^\dagger A} = \sqrt{M}$.

Using your eigendecomposition $M = U^\dagger D U$, the matrix square root is:

$$P = U^\dagger \sqrt{D} U$$

Since you defined your $R$ as $\sqrt{D} U$, we can see that $P = U^\dagger R$. This means your $R$ matrix is exactly the positive-definite polar factor $P$, just rotated by the unitary eigenvector matrix $U$.

[lkdvos]

Suggested change

	@tensoropt MdagM[x2 z z′ x2′] :=
	@tensor MdagM[x2 z z′ x2′] :=

Given that there are only two tensors involved, there is no need to determine the best contraction order here :)

On a separate note, do you think it is possible to write this as a single @tensor call? This makes it easier to play around with allocation stuffafterwards, since the intermediate MdagM isn't required, and the @tensor call will ensure that the intermediate permutations are minimized.

[Yue-Zhengyuan]

Sure, I can write it in a single pass. Then four tensors are involved, which I guess will need @tensoropt?

[lkdvos]

That, or explicit brackets, or even just the @autoopt machinery, since that would fill in the sizes a bit more

[Yue-Zhengyuan]

I changed to @autoopt and added d/D to the contraction labels.