I have a constrained optimization problem of the form:

$$ \min_{Bx=g} \frac{1}{2} x^T A x - x^T f $$

$A \in \mathbb{R}^{n\times n}$ is positive semi-definite (with a tiny null space of dimension <6), sparse (10-20 non-zeros per row) and large ($n=$ thousands or millions).

$B \in \mathbb{R}^{m\times n}$ is dense but short relative to $A$ ($m=$ tens to hundreds).

I'm currently solving this problem by the Lagrange Multiplier method, which reveals the solution via the system of equations:

$$\underbrace{ \begin{bmatrix} A & B^T \\ B & 0 \end{bmatrix}}_{M} \begin{bmatrix} x \\ \lambda \end{bmatrix} = \begin{bmatrix} f\\g \end{bmatrix}. $$ Despite the dense rows and columns from $B$, I currently just treat the entire matrix $M$ as sparse and invoke a sparse direct $LDL^T$ solver.

Is there something more efficient to do in this case?

I considered applying the Schur Complement trick:

$$ x = A^{-1}(f-B^T (B A^{-1} B^T)^{-1} (B A^{-1} f - g)) $$

which would only involve solves against the truly sparse $A$, which I could prefactor, but (unless there's a further reduction I'm not seeing) it would involve $m+2$ back-substitutions (and also a dense solve against $(B A^{-1} B^T) \in \mathbb{R}^{m \times m})$.

Ideally, I'd like to only rely on dense/sparse direct solvers (i.e., avoiding iterative methods or algebraic multigrid etc.).

sparselysolve against $Z^T A Z$. $\endgroup$ – Alec Jacobson May 21 at 16:32