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.
 A: Much of (perhaps most of) the entire field of theoretical condensed matter physics is dedicated to solving this problem.  I don't think you'll be able to find a comprehensive answer on Stack Exchange.
A: 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.
Alternatively, if you're looking to solve for eigenvalues within a band, you can have a look into the Lanczos diagonalization algorithm.
