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][1]. 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.


  [1]: https://arxiv.org/abs/2212.03014