constraint as the $e\mathcal{T}_1$ which has the lowest eigenvalue, [*i.e.*]{} we optimize the following cost function: $$\phi^{e\mathcal{T}_1}(p) = \lambda^{e\mathcal{T}_1}_0(p) + \mu (|p|^2 - C)^2~,$$ as a function of the single-particle state $p$, in which $\mu$ scales the quadratic potential that makes sure the $p$ vector is normalized, and $\lambda_0^{e\mathcal{T}_1}(p)$ is the lowest eigenvalue of the $e\mathcal{T}_1$ matrix for the current value of the state $p$. The gradient of this cost function can be evaluated analytically using the Hellman-Feynman theorem [@hellman_feynman]: $$\frac{\partial \lambda_0^{e\mathcal{T}_1}}{\partial p_{\bar{z}}} = {\langle \Psi_0|}\left(\frac{\partial e\mathcal{T}_1(p)}{\partial p_{\bar{z}}}\right){|\Psi_0\rangle}~,\qquad\text{with}\qquad e\mathcal{T}_1(\Gamma){|\Psi_0\rangle} = \lambda_0^{e\mathcal{T}_1}{|\Psi_0\rangle}~.$$ The cost function can then be optimized using a simple non-linear conjugate gradient algorithm [@jrs]. For the construction of an efficient conjugate gradient algorithm it is essential to have a fast evaluation of the lowest eigenvalue and eigenvector of the constraint with a certain state $p$. This can be achieved because the largest block, the $\mathcal{T}_1$ on the cluster, remains unchanged during the optimization. If we prediagonalize the large block on the cluster, we need a fast method to solve the following problem: $$\left[
\begin{matrix}
\lambda_1 &0&\ldots& 0 & b_{11} & \ldots & b_{1m}\\
0 &\lambda_2& \ldots& 0 & b_{21} & \ldots & b_{2m}\\
0 & 0 &\ddots  &0 