# solve an eigenvalue problem of non-hermitian tensor

Hi,

I want to calculate the eigenvalue of a non-hermitian tensor. After some searching of document and Q&A, I did not find a direct function to solve. I am not familiar with either LAPACK or ITensor Wrapper for LAPACK.

Because the tensor I want to solve is not high dimension (12 indices, each has 2 dims), I can "combiner" 6 of them together to get a 64*64 tensor.

My question is how I convert the tensor type to matrix type so I can use package like in C++.

Or can someone provide me better way to solve this question?

Thanks.

Hi,

I've write a demo to solve arbitrary matrix eigenproblem by using Eigen3, see github page.

you need first installing Eigen3 (it's header files and you can put it anywhere), and edit Makefile, then make and test it yourself.

commented by (170 points)
Thanks. It helps me.

Which eigenvalue are you interested in? Do you want all of the eigenvalues/eigenvectors, or just some dominant ones?

To get all of the eigenvalues/eigenvectors, you could try the ITensor eigen function:

https://github.com/ITensor/ITensor/blob/98fc7022eb812fe68640e662f4051dfd2348e4e9/itensor/decomp.h#L237

Unfortunately it looks like it is not documented yet, I think I was planning on testing it more but you can try it out (just make sure to check that you get proper eigenvalues and eigenvectors).

This function should work as long as your ITensor has pairs of indices i,j,k,...,i',j',k'.

Cheers,
Matt

commented by (170 points)
Thanks.

Is there a function to get the first eigenvalue/eigenvector (or first k eigenvalues/vectors) ?
commented by (4.8k points)
If you are using ITensor V2, you could try the arnoldi function:

https://github.com/ITensor/ITensor/blob/v2/itensor/iterativesolvers.h#L683

It is also not documented yet and not very well tested, so again you should double check to make sure the eigenvalues and eigenvectors you get are correct. If you need it in ITensor V3, let me know and it shouldn't be too hard to port it.

Cheers,
Matt
commented by (170 points)
I need it in ITensor V3. Thank you very much.
commented by (4.8k points)
Ok, I will try to add it next week (I am away at a conference this week).

For now you could use the eigen ITensor function if your tensor is small enough.

I ported arnoldi to ITensor V3 to find dominant eigenvectors and eigenvalues of non-Hermitian sparse matrices. You can take a look at these tests:

https://github.com/ITensor/ITensor/blob/v3/unittest/iterativesolvers_test.cc#L252

to get an idea for how it is used.

Cheers,
Matt

commented by (170 points)

I have read the test code, and run it myself. Because I am not familiar with the LocalMPO class, can I ask some more questions?

From the test code, you first assign an array of Heisenberg Hamiltonian. Then:

LocalMPO PH(H);

psi.position(Nc);
PH.position(Nc,psi);

auto x = psi.A(Nc) * psi.A(Nc+1);
x.randomize();

auto lambda = arnoldi(PH,x,{"MaxIter",20,"ErrGoal",1e-14,"DebugLevel",0,"WhichEig","LargestMagnitude"});
auto PHx = x;
PH.product(x,PHx);

What you have done seems to project the Hamiltonian into the array at sites "Nc" and "Nc+1", then randomize a vector "x" with index of site "Nc" and "Nc+1".

Then the "arnoldi" function returns the first eigenvalue of PH, and assign the eigenvector to x.

Can I apply the "arnoldi" to general tensor "A" with indices i,j,k,i',j',k' other than a LocalMPO class PH, like the "eigen" function.

Or can you recommend me some reference about LocalMPO class?

Thanks.
commented by (170 points)
In my situation, my Matrix is not with high dim, so I think "eigen" is efficient enough for me. BTW, is there a "named arguments" list for "eigen" function? For my case, I want to keep the largest eigenvalue from "eigen". Thanks.
commented ago by (4.8k points)
In order to use arnoldi, you need to define your own class that can be passed into the function. The only requirements are that it defines the functions .size() and .product() which give the size of the matrix (i.e. m for an m x m matrix) and the product of the operator with an ITensor. Here is an example code for getting the eigenvector with the largest magnitude eigenvalue:

#include "itensor/all.h"

using namespace itensor;

class ITensorMap
{
ITensor const& A_;
mutable long size_;
public:

ITensorMap(ITensor const& A)
: A_(A)
{
size_ = 1;
for(auto& I : A_.inds())
{
if(I.primeLevel() == 0)
size_ *= dim(I);
}
}

void
product(ITensor const& x, ITensor& b) const
{
b = A_*x;
b.noPrime();
}

long
size() const
{
return size_;
}

};

int
main()
{
auto i = Index(2,"i");
auto j = Index(2,"j");
auto k = Index(2,"k");
auto A = randomITensor(i,j,k,prime(i),prime(j),prime(k));

// Using the ITensorMap class
// to make an object that can be passed
// to arnoldi
auto AM = ITensorMap(A);

PrintData(AM.size());

// Random starting vector for Arnoldi
auto x = randomITensor(i,j,k);

auto lambda = arnoldi(AM,x,{"ErrGoal=",1E-14,"MaxIter=",20,"MaxRestart=",5});

// Check it is an eigenvector
PrintData(lambda);
PrintData(norm((A*x).noPrime() - lambda*x));

return 0;
}

Unfortunately I have been seeing problems when I try to calculate more than one eigenvalue/eigenvector (you can try it by passing a vector of ITensors as the starting state x).

The eigen function doesn't have any named arguments (other than setting the tag of the newly created index). In order to get only a particular eigenvector, you will have to figure out which one you want (i.e. by searching the eigenvalues and finding the largest one) and projecting the tensor of eigenvectors onto that index (for example by using setElt(d=1) if the new index is d and you want the first eigenvector).

I also seem to be having trouble with the eigen function, we will have to look into it to make sure it is working correctly (I believe it is calculating the correct eigenvalues and eigenvectors, but it may be outputting it in a strange way).

Let me know if the arnoldi code above is sufficient for what you need. You can play around with increasing "MaxIter" and "MaxRestart" to try to get better convergence.

Cheers,
Matt
commented ago by (4.8k points)
I pushed a change to ITensor v3 that should fix the bug I was seeing in eigen, so you can try that as an alternative to arnoldi for finding eigenvectors. To get a single eigenvector, you can try something like the following:

auto i = Index(2,"i"),
j = Index(2,"j"),
k = Index(2,"k");

auto A = randomITensor(i,j,k,prime(i),prime(j),prime(k));

auto [V,D] = eigen(A,{"Tags=","my_tag"});
auto d = commonIndex(V,D);

PrintData(norm(A*V - prime(V)*D));

// Grab the first eigenvector.
//
// Note that if you want the dominant eigenvector,
// search through the diagonal elements
// of D to find the position of the maximum
// eigenvalue. The eigenvectors are not
// necessarily ordered.

auto v1 = V*setElt(d=1);

auto lambda1 = D.eltC(1,1);

PrintData(norm((A*v1).noPrime() - lambda1*v1));

Please let me know if you have any problems with either arnoldi or eigen.

Cheers,
Matt