Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 35 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,39 @@ All notable changes to this project will be documented in this file.
The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).

## [Unreleased]

### Added

- Regression tests for solver and determinant overflow handling [`f763b11`](https://github.com/acgetchell/la-stack/commit/f763b119bcc57276b83f370b0bf7abce654c7eb8)
- Defensive-path test coverage for LU and LDLT solve_vec [`87d426f`](https://github.com/acgetchell/la-stack/commit/87d426fca1627445b804fd26b62fc7d9d4f0ae48)
- Const-ify Lu/Ldlt det + solve_vec and Matrix inf_norm + det_errbound [`81ecb35`](https://github.com/acgetchell/la-stack/commit/81ecb35bdaf159f1f44d1eb24274ecf82c6567d5)

### Changed

- Report infinite vs finite off-diagonal pairs as asymmetric [`1805779`](https://github.com/acgetchell/la-stack/commit/1805779dbca49183fbfa95c68ec00984966aa551)

### Documentation

- Strengthen LDLT symmetry precondition and add is_symmetric API [`1693307`](https://github.com/acgetchell/la-stack/commit/1693307c58e30a804e9b79abc783af99be8dc947)

### Fixed

- Propagate NaN in Matrix::inf_norm [#85](https://github.com/acgetchell/la-stack/pull/85) [`16ffa45`](https://github.com/acgetchell/la-stack/commit/16ffa45ade11b21a179cad8fcecc51d802086a1d)

### Maintenance

- Bump taiki-e/install-action from 2.73.0 to 2.75.7 [`a1d1c1e`](https://github.com/acgetchell/la-stack/commit/a1d1c1edba7d5bf6928de4e8cab86ba853685430)
- Bump actions/download-artifact from 4.3.0 to 8.0.1 [`8772ea2`](https://github.com/acgetchell/la-stack/commit/8772ea2f18cf1b94a21e269b4bc0c8e0a3b650eb)
- Bump rand from 0.9.2 to 0.9.4 [`b91670a`](https://github.com/acgetchell/la-stack/commit/b91670a46f82306dde5433ca7c406a8757e3fc59)
- Bump actions/upload-artifact from 7.0.0 to 7.0.1 [`7579274`](https://github.com/acgetchell/la-stack/commit/7579274d56ab8c46eae6906c07a3c0d0fd41c84f)
- Bump actions/github-script from 8.0.0 to 9.0.0 [`a522abd`](https://github.com/acgetchell/la-stack/commit/a522abdabf211d47f0f52cbfc2a349fa72ce87e2)
- Bump actions/cache from 5.0.4 to 5.0.5 [`b12d95c`](https://github.com/acgetchell/la-stack/commit/b12d95c8ce242ea6014c6a6aaa7845ff501d2e4d)
- Bump actions-rust-lang/setup-rust-toolchain [`c728774`](https://github.com/acgetchell/la-stack/commit/c72877460efba3d598f64ff253d2df8f3695acb1)
- Bump astral-sh/setup-uv from 8.0.0 to 8.1.0 [`7a18733`](https://github.com/acgetchell/la-stack/commit/7a18733ba0f39102678bbe30772fcde10f20e9d8)
- Bump taiki-e/install-action from 2.75.7 to 2.75.18 [`433cfc1`](https://github.com/acgetchell/la-stack/commit/433cfc1b01c9128e93c82cb553aa63d4091bace3)
- Bump MSRV to Rust 1.95 and adopt new stable features [`0ab3c33`](https://github.com/acgetchell/la-stack/commit/0ab3c336074f2b866256fbe5db8a8ec5306d580a)

## [0.4.0] - 2026-04-11

### Added
Expand Down Expand Up @@ -35,6 +68,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
- Bump taiki-e/install-action from 2.70.1 to 2.73.0 [`7de720d`](https://github.com/acgetchell/la-stack/commit/7de720dfb8328d01843cfabd15d23086ee98832b)
- Bump astral-sh/setup-uv from 7.6.0 to 8.0.0 [`af40753`](https://github.com/acgetchell/la-stack/commit/af40753130fac56d48da7fce2f18b11dc391ebe6)
- Add performance regression detection for exact-arithmetic benchmarks [`44bce99`](https://github.com/acgetchell/la-stack/commit/44bce99bcae8f852dbf7500e71fa182730e08bca)
- Release v0.4.0 [`843bb3c`](https://github.com/acgetchell/la-stack/commit/843bb3cfcc88e114bdd5e00d242b605af9b2b434)

### Performance

Expand Down Expand Up @@ -228,6 +262,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0

- Add tarpaulin coverage upload [`7486dfd`](https://github.com/acgetchell/la-stack/commit/7486dfd54e16a6dbde41575c3f35a1acb65f57d2)

[unreleased]: https://github.com/acgetchell/la-stack/compare/v0.4.0..HEAD
[0.4.0]: https://github.com/acgetchell/la-stack/compare/v0.3.0..v0.4.0
[0.3.0]: https://github.com/acgetchell/la-stack/compare/v0.2.2..v0.3.0
[0.2.2]: https://github.com/acgetchell/la-stack/compare/v0.2.1..v0.2.2
Expand Down
53 changes: 37 additions & 16 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -242,6 +242,27 @@ Storage shown above reflects the `f64` implementation.

¹ Requires `features = ["exact"]`.

## 📊 Benchmarks (vs nalgebra/faer)

![LU solve (factor + solve): median time vs dimension](docs/assets/bench/vs_linalg_lu_solve_median.svg)

Raw data: [docs/assets/bench/vs_linalg_lu_solve_median.csv](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):

<!-- BENCH_TABLE:lu_solve:median:new:BEGIN -->
| D | la-stack median (ns) | nalgebra median (ns) | faer median (ns) | la-stack vs nalgebra | la-stack vs faer |
|---:|--------------------:|--------------------:|----------------:|---------------------:|----------------:|
| 2 | 2.309 | 4.365 | 140.156 | +47.1% | +98.4% |
| 3 | 18.331 | 22.706 | 181.074 | +19.3% | +89.9% |
| 4 | 27.430 | 51.372 | 210.451 | +46.6% | +87.0% |
| 5 | 53.819 | 70.722 | 276.064 | +23.9% | +80.5% |
| 8 | 143.611 | 160.309 | 356.960 | +10.4% | +59.8% |
| 16 | 611.393 | 580.793 | 871.704 | -5.3% | +29.9% |
| 32 | 2,631.241 | 2,733.946 | 2,832.816 | +3.8% | +7.1% |
| 64 | 17,233.345 | 14,112.678 | 12,164.571 | -22.1% | -41.7% |
<!-- BENCH_TABLE:lu_solve:median:new:END -->

## 📋 Examples

The `examples/` directory contains small, runnable programs:
Expand Down Expand Up @@ -289,26 +310,26 @@ If you use this library in academic work, please cite it using [CITATION.cff](CI

For canonical references to the algorithms used by this crate, see [REFERENCES.md](REFERENCES.md).

## 📊 Benchmarks (vs nalgebra/faer)
## 🤖 AI Agents

![LU solve (factor + solve): median time vs dimension](docs/assets/bench/vs_linalg_lu_solve_median.svg)
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.

Raw data: [docs/assets/bench/vs_linalg_lu_solve_median.csv](docs/assets/bench/vs_linalg_lu_solve_median.csv)
AI tools (including `ChatGPT`, `Claude`, `CodeRabbit`, `KiloCode`, and `WARP`) are expected to read
and follow `AGENTS.md` when proposing or applying changes.

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):
Portions of this library were developed with the assistance of these AI tools:

<!-- BENCH_TABLE:lu_solve:median:new:BEGIN -->
| D | la-stack median (ns) | nalgebra median (ns) | faer median (ns) | la-stack vs nalgebra | la-stack vs faer |
|---:|--------------------:|--------------------:|----------------:|---------------------:|----------------:|
| 2 | 2.309 | 4.365 | 140.156 | +47.1% | +98.4% |
| 3 | 18.331 | 22.706 | 181.074 | +19.3% | +89.9% |
| 4 | 27.430 | 51.372 | 210.451 | +46.6% | +87.0% |
| 5 | 53.819 | 70.722 | 276.064 | +23.9% | +80.5% |
| 8 | 143.611 | 160.309 | 356.960 | +10.4% | +59.8% |
| 16 | 611.393 | 580.793 | 871.704 | -5.3% | +29.9% |
| 32 | 2,631.241 | 2,733.946 | 2,832.816 | +3.8% | +7.1% |
| 64 | 17,233.345 | 14,112.678 | 12,164.571 | -22.1% | -41.7% |
<!-- BENCH_TABLE:lu_solve:median:new:END -->
- [ChatGPT](https://openai.com/chatgpt)
- [Claude](https://www.anthropic.com/claude)
- [CodeRabbit](https://coderabbit.ai/)
- [KiloCode](https://kilocode.ai/)
- [WARP](https://www.warp.dev)

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

For full tool citation metadata, see the [AI-Assisted Development Tools](REFERENCES.md#ai-assisted-development-tools)
section of `REFERENCES.md`.

## 📄 License

Expand Down
72 changes: 41 additions & 31 deletions REFERENCES.md
Original file line number Diff line number Diff line change
Expand Up @@ -8,14 +8,46 @@ processed by GitHub and other platforms.

A Zenodo DOI will be added for tagged releases.

## AI-Assisted Development Tools

- Anthropic. "Claude." <https://www.anthropic.com/claude>.
- CodeRabbit AI, Inc. "CodeRabbit." <https://coderabbit.ai/>.
- KiloCode. "KiloCode AI Engineering Assistant." <https://kilocode.ai/>.
- OpenAI. "ChatGPT." <https://openai.com/chatgpt>.
- Warp Dev, Inc. "WARP." <https://www.warp.dev/>.

All AI-generated output was reviewed and/or edited by the maintainer.
No generated content was used without human oversight.

## Linear algebra algorithms

### LU decomposition (Gaussian elimination with partial pivoting)
### Absolute error bound for closed-form determinants

The LU implementation in `la-stack` follows the standard Gaussian elimination / LU factorization
approach with partial pivoting for numerical stability.
`Matrix::det_errbound()` returns a conservative Shewchuk-style absolute error bound [8]
for `Matrix::det_direct()` in dimensions 2–4. The bound has the form
`ERR_COEFF_D · p(|A|)`, where `p(|A|)` is the absolute Leibniz sum (the cofactor-expansion
tree with `|·|` at every leaf) and `ERR_COEFF_D ∈ {ERR_COEFF_2, ERR_COEFF_3, ERR_COEFF_4}`
is a dimension-specific constant derived from the rounding-event count of `det_direct`.
The same bound is used internally by `det_sign_exact()`'s fast filter, but
`det_errbound()` itself is available without the `exact` feature, so downstream crates
can build custom adaptive-precision logic with pure f64 arithmetic.

See references [1-3] below.
### Exact determinant sign (adaptive-precision Bareiss)

`det_sign_exact()` uses a Shewchuk-style f64 error-bound filter [8] (the same bound exposed
by `det_errbound()` above) backed by integer-only Bareiss elimination [7] in `BigInt`. Each
f64 entry is decomposed into `mantissa × 2^exponent`, scaled to a common integer base, and
eliminated without any `BigRational` or GCD overhead.
See `src/exact.rs` for the full architecture description.

### f64 → BigRational conversion (`f64_to_bigrational`)

`f64_to_bigrational` converts an f64 to an exact `BigRational` by decomposing the IEEE 754
binary64 bit representation [9] into its sign, exponent, and significand fields. Because
every finite f64 is exactly `±m × 2^e` (where `m` is an integer), the rational can be
constructed directly via `BigRational::new_raw` without GCD normalization — trailing zeros
in the significand are stripped first so the fraction is already in lowest terms. See
Goldberg [10] for background on IEEE 754 representation and exact rational reconstruction.

### LDL^T factorization (symmetric SPD/PSD)

Expand All @@ -26,16 +58,14 @@ pivoting.
For background on the SPD/PSD setting, see [4-5]. For pivoted variants used for symmetric *indefinite*
matrices, see [6].

### Exact determinant sign (adaptive-precision Bareiss)
### LU decomposition (Gaussian elimination with partial pivoting)

`det_sign_exact()` uses a Shewchuk-style f64 error-bound filter [8] backed by integer-only
Bareiss elimination [7] in `BigInt`. Each f64 entry is decomposed into `mantissa × 2^exponent`,
scaled to a common integer base, and eliminated without any `BigRational` or GCD overhead.
See `src/exact.rs` for the full architecture description.
The LU implementation in `la-stack` follows the standard Gaussian elimination / LU factorization
approach with partial pivoting for numerical stability.

## References
See references [1-3] below.

### LU / Gaussian elimination
## References

1. Trefethen, Lloyd N., and Robert S. Schreiber. "Average-case stability of Gaussian elimination."
*SIAM Journal on Matrix Analysis and Applications* 11.3 (1990): 335–360.
Expand All @@ -46,23 +76,14 @@ See `src/exact.rs` for the full architecture description.
3. Huang, Han, and K. Tikhomirov. "Average-case analysis of the Gaussian elimination with partial pivoting."
*Probability Theory and Related Fields* 189 (2024): 501–567.
[Open-access PDF](https://link.springer.com/article/10.1007/s00440-024-01276-2) (also: [arXiv:2206.01726](https://arxiv.org/abs/2206.01726))

### LDL^T / Cholesky (symmetric SPD/PSD)

4. Cholesky, Andre-Louis. "On the numerical solution of systems of linear equations"
(manuscript dated 2 Dec 1910; published 2005).
Scan + English analysis: [BibNum](https://www.bibnum.education.fr/mathematiques/algebre/sur-la-resolution-numerique-des-systemes-d-equations-lineaires)
5. Brezinski, Claude. "La methode de Cholesky." (2005).
[PDF](https://eudml.org/doc/252115)

### Pivoted LDL^T (symmetric indefinite)

6. Bunch, J. R., L. Kaufman, and B. N. Parlett. "Decomposition of a Symmetric Matrix."
*Numerische Mathematik* 27 (1976/1977): 95–110.
[Full text](https://eudml.org/doc/132435)

### Exact determinant (`det_exact`, `det_exact_f64`, `det_sign_exact`)

7. Bareiss, Erwin H. "Sylvester's Identity and Multistep Integer-Preserving Gaussian
Elimination." *Mathematics of Computation* 22.103 (1968): 565–578.
[DOI](https://doi.org/10.1090/S0025-5718-1968-0226829-0) ·
Expand All @@ -72,17 +93,6 @@ See `src/exact.rs` for the full architecture description.
[DOI](https://doi.org/10.1007/PL00009321) ·
[PDF](https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf)
Also: Technical Report CMU-CS-96-140, Carnegie Mellon University, May 1996.

### f64 → BigRational conversion (`f64_to_bigrational`)

`f64_to_bigrational` converts an f64 to an exact `BigRational` by decomposing the IEEE 754
binary64 bit representation into its sign, exponent, and significand fields. Because every
finite f64 is exactly `±m × 2^e` (where `m` is an integer), the rational can be constructed
directly via `BigRational::new_raw` without GCD normalization — trailing zeros in the
significand are stripped first so the fraction is already in lowest terms.

See references [9-10] below.

9. IEEE Computer Society. "IEEE Standard for Floating-Point Arithmetic." *IEEE Std 754-2019*
(Revision of IEEE 754-2008), 2019.
[DOI](https://doi.org/10.1109/IEEESTD.2019.8766229)
Expand Down
7 changes: 2 additions & 5 deletions src/exact.rs
Original file line number Diff line number Diff line change
Expand Up @@ -63,10 +63,7 @@ fn validate_finite<const D: usize>(m: &Matrix<D>) -> Result<(), LaError> {
for c in 0..D {
if !m.rows[r][c].is_finite() {
cold_path();
return Err(LaError::NonFinite {
row: Some(r),
col: c,
});
return Err(LaError::non_finite_cell(r, c));
}
}
}
Expand All @@ -81,7 +78,7 @@ fn validate_finite_vec<const D: usize>(v: &Vector<D>) -> Result<(), LaError> {
for (i, &x) in v.data.iter().enumerate() {
if !x.is_finite() {
cold_path();
return Err(LaError::NonFinite { row: None, col: i });
return Err(LaError::non_finite_at(i));
}
}
Ok(())
Expand Down
Loading
Loading