Fall back to dense eigh in sparse_cholesky for rank-deficient PSD inputs

[Transurgeon]

Description

Please include a short summary of the change.
sparse_cholesky (added in #3070, partially patched in #3270) raises
ValueError: Cholesky factorization failed on simple rank-deficient PSD
inputs such as [[1, 1], [1, 1]]. The root cause is that QDLDL has no
pivoting and bails out as soon as a zero pivot appears mid-elimination —
which the existing rank-deficient handling in #3270 doesn't cover (it
only special-cases all-zero rows and zero pivots that land last in
QDLDL's elimination order).

This adds a dense numpy.linalg.eigh fallback inside
sparse_cholesky's existing except Exception branch. When QDLDL fails,
we eigendecompose A, classify it as PSD/NSD/indefinite using the same
tolerance logic as the QDLDL branch, and return an n × rank(A) sparse
factor built from the nonzero eigenpairs. The function's documented
contract (L[p,:] @ L[p,:].T == sign * A, k = rank(A)) now holds for
all PSD/NSD inputs.

Issue link (if applicable):

Type of change

• New feature (backwards compatible)
• New feature (breaking API changes)
• Bug fix
• Other (Documentation, CI, ...)

Contribution checklist

• Add our license to new files.
• Check that your code adheres to our coding style.
• Write unittests.
• Run the unittests and check that they’re passing.
• Run the benchmarks to make sure your change doesn’t introduce a regression.

[claude]

The except Exception clause (unchanged from before) now triggers an O(n³) dense eigendecomposition + O(n²) memory allocation via A.toarray(), rather than just re-raising as a ValueError. If a large sparse matrix hits this path for an unexpected reason (not rank-deficiency), the dense fallback could OOM or hang.

Consider either:

1. Narrowing the catch to the specific exception type QDLDL raises on zero-pivot (likely a C-level RuntimeError), or
2. Adding a size guard before densifying, e.g.:

if n > 5000:
    raise ValueError(f"{SparseCholeskyMessages.FACTORIZATION_FAILED}: matrix too large for dense fallback")

This way truly unexpected failures on large matrices still surface immediately rather than silently attempting dense computation.

[odow]

This is the same saga as jump-dev/MathOptInterface.jl#2973 ultimately "fixed" jump-dev/MathOptInterface.jl#2967. Even then I don't like it. There's no magic fix. Every option has trade-offs.

We discussed doing this jump-dev/MathOptInterface.jl#2931 but decided not to merge because of the performance issues for larger models.

[PTNobel]

Okay, so I looked at our implementation more; we're using LDL for low-rank factorization in the dense case via Bunch-Kaufman pivoting.

It should be possible to do the same in the sparse case, especially since we don't need symbolic analysis reuse or anything.

Let me try some experiments and see what we can do...

[PTNobel]

Here's a very messy implementation that preserves sparsity and handles low-rank

#3311

[PTNobel]

Here's an easy fix: #3313.

I would like to---at some point---use a Bunch-Kaufman-based sparse-LDL factorization so that we can get low-rank support, but maybe that's a different problem.

[PTNobel]

#3313 fixes this better.