# 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.

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'$.
• 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