solve an eigenvalue problem of non-hermitian tensor

Question

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.

Answer 1

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.

Comments

Thanks. It helps me.

Answer 2

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

Comments

Thanks.

Is there a function to get the first eigenvalue/eigenvector (or first k eigenvalues/vectors) ?

---

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

---

I need it in ITensor V3. Thank you very much.

---

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.

Answer 3

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

Comments

Thanks for your update.

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.

---

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.

---

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

---

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

---

Hi Matt,

I wanted to ask you if the bug that you mentioned in the 'eigen' function is also present in ITensor v2?

I want to use it to compute the exponential of a two-site non-hermitian gate, but at the moment I cannot update to ITensor v3 to use the version you recently updated.

Thanks,
Catalin

---

Hi Catalin,

Unfortunately, yes, it did seem to be a problem in ITensor v2 as well. However, I just added a fix to ITensor v2 as well, please pull the latest "v2" branch of ITensor to try it out.

Cheers,
Matt

---

Hi Matt,

Thank you for you quick reply!

Are there any restrictions of using IQTensors in 'eigen'? Because I think there might be an another small bug in 'eigen' for the case of IQTensors, if my understanding of what the 'Tc' tensor is correct.  

Is the aim to contract all the pairs of site indices from 'T', so you remain with the 2 link indices in 'Tc'? If this is correct, than I think one would get a "Mistmatched IQINdex arrows" error message in computing 'Tc', as the 'prime' function doesn't change the arrow directions of the indices in 'comb' and an extra 'dag' would be necessary.

In the following I copied a piece of my code and the output, with the error message that I get:

   PrintData(Hbond);
  
   auto colinds = std::vector<IQIndex>{};
    for(auto& I : Hbond.inds())
        {
        if(I.primeLevel() == 0) colinds.push_back(I);
        }
    auto comb = combiner(std::move(colinds));
    
    PrintData(comb);
    PrintData(prime(comb));
    PrintData(dag(prime(comb)));
    
    auto Tc = dag(prime(comb)) * Hbond * comb;
    PrintData(Tc);
    
    auto Tc2 = prime(comb) * Hbond * comb;
    PrintData(Tc2);

Output:
"Hbond =
/--------------IQTensor--------------
r=4 div=QN(0) log(scale)=2.15874
IQIndex("Boson site=",3,Site|668) <In>
  ("b0 2",1,Site|479) QN(0)
  ("b1 2",1,Site|699) QN(1)
  ("b2 2",1,Site|441) QN(2)

IQIndex("Boson site=",3,Site|668)' <Out>
  ("b0 2",1,Site|479)' QN(0)
  ("b1 2",1,Site|699)' QN(1)
  ("b2 2",1,Site|441)' QN(2)

IQIndex("Photon site",6,Site|148) <In>
  ("p0 1",1,Site|281) QN(0)
  ("p1 1",1,Site|197) QN(0)
  ("p2 1",1,Site|25) QN(0)
  ("p3 1",1,Site|488) QN(0)
  ("p4 1",1,Site|146) QN(0)
  ("p5 1",1,Site|315) QN(0)

IQIndex("Photon site",6,Site|148)' <Out>
  ("p0 1",1,Site|281)' QN(0)
  ("p1 1",1,Site|197)' QN(0)
  ("p2 1",1,Site|25)' QN(0)
  ("p3 1",1,Site|488)' QN(0)
  ("p4 1",1,Site|146)' QN(0)
  ("p5 1",1,Site|315)' QN(0)

|-- Data -------
QDense Cplx {108 blocks; data size 108}
...I skip the values...
\------------------------------------

comb =
/--------------IQTensor--------------
r=3 div=QN() log(scale)=0
IQIndex("cmb",18,Link|814) <In>
  ("c0",6,Link|603) QN(0)
  ("c1",6,Link|583) QN(1)
  ("c2",6,Link|899) QN(2)

IQIndex("Boson site=",3,Site|668) <Out>
  ("b0 2",1,Site|479) QN(0)
  ("b1 2",1,Site|699) QN(1)
  ("b2 2",1,Site|441) QN(2)

IQIndex("Photon site",6,Site|148) <Out>
  ("p0 1",1,Site|281) QN(0)
  ("p1 1",1,Site|197) QN(0)
  ("p2 1",1,Site|25) QN(0)
  ("p3 1",1,Site|488) QN(0)
  ("p4 1",1,Site|146) QN(0)
  ("p5 1",1,Site|315) QN(0)

|-- Data -------
QCombiner
\------------------------------------


prime(comb) =
/--------------IQTensor--------------
r=3 div=QN() log(scale)=0
IQIndex("cmb",18,Link|814)' <In>
  ("c0",6,Link|603)' QN(0)
  ("c1",6,Link|583)' QN(1)
  ("c2",6,Link|899)' QN(2)

IQIndex("Boson site=",3,Site|668)' <Out>
  ("b0 2",1,Site|479)' QN(0)
  ("b1 2",1,Site|699)' QN(1)
  ("b2 2",1,Site|441)' QN(2)

IQIndex("Photon site",6,Site|148)' <Out>
  ("p0 1",1,Site|281)' QN(0)
  ("p1 1",1,Site|197)' QN(0)
  ("p2 1",1,Site|25)' QN(0)
  ("p3 1",1,Site|488)' QN(0)
  ("p4 1",1,Site|146)' QN(0)
  ("p5 1",1,Site|315)' QN(0)

|-- Data -------
QCombiner
\------------------------------------


dag(prime(comb)) =
/--------------IQTensor--------------
r=3 div=QN() log(scale)=0
IQIndex("cmb",18,Link|814)' <Out>
  ("c0",6,Link|603)' QN(0)
  ("c1",6,Link|583)' QN(1)
  ("c2",6,Link|899)' QN(2)

IQIndex("Boson site=",3,Site|668)' <In>
  ("b0 2",1,Site|479)' QN(0)
  ("b1 2",1,Site|699)' QN(1)
  ("b2 2",1,Site|441)' QN(2)

IQIndex("Photon site",6,Site|148)' <In>
  ("p0 1",1,Site|281)' QN(0)
  ("p1 1",1,Site|197)' QN(0)
  ("p2 1",1,Site|25)' QN(0)
  ("p3 1",1,Site|488)' QN(0)
  ("p4 1",1,Site|146)' QN(0)
  ("p5 1",1,Site|315)' QN(0)

|-- Data -------
QCombiner
\------------------------------------

Tc =
/--------------IQTensor--------------
r=2 div=QN(0) log(scale)=2.15874
IQIndex("cmb",18,Link|814) <In>
  ("c0",6,Link|603) QN(0)
  ("c1",6,Link|583) QN(1)
  ("c2",6,Link|899) QN(2)

IQIndex("cmb",18,Link|814)' <Out>
  ("c0",6,Link|603)' QN(0)
  ("c1",6,Link|583)' QN(1)
  ("c2",6,Link|899)' QN(2)

|-- Data -------
QDense Cplx {3 blocks; data size 108}
...I skip the values...
\------------------------------------

----------------------------------------
IQIndexSet 1 =
IQIndex("cmb",18,Link|814)' <In>
  ("c0",6,Link|603)' QN(0)
  ("c1",6,Link|583)' QN(1)
  ("c2",6,Link|899)' QN(2)

IQIndex("Boson site=",3,Site|668)' <Out>
  ("b0 2",1,Site|479)' QN(0)
  ("b1 2",1,Site|699)' QN(1)
  ("b2 2",1,Site|441)' QN(2)

IQIndex("Photon site",6,Site|148)' <Out>
  ("p0 1",1,Site|281)' QN(0)
  ("p1 1",1,Site|197)' QN(0)
  ("p2 1",1,Site|25)' QN(0)
  ("p3 1",1,Site|488)' QN(0)
  ("p4 1",1,Site|146)' QN(0)
  ("p5 1",1,Site|315)' QN(0)

----------------------------------------
IQIndexSet 2 =
IQIndex("Boson site=",3,Site|668) <In>
  ("b0 2",1,Site|479) QN(0)
  ("b1 2",1,Site|699) QN(1)
  ("b2 2",1,Site|441) QN(2)

IQIndex("Boson site=",3,Site|668)' <Out>
  ("b0 2",1,Site|479)' QN(0)
  ("b1 2",1,Site|699)' QN(1)
  ("b2 2",1,Site|441)' QN(2)

IQIndex("Photon site",6,Site|148) <In>
  ("p0 1",1,Site|281) QN(0)
  ("p1 1",1,Site|197) QN(0)
  ("p2 1",1,Site|25) QN(0)
  ("p3 1",1,Site|488) QN(0)
  ("p4 1",1,Site|146) QN(0)
  ("p5 1",1,Site|315) QN(0)

IQIndex("Photon site",6,Site|148)' <Out>
  ("p0 1",1,Site|281)' QN(0)
  ("p1 1",1,Site|197)' QN(0)
  ("p2 1",1,Site|25)' QN(0)
  ("p3 1",1,Site|488)' QN(0)
  ("p4 1",1,Site|146)' QN(0)
  ("p5 1",1,Site|315)' QN(0)

----------------------------------------
Mismatched IQIndex from set 1 IQIndex("Boson site=",3,Site|668)' <Out>
  ("b0 2",1,Site|479)' QN(0)
  ("b1 2",1,Site|699)' QN(1)
  ("b2 2",1,Site|441)' QN(2)

Mismatched IQIndex from set 2 IQIndex("Boson site=",3,Site|668)' <Out>
  ("b0 2",1,Site|479)' QN(0)
  ("b1 2",1,Site|699)' QN(1)
  ("b2 2",1,Site|441)' QN(2)


Mismatched IQIndex arrows"

I mention that I use a custom SiteSet.

Best,
Catalin

---

Actually, `eigen` is not currently fully implemented for IQTensors (it should fail later on anyway). I can make a better error message about that in v2.

If you need eigenvectors and eigenvalues of a non-Hermitian matrix, you would need to use Arnoldi for now.