Add simple update algorithm

[Yue-Zhengyuan]

This PR adds the core part of the simple update (arXiv 0806.3719) algorithm.

Main changes:

• Added a new struct InfiniteWeightPEPS, which is the collection of tensors on lattice vertices and bond weights. The weights are diagonal Schmidt matrices on the bonds of the iPEPS, and serve as a crude approximation of the environment of each site tensor. They are wrapped in the new struct SUWeight.
• Picked out a sub-function ctmrg_leftmove from ctmrg_iter, which only perform the left move for one column (of the sequential CTMRG (Make :sequential act column-wise #90)). It will be used in the full update (to be added later).
• Added some functions to perform rotation and reflection of the iPEPS, the weights and the LocalOperators.
• Added the Hamiltonian of the Hubbard model (Add Hubbard model #80) and the t-J model in src/operators/models.jl.
• Added simple update test in test/heisenberg.jl, and an example (by @Sander-De-Meyer) on Hubbard model at half-filling at examples/hubbard_su.jl.

Limitations:

• The implementation assumes square lattice geometry.
• For convenience, I made some requirements on whether certain virtual spaces of iPEPS tensors are dual spaces. Hopefully they can be relaxed in the future.
• The current implementation can currently only calculate models with nearest neighbor interaction only. Actually SU can handle next nearest neighbor terms (see arXiv 1908.09359), but that will be implemented in the future.

Things to do before merging:

• Improve docstring.
• Add more constructors for InfiniteWeightPEPS (postponed, when there is absolute need).
• Integrate with Hamiltonians in LocalOperator format.
• Use Julia's more advanced logging system instead of simply printing the information.

[Yue-Zhengyuan]

I found that my simple update test takes abnormally long time to run on Ubuntu (181m58.6s) and Windows (44m09.8s), but seems fine on macOS (14m49.1s). Seems that there are still lots of places for improvement...

[lkdvos]

Actually, the multithreading in the CI seems to have gotten disabled here: 3a9eb40 (the default is just 1 thread on CI)

I fixed the complaints about the required checks, this was a naming thing.

[pbrehmer]

Whoops, okay I wasn't aware that multi-threading was disabled all along. But I agree, that's a separate issue we should fix in a separate PR.

Other than that the PR looks pretty good, it's quite a big step forward! I will quickly add the function signatures in the docstrings.

[pbrehmer]

I added the function signatures as well as some cosmetic changes, and I merged the Heisenberg tests. This is good to go for me once the tests turn green.

Big big thanks to @Yue-Zhengyuan for the great effort!

[lkdvos]

I checked the stacktraces of the failing tests, and I think the issue is a combination of some tolerances and a bug in KrylovKit:
Jutho/KrylovKit.jl#110

The eigensolver in the rrule is failing to converge all requested values, and as a result there is an out of bounds error due to the aforementioned bug in KrylovKit. I think we should be able to avoid running into this altogether, by ensuring that the tolerances are set such that everything converges?

[Yue-Zhengyuan]

Thanks @pbrehmer for updating the docstring! I just found that there may be some mistakes on my part, and will fix them soon.

By the way, I'm not quite familiar with defining the rrule. Could you please also take a look at the updated rrule for sdiag_pow and see if it's correct?

[Yue-Zhengyuan]

When fixing my definition of eltype for InfiniteWeightPEPS, I notice that the definition of eltype is not quite consistent in PEPSKit. For CTMRGEnv it's defined as the number type of each corner/edge tensor:

Base.eltype(env::CTMRGEnv) = eltype(env.corners[1])

However for PEPS/PEPO, it's defined as the type of the tensor

Base.eltype(T::InfinitePEPO) = eltype(T.A)
Base.eltype(T::InfinitePEPS) = eltype(T.A)

For example,

using TensorKit
using PEPSKit
envs = CTMRGEnv(2, 2, 2; unitcell=(1, 2))
println(eltype(envs))
peps = InfinitePEPS(2, 2, 2; unitcell=(1, 2))
println(eltype(peps))

Output:

ComplexF64
TrivialTensorMap{ComplexSpace, 1, 4, Matrix{ComplexF64}}

Should we change the definition for InfinitePEPS/InfinitePEPO as follows?

Base.eltype(T::InfinitePEPO) = eltype(T.A[1])
Base.eltype(T::InfinitePEPS) = eltype(T.A[1])

[pbrehmer]

By the way, I'm not quite familiar with defining the rrule. Could you please also take a look at the updated rrule for sdiag_pow and see if it's correct?

It's almost correct - yesterday I missed a tiny error. When deriving the adjoint for complex variables, one also needs to take the complex conjugate corresponding to spow2'. Apart from that it looks good! (I will push that now as I'm anyway dealing with the eltype thing.)

... I notice that the definition of eltype is not quite consistent in PEPSKit.

Yes, you're right. However, it's not a big issue since the more relevant way to access the scalar element type is using scalartype as defined in VectorInterface. The definitions for scalartype are consistent for InfinitePEPS and CTMRGEnv. (The difference between the two essentially comes down to the fact that eltype is ambiguous for nested structs, whereas scalartype clearly returns the underlying scalar element type of the struct.)

Nonetheless, it makes sense to define eltype consistently, even though it is not needed much.

[Yue-Zhengyuan]

@pbrehmer Thanks for catching the rrule error!

Could you also help add the definition of scalartype for InfiniteWeightPEPS? I want to point out a subtlety that the weight matrices are always real (and there are real types of different precisions), but the vertex tensors can be complex.

By some test I discover that changing the eltype seems to affect the similar function. Maybe you can also help check this? In addition, I think similar should only be used to create a PEPS/PEPO with the same unit cell size and bond dimensions. So maybe we shouldn't define e.g. similar(peps, args...) but only similar(peps)?

[pbrehmer]

Well, my feeling is that InfiniteWeightPEPS are not really vectors since the elements are defined on different fields? So I don't think it can reasonably support the VectorInterface. Also, I don' think we really need it, since we don't do vector operations with that struct, right?

Regarding similar(peps, args...) I think it is fine the way it is. For example, if you have a complex PEPS and want to create a similar real PEPS you can just use similar(peps, Float64). For me similar still works fine like it should, I think.

[Yue-Zhengyuan]

OK, it is then reasonable we don't add scalartype for InfiniteWeightPEPS.

However, the new definition of eltype indeed affects similar for InfinitePEPO. I just pushed a fix for it, but I'm not quite clear what's under the hood. Nevertheless, we can be more explicit that similar for PEPS/PEPO only allows changing the element_type, and changing dims is not allowed.

[pbrehmer]

Ahh thanks for catching that, I didn't realize similar was defined differently for PEPOs.

Regarding the args... of similar which should be allowed: I think it is okay if we leave it up to the user which arguments are passed onto the similar call which acts on the underlying TensorMaps. For example, if the user inputs a domain and codomain such that the tensors are not PEPS/PEPO tensors anymore, the similar call will error.

[lkdvos]

Actually, the fix for eltype should be the other way around: that function is typically reserved for the type of the result of indexing or iterating, which is why it affects similar. At some point we realized that we need to also have access to the number type, so we introduced scalartype, but apparently the CTMRGEnv didn't get updated. In principle, I'd say it isn't even very well-defined for that, since it could be either an edge or a corner, but for the peps and pepo it really is

[pbrehmer]

In that case I would suggest just scrapping eltype for CTMRGEnvs, since we don't use that anyway and perhaps there is no reasonable and consistent way of defining it for CTMRGEnvs. Probably, for all cases where you would want eltype, you could just use scalartype.

[pbrehmer]

So the Heisenberg test does optimize if the derivative of the SVD is performed using GMRES. However, you do see that there are many gauge sensitivity warnings and also warnings about (Heisenberg SU(2)?) degeneracies - as expected in Jutho/KrylovKit.jl#110.

This circumvents Arnoldi failing to converge all values but of course the underlying problem is still there, which leaves the question of how in general such cases should be handled. Also, in principle we would need Lorentzian broadening here (it's still unclear how to implement that together with the KrylovKit adjoints).

Nonetheless, I'd be fine with merging the PR like this - I think we can tackle the Arnoldi issue also in a separate PR :-)

[lkdvos]

Good to for me too!