# Eigenvectors of slightly perturbed matrix

Assume I know eigen-pairs $(\epsilon_i,\vec\phi_i)$ for matrix operator $\hat H$ that $\hat H \vec \phi_i = \epsilon_i \vec \phi_i$

Now I have slightly perturbed $\hat H' = \hat H + \hat R$ were perturbation $\hat R$ is small with respect to $\hat H$ (in some sence, e.g. Frobenius norm ?). Now I would like to find eigen-pairs of the perturbed $\hat H'$. $\hat H' \vec \phi'_i = \epsilon'_i \vec \phi'_i$

I expect there should exist some efficient iterative algorithm which use the knowledge of solution $\hat H \vec \phi_i = \epsilon_i \vec \phi_i$ to obtain solution of $\hat H' \vec \phi'_i = \epsilon'_i \vec \phi'_i$ with modest computational effort ( measured in CPU flops ) especially if I work with sparse matrices.

Is Davidson Method what I search for ? Or, in general, I can use the known solution of un-perturbed Hamiltonian to construct some preconditioner and initial vectors for Arnoldi/Lancozs iteration.

background:

I'm not expert in linear algebra, nor iterative methods. I'm just trying to familiarize myself with the field by googling. I downloaded some papers about Davidson / Arnoldi / Lanczos methods, but at first look I wasn't able to extract clear answer to even the very basic question if it is relevant to my problem.

I also read wikipage Eigenvalue perturbation but it doesn't concern with explicit numerical algorithm and it computational cost.

Mostly I'm interested in case when $\hat H$ and $\hat H'$ are hermitian (ground state Hamiltonian of some quantum system), but for more general understanding I do not want to limit the question on that case.

## 1 Answer

One of the best strategies for finding particular eigenvalues of linear operators is shift-and-invert Lanczos; if you're looking for an eigenvalue of the matrix $A$ close to the value $\sigma$, then you run the Lanczos algorithm on the operator $(A - \sigma I)^{-1}$. The spectrum of this operator will be much better separated than the spectrum of $A$, and since the Lanczos algorithm is great at finding extreme eigenvalues, it will pick out the eigenvalue close to $\sigma$ much faster. Of course, one then has to solve linear systems for $A - \sigma I$, which may be expensive.

In your case, you know the eigenvalue decomposition of $H$, so you can easily solve linear systems $(H - \sigma I)\phi = f$. For a given $f$, the solution $\phi$ should be fairly close to the solution of the system $(H' - \sigma I)\phi' = f$ for the perturbed Hamiltonian, provided that $\sigma$ is not too large. So you can leverage your solution of the unperturbed problem by using $H - \sigma I$ as a preconditioner for $H' - \sigma I$ in a Krylov subspace method such as MINRES, rather than a generic and possibly more expensive method such as $LDL^*$-factorization.

As a blunt approach, you can take $\sigma$ to be the eigenvalue $\epsilon_i$ of $H$; better yet, you can use the usual perturbation theory estimates from quantum mechanics for $\epsilon_i'$.

This strategy is generally most effective if you're looking for the ground state or the first few excited states. I haven't tried the Davidson method myself so I don't know how one might apply the same reasoning as for shift-and-invert Lanczos. This is discussed in the book Numerical Methods for Large Eigenvalue Problems by Yousef Saad, which is (in my opinion) the best resource on this subject. Saad has also done a lot of work on finding interior eigenvalues, which is a substantially harder problem.

• Make sence. Just what I don't like is that when I want to find considerable share of eigenvectors ( let say 25%-30%, consider tight-binding or LCAO with ~2x as much basisfunctions than electrons ), I feel this is not much less costly than full diagonalization. There is double iteration for each eigenvector, outer loop of Lanczos and inner loop of CG ( or MINRES ). I see that CG would converge quickly due to good initial estimate and preconditioning, and Lancozs due to good estimate of $\sigma$. But still, one would expect you exploit known estimate of all, not just one eigenpair. – Prokop Hapala Mar 31 '16 at 15:53
• Ah I see, I forgot to mention restarted Lanczos -- once one eigenpair has converged, say the $k$-th excited state, you can change $\sigma$ to be the current estimate of the $k+1$-th eigenvalue. Using this, you can get a bunch of eigenpairs at once. If you want an implementation, I recommend SLEPc, which has Python bindings if you're so inclined. – Daniel Shapero Mar 31 '16 at 16:47