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 f64 vectors backed by [f64; D]
• Matrix<const D: usize> for fixed-size square f64 matrices backed by [[f64; D]; D]
• 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 release history and docs/roadmap.md for current release planning.

🚫 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
• Alternate floating-point scalar families: la-stack supports f64 and optional exact arithmetic, not f32 / f16 APIs

🔢 Scalar types

The scalar model is intentionally limited to f64 for floating-point work and exact rationals behind the optional "exact" feature. This matches the crate's focus on small, robustness-sensitive numerical and computational geometry workloads. When f64 precision is insufficient (e.g. near-degenerate geometric configurations), the optional "exact" feature provides arbitrary-precision arithmetic via BigRational (see below).

Lower-precision f32 / f16 throughput-oriented workloads are outside the crate's scope; they usually indicate large-matrix or accelerator-oriented use cases better served by broader linear-algebra libraries.

🚀 Quickstart

Add this to your Cargo.toml:

[dependencies]
la-stack = "0.4.1"

Solve a 5×5 system via LU:

use la_stack::prelude::*;

fn main() -> Result<(), LaError> {
    // 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>::try_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>::try_new([14.0, 13.0, 12.0, 11.0, 10.0])?;

    let lu = a.lu(DEFAULT_PIVOT_TOL)?;
    let x = lu.solve_vec(b)?.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);
    }

    Ok(())
}

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::*;

fn main() -> Result<(), LaError> {
    // This matrix is symmetric positive-definite (A = L*L^T) so LDLT works without pivoting.
    let a = Matrix::<5>::try_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)?.det()?;
    assert!((det - 1.0).abs() <= 1e-12);

    Ok(())
}

⚠️ LDLT invariant: The input matrix must be symmetric. Asymmetric inputs passed to Matrix::ldlt return a typed LaError::Asymmetric before factorization starts. Use Matrix::first_asymmetry to locate the offending pair, or fall back to lu() if your matrices may not be symmetric at all. Symmetric inputs with a negative LDLT diagonal return LaError::NotPositiveSemidefinite; zero or too-small non-negative diagonals return LaError::Singular.

⚡ 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 Ok(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: Result<Option<f64>, LaError> = {
    let m = match Matrix::<3>::try_from_rows([
        [2.0, 0.0, 0.0],
        [0.0, 3.0, 0.0],
        [0.0, 0.0, 5.0],
    ]) {
        Ok(matrix) => matrix,
        Err(_) => panic!("matrix entries must be finite"),
    };
    m.det_direct()
};
assert_eq!(DET, Ok(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. Finite inputs return a floating-point determinant estimate in every dimension; det() does not surface LaError::Singular. Tiny nonzero determinants are not flattened by a pivot tolerance. Use lu() directly when you need tolerance-aware singularity detection or the pivot-column diagnostic from the factorization, and use the exact determinant APIs when exact singularity classification matters.

🔬 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.4.1", 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::*;

fn main() -> Result<(), LaError> {
    // Exact determinant
    let m = Matrix::<3>::try_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()?, 0); // exactly singular

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

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

    Ok(())
}

With the exact feature enabled, BigInt and BigRational are re-exported from the crate root and prelude, alongside the most commonly needed num-traits items (FromPrimitive, ToPrimitive, Signed). This lets consumers construct exact values (BigRational::from_f64, from_i64), query sign (is_positive / is_negative), and convert back to f64 (to_f64) with a single use la_stack::prelude::*; — no need to add num-bigint, num-rational, or num-traits to their own Cargo.toml.

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::*;

fn main() -> Result<(), LaError> {
    let m = Matrix::<3>::identity();
    if let Some(bound) = m.det_errbound()? {
        if let Some(det) = m.det_direct()? {
            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()?;
            }
        }
    } else {
        // D ≥ 5: no fast filter, use exact directly (requires `exact` feature)
        let sign = m.det_sign_exact()?;
    }

    Ok(())
}

The error coefficients (ERR_COEFF_2, ERR_COEFF_3, ERR_COEFF_4) are the dimension-specific constants behind that bound. In plain terms, they answer: "how many machine-epsilon-sized rounding mistakes can this closed-form determinant formula accumulate?" To get an absolute error bound, det_errbound() multiplies the coefficient by a size measure of the matrix entries, the absolute Leibniz sum:

p(|A|) = sum over determinant terms of product of absolute values

For a 2×2 matrix [[a, b], [c, d]], that scale is |a*d| + |b*c|, so:

|det_direct(A) - det_exact(A)| <= ERR_COEFF_2 * (|a*d| + |b*c|)

The coefficients are not tolerances and are not meant to be tuned by callers; they are conservative constants derived from the fixed D ≤ 4 formulas and their floating-point rounding chains. They are exposed for advanced users who want to compose the same bound themselves.

🧩 API at a glance

Type	Storage	Purpose	Key methods
Vector<D>	[f64; D]	Fixed-length vector for input and computation	new, zero, dot, norm2_sq
Matrix<D>	[[f64; D]; D]	Square matrix for input and computation	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 intentional f64 scalar model.

Matrix<D> key methods: lu, ldlt, det, det_direct, det_errbound, det_exact¹, det_exact_f64¹, det_sign_exact¹, solve_exact¹, solve_exact_f64¹. Matrix and vector methods validate non-finite inputs at public API boundaries and carry internal proof types through computation so successful factors do not store NaN or infinity.

¹ Requires features = ["exact"].

📊 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.173	4.448	139.923	+51.2%	+98.4%
3	13.989	34.607	180.026	+59.6%	+92.2%
4	27.580	48.435	203.163	+43.1%	+86.4%
5	53.517	75.935	274.375	+29.5%	+80.5%
8	134.859	162.859	371.463	+17.2%	+63.7%
16	635.775	576.171	846.189	-10.3%	+24.9%
32	2,704.570	2,731.740	2,589.494	+1.0%	-4.4%
64	17,381.460	13,744.505	11,276.642	-26.5%	-54.1%

📋 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

The repository uses Rust-native tooling for documentation and config checks: rumdl for Markdown, dprint with pretty_yaml for YAML, taplo for TOML, and typos for spelling. GitHub Actions references are SHA-pinned, restricted to an explicit allowlist, and kept with readable version comments for review.

CI runs just ci on Ubuntu, macOS, and Windows to keep platform coverage aligned with the local comprehensive validation path.

For coverage commands and report locations, see docs/COVERAGE.md. 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.

🤖 AI Agents

This repository contains an AGENTS.md file, which defines the canonical rules and invariants for all AI coding assistants and autonomous agents working on this codebase.

AI tools (including ChatGPT, Claude, CodeRabbit, KiloCode, and WARP) are expected to read and follow AGENTS.md when proposing or applying changes.

Portions of this library were developed with the assistance of these AI tools:

• ChatGPT
• Claude
• CodeRabbit
• KiloCode
• WARP

All code was written and/or reviewed and validated by the author.

For full tool citation metadata, see the AI-Assisted Development Tools section of REFERENCES.md.

📄 License

BSD 3-Clause License. See LICENSE.