la-stack

Fast, stack-allocated linear algebra for fixed dimensions in Rust.

This crate grew from the need to support delaunay with fast, stack-allocated linear algebra primitives and algorithms while keeping the API intentionally small and explicit.

📐 Introduction

la-stack provides a handful of const-generic, stack-backed building blocks:

• Vector<const D: usize> for fixed-length vectors ([f64; D] today)
• Matrix<const D: usize> for fixed-size square matrices ([[f64; D]; D] today)
• Lu<const D: usize> for LU factorization with partial pivoting (solve + det)
• Ldlt<const D: usize> for LDLT factorization without pivoting (solve + det; symmetric SPD/PSD)

✨ Design goals

• ✅ Copy types where possible
• ✅ Const-generic dimensions (no dynamic sizes)
• ✅ const fn where possible (compile-time evaluation of determinants, dot products, etc.)
• ✅ Explicit algorithms (LU, solve, determinant)
• ✅ Robust geometric predicates via optional exact arithmetic (det_sign_exact, det_errbound)
• ✅ Exact linear system solve via optional arbitrary-precision arithmetic (solve_exact, solve_exact_f64)
• ✅ No runtime dependencies by default (optional features may add deps)
• ✅ Stack storage only (no heap allocation in core types)
• ✅ unsafe forbidden

See CHANGELOG.md for details.

🚫 Anti-goals

• Bare-metal performance: see blas-src, lapack-src, openblas-src
• Comprehensive: use nalgebra if you need a full-featured library
• Large matrices/dimensions with parallelism: use faer if you need this

🔢 Scalar types

The core types use f64. When f64 precision is insufficient (e.g. near-degenerate geometric configurations), the optional "exact" feature provides arbitrary-precision arithmetic via BigRational (see below).

🚀 Quickstart

Add this to your Cargo.toml:

[dependencies]
la-stack = "0.3.0"

Solve a 5×5 system via LU:

use la_stack::prelude::*;

// This system requires pivoting (a[0][0] = 0), so it's a good LU demo.
// A = J - I: zeros on diagonal, ones elsewhere.
let a = Matrix::<5>::from_rows([
    [0.0, 1.0, 1.0, 1.0, 1.0],
    [1.0, 0.0, 1.0, 1.0, 1.0],
    [1.0, 1.0, 0.0, 1.0, 1.0],
    [1.0, 1.0, 1.0, 0.0, 1.0],
    [1.0, 1.0, 1.0, 1.0, 0.0],
]);

let b = Vector::<5>::new([14.0, 13.0, 12.0, 11.0, 10.0]);

let lu = a.lu(DEFAULT_PIVOT_TOL).unwrap();
let x = lu.solve_vec(b).unwrap().into_array();

// Floating-point rounding is expected; compare with a tolerance.
let expected = [1.0, 2.0, 3.0, 4.0, 5.0];
for (x_i, e_i) in x.iter().zip(expected.iter()) {
    assert!((*x_i - *e_i).abs() <= 1e-12);
}

Compute a determinant for a symmetric SPD matrix via LDLT (no pivoting).

For symmetric positive-definite matrices, LDL^T is essentially a square-root-free form of the Cholesky decomposition (you can recover a Cholesky factor by absorbing sqrt(D) into L):

use la_stack::prelude::*;

// This matrix is symmetric positive-definite (A = L*L^T) so LDLT works without pivoting.
let a = Matrix::<5>::from_rows([
    [1.0, 1.0, 0.0, 0.0, 0.0],
    [1.0, 2.0, 1.0, 0.0, 0.0],
    [0.0, 1.0, 2.0, 1.0, 0.0],
    [0.0, 0.0, 1.0, 2.0, 1.0],
    [0.0, 0.0, 0.0, 1.0, 2.0],
]);

let det = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap().det();
assert!((det - 1.0).abs() <= 1e-12);

⚡ Compile-time determinants (D ≤ 4)

det_direct() is a const fn providing closed-form determinants for D=0–4, using fused multiply-add where applicable. Matrix::<0>::zero().det_direct() returns Some(1.0) (the empty-product convention). For D=1–4, cofactor expansion bypasses LU factorization entirely. This enables compile-time evaluation when inputs are known:

use la_stack::prelude::*;

// Evaluated entirely at compile time — no runtime cost.
const DET: Option<f64> = {
    let m = Matrix::<3>::from_rows([
        [2.0, 0.0, 0.0],
        [0.0, 3.0, 0.0],
        [0.0, 0.0, 5.0],
    ]);
    m.det_direct()
};
assert_eq!(DET, Some(30.0));

The public det() method automatically dispatches through the closed-form path for D ≤ 4 and falls back to LU for D ≥ 5 — no API change needed.

🔬 Exact arithmetic ("exact" feature)

The default build has zero runtime dependencies. Enable the optional exact Cargo feature to add exact arithmetic methods using arbitrary-precision rationals (this pulls in num-bigint, num-rational, and num-traits for BigRational):

[dependencies]
la-stack = { version = "0.3.0", features = ["exact"] }

Determinants:

• det_exact() — returns the exact determinant as a BigRational
• det_exact_f64() — returns the exact determinant converted to the nearest f64
• det_sign_exact() — returns the provably correct sign (−1, 0, or +1)

Linear system solve:

• solve_exact(b) — solves Ax = b exactly, returning [BigRational; D]
• solve_exact_f64(b) — solves Ax = b exactly, converting the result to Vector<D> (f64)

use la_stack::prelude::*;

// Exact determinant
let m = Matrix::<3>::from_rows([
    [1.0, 2.0, 3.0],
    [4.0, 5.0, 6.0],
    [7.0, 8.0, 9.0],
]);
assert_eq!(m.det_sign_exact().unwrap(), 0); // exactly singular

let det = m.det_exact().unwrap();
assert_eq!(det, BigRational::from_integer(0.into())); // exact zero

// Exact linear system solve
let a = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
let b = Vector::<2>::new([5.0, 11.0]);
let x = a.solve_exact_f64(b).unwrap().into_array();
assert!((x[0] - 1.0).abs() <= f64::EPSILON);
assert!((x[1] - 2.0).abs() <= f64::EPSILON);

BigRational is re-exported from the crate root and prelude when the exact feature is enabled, so consumers don't need to depend on num-rational directly.

For det_sign_exact(), D ≤ 4 matrices use a fast f64 filter (error-bounded det_direct()) that resolves the sign without allocating. Only near-degenerate or large (D ≥ 5) matrices fall through to the exact Bareiss algorithm.

Adaptive precision with det_errbound()

det_errbound() returns the conservative absolute error bound used by the fast filter. This method does NOT require the exact feature — it uses pure f64 arithmetic and is available by default. This enables building custom adaptive-precision logic for geometric predicates:

use la_stack::prelude::*;

let m = Matrix::<3>::identity();
if let Some(bound) = m.det_errbound() {
    let det = m.det_direct().unwrap();
    if det.abs() > bound {
        // f64 sign is guaranteed correct
        let sign = det.signum() as i8;
    } else {
        // Fall back to exact arithmetic (requires `exact` feature)
        let sign = m.det_sign_exact().unwrap();
    }
} else {
    // D ≥ 5: no fast filter, use exact directly (requires `exact` feature)
    let sign = m.det_sign_exact().unwrap();
}

The error coefficients (ERR_COEFF_2, ERR_COEFF_3, ERR_COEFF_4) are also exposed for advanced use cases.

🧩 API at a glance

Type	Storage	Purpose	Key methods
Vector<D>	[f64; D]	Fixed-length vector	new, zero, dot, norm2_sq
Matrix<D>	[[f64; D]; D]	Square matrix	See below
Lu<D>	Matrix<D> + pivot array	Factorization for solves/det	solve_vec, det
Ldlt<D>	Matrix<D>	Factorization for symmetric SPD/PSD solves/det	solve_vec, det

Storage shown above reflects the f64 implementation.

Matrix<D> key methods: lu, ldlt, det, det_direct, det_errbound, det_exact¹, det_exact_f64¹, det_sign_exact¹, solve_exact¹, solve_exact_f64¹.

¹ Requires features = ["exact"].

📋 Examples

The examples/ directory contains small, runnable programs:

• solve_5x5 — solve a 5×5 system via LU with partial pivoting
• det_5x5 — determinant of a 5×5 matrix via LU
• ldlt_solve_3x3 — solve a 3×3 symmetric positive definite system via LDLT
• const_det_4x4 — compile-time 4×4 determinant via det_direct()
• exact_det_3x3 — exact determinant value of a near-singular 3×3 matrix (requires exact feature)
• exact_sign_3x3 — exact determinant sign of a near-singular 3×3 matrix (requires exact feature)
• exact_solve_3x3 — exact solve of a near-singular 3×3 system vs f64 LU (requires exact feature)

just examples
# or individually:
cargo run --example solve_5x5
cargo run --example det_5x5
cargo run --example ldlt_solve_3x3
cargo run --example const_det_4x4
cargo run --features exact --example exact_det_3x3
cargo run --features exact --example exact_sign_3x3
cargo run --features exact --example exact_solve_3x3

🤝 Contributing

A short contributor workflow:

cargo install just
just setup        # install/verify dev tools + sync Python deps
just check        # lint/validate (non-mutating)
just fix          # apply auto-fixes (mutating)
just ci           # lint + tests + examples + bench compile

For the full set of developer commands, see just --list and AGENTS.md.

📝 Citation

If you use this library in academic work, please cite it using CITATION.cff (or GitHub's "Cite this repository" feature). A Zenodo DOI will be added for tagged releases.

📚 References

For canonical references to the algorithms used by this crate, see REFERENCES.md.

📊 Benchmarks (vs nalgebra/faer)

Raw data: docs/assets/bench/vs_linalg_lu_solve_median.csv

Summary (median time; lower is better). The “la-stack vs nalgebra/faer” columns show the % time reduction relative to each baseline (positive = la-stack faster):

D	la-stack median (ns)	nalgebra median (ns)	faer median (ns)	la-stack vs nalgebra	la-stack vs faer
2	2.026	4.476	142.364	+54.7%	+98.6%
3	15.718	23.857	191.028	+34.1%	+91.8%
4	28.171	53.516	213.492	+47.4%	+86.8%
5	47.595	72.861	287.763	+34.7%	+83.5%
8	137.876	163.720	365.792	+15.8%	+62.3%
16	609.456	594.194	910.985	-2.6%	+33.1%
32	2,719.556	2,812.766	2,921.820	+3.3%	+6.9%
64	17,776.557	14,083.938	12,541.345	-26.2%	-41.7%

📄 License

BSD 3-Clause License. See LICENSE.