Crate la_stack

Expand description

§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.4.1"

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);

⚠️ LDLT precondition: The input matrix must be symmetric. Symmetry is verified by a debug_assert! in debug builds only; release builds silently accept asymmetric inputs and produce a meaningless factorization. Pre-validate with Matrix::is_symmetric (or locate the offending pair with Matrix::first_asymmetry) when you cannot statically guarantee symmetry, or fall back to lu() if your matrices may not be symmetric at all. See Matrix::ldlt for details.

§⚡ 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.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::*;

// 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);

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

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"].

§📊 Benchmarks (vs nalgebra/faer)

LU solve (factor + solve): median time vs dimension

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

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:

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.

Modules§

prelude: Common imports for ergonomic usage.

Structs§

BigInt: A big signed integer type.
Ldlt: LDLT factorization (A = L D Lᵀ) for symmetric positive (semi)definite matrices.
Lu: LU decomposition (PA = LU) with partial pivoting.
Matrix: Fixed-size square matrix D×D, stored inline.
Vector: Fixed-size vector of length D, stored inline.

Enums§

LaError: Linear algebra errors.

Constants§

DEFAULT_PIVOT_TOL: Default absolute pivot magnitude threshold used for LU pivot selection / singularity detection.
DEFAULT_SINGULAR_TOL: Default absolute threshold used for singularity/degeneracy detection.
ERR_COEFF_2: Absolute error coefficient for Matrix::<2>::det_direct.
ERR_COEFF_3: Absolute error coefficient for Matrix::<3>::det_direct.
ERR_COEFF_4: Absolute error coefficient for Matrix::<4>::det_direct.

Traits§

FromPrimitive: A generic trait for converting a number to a value.
Signed: Useful functions for signed numbers (i.e. numbers that can be negative).
ToPrimitive: A generic trait for converting a value to a number.

Type Aliases§

BigRational: Alias for arbitrary precision rationals.

Crate la_stack

Crate la_stack Copy item path

§la-stack

§📐 Introduction

§✨ Design goals

§🚫 Anti-goals

§🔢 Scalar types

§🚀 Quickstart

§⚡ Compile-time determinants (D ≤ 4)

§🔬 Exact arithmetic ("exact" feature)

§Adaptive precision with det_errbound()

§🧩 API at a glance

§📊 Benchmarks (vs nalgebra/faer)

§📋 Examples

§🤝 Contributing

§📝 Citation

§📚 References

§🤖 AI Agents

§📄 License

Modules§

Structs§

Enums§

Constants§

Traits§

Type Aliases§

Crate la_stack

§🔬 Exact arithmetic (`"exact"` feature)

§Adaptive precision with `det_errbound()`