Suppose we have a system of equations:
$$ \begin{bmatrix} \boldsymbol{A} & \boldsymbol{1} \\ \boldsymbol{1}^\top & 0 \end{bmatrix} \begin{bmatrix} \boldsymbol{x} \\ x_0 \end{bmatrix} = \begin{bmatrix} \boldsymbol{b} \\ 1 \end{bmatrix} $$
where $\boldsymbol{A} \succ 0$ is symmetric positive definitive. How would you exploit the Cholesky factorization of $\boldsymbol{A} = \boldsymbol{L}\boldsymbol{L}^\top$ to efficiently solve this system for various right-hand-side $\boldsymbol{b}$?
If there is no efficient trick with the Cholesky factorization, what would you recommend to guarantee decent numerical stability and speed?
We have $$ \mathbf{A}\mathbf{x}+x_0\mathbf{1}=\mathbf{b}, \quad \mathbf{1}^T\mathbf{x}=1. $$ Eliminate $\mathbf{x}$ from the first equation and substitute to the second: $$ \mathbf{x}=\mathbf{A}^{-1}(\mathbf{b}-x_0\mathbf{1}) \quad\implies\quad \mathbf{1}^T\mathbf{A}^{-1}(\mathbf{b}-x_0\mathbf{1})=1. $$ This gives $$ x_0=\frac{\mathbf{1}^T\mathbf{A}^{-1}\mathbf{b}-1}{\mathbf{1}^T\mathbf{A}^{-1}\mathbf{1}}. $$
So you can first solve using Cholesky the systems $$ \mathbf{y}=\mathbf{A}^{-1}\mathbf{b}, \quad \mathbf{f}=\mathbf{A}^{-1}\mathbf{1}, $$ and then form $x_0$ and $\mathbf{x}$ as $$ x_0=\frac{\mathbf{1}^T\mathbf{y}-1}{\mathbf{1}^T\mathbf{f}}, \quad \mathbf{x}=\mathbf{y}-x_0\mathbf{f}. $$