fast diagonalisation in itensor

ITensor C++ Questions
1

Dear all,

I recently noticed that diagPosSemiDef it is quite slow for sizes >= order 2^11, even with truncation of order 10^-4 (I am diagonalizing a reduced density matrix, written as an ITensor with (combined) indices s, prime (s)).
I guess one of the reason is that it also compute the eigenvectors, but what to do if we only care about the leading eigenvalues of the spectrum? In principle, in this situation, a standard diagonalization routine should be easily able to access matrices of size 2^18 - 2^20. Do you have any suggestions, besides moving the ITensor to a matrix and then using lapack routines or something (which I guess it would be very expensive too).
Thanks a lot!
Jacopo.

2

Hi Jacopo,
If you only require the leading eigenvalues of a large, Hermitian matrix, then a good option can be to use iterative linear algebra methods, based on Krylov-subspace ideas. This is how DMRG solves for the dominant or leading eigenvector in the “core” step of the algorithm.

The relevant method implemented in ITensor is the Davidson algorithm. (It solves a similar problem to the more famous Lanczos algorithm but works internally in a slightly different way.)

ITensor has an implementation of Davidson you can use, and it has this signature:

template <class BigMatrixT>
std::vector<Real>
davidson(BigMatrixT const& A,
         std::vector<ITensor>& phi,
         Args const& args)

Steps to use it are:

  1. define a class (name of your own choosing) which holds your matrix, or pieces of your matrix, internally and then implements the method:
ITensor product(ITensor const& input, ITensor & output);

where conceptually it is doing output = A x input where “x” means matrix multiplication. Also this class should define a .size() method which returns the linear size of the matrix it represents.

  1. input a std::vector of initial guesses for the eigenvectors, as ITensors

  2. call auto eigs = davidson(A,phis,{"MaxIter=",30}); to get the computed eigenvalues. The variable phis is overwritten with the computed eigenvectors. Here “MaxIter” is a parameter you can choose and the larger you make it the harder the algorithm will work to get you a better accuracy.

You can see more of the code here:
https://github.com/ITensor/ITensor/blob/45bf86103b232f6b0d9992137253945f44e2a2f3/itensor/iterativesolvers.h#L101

Lastly, diagPosSemiDef does call a LAPACK routine internally, so I don’t think you’d be able to gain any speed by trying to call one yourself.