sparse-ldlt

Pure-Rust, dependency-free sparse symmetric-indefinite LDLᵀ factorization and solver.

Factors a symmetric sparse matrix A = L·D·Lᵀ (with L unit-lower-triangular and D a signed diagonal), then solves A x = b. Because D may hold negative entries, it handles symmetric indefinite systems - not just positive-definite ones - and it exposes D, so you can read the matrix inertia (the number of negative eigenvalues, by Sylvester's law of inertia) directly from the pivots.

Why this exists

Most sparse direct solvers available in pure Rust only provide positive-definite Cholesky and never expose the signed pivots. That leaves a gap for the problems where the sign of the pivots is the whole point:

• saddle-point / KKT systems from constrained optimisation and mixed finite elements,
• shifted eigenvalue matrices K - σM, which are indefinite for shifts inside the spectrum - factoring them and counting negative pivots gives a Sturm eigenvalue count,
• quasi-definite systems (interior-point methods).

sparse-ldlt is a small, self-contained implementation of the standard up-looking sparse LDLᵀ (elimination-tree) method - see T. A. Davis, Direct Methods for Sparse Linear Systems (SIAM, 2006) - with zero runtime dependencies, no unsafe, and stable-Rust only. (Its own tests and benchmarks use dev-only crates - criterion - which never ship to consumers.)

No pivoting - read this before relying on the inertia count

The factorization is un-pivoted. A pivot D[k] that reaches exactly zero fails loudly with LdltError::ZeroPivot(k). A pivot whose magnitude has been destroyed by cancellation is just as dangerous and is now reported too, as LdltError::NearZeroPivot { column, pivot, scale, suggested_shift }, whenever |D[k]| < NEAR_ZERO_PIVOT_REL * scale (1e-13 relative to the largest absolute diagonal entry of the input, about 1000x f64::EPSILON). That case is the one worth naming: such a pivot still carries a sign, but the sign is rounding noise, and the sign pattern of D IS the matrix inertia - so on a matrix near a singular point (a shift σ landing on an eigenvalue, a mechanism in a structure) the old behaviour was a WRONG inertia count with no error raised. It is no longer returned silently.

Recovering from it is no longer the caller's problem either. factor_shifted (and factor_perm_shifted) try the unshifted factorization first, and on a breakdown retry with a positive diagonal shift - starting at suggested_shift, multiplying by 8, at most 8 attempts. shift() reports how far the matrix was moved:

let f = SparseLdlt::factor_shifted(n, &col_ptr, &row_idx, &values)?;
if f.shift() != 0.0 {
    // This is an exact factorization of A + shift*I, NOT of A. Its inertia is the shifted
    // matrix's inertia, so a Sturm count from it is a count at (sigma - shift): correct for
    // it. Ignoring shift() is a bug.
}

If your use case needs certified pivots without any perturbation, you need a pivoting solver (Bunch-Kaufman / multifrontal); this crate trades that machinery for ~800 dependency-free lines, and reports honestly where the trade bites.

Fill-reducing ordering (AMD)

use sparse_ldlt::{amd, SparseLdlt};

let order = amd(n, &col_ptr, &row_idx);                 // Amestoy-Davis-Duff ordering
let f = SparseLdlt::factor_perm(n, &col_ptr, &row_idx, &values, &order).unwrap();
let x = f.solve(&b).unwrap();                           // permutation handled for you

Without an ordering, fill-in on an irregular sparsity can cost orders of magnitude (measured, cargo bench: random 2%-dense n=1024 - unordered factor ~0.27 s vs ~30 µs on a banded matrix of the same order; with amd + factor_perm the same random matrix factors in ~70 ms). amd is a quotient-graph approximate minimum degree (Amestoy-Davis-Duff 1996) implemented in this crate with the same zero-dependency rules: eliminated nodes become elements, degrees are the AMD external degrees recomputed over the neighbourhood only. Ordering NEVER changes inertia (a symmetric permutation is a congruence - Sylvester's law), so Sturm counts are identical with or without it; the tests/ordering.rs gate asserts exactly that, alongside measured fill reduction.

Usage

Supply the matrix in compressed-sparse-column (CSC) form. Only the upper triangle (entries with row ≤ col in each column) is read, so a fully-populated symmetric matrix is also fine.

use sparse_ldlt::SparseLdlt;

// Symmetric indefinite 3x3 matrix (full storage), in CSC:
//   [ 2  1  0 ]
//   [ 1 -3  1 ]
//   [ 0  1  2 ]
let col_ptr = vec![0, 2, 5, 7];
let row_idx = vec![0, 1,  0, 1, 2,  1, 2];
let values  = vec![2.0, 1.0,  1.0, -3.0, 1.0,  1.0, 2.0];

let f = SparseLdlt::factor(3, &col_ptr, &row_idx, &values).unwrap();

// Solve A x = b
let x = f.solve(&[1.0, 2.0, 3.0]).unwrap();

// Inertia: number of negative eigenvalues == number of negative pivots
let negative_eigenvalues = f.d().iter().filter(|&&v| v < 0.0).count();
assert_eq!(negative_eigenvalues, 1);

Notes

• Ordering: amd + factor_perm are built in (see above). The plain factor still applies none - deterministic and unchanged since v0.1.0.
• No pivoting (see the section above): breakdown is loud - LdltError::ZeroPivot for an exact zero, LdltError::NearZeroPivot for a pivot whose sign has become rounding noise. factor_shifted does the nudging and shift() says how much it nudged. Non-finite input values (NaN / ±inf) are rejected.
• solve returns Result<Vec<f64>, LdltError> - a right-hand side that does not match the factored order is LdltError::SizeMismatch, never a panic.
• Correctness is gated by an inertia oracle (tests/inertia_oracle.rs) in the spirit of feral's consensus validation, scoped to what a dependency-free crate can run anywhere: wherever a factorization succeeds, the pivot-sign inertia must be exactly correct - no tolerance. The oracle families are matrices whose inertia is known by construction (congruence A = XᵀSX, quasi-definite KKT blocks, Sturm shifts with exact endpoints and a monotonicity sweep), plus dense-residual checks at machine precision. A fourth, adversarial family targets the near-zero pivot directly - shifts landing 1e-15 from an eigenvalue, quasi-definite blocks with a diagonal entry driven to 1e-18, and KKT saddle points with a zero (2,2) block and a rank-deficient constraint - and is oracled against a dependency-free dense cyclic Jacobi eigensolver that is itself checked against closed-form spectra. Each fixture must either be refused or produce the exact inertia, and every refusal must then be recovered by factor_shifted with a small residual against A + shift*I.
• Property tests (tests/property.rs) pin the adversarial-CSC contract: duplicate entries are summed, explicit zeros are harmless, any row order within a column is accepted, degenerate shapes (n = 0, empty columns) never panic, malformed arrays return InvalidInput, and ~330 random valid-shape CSCs (wild magnitudes included) produce either a correct factorization or an honest pivot breakdown - never a panic or a NaN pivot.
• Real-matrix corpus (tests/corpus.rs): real structural stiffness matrices from the SuiteSparse (Harwell-Boeing) collection are bundled as fixtures and gated on external metadata - SPD by the collection, so inertia must be exactly 0 - plus a dense Jacobi cross-check and a corpus-tests feature that sweeps any directory of .mtx files (CK_LDLT_CORPUS_DIR) for CI-scale validation. No network, no dependencies.
• Benchmarks (cargo bench, criterion): factor/solve vs n on banded (structural) and random-sparse patterns, with and without AMD. The measured fill wall - n = 1024, banded ~30 µs vs random-2% ~0.27 s - is what amd addresses: the same random matrix factors in ~70 ms ordered (~0.27 s unordered). On a banded matrix AMD neither helps nor much hurts (n = 4096: ~121 µs unordered, ~187 µs ordered). Numbers from one machine; run cargo bench for yours.

Provenance

Written from the algorithm's published description - T. A. Davis, Direct Methods for Sparse Linear Systems (SIAM, 2006) - as an independent implementation. No source code from Tim Davis's LDL or from sprs-ldl was copied, translated, or consulted while writing it; any resemblance is the algorithm itself, which is published mathematics.

Created by MAGE Engineering for its FEM Analysis Studio, where it replaces an LGPL sparse LDLᵀ dependency in the structural analysis engine. Released under the MIT licence so the wider Rust community can use it too.

Licence

MIT - see LICENSE.