error in calling arnoldi function

[shencebebetterme]

'rom this post I can see how to use the arnoldi function to get the largest eigenvalue of a square matrix. However, if I attempt to calculate the first k eigenvalues (k>1) and feed the arnoldi function with a vector of ITensor rather than a single ITensor, the compiler throws a long template error.

The minimal codes to reproduce the error is

#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");
  ITensor A =randomITensor(i,prime(i));

  auto AM = ITensorMap(A);

  ITensor x = randomITensor(i);
  ITensor y = randomITensor(i);
  std::vector<ITensor> vec = {x,y};

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

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

  return 0;
  }

In the post MattFishman said he also encountered a problem when calculating more than one eigenvalues/eigenvectors, it seems this issue hasn't been fully solved.

[mtfishman]

Hi,

It seems like the error you are seeing is when you are checking the eigenvector, since you are doing lambda * x but lambda is now a vector of the eigenvalues. This code runs:

#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 A = randomITensor(i,prime(i));

  auto AM = ITensorMap(A);

  ITensor x = randomITensor(i);
  ITensor y = randomITensor(i);
  std::vector<ITensor> vec = {x,y};

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

  PrintData(norm((A*vec[0]).noPrime() - lambda[0]*vec[0]));
  PrintData(norm((A*vec[1]).noPrime() - lambda[1]*vec[1]));

  // Compare to ED
  auto [U,D] = eigen(A);

  auto l = commonIndex(U, D);

  PrintData(lambda[0]);
  PrintData(lambda[1]);

  PrintData(D.elt(1,1));
  PrintData(D.elt(2,2));

  PrintData(vec[0]);
  PrintData(vec[1]);

  PrintData(U*setElt(l=1));  
  PrintData(U*setElt(l=2));

  return 0;
  }

However I see something like:

norm((A*vec[0]).noPrime() - lambda[0]*vec[0]) = 1.11022e-16
norm((A*vec[1]).noPrime() - lambda[1]*vec[1]) = 0.358771
lambda[0] = (0.598975,0)
lambda[1] = (0.240204,0)
D.elt(1,1) = 0.598975
D.elt(2,2) = -0.198257
vec[0] = 
ITensor ord=1: (dim=2|id=283|"i") 
{norm=1.00 (Dense Real)}
(1) 0.6263638
(2) 0.7795309

vec[1] = 
ITensor ord=1: (dim=2|id=283|"i") 
{norm=1.00 (Dense Real)}
(1) 0.6263638
(2) 0.7795309

U*setElt(l=1) = 
ITensor ord=1: (dim=2|id=283|"i") 
{norm=1.00 (Dense Real)}
(1) 0.6263638
(2) 0.7795309

U*setElt(l=2) = 
ITensor ord=1: (dim=2|id=283|"i") 
{norm=1.00 (Dense Real)}
(1) -0.571766
(2) 0.8204164

so the second eigenvector/eigenvalue is not being calculated correctly.

I won't have much time to look into this problem. However, if someone really needs multiple dominant eigenvalues/eigenvectors of a non-Hermitian operator, they could try first calculating the dominant left and right eigenvalues/eigenvectors, and then project out the dominant eigenspace and calculate the next eigenvalue/eigenvector, so schematically:

  auto i = Index(2,"i");
  ITensor A =randomITensor(prime(i), i);

  auto v1r = randomITensor(i);
  auto lambda1r = arnoldi(ITensorMap(A),v1r,{"ErrGoal=",1E-14,"MaxIter=",20,"MaxRestart=",5});

  auto v1l = randomITensor(i);
  auto lambda1l = arnoldi(ITensorMap(swapPrime(A,0,1)),v1l,{"ErrGoal=",1E-14,"MaxIter=",20,"MaxRestart=",5});

  auto v2 = randomITensor(i);
  auto lambda2 = arnoldi(ITensorMap(A-lambda1r*prime(v1r)*v1l),v2,{"ErrGoal=",1E-14,"MaxIter=",20,"MaxRestart=",5});

  PrintData(lambda1l);
  PrintData(lambda1r);
  PrintData(lambda2);

  PrintData(norm((A*v1r).noPrime() - lambda1r*v1r));
  PrintData(norm((A*v2).noPrime() - lambda2*v2));

In a real application, you would want to make a special ITensor map object that projects out the dominant eigenspace more efficiently.

Cheers,
Matt