# Finding the ground state of a Hamiltonian Matrix

I have numerically constructed a Hamiltonian matrix. I am currently finding the ground state by full diagonalisation of the matrix (with the GSL library) and finding the most negative eigenvalue and its associated eigenvector. This is very slow and can surely be done better. I thought about just finding the eigenvalues and the solving the linear equation associated with the most negative eigenvalue for the eigenvector but this is probably not much quicker?

I am wondering if any one knows of any faster numerical routines for finding the ground state associated with a Hamiltonian matrix. If it makes any difference this is for the unit filled 1D Bose Hubbard Model in second quantised form.

i.e. $\hat{H} = -J\sum_{<ij>}\hat{b}_{i}^{\dagger}\hat{b}_{j} + \frac{U}{2}\sum_{i}\hat{n}_{i}(\hat{n}_{i} - 1)$

Also I have, already, fairly good estimates for the ground state energy (within 5%) and wavefunction if that may help any possible routines.

Many thanks.

• But I was under the assumption the most negative eigenvalue of this matrix does not necessarily have the largest magnitude so how would this work? Unless my assumption is provably false. – EducationalFerret Dec 1 '16 at 23:16
• Yes I believe I could produce an estimate to within 5% for it. If any other eigenvalues are within 5% of the smallest one then I'm guessing it would fail? I'm not sure if that's a likely scenario or not. – EducationalFerret Dec 1 '16 at 23:41
• Would Computational Science be a better home for this question? – Emilio Pisanty Dec 1 '16 at 23:53
• yes you're probably right, Rayleigh quotient iteration may well work co-incidentally because I also have a strong estimate for the ground state wavefunction. The need for multiple matrix inverses may be problematic. Thank you for your help – EducationalFerret Dec 1 '16 at 23:55

If all you want is the ground state, then you should look into imaginary-time propagation. In this method, you propagate the Schrödinger equation along time $t=-i\tau$, so you turn what used to be phases $e^{iEt}$ into decay factors $e^{-E\tau}$, with a higher decay constant for higher energies.
Thus, if you start with some random initial wavefunction $\psi_0$ and then propagate it along $$\hbar\partial_t\psi = -H\psi,$$ the formal solution for $\psi$ as a superposition of eigenstates, $\psi = \sum_n e^{-E_n\tau}\psi_{n,0}$, will reduce to just the ground-state component of $\psi_0$. (Normally, it is useful to re-set the norm of $\psi$ to $1$ at each timestep, to help with numerical stability.) This tends to be more efficient than a full numerical diagonalization, particularly since the decay of the unwanted components tends to be quite fast (and more so if the initial guess $\psi_0$ is close to the true ground state), but on the other hand I don't know how it will behave for a full-borne Hubbard model.