More general `InfiniteMPO` constructors and compression, iDMRG, and subspace expansion

[mtfishman]

Hey Matt,

On my own fork, I have finished the implementations of a few useful features (paper coming out soon using it). This includes, from trivial to less trivial:

• SU(3) and SU(4) irrep using color conservation (nothing non-Abelian)
• some small fixes in general
• implementation of a few convenient functions for data analysis (QN-resolved analysis of ent spectrum, transfer matrix spectrum + some correlation functions)
• a working iDMRG for InfiniteMPOMatrix
• a working InfiniteMPO constructed from an InfiniteMPOMatrix
• a working iDMRG for InfiniteMPO
• compression of iMPOs using InfiniteMPOMatrix and the approach by Parker, Cao et Zaletel in Phys.Rev.B 102 . It is not yet fully benchmarked so that should clearly be experimental , especially as it is the simplest, non-fail proof version.
• I also plan to do this week some modifications on the subspace expansion as I cannot afford the default version in my new project (the one that required the compression of iMPOs)

It might be a lot for a single merge request. So how would you like to proceed? It will take me one or two weeks anyway to make things a bit cleaner.

Originally posted by @LHerviou in #69 (comment)

[mtfishman]

It sounds like an impressive list of features you've developed. My main reaction is that I am trying to minimize my maintenance burden from this package and others so I can focus on the development of ITensorNetworks.jl, which has our next generation MPS, TTN, and PEPS/general TN solver and evolution code, using a graph-inspired interface for tensor networks. My goal there is to take lessons from this package like the CelledVector design and design a general infinite tensor network type, with infinite MPS/VUMPS just as a special case.

So with that in mind, let me go through your list:

• SU(3) and SU(4) irrep using color conservation (nothing non-Abelian)

Can this just be shared as an example code for others interested in using that? It doesn't sound like it requires modifications to VUMPS/InfiniteMPS itself.

• some small fixes in general

General small fixes would definitely be appreciated (maybe split up into PRs with fixes that are thematically similar).

• implementation of a few convenient functions for data analysis (QN-resolved analysis of ent spectrum, transfer matrix spectrum + some correlation functions)

I was hoping to provide this kind of functionality as examples at first, and think of this package as more bare-bones. For example, my goal was to provide an InfiniteMPS type and VUMPS but not much else, but provide examples and a nice interface so people can do that kind of thing (measurements and analysis of the infinite MPS) themselves. Also, I was hoping that by making it easy to slice an infinite MPS into a finite MPS then users can leverage finite MPS tools from ITensors.jl to perform many useful operations for analyzing the state. The developer burden of maintaining and improving the MPS/MPO code in ITensors.jl was immense. Things are more stable in that package now, but since I am trying to push the development of ITensorNetworks.jl as a replacement for the existing MPS/MPO code in ITensors.jl (and ultimately a replacement for this package), and also have the responsibility of maintaining the other packages in the ITensors.jl ecosystem, I don't want to double that effort by providing and maintaining too many features in this package. So I think it is best to share specialized tools like that in examples unless there is something that crosses the threshold of being very general and that has multiple use cases. If anything, my goal would be to try to compress the code of this package where possible, as I discuss below.

• a working iDMRG for InfiniteMPOMatrix

That sounds great. I would appreciate a PR for that. If the code is simple and it is well tested it would make sense to include in this package. However, I think we should either choose InfiniteMPO or InfiniteMPOMatrix and stick with that, or make the codes only very minimally different, in order to minimize maintenance burden.

• a working InfiniteMPO constructed from an InfiniteMPOMatrix

That sounds useful. I think it is unfortunate how many different MPO formats there are in this package, and makes maintaining the package challenging. I think it would be best if we stuck with just one external format if possible, and a minimal number of internal formats, and try to design it so that we just have one version of VUMPS/iDMRG to maintain, optimize, etc. For example, perhaps you can appreciate the challenges involved given the number of performance PRs I had to make: #56, #64, #65, #67. I think this would have been largely unnecessary if we had one VUMPS implementation with a single infinite MPO format. Do you think at this point we can just go with one format? Is there a reason to keep around both InfiniteMPO and InfiniteMPOMatrix? Additionally, I think we should drop support for InfiniteSum{MPO}/InfiniteSum{ITensor} since I would hope those become largely unnecessary with good support for a single infinite MPO of some form.

• a working iDMRG for InfiniteMPO

Same comments as above.

• compression of iMPOs using InfiniteMPOMatrix and the approach by Parker, Cao et Zaletel in Phys.Rev.B 102 . It is not yet fully benchmarked so that should clearly be experimental , especially as it is the simplest, non-fail proof version.

I think @JanReimers has a version of that algorithm working in ITensorMPOCompression.jl. Can you join forces on that? @JanReimers, is the infinite MPO version working?

• I also plan to do this week some modifications on the subspace expansion as I cannot afford the default version in my new project (the one that required the compression of iMPOs)

That sounds great. Improved subspace expansion tools would be very appreciated, hopefully they can also be shared with the solvers in ITensorTDVP.jl/ITensorNetworks.jl since those would also benefit from better subspace expansion methods. @b-kloss and @emstoudenmire have worked on some subspace expansion methods in https://github.com/ITensor/ITensorTDVP.jl/pulls so I'm curous how it relates and if it has some analog for finite MPS.

Originally posted by @mtfishman in #69 (comment)

[mtfishman]

On Thu, Feb 23, 2023 at 10:11 PM Matt Fishman @.***>
wrote:

• compression of iMPOs using InfiniteMPOMatrix and the approach by
  Parker, Cao et Zaletel in Phys.Rev.B 102
  https://doi.org/10.1103/PhysRevB.102.035147 . It is not yet fully
  benchmarked so that should clearly be experimental , especially as
  it is the simplest, non-fail proof version.

I think @JanReimers https://github.com/JanReimers has a version of that
algorithm working in ITensorMPOCompression.jl
https://github.com/JanReimers/ITensorMPOCompression.jl. Can you join
forces on that? @JanReimers https://github.com/JanReimers, is the
infinite MPO version working?

Yes it is working for InfiniteMPO with dense storage and for some cases
(Trans. Ising) with block sparse storage. The caveats are:

1. It assumes all the new QR/RQ code is in ITensors. But I think those
   changes are still in the PR queue.
2. For block sparse storage the combiners used for SVD/QR/RQ will in
   general re-order and combine the QN spaces, breaking regular form. So
   right now I am working on gauge transforming back to regular form after the
   decomp operations. For simple systems like transverse Ising this problem
   does not arise. On the other hand, for a Hubbard model the spaces are much
   more complicated so I am using that as a benchmark for success.

In any case I am happy to coordinate with @LHerviou
https://github.com/LHerviou! @nbaldelli https://github.com/nbaldelli .
Perhaps their work has eclipsed what I have been doing, in which case I can
possibly contribute with documentation work.

Regards
Jan

Message ID: @.***>

Originally posted by @JanReimers in #69 (comment)

[LHerviou]

Again, sorry for the delays, everything went into spam.

Just want to preface that I am very grateful for your work and help. I am also maintaining some much smaller librairies for our group, and it is quite the hassle already.

I am happy to discuss with you @JanReimers, and if your version is better, so be it! And given that I selfishly implemented it only for symmetric tensors, and using my weird format like the InfiniteMPO matrix (which simplifies a lot some of the discussions, it surely is the case).
I am still having a minor issue on mine (some imprecision when the MPO becomes very large and non-orthonormal), so it would be useful for me.

For the rest, I will try to submit some independent PRs in the month to come. I am again in the applying season, so it might get complicated before mid april.

I am looking forward to the new ITensors in any case.

[JanReimers]

Hi @LHerviou if you have some code to share for generating your MPO, I would be interested in putting it through the iMPO compression module to see what it can do.

[LHerviou]

Well, it is a bit complicated as it is related to the tensor approach.
Here is a snippet that should work (it uses that last fix for the generation of the indices).

I solved my problem, though. It was a simple orthogonalization issue.

N = 6
model = Model"fqhe_2b_pot"()
model_params = (Vs= [1.0, 0.0, 1.0, 0.0, 0.1], Ly = 6.0, prec = 1e-8)
function initstate(n)
	mod1(n, 3) == 1 && return 2
	return 1
end
p = 1
q = 3
conserve_momentum = true

function ITensors.space(::SiteType"FermionK", pos::Int; p=1, q=1, conserve_momentum=true)
  if !conserve_momentum
    return [QN("Nf", -p) => 1, QN("Nf", q - p) => 1]
  else
    return [
      QN(("Nf", -p), ("NfMom", -p * pos)) => 1,
      QN(("Nf", q - p), ("NfMom", (q - p) * pos)) => 1,
    ]
  end
end

function fermion_momentum_translator(i::Index, n::Integer; N=N)
      ts = tags(i)
      translated_ts = ITensorInfiniteMPS.translatecelltags(ts, n)
      new_i = replacetags(i, ts => translated_ts)
      for j in 1:length(new_i.space)
        ch = new_i.space[j][1][1].val
        mom = new_i.space[j][1][2].val
        new_i.space[j] = Pair(
          QN(("Nf", ch), ("NfMom", mom + n * N * ch)), new_i.space[j][2]
        )
      end
      return new_i
    end

function ITensors.op!(Op::ITensor, opname::OpName, ::SiteType"FermionK", s::Index...)
  return ITensors.op!(Op, opname, SiteType("Fermion"), s...)
end

s = infsiteinds("FermionK", N; translator = fermion_momentum_translator, initstate, conserve_momentum, p, q);
ψ = InfMPS(s, initstate);

H = InfiniteMPOMatrix(model, s; model_params...);

[LHerviou]

The update on the automatic shifting actually broke this part of my code.
EDIT: Solved

[LHerviou]

I finished pretty much everything I wanted to implement on my branch, and I am ready to bring it to the main branch.
Related to that, I had a question or two on some conventions we might want to change/simplify, but which would be a breaking change.
By the way, thanks to @JanReimers for the updates. Your code is nicer than mine. Is your QR working now?

Currently the InfiniteMPOMatrix has many empty/undef elements, and a non-uniform index pattern.
This is entirely my fault. I would like to push for the following changes, with two options

Option 1, a minima imho

Hm should have the form below. Only link indices are shown.
I       l=n-->0
M-->l=n   l=n-->0-->l=n
0-->l=n-1   l=n-->M-->l=n-1 ...
:       :
0-->l=1     l=n-->0-->l=1
0        l=n-->0

That is to say, each line and column has a well defined index (or no index for first and last).
Each tensor is initialized to an empty tensor with the correct indices, if it is not occupied.
We no longer need to assume an underdiagonal form of Hm, which allows us to easily add InfiniteMPOs
Perks: better structure, applying the MPO on the environments for example, can be implemented as a simple multiplication.

Option 2, same, but we add the dummy tensors Jan needed on the first and last lines and columns, so that everything is a four tensor. This simplifies marginally the code.

[JanReimers]

Hi Loic, good to hear from you. Yes the Ac/bA block orthogonalization using column pivoting QR is working well. I use your fqhe model as a tough unit test case. In case you are interested here is what I am finding for that system. The iMPO construction goes as

Model"fqhe_2b_pot"() --> Opsum --> InfiniteSum{MPO} -- > InfiniteMPOMatrix --> InfiniteMPO

The bond dimension of the final iMPO is Dw=97 with 97 QN blocks ... i.e. no blocks are merged. So this is quite and expensive MPO to work with either in the InfiniteMPOMatrix or InfiniteMPO representation. If we simply combine blocks (not possible for InfiniteMPOMatrix?) the QN spaces drops to 25. Furthermore if we truncate/compress the InfiniteMPO with cutoff=1e-15 then Dw drops to 37. The resulting bond spectra is very spread out, so if you raise the cutoff you will probably start eating into the real Hamiltonian, rather just removing numerical noise.

I would like to benchmark a GS calculation for InfiniteMPOMatrix vs. compressed InfiniteMPO. But as you already know I am running into the zero-norm expansion issue when trying to expand the subspace beyond D=1. Is that why you switched to iDMRG in you research?
I don't know if the "long range" explanation in the warning message is correct. The reason is that I can do Heisenberg and Hubbard models with arbitrary distant interactions and not the get that message. Do you have any ideas about how to get the subspace expansion working for fqhe? (maybe I accidentally got a default translator hiding in H ? I will check)

With regards to code complexity: For InfiniteMPOMatrix one must do a lot of manually coded contractions. For InfiniteMPO we can just leverage the contraction algos built into ITensors. For example the generate_twobody_nullspace function is a one liner:

return noprime((L[k-1]*H[k])*ψ.AL[k]*ψ.C[k]*ψ.AR[k+1]*(H[k+1]*R[k+1]))

I have three branches ready to pull into the main branch, to support VUMPS for InfiniteMPO and compressed InfiniteMPO. But if you prefer we can wait until your changes are pulled in and then I can deal with any conflicts. If you want to change the structure of InfiniteMPOMatrix that is fine with me, it is no hard to adapt the conversion code. I don't need to receive order 4 tensors, my code will convert to order 4 as needed. So again you can do whatever works best for your code.

[LHerviou]

Hi Jan,
Yeah, sorry for the long silence. Hopefully I have a bit more time now to make everything public.

I have also implemented my own MPO compression following the same method, and I get very similar results. I tried a bit using iDMRG and it seemed cutting away the smaller singular values do not affect significantly the wavefunctions, though indeed there is no clear gap and the cut off is a bit arbitrary.
I get good results with iDMRG for MPO bond dimensions of up to 500-600 without too much problem.
It does seem like I eat a wall with the ITensors MPO construction though when the number of operators reach a few thousand. It is not a big problem as I can split the Hamiltonian in more blocks, but I may try to have a look at that at some point.

I switched to iDMRG mainly due to some problems with convergence using VUMPS in frustrated models (where several groundstates exist). I found it was better/more stable to start with iDMRG, and then switch to VUMPS once converged close enough to a specific state.
The long-range is probably only partially correct. What matters is more the number of symmetry. The typical example are pair hopping style models (or FQHE equivalently) of the form:
c^d_i c_{i+1} c_{i+2} c^d_{i+3}.
Two site idmrg also cannot explore the phase space by itself (without the subspace expansion). And similarly the default subspace expansion that we have which only takes into account two sites cannot work (it is not an error from your side in principle).
One needs either to have subspace expansion on more than 2 sites (I implemented it at some point, but it is slow and had some issues) I have a working subspace expansion for InfiniteMPO, is it what you are looking for? I can also incorporate it.

Finally, concerning the changes I wanted: if you can wait a few days, I will try to submit something tomorrow or Thursday. If I cannot, feel free to proceed, and I will manage the integration.
If you have no strong opinion, I will think a bit about what I believe is the best here.

[JanReimers]

I will experiment with adding a four site term to the Hubbard model (if opsum supports that?). Good thing to have in the unit tests anyway. I would be interested to see what your subspace expansion code for handling >2-site interactions looks like. If you give me branch name in your repo I can probably find it. It would be great to generalize that code accordingly. For efficiency we should find a way to detect and store the max_interaction_order in the MPO struct.

I am also interested in how you calculate the L/R environments on a compressed InifniteMPO. The reason is that in general compression results in a non-triangular A-block which interferes with triangular algo that avoids eigen solvers. My current workaround is to store the gauge transform from the uncompressed MPO to the compressed one. Calculate L/R on the lower triangular uncompressed MPO, and then transform L/R to the new compressed basis. The code is small, but as you can imagine this defeats some of the performance benifits you might expect from compression.

Thanks
Jan

[LHerviou]

OpSum will have no problems (it s for odd number of fermionic terms that one has problems, and in any case, one could work with hardcore bosons and have the same ergodicity issues).
You can use the following model:
function unit_cell_terms(::Model"correlated_hopping"; J=1.0)
opsum = OpSum()
opsum += -J, "Cdag", 1, "C", 2, "C", 3, "Cdag", 4
opsum += -J, "C", 1, "Cdag", 2, "Cdag", 3, "C", 4
return opsum
end
Starting with a 4 site unitcell and and the initial state 1 2 2 1, the target energy per site is half the one of the XX model.

For the subspace expansion, I follow the method proposed by Hubig. The code can be found in the branch [idmrginfinitempo]. You can look at subspace_expansion.jl, l437 for the definition of the added term, and for example new_idmrg_infinitempo l49 for the insertion into iDMRG. I can send you a drawing of the procedure if you want.

For the environments: I have not yet thought about what to do for the VUMPS environment. It is not really a problem in iDMRG as you update an environment (though it does require some careful treatment given that the outer bonds have to be updated).
I guess one could try to perform some form of LU/LQ decomposition on the center ITensor.
There is an alternative (which I am not necessarily recommending), which is to perform the QR transformation while keeping the starting structure of the InfiniteMPOMatrix. In principle, one can sequentially (by column or line depending whether we go left or right) do the QR decomposition, which guarantees that the resulting MPO is triangular inferior. I wanted to explore that possibility at some point.

[JanReimers]

OK great, A number of things for me to think about.
Thanks
Jan