Obtaining transfer matrix using diagHermitian in iDMRG

ITensor C++ Questions
1

Hello,

I am studying a spin system using the idmrg.h file and trying to construct the transfer matrix to determine the correlation length of the system. I can construct the transfer matrix (bond dimension is around 50-80) with no issues. Once I pass the order 4 tensor to diagHermitian, the 1st eigenvalue is not 1 (In my case =0.403) but the state is normalized. Below is a snippet of how I construct the transfer matrix. Any help is appreciated, I can provide more details if needed.

`
auto wf_lower = prime(wf(0)*wf(1),“Site”);
auto wf_upper = dag(prime(wf(0)*wf(1)));

//Loop over wf sites to make lower and upper part of contraction diagram:
for(int i = 2;i <= N; i+=1)
{
    wf_lower *= prime(wf(i),"Site");
}
for(int i = 2;i <= N; i +=1)
{
    wf_upper *= dag(prime(wf(i)));
}

//Make Transfer Matrix

auto T = wf_lower*wf_upper;
print(T);

`

1 Like
2

Thanks for the question. One issue here may be that the transfer matrix is not Hermitian. Typically for this reason its eigenvalues are found using the Arnoldi algorithm for this reason.

We have an implementation of arnoldi in ITensor which you can find here:
<https://github.com/ITensor/ITensor/blob/45bf86103b232f6b0d9992137253945f44e2a2f3/itensor/iterativesolvers.h#L715>
It’s not documented but basically you pass an object A which implements a method A.product(ITensor const& input, ITensor & output) that internally multiplies an ITensor “input” by the transfer matrix and writes the result to the ITensor “output”. (This approach actually lets you avoid bond dimension to the power of 4 scaling, because you can contract the MPS tensors one by one onto the eigenvector.) You also pass a std::vector of ITensors corresponding to the initial guesses of the eigenvectors you want to find. The returned output are the largest magnitude eigenvalues and I believe the eigenvectors also get overwritten with the converged ones.

3

Hi Miles, thanks for the fast response.

I tried calculating norm(T - dag(swapTags(T,“0”,“1”))), where the term in the parentheses I believe is the equivalent of doing T - T^dag. The output is 0 so I assumed that the matrix was Hermitian.

Any tips on writing the object that can be passed to the Arnoldi algorithm/ any examples? I am not the best at C++ so anything helps!

Thanks again,
Sebastien

4

It may well be that the transfer matrix is coming out Hermitian in your case and I guess there’s nothing to prevent it from being. It’s more that it’s just not guaranteed in general so it could be a bit “brittle” of an assumption. If you do keep the norm(T - dag(swapTags(T,“0”,“1”))) check in your code then it could be a good way to catch if T ever failed to be Hermitian.

Looking back over your code, which I should have done more carefully, I think I may have spotted the issue. You are including wf(0) into your transfer matrix but wf(0) as output by the iDMRG code is the “center” tensor of the MPS. It is like the \Lambda in the \Gamma-\Lambda or Vidal canonical gauge of an MPS if you know that form. The key point is that wf(1) and wf(2) and so on are right-orthogonal and if I’m thinking about it correctly the transfer matrix should only be made of those.

Can you please try:

  • just leaving out wf(0)
  • checking again if T is Hermitian, in which case for now I think it’s fine to use diagHermitian

If those steps work, you may still want to eventually switch to the Arnoldi approach especially if you plan to study larger bond dimensions.

For a good example of how to implement a class that offers the .product interface, you can see
https://github.com/ITensor/ITensor/blob/v3/itensor/mps/localmpo.h
Basically you could implement a class like that with a new name (any name you like) and just leave out most of the code except a minimal constructor and the void product(ITensor const& input, ITensor & output) method with exactly that signature and it should work. It might be necessary to implement one or two other methods like .size() but you could try just .product first to see if that’s enough.

5

Here is what I ended up doing finally: I left out psi(0) and relabeled the last index to avoid the network contracting to a scalar. The final transfer matrix is Hermitian, norm(T - T^dag ) = 0 but the first and second eigenvalues are greater than 1. Here is the bit of code I used:

`
auto wf = psi;
auto lb1 = uniqueIndex(psi(1),psi(2),“Link”);
auto rb1 = uniqueIndex(psi(N),psi(N-1),“Link”);

//Relabel the last index do the netowrk doesnt contract:
auto new_ind = Index(dim(rb1),"l=1 to 8");
auto kron = delta(new_ind,rb1);

auto wf_lower = prime(wf(1),"Site");
auto wf_upper = dag(prime(wf(1)));

//Loop over wf sites to make lower and upper part of contraction diagram:
for(int i = 2;i < N; i+=1)
{
    wf_lower *= prime(wf(i),"Site");
}
for(int i = 2;i < N; i +=1)
{
    wf_upper *= dag(prime(wf(i)));
}

auto A8 = kron*wf(N);

wf_lower *= prime(A8,"Site");
print(wf_lower);
wf_upper *= dag(prime(A8));

//Make Transfer Matrix

auto T = wf_lower*wf_upper;
print(T);
printfln("\nNorm of T - T^dag:\n");
print(norm(T - dag(swapTags(T,"0","1")))); \\returns zero
printfln("\n");

`
I could try using Arnoldi if this wont work, I am unsure however why this is not working still.

6

Here are some key things you can do to debug this code:

  • a very important one is to “diagram it out”, drawing pictures of the contractions step-by-step and then pausing the code at each step (you can put PAUSE; to do this) then printing out wf_lower and wf_upper as they are being made.
  • one thing that looks suspicious to me is that you are never priming any “Link” indices. I wonder if some of the link indices are contracting between wf_upper and wf_lower in a way they are not supposed to
  • lastly, because T is supposed to be made of purely right-orthogonal MPS tensors, one property it rigorously must have is that the identity operator or matrix is a right eigenvalue of it with eigenvalue 1. So instead of finding this eigenvalue numerically, you can form this (supposed) right eigenvector and see if it is actually an eigenvector and with what eigenvalue

Best,
Miles

7

Sorry for the late reply,

I did draw the diagram out, but since the psi(i) matrices are periodic in the sense that the left bond of psi(1) is equal to the right bond of psi(N) I decided to relabel the indices which is why there are no primes. Since my unit cell is small (N=8) I decided to construct the transfer matrix by leaving the left and right bonds free with the site indices primed as I looped through the sites. I thought this was better as the tensor size would be of order 2^N * chi^2 versus the 2chi^4. No accidental contracting was done, I checked explicitly by printing everything out.
I checked quickly if delta(lb1,prime(lb1)) is an eigenvector of the transfer matrix. For N=8 and a bond dim of 60, the norm of Tdelta is ~7.9 while the norm of delta is 7.75. I tried increasing the bond dimension but it only makes it worse (bond_dim =70 , 80 makes makes the difference in norm around 4,5). If this applicable, dividing the norm of T*delta by norm(delta) gives 1.02 in the spirit of a Raleigh quotient which is close but larger than 1.
I am unsure if this check should always give me a close enough eigenvalue of 1 or it should be roughly close to 1 (here I assumed the norm would reflect the eigenvalue of ).

If the case is that for low bond dimension the approximation of the transfer matrix is bad, I will just move on to Arnoldi. Otherwise, it might be possible that I am not constructing T properly.

Best,
Seb

8

It’s good you tried all of those things.

I think as you know, one can show rigorously that if the MPS tensors making up the transfer matrix T are right-orthogonal, then T has the identity as an exact right eigenvector with eigenvalue exactly 1.0. So if you found that T*delta had a different norm than delta, that seems like a “red flag” right there and my bet would be on something being wrong with the transfer matrix itself. Note that the property about the identity being an exact right eigenvector would be exactly true for any bond dimension MPS, so if it’s not true there’s a more fundamental error happening.

Can you please numerically check one by one that the tensors making up the transfer matrix are right orthogonal? You could contract them together pairwise on their site and right link index and then print out the result to confirm that it’s the identity. If this works, then almost literally the same code would be the most efficient way to apply T to the identity on the right, so it would have to then work. Otherwise you will see that one of the pair contractions gives something not equal to the identity and you will have found your bug.

9

Hello,

I tried checking the orthogonality of the matrices and confirmed they are indeed right orthogonal. I increased the sweep count and all have right eigenvectors of 1. It took sometime to figure out why I was not and it turns out that the bare iDMRG file available on GitHub doubled the bond dimensions halfway through the chain on one of the bonds. So, when I was comparing norms I was slightly off by a factor of sqrt(2) somewhere. Now that the transfer matrix is confirmed to have its first right eigenvalue of 1 I will start diagonalizing this. Thanks so much for the help!

10

Glad to hear you made some definite progress with it!

1 Like
11

Just wanted to follow up quickly: I was able to diagonalize the matrix by using eigen(), the tricky part was the priming and naming convention of the indices in the documentation.

Thanks a ton!

12

Really glad to hear it - thakns for the update!

1 Like
13

Sorry to ask something more. I find the diagHermitian decomposing Transfer Matrix along the direction in default as the upper of my picture tells. But what eigen values we need is the way in lower picture.

20230106033146845
20230106033146845698×395 15.3 KB

I met the same problem, so I wanna ask how to solve it.
Any help is welcome. Thanks

14

According to the eigen documentation on the github: ITensor/decomp.h at 98fc7022eb812fe68640e662f4051dfd2348e4e9 · ITensor/ITensor · GitHub. You just need to change how you label the indices of your transfer matrix, so in the digram in the doc it would just be changing I and I’ prime indices to be I and J instead for example.

Hope that helps!

15

Thank you!
You answer is very useful.
I have some questions in detail, What do you mean by changing indices? Is it means changing the tags? I have tried but it comes to an error.
Tranfer Matrix

findInds(TM, “”) =
(100|id=137|Lup,Link)
1: 1 QN({“L8”,-6},{“Sz”,6})
2: 10 QN({“L8”,0},{“Sz”,4})
3: 1 QN({“L8”,6},{“Sz”,6})
4: 10 QN({“L8”,6},{“Sz”,2})
5: 2 QN({“L8”,12},{“Sz”,0})
6: 10 QN({“L8”,-6},{“Sz”,2})
7: 23 QN({“L8”,0},{“Sz”,0})
8: 12 QN({“L8”,6},{“Sz”,-2})
9: 1 QN({“L8”,12},{“Sz”,-4})
10: 1 QN({“L8”,-12},{“Sz”,0})
11: 11 QN({“L8”,-6},{“Sz”,-2})
12: 12 QN({“L8”,0},{“Sz”,-4})
13: 2 QN({“L8”,6},{“Sz”,-6})
14: 1 QN({“L8”,-12},{“Sz”,-4})
15: 2 QN({“L8”,-6},{“Sz”,-6})
16: 1 QN({“L8”,0},{“Sz”,-8})
(100|id=137|Ldn,Link)
1: 1 QN({“L8”,-6},{“Sz”,6})
2: 10 QN({“L8”,0},{“Sz”,4})
3: 1 QN({“L8”,6},{“Sz”,6})
4: 10 QN({“L8”,6},{“Sz”,2})
5: 2 QN({“L8”,12},{“Sz”,0})
6: 10 QN({“L8”,-6},{“Sz”,2})
7: 23 QN({“L8”,0},{“Sz”,0})
8: 12 QN({“L8”,6},{“Sz”,-2})
9: 1 QN({“L8”,12},{“Sz”,-4})
10: 1 QN({“L8”,-12},{“Sz”,0})
11: 11 QN({“L8”,-6},{“Sz”,-2})
12: 12 QN({“L8”,0},{“Sz”,-4})
13: 2 QN({“L8”,6},{“Sz”,-6})
14: 1 QN({“L8”,-12},{“Sz”,-4})
15: 2 QN({“L8”,-6},{“Sz”,-6})
16: 1 QN({“L8”,0},{“Sz”,-8})
(100|id=415|Lup,Link)’
1: 10 QN({“L8”,6},{“Sz”,2})
2: 1 QN({“L8”,12},{“Sz”,0})
3: 24 QN({“L8”,0},{“Sz”,0})
4: 12 QN({“L8”,6},{“Sz”,-2})
5: 1 QN({“L8”,12},{“Sz”,-4})
6: 11 QN({“L8”,-6},{“Sz”,-2})
7: 12 QN({“L8”,0},{“Sz”,-4})
8: 2 QN({“L8”,6},{“Sz”,-6})
9: 1 QN({“L8”,-12},{“Sz”,-4})
10: 2 QN({“L8”,-6},{“Sz”,-6})
11: 1 QN({“L8”,0},{“Sz”,-8})
12: 1 QN({“L8”,6},{“Sz”,6})
13: 10 QN({“L8”,0},{“Sz”,4})
14: 10 QN({“L8”,-6},{“Sz”,2})
15: 1 QN({“L8”,-12},{“Sz”,0})
16: 1 QN({“L8”,-6},{“Sz”,6})
(100|id=415|Ldn,Link)’
1: 10 QN({“L8”,6},{“Sz”,2})
2: 1 QN({“L8”,12},{“Sz”,0})
3: 24 QN({“L8”,0},{“Sz”,0})
4: 12 QN({“L8”,6},{“Sz”,-2})
5: 1 QN({“L8”,12},{“Sz”,-4})
6: 11 QN({“L8”,-6},{“Sz”,-2})
7: 12 QN({“L8”,0},{“Sz”,-4})
8: 2 QN({“L8”,6},{“Sz”,-6})
9: 1 QN({“L8”,-12},{“Sz”,-4})
10: 2 QN({“L8”,-6},{“Sz”,-6})
11: 1 QN({“L8”,0},{“Sz”,-8})
12: 1 QN({“L8”,6},{“Sz”,6})
13: 10 QN({“L8”,0},{“Sz”,4})
14: 10 QN({“L8”,-6},{“Sz”,2})
15: 1 QN({“L8”,-12},{“Sz”,0})
16: 1 QN({“L8”,-6},{“Sz”,6})

when I use eigen, it shows Segmentation fault.
When I use diagHermitian, it tells

From line 743, file decomp.cc

Input tensor to diagHermitian should have pairs of indices with equally spaced prime levels

16

The link I posted above, if you scroll slightly up, contains a diagram that shows how the transfer matrix contracts with the eigenvector to give the eigenvalue equation.
What you need to do is match this convention by either relabeling the indices, which can be done by creating new indices and using a Kronecker delta to swap the indices around.

So for example, if your transfer matrix is T and has indices i,j,ik and l, the indices need to be changed to i,j,i’ and j’ due to this convention. The choice of what k and l become, either i’ or j’, is what determines the correct eigenvalue equation.