Add PEPS initialization based on fidelity maximization

[pbrehmer]

This PR adds tools to initialize a PEPS based on maximizing the fidelity with respect to another PEPS. As we discussed on Tuesday, this would be useful e.g. to extract a single-site approximation of a SU evolved PEPS on a $2\times2$ unit cell. This would be great to have since SU/FU evolved initial PEPS can then be used for single-site variational optimization.

Note that I use TensorKit's embed branch in order to use the embed! function:

] add TensorKit#ld-embed

[codecov]

Codecov Report

❌ Patch coverage is 0% with 74 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/algorithms/optimization/fidelity_initialize.jl	0.00%	74 Missing ⚠️

Files with missing lines	Coverage Δ
src/Defaults.jl	85.71% <ø> (ø)
src/PEPSKit.jl	100.00% <ø> (ø)
src/algorithms/optimization/fidelity_initialize.jl	0.00% <0.00%> (ø)

... and 2 files with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[pbrehmer]

@leburgel, @lkdvos I tried to recreate what you explained on Tuesday but I'm still a bit unsure about some things here:

• It's not quite clear how to treat unit cells in the fidelity iteration in maximize_fidelity!
• To compute the fidelity, should I use network_value which takes care of environment normalization or do we really just need to contract the ket and bra with the environment (as in contract_local_norm)?
• The procedure is not yet properly normalized since the fidelity will get larger than 1 sometimes

I'll figure this out over the next days, but I'm happy about any hints/comments.

[lkdvos]

I'll try and go over this in a bit more detail later, just wanted to leave these small comments before I forget

[Yue-Zhengyuan]

I think it will be straightforward to do this for mixed states (PEPOs) as well? One particular application is first obtaining exp(-βH) from simple/full update, tracing out the physical axes, and putting it on a periodic lattice and contracting using TNR (when we want to extract CFT information). For TNR, afaik it is often cumbersome to work with non-trivial unit cells. Anyway, there are other ways to directly obtain the exp(-βH) PEPO with trivial unit cell (e.g. HOTRG, cluster expansion) and we don't need to do this right now.

[pbrehmer]

The new tests should pass once a new TensorKit version is released since we need absorb!. Other than that I think the fidelity maximization is now somewhat working. I will test this with SU evolved 2x2 PEPS as initial states for optimization in the next days.

[pbrehmer]

Small update on applying this to approximate a $2\times2$ SU evolved PEPS by a single-site PEPS: This doesn't seem work. Note that the SU PEPS has inhomogeneous virtual spaces across the unit cell (here for $D=4$):

julia> domain(peps_SU[1])
(Rep[U₁](0=>2, 1=>1, -1=>1) ⊗ Rep[U₁](0=>2, 1=>1, -1=>1) ⊗ Rep[U₁](0=>2, 1=>1, -1=>1)' ⊗ Rep[U₁](1/2=>2, -1/2=>1, 3/2=>1)')

With that there are problems: One is that _spacemax will actually create a U1Space with bond dimension 9 which is not really what we want. If one takes e.g. U1Space(0 => 2, 1 => 1, -1 => 1) as the approximate virtual space, then there will be problems with having no fusion channels. In any case the ultimate problem is that the contraction in

PEPSKit.jl/src/algorithms/optimization/fidelity_initialize.jl

Line 78 in 2fc9f5e

peps′ = _∂local_norm(peps₀, envnw)

will yield mostly zero entries in this symmetric case. If we then try to contract that PEPS, the network_value gets super small. This process is repeated a few times so that the algorithm becomes numerically unstable and the final PEPS is mostly zeros and cannot be contracted anymore.

For dense PEPS the algorithm does seem fine and maybe for symmetric PEPS with compatible virtual spaces it also works. But for the $2\times2$ SU PEPS with varying symmetry sectors I need to try something else I think.

[lkdvos]

@pbrehmer, the vanishing value and no fusion channels sounds like you are trying a case where a 1x1 unitcell is actually forbidden by the symmetry, in which case we are not expecting this to work indeed. I think the main idea was to try this in the case where we know the state has a 1x1 unitcell, and we had to artificially enlarge it to 2x2 because the optimization algorithm (typically SU) breaks that, but using a different optimization algorithm would have produced a 1x1 unitcell without issues.

Given this, we would typically expect the virtual spaces to be somewhat compatible, and the variations to be mostly due to numerical noise, rather than structural differences, such that all of the corners we are cutting are warranted.

[leburgel]

Small update on applying this to approximate a 2 × 2 SU evolved PEPS by a single-site PEPS: This doesn't seem work. Note that the SU PEPS has inhomogeneous virtual spaces across the unit cell (here for D = 4 ):

julia> domain(peps_SU[1])
(Rep[U₁](0=>2, 1=>1, -1=>1) ⊗ Rep[U₁](0=>2, 1=>1, -1=>1) ⊗ Rep[U₁](0=>2, 1=>1, -1=>1)' ⊗ Rep[U₁](1/2=>2, -1/2=>1, 3/2=>1)')

Ah, I see. I thought the difference in virtual spaces between different directions would just be the distribution of degeneracies, I didn't know there was just one space with halfinteger charges compared to all the rest. It's quite curious that there's just the one that is all halfintegers and the rest are integers, but they're never mixed anywhere. Just from looking at this distribution, to me this looks like something that could be better addressed by shifting the physical space through adding a uniform U1Irrep(1/2) auxiliary charge everywhere.

Would this make sense, or would you not expect a nonzero charge density for the ground state here?

EDIT: missed Lukas' reply. Do we expect the ground state to have a 1x1 unit cell here physically, or does it really break translation invariance?

[pbrehmer]

Just from looking at this distribution, to me this looks like something that could be better addressed by shifting the physical space through adding a uniform U1Irrep(1/2) auxiliary charge everywhere.

Would this make sense, or would you not expect a nonzero charge density for the ground state here?

Not quite sure what you mean, could you elaborate a bit more on that? How would adding a uniform U1Irrep(1/2) charge help here?

Do we expect the ground state to have a 1x1 unit cell here physically, or does it really break translation invariance?

I'm not super sure about the ground state properties at $J_2=0.5$ but it seems from other simulations that the trick of applying a sublattice rotation to create a translationally symmetric PEPS on a single-site unit cell does work, see e.g. Juraj's repository.

[leburgel]

Not quite sure what you mean, could you elaborate a bit more on that? How would adding a uniform U1Irrep(1/2) charge help here?

Just for the example of the tensor you printed, if you split off an auxiliary physical leg with space U1Space(1/2 => 1) from the one virtual leg that has halfinteger charges, then you can actually use the same virtual charges everywhere. A similar trick is used here to use uniform virtual spaces for the U(1)-symmetric XXZ model in the region where the ground state has Néel order, by adding a staggered charge shift which essentially removes the need to have virtual spaces with halfinteger charges.

I'm not super sure about the ground state properties at J 2 = 0.5 but it seems from other simulations that the trick of applying a sublattice rotation to create a translationally symmetric PEPS on a single-site unit cell does work, see e.g. Juraj's repository.

I can't say I'm entirely sure about this, but from reading the explanation here it seems like the local tensors here do have a non-trivial total U(1) charge. I think our labels are half of the ones used there (the physical space there has +1 and -1 charges, whereas we use +1/2 and -1/2), so a total U(1) charge of +1 as indicated there would be a charge shift by U1Irrep(1/2) for us. Of course we can only do this in the regime where there is no real inherent spatial symmetry breaking in the ground state. If there is some other emergent magnetization pattern (e.g. Néel order) we would have to use a non-trivial unit cell with a non-uniform charge shift.

[Yue-Zhengyuan]

If there is some other emergent magnetization pattern (e.g. Néel order) we would have to use a non-trivial unit cell with a non-uniform charge shift.

The simple update states I provided to Paul do have a nonzero staggered magnetization. This might indeed be the physical reason why 1x1 unit cell won't work.

@pbrehmer Is it possible to perform the sub-lattice spin rotation on my PEPS before optimization? If not possible with U(1) (which is kind of expected if this is the case), we may try not to use any spin rotation symmetry.

[leburgel]

The simple update states I provided to Paul do have a nonzero staggered magnetization. This might indeed be the physical reason why 1x1 unit cell won't work.

I see, and is it staggered because we didn't do a sublattice rotation to map antiferromagnetic order to ferromagnetic order at the very start?

But even in that case, I think we can do a non-uniform charge shift to account for the antiferromagnetism (exactly like is done here) while still using the same virtual space everywhere. In that case, a fidelity initialization would still be possible, it just wouldn't be single-site anymore

Maybe this is going a bit beyond the scope of the actual PR here, but I do think the issue can be solved with the approach here. How about we break off a discussion about fixing up the J1-J2 example somewhere, and return to this once we have figured that out?

[pbrehmer]

Is it possible to perform the sub-lattice spin rotation on my PEPS before optimization? If not possible with U(1) (which is kind of expected if this is the case), we may try not to use any spin rotation symmetry.

I think it would be possible to sublattice rotate it but applying a $U(1)$-symmetric $\sigma^y$ to the PEPS will result in an auxiliary charge on the physical leg; not sure if that'd be a problem.

I see, and is it staggered because we didn't do a sublattice rotation to map antiferromagnetic order to ferromagnetic order at the very start?

Yes, precisely. (Doing a sublattice rotation at the start corresponds to changing the signs in front of the next-nearest neighbor Hamiltonian terms such that we get a ferromagnetic single-site ground state.)

Maybe this is going a bit beyond the scope of the actual PR here, but I do think the issue can be solved with the approach here. How about we break off a discussion about fixing up the J1-J2 example somewhere, and return to this once we have figured that out?

Yes, very much agree, this specific test case is beyond the point of the PR. Maybe we can discuss this a bit further tomorrow. Actually I think the approximate method is kind of useful in general for PEPS optimizations. For example, when one wants to initialize a higher $D$ PEPS from a smaller $D$ optimized one this could be a good tool.

[pbrehmer]

As discussed earlier, I added an algorithm struct for approximate. To that end I also implemented the whole default/kwarg algorithm struct mechanism such that the approximate algorithm can be specified either as a struct or just using keyword arguments. Of course the names ApproximateAlgorithm (for the abstract type) and FidelityMaxCrude (for the specific fidelity maximization algorithm) are up for discussion and I'm happy to change them.

I changed the method signature to match the one from MPSKit. In that process I decided that specifying the envspace isn't really necessary and would make the interface weird, so I dropped that and instead now one can specify the truncdim via the boundary algorithm. Note that I also got rid of the mutating approximate! as it wasn't actually doing anything in-place anymore since absorb! was moved out.

[pbrehmer]

@leburgel @Yue-Zhengyuan Quick update on running SU with staggered charges on the physical spaces (as discussed yesterday): This does seem to work indeed and I obtain SU energies which are quite similar to the previous runs. Note that I first had to fix the index rotation functions to make this work, see #267. I haven't run variational optimizations on these states yet, but I'll keep you posted.