It depends a lot on the specifics of the Hamiltonian and system being examined. One approach that works sometimes is to cast the computation as a semidefinite program and use the duality theory there to obtain bounds.
One can write \begin{align} \lambda = \min_\rho \operatorname{tr}(\rho H), \end{align} where the minimization is over density operators $\rho$ (that is positive semidefinite, trace $1$) operators. Weak duality then gives you that \begin{align} \lambda \geq &\max y\\ &\ \ \text{s.t. } yI \preceq H, \end{align} which is almost trivial if you unpack it, the ground state energy is the largest constant you can subtract from your Hamiltonian such that it is a positive operator.
You can take this approach and make it essentially arbitrarily complicated. For example see this arXiv paper. In this work they rewrite the problem as instead of minimizing over some global state they minimize over some local state (e.g. $2$ body) state with compatibility conditions added to ensure all the little local states consistently add up to some global state. Then by adding more and more of these compatibility conditions to their optimization problem they obtain tighter and tighter lower bounds.
Edit: This is the paper I was originally searching for to find the answer to this question. Here what they do is observe that the hard part of optimizing over all density matrices is ensuring positivity, so they relax the assumption of positivity. A matrix $\rho$ is positive semidefinite if and only if \begin{align} \operatorname{tr}(A A^\dagger\rho) \geq 0, \end{align} for all choices of $A$. So what they do is just pick a "smart" set of operators $A_j$ and optimize the energy over all matrices $\rho$ satisfying $\operatorname{tr}(A_j A^\dagger_j\rho) \geq 0$ for all $j$, rather than those for which the trace is positive for all $A$ *. This is a bigger set than positive-semidefinite $\rho$ so you get a lower bound.