std::bad_alloc in Transfer Matrix calculation

Question

Hi all,

I am studying the ground state of a Heisenberg antiferromagnets on a 2D lattice (one periodic and one open), and I want to calculate the transfer matrix which is defined like this:

//got the ground state wave function psi through dmrg;
psi.position(N/2);
ITensor wf = psi(N/2 - Ny + 1);
for (int i = 2; i <= Ny; i++){
wf = wf * psi(N/2 - Ny + i);
}

then contract the wf with its own copy but keep the "Link" indices open,

wf *= prime(wf, "Link");

and I should get the transfer matrix where I want to calculate its largest eigenvalues and eigenvectors. However the second step throws an error of "std::bad_alloc" because the "Link" index dimension is too large ( around 3000), if I try 300 as "Link" index dimension it works well.

I have found a similar post http://itensor.org/support/1254/std-bad_alloc-in-spin-polarized-system where the solution is to avoid such large dimensions, but I can not think of a way to do that in my problem. Will it be better if I use a larger memory (like 512GB)? Any help would be much appreciated.

Best,

Yixuan

Answer 1

Hi Yixuan,
Yes, the reason for this error message is that your system has run out of memory / RAM. So one solution is definitely to try a system with more RAM if possible.

But another solution could be to write any large tensors you are not using at a particular step to disk. You can estimate the size of various tensors by multiplying the dimension of their indices together. Using this estimate, you might find that individual MPS tensors themselves aren't using a lot of memory, but tensors formed from the MPS can use a very large amount. (Or you might find that the MPS tensors alone are using a lot of memory; it depends on the details.)

We have some automatic facilities to write parts of an MPS to disk. But these may not work ideally for a transfer matrix calculation and are designed mostly to be used inside a DMRG calculation.

But you can even do writing to disk manually by just writing a tensor to a particular file with a good file name, then overwriting the tensor with an empty tensor to 'free' it from memory. Later, when you need the tensor again you can read it back in from the file and overwrite the read-in tensor onto your empty tensor to restore it.

If you decide to follow one of the three approaches above and run into issues, please post a comment below or a new question.

Best regards,
Miles

Comments

Dear Miles,

Thank you for the prompt reply. Now that I have estimated roughly the size of that tensor, it seems like an impossible task because 3000^4 numbers * 64 Bytes per number ~ 4000 TB, although I am using conserved quantum number so the matrix should be sparse. I will let you know if I have questions about writing a tensor to disk.

Best,

Yixuan

---

Yes, if you are trying to make a tensor with 4 indices of size 3000 I'm afraid it is rather hopeless on current machines.

*However*, if your goal is to study the MPS transfer matrix, then what people usually do is not to actually make this matrix. Instead, they try to find its dominant eigenvectors and eigenvalues (just two of these are enough to estimate the correlation length, for example). This is done by using an iterative eigensolver method such as Arnoldi where at each iteration, you multiply the current guess for the eigenvector by the transfer matrix. This multiplication can be done by contracting one tensor at a time such that no intermediate tensor ever has more than 3 large indices. In this way you can reduce the memory usage by a factor of 3000 in your case.

Maybe you already knew all that, though, but just wanted to make sure.

Best regards,
Miles

---

That is exactly what I want! Thank you so much!! Actually I was looking at the Arnoldi part added in ITensor v3 recently (https://github.com/ITensor/ITensor/blob/v3/unittest/iterativesolvers_test.cc#L252) and thought it was meant for the non-Hermitian matrix, but you are right this should be the best way to get the dominant eigenvector.

Best regards,
Yixuan

---

Hi Yixuan, great - glad that's what you were looking for. It's definitely the best way to go about extracting transfer matrix information.

I believe that in general the transfer matrix is not Hermitian, though it can be in certain cases. Best not to assume it is unless you can prove that it must be.

Miles

---

Dear Miles,

Thanks for the explanation. I have tried and successfully run the Arnoldi to find the dominant eigenvector of a simple random ITensor with (i,j,k...) and (i',j',k'...) but I run into an error when I deal with the transfer matrix which is an ITensor pulled from the ground state MPS (with quantum number). I need to change the "Link" index of it in order to look like (i,j) and (i',j') to use Arnoldi, so I use the replaceInds(ITensor T, IndexSet is1, IndexSet is2) function:

auto lf = leftLinkIndex(psi, XX - Ny + 1);
auto rt = rightLinkIndex(psi, XX);
wf.replaceInds({rt}, {dag(prime(lf))});

and the error is "Diagonal elements of QDiag ITensor would have inconsistent divergence". I  know I could also replace the index by the delta function

wf *= dag(delta(rt, prime(lf)));

but unfortunately that gives out the same error.  I thought the diagonal tensor made from delta() should always have divergence of 0 because it only has non-zero diagonal elements. Why would it have inconsistent divergence with the original ITensor?

---

Hi, so I'm not sure why you are getting this error. I think I would need to see a printout of the wf tensor, its indices, and its flux to know what is happening.

It could be a bug in the replaceInds function. If so, one alternative which could be easier for you and more transparent to debug is to make delta tensors yourself, and contract them with wf to replace the indices (I think this is all replaceInds does anyway). Then you could be sure the delta tensors are properly defined, *especially* regarding the index arrow directions. I suspect that could be the source of the problem, namely that replaceInds might not be handling arrows properly either in general or for your case.

---

Dear Miles,

Thanks for the suggestions, I have tested it using the delta function and tried to replace the indices but there is still the error saying "Diagonal elements of QDiag ITensor would have inconsistent divergence". I should have pointed out before that the error message only appears in the debug mode, but that is still a little bit concerning. Here is the printout of wf,

 ITensor ord=14:
(2|id=355|n=97,Site,S=1/2) <Out>
  1: 1 QN({"Sz",1})
  2: 1 QN({"Sz",-1})
(100|id=705|l=96,Link) <Out>
  1: 5 QN({"Sz",4})
  2: 26 QN({"Sz",2})
  3: 37 QN({"Sz",0})
  4: 26 QN({"Sz",-2})
  5: 6 QN({"Sz",-4})
(2|id=984|n=98,Site,S=1/2) <Out>
  1: 1 QN({"Sz",1})
  2: 1 QN({"Sz",-1})
(2|id=494|n=99,Site,S=1/2) <Out>
  1: 1 QN({"Sz",1})
  2: 1 QN({"Sz",-1})
(2|id=351|Site,n=100,S=1/2) <Out>
  1: 1 QN({"Sz",1})
  2: 1 QN({"Sz",-1})
(2|id=39|Site,n=101,S=1/2) <Out>
  1: 1 QN({"Sz",1})
  2: 1 QN({"Sz",-1})
(2|id=302|Site,S=1/2,n=102) <Out>
  1: 1 QN({"Sz",1})
  2: 1 QN({"Sz",-1})
(2|id=902|Site,S=1/2,n=103) <Out>
  1: 1 QN({"Sz",1})
  2: 1 QN({"Sz",-1})
(2|id=677|Site,S=1/2,n=104) <Out>
  1: 1 QN({"Sz",1})
  2: 1 QN({"Sz",-1})
(2|id=434|Site,S=1/2,n=105) <Out>
  1: 1 QN({"Sz",1})
  2: 1 QN({"Sz",-1})
(2|id=693|Site,S=1/2,n=106) <Out>
  1: 1 QN({"Sz",1})
  2: 1 QN({"Sz",-1})
(2|id=22|Site,S=1/2,n=107) <Out>
  1: 1 QN({"Sz",1})
  2: 1 QN({"Sz",-1})
(100|id=525|Link,l=108) <Out>
  1: 6 QN({"Sz",4})
  2: 26 QN({"Sz",2})
  3: 37 QN({"Sz",0})
  4: 25 QN({"Sz",-2})
  5: 6 QN({"Sz",-4})
(2|id=568|Site,S=1/2,n=108) <Out>
  1: 1 QN({"Sz",1})
  2: 1 QN({"Sz",-1})
{norm=1.00 (QDense Cplx)}


where lf is the Link index with l=96, and rt is the Index with l=108. I think the error is because that the bond dimensions of certain quantum number sector does not match with the one in the other index, if I do wf.replaceInds({rt}, {prime(rt)}); then there is no error.  I was wondering if there is any way to fix the bond dimensions of every quantum number sector for those two Indices the same during the DMRG process as a constrained sweeps.

Best,

Yixuan

---

Hi Yixuan,
Is your transfer matrix from the unit cell of an infinite system calculation? Or is coming from a finite-size system? Thanks -

Miles

---

It is from a finite-size 2D system.

Best,

Yixuan

---

I see, so in general I'm not sure if this idea of taking the transfer matrix from a finite system is viable, especially when conserving quantum numbers. The issue is that while the index on the left-hand-side of the transfer matrix *might* be the same size as on the right-hand-side, there's no generic reason it has to be so it's kind of a brittle assumption. Even more brittle is assuming the quantum numbers will match on the left and right hand sides in terms of which sectors are present and which sizes they each have. It's hard to ensure that in general, unless one does special steps to try to ensure it (I'm not sure what these steps are, just saying they could exist maybe).

So I would mainly recommend using transfer matrix methods only for infinite systems, unless one has a way to ensure they will also work for finite systems.

---

Hi Miles,

Thank you for the advice, now that I think more about it, it is hard to keep those two "Link" indices of the transfer matrix the same (differ by a prime) because they changed during the sweeps, and because I have an orthogonal center at a certain site, the transfer matrix I take from a finite system can not be considered as transnational invariant (even ignoring the finite size effect). I will try the infinite system for transfer matrix calculation.

Best,

Yixuan