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);

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.

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!

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.

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.

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

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.

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.

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!

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.