Add 3-site simple update (aka 3-site cluster update)

[Yue-Zhengyuan]

This PR adds the 3-site simple update (aka 3-site cluster update; see discussion in Issue #164), which can handle Hamiltonians containing up to next nearest neighbor terms (e.g. the J1-J2 model).

• It convert the Hamiltonian from LocalOperator format to a bunch of 3-site MPOs acting on the clusters. Relevant functions are added in src/algorithms/time_evolution/evoltools.jl.
• Truncation of bonds in the 3-site cluster after applying the time evolution gate is done using some "projectors" (I didn't bother to customize OBC MPS algorithms from MPSKit to suit my needs). Their construction is explained in detail in the beginning of the file src/algorithms/time_evolution/simpleupdate3site.jl. It is tested in the new test file test/timeevol/cluster_projectors.jl.
• Added an J1-J2 model SU example at examples/heisenberg_evol/j1j2.jl. We may also consider using the 3-site SU to provide initialization in the J1-J2 AD test. But note that the 3-site SU also requires a minimal unit cell size (2, 2) (in addition, bipartite structure is now not allowed), while the AD test currently uses 1-site unit cell with C4v symmetry.

This is kind of a large unseparated PR. Suggestions are welcome. I'm still testing fermion models, so the code is not finalized yet.

[codecov]

Codecov Report

Attention: Patch coverage is 96.71362% with 7 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/operators/localoperator.jl	45.45%	6 Missing ⚠️
src/algorithms/time_evolution/evoltools.jl	97.29%	1 Missing ⚠️

Files with missing lines	Coverage Δ
src/PEPSKit.jl	100.00% <ø> (ø)
src/algorithms/time_evolution/simpleupdate.jl	100.00% <100.00%> (+13.04%)	⬆️
src/algorithms/time_evolution/simpleupdate3site.jl	100.00% <100.00%> (ø)
src/algorithms/time_evolution/evoltools.jl	98.33% <97.29%> (-1.67%)	⬇️
src/operators/localoperator.jl	69.09% <45.45%> (+7.88%)	⬆️

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

[Yue-Zhengyuan]

I further added a test to compare the SU and 3-site SU energy for the Hubbard model without next-nearest neighbor terms (they should be close). I think the code is now ready for fermions.

[lkdvos]

I left some general comments throughout, keep in mind that I do not have much experience with the algorithm itself so I definitely might be missing some things or getting things wrong.
As always though, this is an impressive amount of work, and a great addition.

[lkdvos]

I think there is functionality in TensorKit to "flip" the arrow of a tensor: flip(t, i) should do this. This might be more efficient and readable and avoid having to reconstruct the diagonal tensor.

[Yue-Zhengyuan]

Thanks for informing me! But I find that flip acting on s will return a ordinary TensorMap, so I still need to add the conversion to DiagonalTensorMap.

[lkdvos]

Jutho/TensorKit.jl#243
Will look into fixing this soon ™️

[lkdvos]

So I didn't really think this through: defining flip(::DiagonalTensorMap) is a little bit more involved... In this case, I think you can avoid the entire flip by doing the svd in the other direction, ie tsvd!(transpose(rl); trunc) and then transposing the results. Why does this show up though, i.e. why do we need the rev option?

[Yue-Zhengyuan]

I did consider doing tsvd in the other direction, but didn't use this approach since I'm exhausted getting the twists right... I want to keep the current implementation and wait for the upstream improvement in TensorKit.

I need rev because the arrows in the cluster do not always point in the same direction. Under the usual choice of arrows in an InfiniteWeightPEPS, each :sw cluster looks like this:

    M1
    ↓
    y       
    ↓
    M2 ← x ← M3

Thus if I don't transpose rl for tsvd, I need to reverse the arrow direction between M1 and M2 afterwards.

In addition, After I apply the 3-site gate MPO, I fuse the internal virtual legs of the gate MPO and the cluster such that the result is M1 ← M2 ← M3. So I need to convert the result back to M1 → M2 ← M3. In principle I could construct the fusers that preserve the arrow direction in the cluster, so I don't need to use an additional rev argument to specify on which bonds should I reverse the arrows. But I want to be more flexible and enable manual tweaking of the arrows on each bond.

BTW, I found that transpose(rl) behave differently from permute(rl, ((2,), (1,))). Is it because transpose will not produce fermion parity signs?

using TensorKit
V = Vect[FermionParity](0 => 2, 1 => 2)
t = rand(V ← V);
transpose(t) ≈ permute(t, ((2,), (1,))) # false
transpose(t) ≈ twist(permute(t, ((2,), (1,))), 1) # true

[lkdvos]

The difference between the two is that transpose is planar, while permute is not. In other words, transpose will bend the legs either clockwise or anticlockwise (does not alter the result), while permute will cross them.

I understand the feeling about the twists, it really is not the nicest thing to work around...

[lkdvos]

Is there a way to detect the kind of hamiltonian that is supplied, and automatically map this to either the nearest-neighbour or the 3-site version?

[Yue-Zhengyuan]

This can surely be done: just check that for each term in the Hamiltonian, the distance between the sites of this term is 1. But this can be added in another PR. For test purposes I may deliberately to pass a Hamiltonian with only nearest neighbor terms to the 3-site SU.

[lkdvos]

The reason I'm asking is mostly since I was wondering if we really want to export simpleupdate3site as an external-facing function. I would prefer it if you could just use simpleupdate, and this would then either tell you that it is not supported or dispatch to a lower level function to implement the specific case at hand.
This has the benefit that if we ever extend or change this, we don't necessarily have to perform breaking changes for this.

[Yue-Zhengyuan]

Now simpleupdate can automatically decide whether to use the 2-site or the 3-site version depending on whether the input ham only has nearest neighbor terms. But I also added a keyword argument force_3site to use the 3-site version for any ham.

[Yue-Zhengyuan]

@pbrehmer I have added an example to do SU for J1-J2 model (with U(1) spin symmetry), but haven't generate the docs yet. Could you please try if the obtained SU state can be a good starting point for AD optimization?

[pbrehmer]

I'll take a look tomorrow and try out the SU J1-J2 state. I you want, I can then add a section to the example which uses the SU state to initialize an AD optimization (if it works well).

[pbrehmer]

AD optimization seems to work but it has some hickups at some point due to linesearching. I will try to improve that tomorrow. I also extended the commentary in the J1-J2 example a bit to make it similar style-wise to other examples.

[pbrehmer]

I decided to just restrict the maxiter of the optimization since the AD-optimized estimate already improves the SU estimate quite a bit and a proper full optimization would probably take long and be more complicated. I added the rendered example as well, so it should show up in the docs now. Ready to merge for me!

[Yue-Zhengyuan]

@pbrehmer I'm glad to see that with SU initialization the AD energy is lower than YASTN result, so maybe a good initialization does help. How does the energy compare with the AD optimization with just 1x1 unit cell starting from random initialization?

[pbrehmer]

I was also happy to see that! Actually, on a 1x1 unit cell from random initialization, I'm having a lot of trouble to optimize the PEPS. There, I'm encountering degenerate singular values sometimes and CTMRG sometimes struggles to converge during optimization such that the linesearch fails. I also tried Lander's backtracking linesearch algorithm and also there I'm having problems.

So my feeling is that the SU-evolved PEPS is actually quite a nice starting point for the AD optimization. I'll try to optimize that fully (up to higher accuracy) but I suspect that the environment dimension might be a little low for that - for small chi, PEPS optimization runs tend to get stuck at some point.

[Yue-Zhengyuan]

Looking at YASTN data, the first line is

# uuid CTM_BFGS100LS_U1B_D4-chi96-j20.5-run0-c2U1BD4j2045chi96n0

So he is already using $\chi = 96$. Then he measure the optmized state with various $\chi$'s, and the energy converges to 1e-6 precision for $\chi \gtrsim 128$. Indeed a large $\chi$ is necessary. Maybe now IterSVD can be used to speed up the optimization?

In the calculation he also imposes the geometric $C_{4v}$ symmetry to the 1x1 unit cell in addition to a sub-lattice spin rotation, which further reduces the number of parameters. Did you do this for the 1x1 optimization before?

[pbrehmer]

Indeed! I'm not yet sure if IterSVD already pays off here with D=4 but in any case we can push the environment dimension quite a bit. I didn't impose $C_{4v}$ symmetry but I'll try that next. Of course we don't have proper $C_{4v}$ CTMRG yet because the environment should consist only of one corner and edge tensor (per unit cell entry). But I'm thinking that we should implement that anyway at some point (also if we want to run proper benchmarks).