From 31b7b1e7c547fc69edbb601e68b2f29e83d5cd24 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Sun, 12 Jul 2026 09:32:38 -0700 Subject: [PATCH 1/2] docs(science): add mathematical basis and clarify numerical guarantees - define the binary64, factorization, determinant-filter, exact-arithmetic, tolerance, and error contracts - align README, API documentation, examples, and LDLT diagnostics with guarantees over stored values - audit academic references and extend Criterion coverage for determinant error-bound paths Closes #159 Closes #164 --- CONTRIBUTING.md | 1 + README.md | 153 +++++++++++-------- REFERENCES.md | 76 ++++++---- benches/common/exact.rs | 33 +++-- benches/exact.rs | 55 ++++++- benches/vs_linalg.rs | 3 +- docs/BENCHMARKING.md | 18 ++- docs/mathematical_basis.md | 220 ++++++++++++++++++++++++++++ examples/exact_det_3x3.rs | 4 +- examples/exact_sign_3x3.rs | 4 +- examples/exact_solve_3x3.rs | 10 +- justfile | 2 +- scripts/bench_compare.py | 40 ++++- scripts/tests/test_bench_compare.py | 37 ++++- src/error.rs | 30 ++-- src/exact.rs | 50 ++++--- src/ldlt.rs | 23 ++- src/lib.rs | 18 --- src/lu.rs | 7 + src/matrix.rs | 43 ++++-- tests/exact_bench_config.rs | 10 ++ tests/vs_linalg_inputs.rs | 25 +++- typos.toml | 2 +- 23 files changed, 670 insertions(+), 194 deletions(-) create mode 100644 docs/mathematical_basis.md diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 81f0844..8fa399c 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -69,6 +69,7 @@ Use the existing canonical documents instead of duplicating their guidance: |-------|---------------------| | Agent rules and repository invariants | [`AGENTS.md`](AGENTS.md) | | User-facing API, examples, and project scope | [`README.md`](README.md) | +| Mathematical basis and numerical validity | [`docs/mathematical_basis.md`](docs/mathematical_basis.md) | | Package metadata, features, and dependencies | [`Cargo.toml`](Cargo.toml) | | Commands and validation workflow | [`justfile`](justfile), `just --list` | | Python support tooling | [`scripts/README.md`](scripts/README.md) | diff --git a/README.md b/README.md index 38f162c..89f44b0 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,9 @@ # la-stack [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.18158926.svg)](https://doi.org/10.5281/zenodo.18158926) -[![Crates.io](https://img.shields.io/crates/v/la-stack.svg)](https://crates.io/crates/la-stack) -[![Downloads](https://img.shields.io/crates/d/la-stack.svg)](https://crates.io/crates/la-stack) -[![License](https://img.shields.io/crates/l/la-stack.svg)](https://github.com/acgetchell/la-stack/blob/v0.4.3/LICENSE) +[![Crates.io](https://badgen.net/crates/v/la-stack)](https://crates.io/crates/la-stack) +[![Downloads](https://badgen.net/crates/d/la-stack)](https://crates.io/crates/la-stack) +[![License](https://badgen.net/github/license/acgetchell/la-stack)](https://github.com/acgetchell/la-stack/blob/v0.4.3/LICENSE) [![Docs.rs](https://docs.rs/la-stack/badge.svg)](https://docs.rs/la-stack) [![CI](https://github.com/acgetchell/la-stack/actions/workflows/ci.yml/badge.svg)](https://github.com/acgetchell/la-stack/actions/workflows/ci.yml) [![rust-clippy analyze][clippy-badge]][clippy-workflow] @@ -24,18 +24,40 @@ while keeping the API intentionally small and explicit. - `Vector` for fixed-length `f64` vectors backed by `[f64; D]` - `Matrix` for fixed-size square `f64` matrices backed by `[[f64; D]; D]` - `Lu` for LU factorization with partial pivoting (solve + det) -- `Ldlt` for LDLT factorization without pivoting (solve + det; symmetric SPD/PSD) +- `Ldlt` for no-pivot factorization intended for exactly + symmetric positive-definite matrices (solve + det; typed pivot diagnostics) + +## 🧮 Mathematical basis + +`la-stack` operates on finite IEEE 754 binary64 values in small, fixed +dimensions. Its floating-point paths use LU with partial pivoting, LDLT without +pivoting for exactly symmetric positive-definite matrices, and closed-form +determinants through D=4. These results remain subject to conditioning and +binary64 rounding; +factorization tolerances are rejection thresholds, not accuracy guarantees. For +D≤4, direct determinants can be paired with a conservative absolute roundoff +bound when its range preconditions hold. + +With `features = ["exact"]`, stored binary64 inputs are lifted losslessly to +rationals for exact determinant signs, determinant values, and solves. Exactness +starts at the stored values and cannot recover information rounded away before +construction. See the +[mathematical basis](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md) +for the algorithms, validity boundaries, and supporting references. ## ✨ Design goals - ✅ `Copy` types where possible -- ✅ Const-generic dimensions (no dynamic sizes) +- ✅ Const-generic storage (no dynamically sized matrix or vector representation) - ✅ `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`, strict/rounded f64 conversions) +- ✅ Error-bounded f64 determinant filtering plus optional exact signs + (`det_errbound`, `det_sign_exact`) +- ✅ Exact determinant values and linear solves via optional arbitrary-precision + arithmetic (`det_exact`, `solve_exact`, strict/rounded f64 conversions) - ✅ No runtime dependencies by default (optional features may add deps) -- ✅ Stack storage only (no heap allocation in core types) +- ✅ Inline, stack-backed storage for core types; optional arbitrary-precision + exact values allocate as required - ✅ `unsafe` forbidden See [CHANGELOG.md](https://github.com/acgetchell/la-stack/blob/v0.4.3/CHANGELOG.md) @@ -47,8 +69,8 @@ for current release planning. - Bare-metal performance: see [`blas-src`](https://crates.io/crates/blas-src), [`lapack-src`](https://crates.io/crates/lapack-src), or [`openblas-src`](https://crates.io/crates/openblas-src) -- Comprehensive: use [`nalgebra`](https://crates.io/crates/nalgebra) if you need a full-featured library -- Large matrices/dimensions with parallelism: use [`faer`](https://crates.io/crates/faer) if you need this +- Broad general-purpose linear algebra: use [`nalgebra`](https://crates.io/crates/nalgebra) +- Large matrices/dimensions with parallelism: use [`faer`](https://crates.io/crates/faer) - Alternate floating-point scalar families: `la-stack` supports `f64` and optional exact arithmetic, not `f32` / `f16` APIs ## ✅ Use this crate when @@ -87,8 +109,10 @@ la-stack = "0.4.3" - `default`: no runtime dependencies - `exact`: `BigRational` exact determinant and solve APIs -- `bench`: cfg-only gate for benchmark fixtures and benchmark-input tests; - benchmark libraries remain development dependencies +- `bench`: repository-development gate used only by benchmark targets and + benchmark-input tests; application crates should not enable it + +### LU solve Solve a 5×5 system via LU: @@ -121,10 +145,14 @@ fn main() -> Result<(), LaError> { } ``` -Compute a determinant for a symmetric SPD matrix via LDLT (no pivoting). +### LDLT determinant + +Compute a determinant for a symmetric positive-definite 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`): +For these matrices, `LDLᵀ` is a square-root-free Cholesky form. Multiplying each +column of `L` by the square root of the corresponding diagonal entry yields a +Cholesky factor: ```rust use la_stack::prelude::*; @@ -172,20 +200,24 @@ fn main() -> Result<(), LaError> { > required allowed difference of zero. The tolerance-based > [`Matrix::first_asymmetry`](https://docs.rs/la-stack/latest/la_stack/struct.Matrix.html#method.first_asymmetry) > and `Matrix::is_symmetric` methods remain useful diagnostics, but do not prove -> the exact precondition required by LDLT. Fall back to `lu()` if your matrices -> may not be symmetric at all. A negative LDLT diagonal or a zero diagonal with nonzero -> remaining coupling returns `LaError::NotPositiveSemidefinite` with a typed -> `PositiveSemidefiniteViolation`. An uncoupled zero or other non-negative pivot +> the exact precondition required by LDLT. Use `lu()` when exact symmetry or +> positive definiteness is not guaranteed. A negative LDLT diagonal or a zero +> diagonal with nonzero remaining coupling returns +> `LaError::NotPositiveSemidefinite` with a typed +> `PositiveSemidefiniteViolation`. An uncoupled zero or positive pivot > at or below the caller's tolerance returns `LaError::Singular` with a -> numerical `SingularityReason`. +> numerical `SingularityReason`. Because these pivots are computed in binary64, +> success is not an exact positive-definiteness certificate for the stored +> matrix. ## ⚡ 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: +using fused multiply-add where applicable. It returns `Ok(Some(det))` for those +dimensions and `Ok(None)` for D ≥ 5. `Matrix::<0>::zero().det_direct()` returns +`Ok(Some(1.0))` (the empty-product convention). For D=1–4, direct formulas +bypass LU factorization entirely. This enables compile-time evaluation when +inputs are known: ```rust use la_stack::prelude::*; @@ -227,6 +259,12 @@ rationals (this pulls in `num-bigint`, `num-rational`, and `num-traits` for la-stack = { version = "0.4.3", features = ["exact"] } ``` +These routines are exact with respect to the finite binary64 values stored in +`Matrix` and `Vector`. They treat each stored value as the exact rational number +represented by its bits, so the exact determinant or solve stage introduces no +further roundoff. They cannot recover information already lost when source +values were rounded to `f64` before construction. + **Determinants:** - **`det_exact()`** — returns the exact determinant as a `BigRational` @@ -248,6 +286,12 @@ la-stack = { version = "0.4.3", features = ["exact"] } - **`ExactF64Conversion`** — converts an existing exact determinant or solution under the strict or rounded contract without repeating exact elimination +For exact-to-f64 output, strict conversions use +`UnrepresentableReason::RequiresRounding` when explicit rounding can produce a +finite value and `UnrepresentableReason::NotFinite` otherwise. Rounded +conversions opt into nearest-even rounding but still report `NotFinite` when no +finite `f64` exists. + ```rust,ignore use la_stack::prelude::*; @@ -328,12 +372,13 @@ filter and uses fraction-free Bareiss elimination in `BigInt`. Because `Matrix` stores only finite entries, arithmetic range failures in the filter are inconclusive rather than errors and the exact fallback is total. -### Adaptive precision with `det_direct_with_errbound()` +## 🛡️ Adaptive determinant filtering (D ≤ 4) `det_direct_with_errbound()` returns a closed-form determinant together with the conservative absolute error bound used by the fast filter, computed from -one shared traversal. It returns `None` when a D ≤ 4 computation may be -affected by gradual underflow, as well as for unsupported D ≥ 5 dimensions. +one call that evaluates the determinant once and computes its matching bound. +It returns `None` when a D ≤ 4 computation may be affected by gradual +underflow, as well as for unsupported D ≥ 5 dimensions. This method does NOT require the `exact` feature — it uses pure f64 arithmetic and is available by default. Use `det_errbound()` when only the bound is needed. The paired API enables custom adaptive-precision logic for geometric predicates: @@ -390,28 +435,11 @@ fn main() -> Result<(), LaError> { } ``` -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**, equivalently the permanent of `|A|`: - -```text -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: - -```text -|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 when gradual underflow is absent. They are -explicit crate-root exports for -advanced users who want to compose the same bound themselves: +The error coefficients (`ERR_COEFF_2`, `ERR_COEFF_3`, `ERR_COEFF_4`) are +conservative, dimension-specific constants, not caller-tunable tolerances. The +[mathematical basis](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#determinants-and-certified-sign-filtering) +documents the bound and states its range preconditions. The constants are explicit +crate-root exports for advanced users who want to compose the same bound: `use la_stack::{ERR_COEFF_2, ERR_COEFF_3, ERR_COEFF_4};`. They intentionally stay out of the common prelude. @@ -421,15 +449,20 @@ out of the common prelude. |---|---|---|---| | `Vector` | `[f64; D]` | Finite fixed-length vector for input and computation | `try_new`, `as_array`, `into_array`, `dot`, `norm2_sq` | | `Matrix` | `[[f64; D]; D]` | Finite square matrix for input and computation | See below | -| `DeterminantWithErrorBound` | two private `f64` fields | Paired direct determinant and certified absolute bound | `determinant`, `absolute_error_bound` | -| `Lu` | `Matrix` + pivot array | Factorization for solves/det | `solve`, `det` | -| `Ldlt` | `Matrix` | Factorization for symmetric SPD/PSD solves/det | `solve`, `det` | +| `DeterminantWithErrorBound` | Opaque validated pair | Paired direct determinant and certified absolute bound | `determinant`, `absolute_error_bound` | +| `Lu` | Inline factors + permutation | Factorization for solves/det | `solve`, `det` | +| `Ldlt` | Inline factors | No-pivot SPD factorization for solves/det | `solve`, `det` | | `Tolerance` | finite non-negative `f64` | Validated numerical threshold | `try_new`, `get` | -| `LaError` | typed variants and reasons | Structured, actionable failure reporting | See error enums below | +| `LaError` | typed variants and reasons | Structured, actionable failure reporting | See error semantics below | | `DeterminantSign`¹ | enum | Exact determinant sign | `as_i8` | Storage shown above reflects the intentional `f64` scalar model. +For a runtime dimension from 0 through `MAX_STACK_MATRIX_DISPATCH_DIM` (7), +`try_with_stack_matrix!` dispatches to a concrete `Matrix` while preserving +inline stack storage. Larger dimensions return `LaError::UnsupportedDimension`; +the macro does not introduce a dynamically sized matrix representation. + `Matrix` key methods: `as_rows`, `into_rows`, `lu`, `ldlt`, `det`, `det_direct`, `det_direct_with_errbound`, `det_errbound`, `det_exact`¹, `det_exact_f64`¹, `det_exact_rounded_f64`¹, `det_sign_exact`¹, @@ -444,12 +477,13 @@ observable results. `Matrix::into_rows` and `Vector::into_array` consume the value and return the owned fixed-size arrays. -`Matrix::get` returns `Option` for bounds-only access; `Matrix::try_get` -preserves invalid coordinates in `LaError`. The single fallible `Matrix::set` -validates both bounds and finiteness before mutating the matrix. +`Matrix::get(row, col)` returns `None` for out-of-bounds coordinates; +`Matrix::try_get` instead returns a structured `LaError` preserving those +coordinates. The single fallible `Matrix::set` validates both coordinates and +finiteness before mutating the matrix. `LaError` and its reason/location enums are non-exhaustive. Numerical -singularity records the [`FactorizationKind`](https://docs.rs/la-stack/latest/la_stack/enum.FactorizationKind.html), +singularity records the `FactorizationKind`, observed pivot magnitude, and tolerance, while exact-arithmetic singularity is identified separately. `LaError::NonFinite` retains the crate-wide non-finite contract but uses `NonFiniteOrigin`, `NonFiniteLocation`, and @@ -572,7 +606,8 @@ For the full contributor workflow, see If you use this library in academic work, please cite it using [CITATION.cff](https://github.com/acgetchell/la-stack/blob/v0.4.3/CITATION.cff) (or GitHub's "Cite this repository" feature). Tagged releases are archived on -Zenodo. +Zenodo under the +[all-versions concept DOI](https://doi.org/10.5281/zenodo.18158926). ## 📚 References @@ -593,7 +628,7 @@ BSD 3-Clause License. See [LICENSE](https://github.com/acgetchell/la-stack/blob/ [audit-badge]: https://github.com/acgetchell/la-stack/actions/workflows/audit.yml/badge.svg [audit-workflow]: https://github.com/acgetchell/la-stack/actions/workflows/audit.yml -[benchmark-provenance]: https://github.com/acgetchell/la-stack/blob/v0.4.3/docs/assets/bench/vs_linalg_lu_solve_median.provenance.json +[benchmark-provenance]: https://github.com/acgetchell/la-stack/blob/668daed6/docs/assets/bench/vs_linalg_lu_solve_median.provenance.json [clippy-badge]: https://github.com/acgetchell/la-stack/actions/workflows/rust-clippy.yml/badge.svg [clippy-workflow]: https://github.com/acgetchell/la-stack/actions/workflows/rust-clippy.yml [lu-solve-benchmark]: https://raw.githubusercontent.com/acgetchell/la-stack/v0.4.3/docs/assets/bench/vs_linalg_lu_solve_median.svg diff --git a/REFERENCES.md b/REFERENCES.md index e54134d..0733cbe 100644 --- a/REFERENCES.md +++ b/REFERENCES.md @@ -15,6 +15,7 @@ Tagged releases are archived on Zenodo under the all-versions concept DOI - CodeRabbit AI, Inc. "CodeRabbit." . - KiloCode. "KiloCode AI Engineering Assistant." . - OpenAI. "ChatGPT." . +- OpenAI. "Codex." . - Warp Dev, Inc. "WARP." . All AI-generated output was reviewed and/or edited by the maintainer. @@ -62,50 +63,62 @@ matrix and RHS. After back-substitution, multiplying by the exact power-of-two s Both the determinant and solve paths convert their finite-by-construction entries via `decompose_proven_finite_f64`, which extracts the IEEE 754 binary64 sign, unbiased exponent, -and significand [9] and strips trailing zeros from the significand so `|x| = m · 2^e` with -`m` odd. The integer matrix is then assembled by shifting each mantissa left by +and significand [9]. For nonzero `x`, it strips trailing zeros from the +significand so `|x| = m · 2^e` with `m` odd; signed zeros use a separate zero +component. The integer matrix is then assembled by shifting each mantissa left by `exp − e_min`, giving a GCD-free exact-integer starting point. Solves and D ≥ 5 determinants then apply Bareiss elimination; D ≤ 4 determinants use direct expansions. The test-only fallible wrapper `decompose_f64` verifies rejection of non-finite raw scalars, while the test-only `f64_to_big_rational` helper packages the same decomposition into a single -`BigRational`. See Goldberg [10] for background on IEEE 754 representation and exact rational -reconstruction. +`BigRational`. See Goldberg [10] for background on floating-point representation and +conversion. -### LDL^T factorization (symmetric SPD/PSD) +### LDLᵀ factorization (exactly symmetric positive-definite inputs) -The LDL^T (often abbreviated "LDLT") implementation in `la-stack` is intended for symmetric positive -definite (SPD) and positive semi-definite (PSD) matrices (e.g. Gram matrices), and does not perform -pivoting. +The no-pivot LDLT implementation targets `A = L D Lᵀ` for exactly symmetric +positive-definite inputs [4-5, 11-12]. Successful construction requires every +computed diagonal pivot to be positive and greater than the caller's tolerance. +Computed zero and tolerance-small positive pivots are therefore part of the +typed diagnostic domain, not returned in a usable factorization. Because the +pivots are computed in binary64, a successful factorization is not an exact +certificate that the represented matrix is positive definite. -For background on the SPD/PSD setting, see [4-5]. For pivoted variants used for symmetric *indefinite* -matrices, see [6]. +For pivoted variants used for symmetric *indefinite* matrices, see [6, 11-12]. ### LU decomposition (Gaussian elimination with partial pivoting) -The LU implementation in `la-stack` follows the standard Gaussian elimination / LU factorization -approach with partial pivoting for numerical stability. +The LU implementation targets `P A = L U` and uses partial pivoting: each step +selects the remaining entry of largest magnitude in the active column. Partial +pivoting is a practical stability strategy, not an unconditional accuracy +guarantee; worst-case growth and average-case behavior are distinct concerns. -See references [1-3] below. +See [1-3, 11-12] for stability analysis, finite-precision behavior, and standard +algorithmic background. ## 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. + [DOI](https://doi.org/10.1137/0611023) · [PDF](https://people.maths.ox.ac.uk/trefethen/publication/PDF/1990_44.pdf) 2. Businger, P. A. "Monitoring the Numerical Stability of Gaussian Elimination." - *Numerische Mathematik* 16 (1970/71): 360–361. - [Full text](https://eudml.org/doc/132040) + *Numerische Mathematik* 16.4 (1971): 360–361. + [DOI](https://doi.org/10.1007/BF02165006) · [Full text](https://eudml.org/doc/132040) 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)) -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) -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) + [DOI](https://doi.org/10.1007/s00440-024-01276-2) · + [Open-access article](https://link.springer.com/article/10.1007/s00440-024-01276-2) · + [arXiv:2206.01726](https://arxiv.org/abs/2206.01726) +4. Cholesky, André-Louis. "Sur la résolution numérique des systèmes d'équations linéaires." + *Bulletin de la Sabix* 39 (2005): 81–95. Manuscript dated 2 December 1910. + [DOI](https://doi.org/10.4000/sabix.529) +5. Brezinski, Claude. "La méthode de Cholesky." + *Revue d'histoire des mathématiques* 11.2 (2005): 205–238. + [DOI](https://doi.org/10.24033/rhm.30) · + [Full text](https://www.numdam.org/articles/10.24033/rhm.30/) +6. Bunch, James R., Linda Kaufman, and Beresford N. Parlett. "Decomposition of a Symmetric Matrix." + *Numerische Mathematik* 27 (1976): 95–109. + [DOI](https://doi.org/10.1007/BF01399088) · [Full text](https://eudml.org/doc/132435) 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) · @@ -123,6 +136,15 @@ See references [1-3] below. 10. Goldberg, David. "What Every Computer Scientist Should Know About Floating-Point Arithmetic." *ACM Computing Surveys* 23.1 (1991): 5–48. [DOI](https://doi.org/10.1145/103162.103163) · - [PDF](https://www.validlab.com/goldberg/paper.pdf) - Comprehensive survey of IEEE 754 representation, rounding, and exact rational - reconstruction of floating-point values. + [Authorized HTML reprint](https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html) + Comprehensive survey of floating-point representation, rounding, and conversion. +11. Higham, Nicholas J. *Accuracy and Stability of Numerical Algorithms*. 2nd ed. + Society for Industrial and Applied Mathematics, 2002. + [DOI](https://doi.org/10.1137/1.9780898718027) +12. Golub, Gene H., and Charles F. Van Loan. *Matrix Computations*. 4th ed. + Johns Hopkins University Press, 2013. + [DOI](https://doi.org/10.56021/9781421407944) · + [Publisher record](https://www.press.jhu.edu/books/title/10678/matrix-computations) +13. Kalibera, Tomas, and Richard Jones. "Rigorous Benchmarking in Reasonable Time." + *Proceedings of the 2013 International Symposium on Memory Management* (ISMM '13), + 2013: 63–74. [DOI](https://doi.org/10.1145/2464157.2464160) diff --git a/benches/common/exact.rs b/benches/common/exact.rs index 328bddb..8b58bba 100644 --- a/benches/common/exact.rs +++ b/benches/common/exact.rs @@ -493,9 +493,9 @@ fn assert_approximate_determinant(actual: f64, exact: &BigRational, operation: & /// Validate the floating-point determinant operations used by Criterion. /// /// This runs during setup, outside timed closures. The deterministic `det` and -/// `det_direct` results are compared with the independent Leibniz oracle; on -/// current revisions the combined direct result is additionally checked -/// against its certified absolute bound. +/// `det_direct` results are compared with the independent Leibniz oracle, and +/// `det_errbound` must certify the observed direct error. On current revisions, +/// the paired result must also match the separate determinant and bound calls. /// /// # Panics /// @@ -513,11 +513,24 @@ pub fn validate_f64_determinant_benchmarks(input: &ValidatedExac .matrix() .det_direct() .or_abort("direct f64 determinant oracle check"); + let standalone_bound = input + .matrix() + .det_errbound() + .or_abort("standalone determinant error-bound oracle check"); if D <= 4 { let Some(direct) = direct else { panic!("det_direct must support benchmark dimension {D}"); }; + let Some(standalone_bound) = standalone_bound else { + panic!("det_errbound must support benchmark dimension {D}"); + }; assert_approximate_determinant(direct, &exact, "direct f64 determinant"); + let observed_error = (rational_from_f64(direct) - &exact).abs(); + let certified_bound = rational_from_f64(standalone_bound); + assert!( + observed_error <= certified_bound, + "direct determinant error {observed_error} exceeds standalone certified bound {certified_bound}", + ); #[cfg(not(la_stack_v0_4_3_api))] { @@ -529,19 +542,21 @@ pub fn validate_f64_determinant_benchmarks(input: &ValidatedExac panic!("the baseline fixture must have a certified D={D} determinant bound"); }; assert_eq!(estimate.determinant().to_bits(), direct.to_bits()); + assert_eq!( + estimate.absolute_error_bound().to_bits(), + standalone_bound.to_bits(), + ); assert!( estimate.determinant().abs() > estimate.absolute_error_bound(), "the headline D={D} det_sign_exact benchmark must exercise the fast filter", ); - let observed_error = (rational_from_f64(direct) - &exact).abs(); - let certified_bound = rational_from_f64(estimate.absolute_error_bound()); - assert!( - observed_error <= certified_bound, - "direct determinant error {observed_error} exceeds certified bound {certified_bound}", - ); } } else { assert!(direct.is_none(), "det_direct unexpectedly supports D={D}"); + assert!( + standalone_bound.is_none(), + "det_errbound unexpectedly supports D={D}", + ); #[cfg(not(la_stack_v0_4_3_api))] assert!( diff --git a/benches/exact.rs b/benches/exact.rs index 3f71646..a004754 100644 --- a/benches/exact.rs +++ b/benches/exact.rs @@ -8,7 +8,8 @@ //! 1. **General-case benches** (`exact_d{2..5}`) — a single //! well-conditioned diagonally-dominant matrix per dimension. These //! measure typical-case performance and track regressions against a -//! reproducible input. +//! reproducible input. D=2..=4 also include the direct determinant and +//! certified error-bound f64 baselines. //! 2. **Adversarial / extreme-input benches** — matrices chosen to //! stress specific corners of the exact-arithmetic pipeline: //! near-singularity (forces the exact integer fallback), large f64 entries @@ -200,9 +201,57 @@ fn bench_det_direct( }); } -macro_rules! register_det_direct_benchmark { +/// Add the standalone certified determinant-bound baseline. +fn bench_det_errbound( + group: &mut BenchmarkGroup<'_, WallTime>, + input: &ValidatedExactInput, +) { + let bound = input + .matrix() + .det_errbound() + .or_abort("determinant error-bound setup") + .or_abort("determinant error-bound setup"); + black_box(bound); + group.bench_function("det_errbound", |bencher| { + bencher.iter(|| { + let bound = black_box(input.matrix()) + .det_errbound() + .or_abort("f64 determinant error bound") + .or_abort("f64 determinant error bound"); + black_box(bound); + }); + }); +} + +/// Add the paired direct-determinant and certified-bound baseline. +#[cfg(not(la_stack_v0_4_3_api))] +fn bench_det_direct_with_errbound( + group: &mut BenchmarkGroup<'_, WallTime>, + input: &ValidatedExactInput, +) { + let estimate = input + .matrix() + .det_direct_with_errbound() + .or_abort("paired determinant-bound setup") + .or_abort("paired determinant-bound setup"); + black_box((estimate.determinant(), estimate.absolute_error_bound())); + group.bench_function("det_direct_with_errbound", |bencher| { + bencher.iter(|| { + let estimate = black_box(input.matrix()) + .det_direct_with_errbound() + .or_abort("paired f64 determinant and error bound") + .or_abort("paired f64 determinant and error bound"); + black_box((estimate.determinant(), estimate.absolute_error_bound())); + }); + }); +} + +macro_rules! register_det_filter_benchmarks { ($group:expr, $input:expr, supported) => {{ bench_det_direct(&mut $group, &$input); + #[cfg(not(la_stack_v0_4_3_api))] + bench_det_direct_with_errbound(&mut $group, &$input); + bench_det_errbound(&mut $group, &$input); }}; ($group:expr, $matrix:expr, unsupported) => {}; } @@ -227,7 +276,7 @@ macro_rules! gen_exact_benches_for_dim { }); }); - register_det_direct_benchmark!(group, input, $direct); + register_det_filter_benchmarks!(group, input, $direct); for &operation in GENERAL_OPERATIONS { bench_exact_operation(&mut group, operation, &input); diff --git a/benches/vs_linalg.rs b/benches/vs_linalg.rs index bc74007..da484f2 100644 --- a/benches/vs_linalg.rs +++ b/benches/vs_linalg.rs @@ -6,7 +6,8 @@ //! fixed dimensions. //! //! Notes: -//! - Determinant is benchmarked via LU on all sides (nalgebra uses closed-forms for 1×1/2×2/3×3). +//! - Determinant groups distinguish factorization-inclusive LU, `Matrix::det`, +//! precomputed LU determinant queries, and precomputed LDLT/Cholesky queries. //! - Matrix infinity norm is the maximum absolute row sum on all sides. use std::hint::black_box; diff --git a/docs/BENCHMARKING.md b/docs/BENCHMARKING.md index d1db7d5..162852d 100644 --- a/docs/BENCHMARKING.md +++ b/docs/BENCHMARKING.md @@ -56,8 +56,9 @@ factorization in the dependency version used here. **`exact`** (`benches/exact.rs`) measures exact-arithmetic methods (`det_exact`, `solve_exact`, `det_sign_exact`, strict `*_result` conversions, and lossy `*_rounded_f64` conversions) alongside the f64 `det` baseline across -D=2-5 and `det_direct` across its supported D=2-4 range. Use this suite to -understand exact-arithmetic cost and track optimization progress. +D=2-5. Its supported D=2-4 range also includes `det_direct`, the paired +`det_direct_with_errbound`, and the bound-only `det_errbound`. Use this suite +to understand exact-arithmetic cost and track optimization progress. ## Common Workflows @@ -91,6 +92,11 @@ at warning for both revisions because the current manifest's lint policy may reject historical source that predates a lint, even though that source remains valid benchmark input. +The paired `det_direct_with_errbound` API postdates v0.4.3, so its D=2-4 +baseline rows remain explicitly unavailable in that comparison. The older +bound-only `det_errbound` API is source-compatible and remains a required +shared-harness baseline. + The v0.4.3 LU/LDLT balanced-range determinant paths return an incorrect zero, so their two D=8 stress rows are deliberately not timed as baselines. Reports leave those baselines explicitly unavailable rather than presenting invalid @@ -179,7 +185,7 @@ coverage or provenance aborts publication. Use `--allow-partial` only for exploratory CSV/SVG output; it cannot update README and its sidecar explicitly marks measurement provenance unavailable. -See `scripts/criterion_dim_plot.py --help` for plotting options. +See `uv run --locked criterion-dim-plot --help` for plotting options. ### Create The Release Performance Report @@ -247,6 +253,9 @@ The README table uses `median.point_estimate` in nanoseconds. Lower is better, but point-estimate ratios alone are descriptive and do not establish a statistically supported performance difference. Preserve Criterion confidence intervals or repeat controlled runs when making a stronger claim. +For experimental-design background on controlled repetitions and uncertainty, +see [REFERENCES.md](../REFERENCES.md) \[13\]; these workflows do not claim to +implement every recommendation in that study. All three crates receive equivalent deterministic inputs for a given dimension: @@ -341,7 +350,8 @@ Before timing begins, every fixed, adversarial, and corpus input is consumed int a private-field `ValidatedExactInput` after checks by an independent exact oracle. Timed and registration helpers accept only that proof-bearing wrapper. A factorial-time Leibniz determinant over exact rational reconstructions verifies -determinant values and signs; exact residuals verify `A x = b`; and +determinant values and signs, including each direct determinant's certified +absolute bound; exact residuals verify `A x = b`; and strict/rounded binary64 results are checked for their exact bits, typed reason, and first failing component. These checks run outside timed Criterion closures. Any disagreement or unexpected error fails setup instead of becoming an diff --git a/docs/mathematical_basis.md b/docs/mathematical_basis.md new file mode 100644 index 0000000..da24132 --- /dev/null +++ b/docs/mathematical_basis.md @@ -0,0 +1,220 @@ +# Mathematical basis + +`la-stack` provides fixed-dimension numerical linear algebra over finite IEEE +754 binary64 values. For a compile-time dimension `D`, `Matrix` is a dense +`D × D` square matrix and `Vector` is a length-`D` vector over the finite +binary64 set `F`. The default algorithms operate in binary64 and are therefore +approximate. The optional `exact` feature instead lifts each stored binary64 +value to the exact rational number it represents. + +This document separates three questions that are easy to conflate: + +1. Which mathematical factorization or determinant identity is being used? +2. Which part of the computation is rounded in binary64? +3. What does a tolerance or typed error prove about the supplied values? + +Reference numbers point to [REFERENCES.md](../REFERENCES.md). + +## Represented values and arithmetic model + +Public matrix and vector construction rejects NaN and infinity. Because the +backing arrays are private and mutation is validated, every stored entry is +finite. Arithmetic can still overflow, underflow, or round; computed non-finite +values are reported through `LaError::NonFinite` rather than stored in a public +matrix, vector, or factorization. + +Every finite binary64 value is a dyadic rational with an exact representation + +```text +x = (-1)^s m × 2^e, +``` + +for integer `m` and exponent `e` \[9-11\]. Both signed-zero bit patterns +represent rational zero. Exact APIs exploit this representation, but exactness +begins only after construction: they cannot recover information lost when a +decimal value or an earlier computation was rounded to `f64`. + +`Matrix`, `Vector`, `Lu`, and `Ldlt` use inline fixed-size storage. +Arbitrary-precision `BigInt` and `BigRational` values allocate when the `exact` feature is +used. Const-generic `D` is not itself a mathematical dimension limit; “small, +fixed dimension” is the intended performance scope. `try_with_stack_matrix!` is +separately limited to `D = 0..=MAX_STACK_MATRIX_DISPATCH_DIM` (currently 7) +because it enumerates concrete stack types. + +Except for the determinant filter described below, the floating-point APIs do +not provide certified forward, backward, or absolute error bounds. This includes +dot products, squared norms, matrix norms, factorizations, and solves. Some +kernels use FMA to reduce rounding steps, but that does not make them exact. + +## Floating-point factorizations + +### LU with partial pivoting + +LU factorization targets + +```text +P A = L U, +``` + +where `P` is a row permutation, `L` is unit lower triangular, and `U` is upper +triangular. At column `k`, the implementation selects the largest-magnitude +remaining entry in the active column. It rejects the factorization when that +magnitude is less than or equal to the caller's tolerance. + +The tolerance is an absolute finite, non-negative threshold. It is not divided +by a matrix norm, so rescaling a system can change whether a pivot is accepted. +Partial pivoting is a practical stability strategy, not an unconditional +accuracy guarantee: worst-case element growth and typical behavior are distinct +questions \[1-3, 11-12\]. `Lu::solve` does not estimate conditioning or refine +the result. + +`Lu::det` combines permutation parity with the product of the diagonal of `U`. +Scaled product accumulation avoids some premature overflow and underflow, but +the final binary64 determinant is still rounded. A returned zero or a numerical +`LaError::Singular` is therefore not proof that the represented matrix is +exactly singular. + +### LDLT without pivoting + +LDLT factorization targets + +```text +A = L D Lᵀ, +``` + +with unit lower-triangular `L` and diagonal `D` \[4-5, 11-12\]. The input +must be exactly symmetric under binary64 comparison: every mirrored pair must +satisfy `A[i][j] == A[j][i]`. Signed zeros compare equal and are accepted. + +A successful `Ldlt` requires every computed diagonal pivot to be positive and +greater than the caller's tolerance. An uncoupled computed zero or a positive +pivot at or below tolerance returns `LaError::Singular`; a negative pivot or +zero pivot with remaining coupling returns `LaError::NotPositiveSemidefinite` +with a typed violation. + +This is not a pivoted symmetric-indefinite factorization such as Bunch-Kaufman +\[6, 11-12\]. Pivots are computed in binary64, so a singular represented matrix +can produce a small positive pivot above a low tolerance. Successful +factorization is therefore not an exact positive-definiteness certificate for +the stored matrix, much less for ideal values before binary64 input conversion. + +## Determinants and certified sign filtering + +`det_direct()` evaluates closed forms for `D = 0..=4`, with the empty-product +convention `det(Matrix::<0>) = 1`. `Matrix::det()` uses that path through `D = 4` +and zero-tolerance LU for `D ≥ 5`. The general `det()` result has no certified +roundoff bound, and its LU fallback can report numerical singularity even when +exact arithmetic over the stored entries would find a nonzero determinant. + +For `D = 0` and `D = 1`, `det_direct_with_errbound()` returns the exact direct +determinant with a zero bound. For `D = 2..=4`, it can return a determinant +estimate and a conservative absolute bound. Let `ι(A)` denote the exact rational +lift of the stored entries and let + +```text +p(|A|) = Σ_(σ ∈ S_D) Π_i |a_(i,σ(i))| = perm(|A|). +``` + +When every rounded intermediate in the implemented determinant and `p(|A|)` +evaluation trees is normal or an exact structural zero, the implementation uses + +```text +|det_direct(A) - det(ι(A))| ≤ c_D × p(|A|), +``` + +with `ε = f64::EPSILON` and project-specific coefficients + +```text +c_2 = 3ε + 16ε² +c_3 = 8ε + 64ε² +c_4 = 12ε + 128ε². +``` + +The analysis follows Shewchuk's adaptive-filter framework \[8\], while the +FMA evaluation trees and constants are derived for this crate rather than copied +from that source. IEEE 754 and standard floating-point error analysis provide +the arithmetic model \[9-11\]. + +If `|determinant| > absolute_error_bound`, the sign is certified. Otherwise the +filter is inconclusive; it does not prove that the matrix is singular. The API +returns `Ok(None)` for `D ≥ 5` or when gradual underflow could invalidate the +relative-error model. A non-finite computed determinant or bound returns +`LaError::NonFinite`. + +## Exact arithmetic over binary64 inputs + +The `exact` feature decomposes each stored entry into an integer mantissa and a +power of two. The entries are scaled to integer matrices without changing their +represented rational values \[9-10\]. + +Exact determinants use direct `BigInt` expansions for `D ≤ 4` and fraction-free +Bareiss elimination for `D ≥ 5` \[7\]. Exact solves apply Bareiss updates to an +integer augmented system, then use `BigRational` for back-substitution. Matrix +and right-hand-side scaling are tracked separately and reconciled with an exact +power-of-two factor. First-nonzero pivoting is sufficient for correctness in +exact arithmetic, although pivot choice can still affect computational cost. + +`det_sign_exact()` first attempts the certified binary64 filter for `D ≤ 4` and +falls back to exact integer arithmetic when the filter is inconclusive. It +returns the exact determinant sign for every finite stored matrix. `det_exact()` +and `solve_exact()` return values exact for the stored `A` and `b`, subject to +their documented scale and singularity errors. + +Strict `*_exact_f64` conversions succeed only when the exact result already has +an exactly representable finite binary64 value. `RequiresRounding` means a +finite output exists only after rounding; `NotFinite` means even the rounded +result cannot be finite. The explicit rounded conversions use round-to-nearest, +ties-to-even \[9-10\]. A nonzero exact value may consequently round to zero. + +## Tolerances and typed errors + +`Tolerance::try_new` accepts finite values greater than or equal to zero. +`DEFAULT_SINGULAR_TOL` is the absolute value `1e-12`. A tolerance is a rejection +policy for numerical pivots or diagnostics, not an error estimate and not a +condition-number threshold. + +The error model keeps distinct mathematical conclusions separate: + +- `LaError::Singular` distinguishes numerical rejection from exact singularity. + Numerical context retains the factorization kind, observed pivot magnitude, + and tolerance. +- `LaError::Asymmetric` reports the mirrored values that violate LDLT's exact + symmetry precondition. +- `LaError::NotPositiveSemidefinite` records a negative pivot or a zero pivot + with remaining coupling. +- `LaError::NonFinite` distinguishes invalid input locations from arithmetic + operations that overflowed during computation. +- `LaError::Unrepresentable` distinguishes required rounding from the absence + of any finite binary64 output. + +`Matrix::is_symmetric` and `Matrix::first_asymmetry` are tolerance-based, +scale-aware diagnostics. They do not establish the exact mirrored equality +required by `Matrix::ldlt`. + +## Choosing an API + +| Need | API | Important boundary | +|------|-----|--------------------| +| General floating solve | `lu(tol)` then `solve` | Approximate; absolute pivot policy | +| Positive-definite floating solve | `ldlt(tol)` then `solve` | Exact symmetry; computed pivots must exceed tolerance; success is not a certificate | +| Floating determinant, any `D` | `det` | No certified bound; zero is not exact singularity | +| `D ≤ 4` error-bounded determinant/sign test | `det_direct_with_errbound` | Sign is certified only when `|estimate| > bound`; otherwise inconclusive | +| Exact determinant sign | `det_sign_exact` | Exact for stored binary64 entries | +| Exact determinant value or solve | `det_exact`, `solve_exact` | Exact for represented inputs | +| Binary64 output from an exact result | Strict or rounded conversions | Strict conversion forbids rounding | + +## Geometry relationship and scope + +Orientation, in-sphere, and related geometric predicates can be reduced to +determinant signs, which is why an adaptive exact sign is useful near degeneracy +\[8\]. `la-stack` supplies the determinant primitive; callers construct the +problem-specific predicate matrix and remain responsible for any rounding that +occurs during that construction. The crate originated to support +[`delaunay`](https://crates.io/crates/delaunay), but its matrix, factorization, +and exact-arithmetic APIs are general numerical infrastructure. + +The deliberate anti-goals are dynamically sized or rectangular matrices, +sparse storage, broad decomposition coverage, alternate floating scalar +families, and GPU-, parallel-, BLAS-, or LAPACK-scale throughput. Those problems +need different storage models and numerical policies rather than extensions to +this small fixed-dimension design. diff --git a/examples/exact_det_3x3.rs b/examples/exact_det_3x3.rs index 70fb897..a4bfadc 100644 --- a/examples/exact_det_3x3.rs +++ b/examples/exact_det_3x3.rs @@ -3,8 +3,8 @@ //! Exact determinant value for a near-singular 3×3 matrix. //! //! This example demonstrates `det_exact()` and [`ExactF64Conversion`], retaining -//! the provably correct rational determinant before converting that same value -//! to binary64 without repeating exact evaluation. +//! the rational determinant exact for the stored binary64 entries before +//! converting that same value without repeating exact evaluation. //! //! Run with: `cargo run --features exact --example exact_det_3x3` diff --git a/examples/exact_sign_3x3.rs b/examples/exact_sign_3x3.rs index 698e1e7..896f94b 100644 --- a/examples/exact_sign_3x3.rs +++ b/examples/exact_sign_3x3.rs @@ -3,8 +3,8 @@ //! Exact determinant sign for a near-singular 3×3 matrix. //! //! This example demonstrates `det_sign_exact()`, which uses adaptive-precision -//! arithmetic to return the provably correct sign of the determinant — even when -//! the matrix is so close to singular that f64 rounding could give the wrong answer. +//! arithmetic to return the determinant sign exact for the stored binary64 +//! entries, even when f64 evaluation could give the wrong answer. //! //! Run with: `cargo run --features exact --example exact_sign_3x3` diff --git a/examples/exact_solve_3x3.rs b/examples/exact_solve_3x3.rs index 58b806b..7d1d302 100644 --- a/examples/exact_solve_3x3.rs +++ b/examples/exact_solve_3x3.rs @@ -2,11 +2,11 @@ //! Exact linear system solve for a near-singular 3×3 system. //! -//! This example demonstrates `solve_exact()` and [`ExactF64Conversion`]. The exact -//! solve uses arbitrary-precision rational arithmetic to compute a provably -//! correct solution — even when the matrix is so close to singular that the f64 -//! LU solve produces a wildly inaccurate result. `try_to_f64()` only succeeds -//! when every exact component is exactly representable as `f64`, while +//! This example demonstrates `solve_exact()` and [`ExactF64Conversion`]. The +//! solve uses arbitrary-precision rational arithmetic to compute a solution +//! exact for the stored binary64 inputs, even when the f64 LU solve is wildly +//! inaccurate. `try_to_f64()` only succeeds when every exact component is +//! exactly representable as `f64`, while //! `to_rounded_f64()` explicitly opts into nearest-even rounding. //! //! Run with: `cargo run --features exact --example exact_solve_3x3` diff --git a/justfile b/justfile index ebe763b..7793571 100644 --- a/justfile +++ b/justfile @@ -17,7 +17,7 @@ cargo_llvm_cov_version := "0.8.7" dprint_version := "0.55.1" git_cliff_version := "2.13.1" just_version := "1.56.0" -rumdl_version := "0.2.30" +rumdl_version := "0.2.31" taplo_version := "0.10.0" typos_version := "1.48.0" uv_version := "0.11.28" diff --git a/scripts/bench_compare.py b/scripts/bench_compare.py index 8c8b5b9..44975a5 100644 --- a/scripts/bench_compare.py +++ b/scripts/bench_compare.py @@ -70,6 +70,8 @@ _EXACT_DIMENSION_BENCHES_WITH_DIRECT: list[str] = [ "det", "det_direct", + "det_direct_with_errbound", + "det_errbound", *_EXACT_DIMENSION_BENCHES[1:], ] @@ -140,6 +142,9 @@ _V0_4_3_API_COMPATIBILITY = "la_stack_v0_4_3_api" _V0_4_3_UNAVAILABLE_BASELINE_ROWS: frozenset[tuple[str, str]] = frozenset( { + ("exact_d2", "det_direct_with_errbound"), + ("exact_d3", "det_direct_with_errbound"), + ("exact_d4", "det_direct_with_errbound"), ("d8", "la_stack_det_from_lu_balanced_range"), ("d8", "la_stack_det_from_ldlt_balanced_range"), } @@ -772,14 +777,18 @@ def _collect_exact_comparisons( criterion_dir: Path, baseline_name: str, stat: str, - scope: str, + policy: ComparisonPolicy, ) -> ComparisonCollection: """Compare exact results while retaining every missing expected row.""" comparisons: list[Comparison] = [] gaps: list[CoverageGap] = [] + unavailable_baseline_rows = _UNAVAILABLE_BASELINE_ROWS_BY_COMPATIBILITY.get( + policy.baseline_api_compatibility or "", + frozenset(), + ) for group, benches in EXACT_GROUPS.items(): - if scope == "release-signal" and group not in EXACT_RELEASE_SIGNAL_GROUPS: + if policy.scope == "release-signal" and group not in EXACT_RELEASE_SIGNAL_GROUPS: continue group_dir = criterion_dir / group @@ -788,6 +797,23 @@ def _collect_exact_comparisons( baseline_bench, base_path = _exact_baseline_path(group_dir, bench, baseline_name) missing_current = not new_path.exists() missing_baseline = not base_path.exists() + baseline_unavailable = (group, bench) in unavailable_baseline_rows + + if baseline_unavailable: + if missing_current: + gaps.append( + CoverageGap( + suite="exact", + group=group, + bench=bench, + baseline_bench=baseline_bench, + missing_current=True, + missing_baseline=False, + ) + ) + else: + _read_estimate(new_path, stat) + continue if missing_current or missing_baseline: gaps.append( @@ -817,7 +843,7 @@ def _collect_exact_comparisons( ) ) - expected_groups = [group for group in EXACT_GROUPS if scope != "release-signal" or group in EXACT_RELEASE_SIGNAL_GROUPS] + expected_groups = [group for group in EXACT_GROUPS if policy.scope != "release-signal" or group in EXACT_RELEASE_SIGNAL_GROUPS] if not any((criterion_dir / group).is_dir() for group in expected_groups): gaps.append(_entire_suite_gap("exact")) @@ -1047,7 +1073,7 @@ def _collect_comparisons( comparisons: list[Comparison] = [] gaps: list[CoverageGap] = [] if suite in ("all", "exact"): - exact = _collect_exact_comparisons(criterion_dir, baseline_name, stat, policy.scope) + exact = _collect_exact_comparisons(criterion_dir, baseline_name, stat, policy) comparisons.extend(exact.comparisons) gaps.extend(exact.gaps) if suite in ("all", "vs_linalg"): @@ -1501,6 +1527,12 @@ def _provenance_markdown(provenance: HarnessProvenance) -> list[str]: "`d8/la_stack_det_from_ldlt_balanced_range` were not timed because v0.4.3 returns zero for a fixture " "whose exact determinant is one; current samples remain required, but no speedup is claimed." ) + if compatibility == _V0_4_3_API_COMPATIBILITY and criterion.suite in {"all", "exact"}: + lines.append( + "- Baseline-unavailable rows: `exact_d2/det_direct_with_errbound`, " + "`exact_d3/det_direct_with_errbound`, and `exact_d4/det_direct_with_errbound` were not timed because " + "v0.4.3 predates the paired API; the comparable `det_errbound` baselines remain required." + ) return lines diff --git a/scripts/tests/test_bench_compare.py b/scripts/tests/test_bench_compare.py index 55264f8..e5a56f1 100644 --- a/scripts/tests/test_bench_compare.py +++ b/scripts/tests/test_bench_compare.py @@ -236,10 +236,13 @@ def test_unknown_passthrough(self) -> None: assert bench_compare._group_heading("something_else") == "something_else" -def test_exact_registry_only_tracks_supported_direct_determinants() -> None: +def test_exact_registry_only_tracks_supported_direct_determinant_filters() -> None: + filter_benches = ["det_direct", "det_direct_with_errbound", "det_errbound"] for dimension in (2, 3, 4): - assert "det_direct" in bench_compare.EXACT_GROUPS[f"exact_d{dimension}"] - assert "det_direct" not in bench_compare.EXACT_GROUPS["exact_d5"] + benches = bench_compare.EXACT_GROUPS[f"exact_d{dimension}"] + start = benches.index("det_direct") + assert benches[start : start + len(filter_benches)] == filter_benches + assert not set(filter_benches) & set(bench_compare.EXACT_GROUPS["exact_d5"]) # --------------------------------------------------------------------------- @@ -447,12 +450,34 @@ def test_collect_comparisons_records_each_missing_side_in_registry_order(tmp_pat collection = bench_compare._collect_comparisons(tmp_path, "last", "median", suite="exact") first_three = collection.gaps[:3] - assert [gap.bench for gap in first_three] == ["det", "det_direct", "det_exact"] + assert [gap.bench for gap in first_three] == ["det", "det_direct", "det_direct_with_errbound"] assert (first_three[0].missing_current, first_three[0].missing_baseline) == (False, True) assert (first_three[1].missing_current, first_three[1].missing_baseline) == (True, False) assert (first_three[2].missing_current, first_three[2].missing_baseline) == (True, True) +def test_v043_adapter_requires_bound_only_baselines_but_allows_missing_paired_rows(tmp_path: Path) -> None: + unavailable = bench_compare._V0_4_3_UNAVAILABLE_BASELINE_ROWS + for group, benches in bench_compare.EXACT_GROUPS.items(): + for bench in benches: + _write_estimates(tmp_path / group / bench / "new" / "estimates.json", "median", 10.0) + if (group, bench) not in unavailable: + _write_estimates(tmp_path / group / bench / "v0.4.3" / "estimates.json", "median", 20.0) + + collection = bench_compare._collect_comparisons( + tmp_path, + "v0.4.3", + "median", + suite="exact", + policy=bench_compare.ComparisonPolicy(baseline_api_compatibility="la_stack_v0_4_3_api"), + ) + + assert collection.gaps == [] + comparison_rows = {(comparison.group, comparison.bench) for comparison in collection.comparisons} + assert all((f"exact_d{dimension}", "det_errbound") in comparison_rows for dimension in (2, 3, 4)) + assert all((f"exact_d{dimension}", "det_direct_with_errbound") not in comparison_rows for dimension in (2, 3, 4)) + + def test_collect_comparisons_reports_wholly_absent_selected_suite(tmp_path: Path) -> None: group = tmp_path / "exact_d2" _write_estimates(group / "det" / "new" / "estimates.json", "median", 10.0) @@ -762,6 +787,10 @@ def test_read_schema2_provenance_records_versions_dirty_source_and_both_gates(tm assert "d8/la_stack_det_from_lu_balanced_range" in rendered assert "d8/la_stack_det_from_ldlt_balanced_range" in rendered assert "exact determinant is one" in rendered + assert "exact_d2/det_direct_with_errbound" in rendered + assert "exact_d3/det_direct_with_errbound" in rendered + assert "exact_d4/det_direct_with_errbound" in rendered + assert "the comparable `det_errbound` baselines remain required" in rendered assert bench_compare._comparison_policy("release-signal", provenance) == bench_compare.ComparisonPolicy( scope="release-signal", baseline_api_compatibility="la_stack_v0_4_3_api", diff --git a/src/error.rs b/src/error.rs index 90e6144..f0d0e5d 100644 --- a/src/error.rs +++ b/src/error.rs @@ -172,7 +172,11 @@ pub enum NonFiniteOrigin { }, } -/// Reason a symmetric matrix is outside the positive-semidefinite LDLT domain. +/// Computed LDLT condition that violates the no-pivot positive-semidefinite +/// factorization requirements. +/// +/// These values are computed in binary64. They explain why LDLT rejected the +/// matrix, but do not prove that the stored matrix is exactly indefinite. /// /// # Examples /// ``` @@ -359,7 +363,11 @@ pub enum LaError { /// Maximum absolute difference allowed by the symmetry check. allowed_abs_diff: f64, }, - /// A symmetric matrix is outside the positive-semidefinite LDLT domain. + /// A computed LDLT pivot violates the no-pivot positive-semidefinite + /// factorization requirements. + /// + /// This diagnoses a binary64 factorization rejection; it is not an exact + /// certificate that the stored matrix is indefinite. #[non_exhaustive] NotPositiveSemidefinite { /// LDLT pivot column or step where the violation was detected. @@ -570,8 +578,8 @@ impl LaError { } } - /// Construct a [`LaError::NotPositiveSemidefinite`] error for a negative - /// LDLT diagonal pivot. + /// Construct a [`LaError::NotPositiveSemidefinite`] error for a computed + /// negative LDLT diagonal pivot. #[inline] #[must_use] pub const fn not_positive_semidefinite_negative(pivot_col: usize, value: f64) -> Self { @@ -581,9 +589,9 @@ impl LaError { } } - /// Construct a [`LaError::NotPositiveSemidefinite`] error for a zero LDLT - /// diagonal with a non-zero coupling at `row`, distinguishing the observed - /// positive-semidefinite violation from an uncoupled singular pivot. + /// Construct a [`LaError::NotPositiveSemidefinite`] error for a computed + /// zero LDLT diagonal with a non-zero coupling at `row`, distinguishing the + /// observed factorization violation from an uncoupled singular pivot. #[inline] #[must_use] pub const fn not_positive_semidefinite_zero_coupling( @@ -718,14 +726,14 @@ impl fmt::Display for LaError { violation: PositiveSemidefiniteViolation::NegativePivot { value }, } => write!( f, - "matrix is not positive semidefinite at LDLT pivot column {pivot_col}: diagonal value {value} < 0" + "LDLT rejected the matrix at pivot column {pivot_col}: computed diagonal value {value} < 0" ), Self::NotPositiveSemidefinite { pivot_col, violation: PositiveSemidefiniteViolation::ZeroPivotCoupling { row, value }, } => write!( f, - "matrix is not positive semidefinite at LDLT pivot column {pivot_col}: zero diagonal has non-zero coupling at row {row} with value {value}" + "LDLT rejected the matrix at pivot column {pivot_col}: computed zero diagonal has non-zero coupling at row {row} with value {value}" ), } } @@ -925,11 +933,11 @@ mod tests { fn positive_semidefinite_errors_preserve_distinct_violations() { assert_eq!( LaError::not_positive_semidefinite_negative(1, -3.0).to_string(), - "matrix is not positive semidefinite at LDLT pivot column 1: diagonal value -3 < 0" + "LDLT rejected the matrix at pivot column 1: computed diagonal value -3 < 0" ); assert_eq!( LaError::not_positive_semidefinite_zero_coupling(0, 1, 2.0).to_string(), - "matrix is not positive semidefinite at LDLT pivot column 0: zero diagonal has non-zero coupling at row 1 with value 2" + "LDLT rejected the matrix at pivot column 0: computed zero diagonal has non-zero coupling at row 1 with value 2" ); } diff --git a/src/exact.rs b/src/exact.rs index ae21675..2f0a5ff 100644 --- a/src/exact.rs +++ b/src/exact.rs @@ -3,6 +3,9 @@ //! Exact arithmetic operations via arbitrary-precision rational numbers. //! //! This module is only compiled when the `"exact"` Cargo feature is enabled. +//! Exactness begins with the finite binary64 values already stored in +//! [`Matrix`] and [`Vector`]: each value is lifted losslessly to a rational. +//! These APIs cannot recover information rounded away before construction. //! //! # Architecture //! @@ -14,7 +17,7 @@ //! `decompose_proven_finite_f64` into `mantissa × 2^exponent`, then all entries //! are scaled to a common `BigInt` //! matrix (shifting by `e - e_min`). D≤4 uses direct integer expansions; larger -//! matrices use fraction-free Bareiss elimination entirely in `BigInt` +//! matrices use fraction-free Bareiss elimination \[7\] entirely in `BigInt` //! arithmetic — no `BigRational`, no GCD, no denominator tracking. The result //! is `(det_int, total_exp)` where `det = det_int × 2^(D × e_min)`. `det_exact` //! wraps this with `big_int_exp_to_big_rational` to reconstruct a reduced @@ -24,7 +27,7 @@ //! directly from `det_int` (the scale factor is always positive). //! //! `det_sign_exact` adds a two-stage adaptive-precision optimisation inspired -//! by Shewchuk's robust geometric predicates: +//! by Shewchuk's robust geometric predicates \[8\]: //! //! 1. **Fast filter (D ≤ 4)**: compute `det_direct()` and a conservative error //! bound. If `|det| > bound`, the f64 sign is provably correct — return @@ -45,14 +48,14 @@ //! the exact power-of-two ratio between those scales. This avoids inflating one //! side's integers merely because the other side has much smaller entries. //! Forward elimination runs entirely in `BigInt` with -//! fraction-free Bareiss updates — no `BigRational`, no GCD +//! fraction-free Bareiss updates \[7\] — no `BigRational`, no GCD //! normalisation in the `O(D³)` phase. Once the system is upper //! triangular, back-substitution is performed in `BigRational`, where //! fractions are inherent; this phase is only `O(D²)` so the rational //! overhead is modest. First-non-zero pivoting is used throughout; //! since all arithmetic is exact, any non-zero pivot gives the correct -//! result (no numerical stability concern). Every finite `f64` is -//! exactly representable as a rational, so the result is provably correct. +//! result (no numerical stability concern). Every finite `f64` is exactly +//! representable as a rational, so the result is exact for the stored inputs. //! `solve_exact_f64` returns `Vector` only when every exact component is //! exactly representable as finite binary64; `solve_exact_rounded_f64` returns //! the exact components rounded to finite binary64. @@ -68,7 +71,7 @@ //! path builds a `BigInt` augmented system and lifts the //! upper-triangular result into `BigRational` for back-substitution. //! See Goldberg \[10\] for background on floating-point representation -//! and exact rational reconstruction. Reference numbers refer to +//! and conversion. Reference numbers refer to //! `REFERENCES.md`. //! //! ## Validation @@ -1368,13 +1371,15 @@ impl Matrix { /// /// Returns the determinant as an exact [`BigRational`] value. Every finite /// `f64` is exactly representable as a rational, so the conversion is - /// lossless and the result is provably correct. + /// lossless and the result is exact for the stored binary64 entries. It + /// cannot recover precision lost before matrix construction. /// /// # When to use /// /// Use this when you need the exact determinant *value* — for example, - /// exact volume computation or distinguishing truly-degenerate simplices - /// from near-degenerate ones. If you only need the *sign*, prefer + /// volume computation over stored coordinates or distinguishing simplices + /// that are exactly degenerate at those coordinates from near-degenerate + /// ones. If you only need the *sign*, prefer /// [`det_sign_exact`](Self::det_sign_exact) which has a fast f64 filter. /// /// # Examples @@ -1485,15 +1490,17 @@ impl Matrix { /// Requires the `exact` Cargo feature. /// /// Solves `A x = b` where `A` is `self` and `b` is the given vector. - /// Returns the exact solution as `[BigRational; D]`. Every finite `f64` - /// is exactly representable as a rational, so the conversion is lossless - /// and the result is provably correct. + /// Returns the exact solution as `[BigRational; D]`. Every finite `f64` is + /// exactly representable as a rational, so the conversion is lossless and + /// the result is exact for the stored binary64 entries. It cannot recover + /// precision lost before matrix or vector construction. /// /// # When to use /// - /// Use this when you need a provably correct solution — for example, - /// exact circumcenter computation for near-degenerate simplices where - /// f64 arithmetic may produce wildly wrong results. + /// Use this when you need a solution exact for the stored inputs — for + /// example, circumcenter computation over stored coordinates for + /// near-degenerate simplices where f64 arithmetic may produce wildly wrong + /// results. /// /// # Algorithm /// @@ -1615,7 +1622,9 @@ impl Matrix { /// Requires the `exact` Cargo feature. /// /// Returns [`DeterminantSign::Positive`], [`DeterminantSign::Negative`], or - /// [`DeterminantSign::Zero`] according to the exact determinant. + /// [`DeterminantSign::Zero`] according to the determinant that is exact for + /// the stored binary64 entries. This cannot recover precision lost before + /// matrix construction. /// /// For D ≤ 4, a fast f64 filter is tried first: `det_direct()` is compared /// against a conservative error bound derived from the matrix permanent. @@ -1626,10 +1635,11 @@ impl Matrix { /// /// # When to use /// - /// Use this when the sign of the determinant must be correct regardless of - /// floating-point conditioning (e.g. geometric predicates on near-degenerate - /// configurations). For well-conditioned matrices the fast filter resolves - /// the sign without touching `BigRational`, so the overhead is minimal. + /// Use this when the sign of the determinant over the stored entries must be + /// correct regardless of floating-point conditioning (e.g. geometric + /// predicates on near-degenerate stored coordinates). For well-conditioned + /// matrices the fast filter resolves the sign without touching + /// `BigRational`, so the overhead is minimal. /// /// # Examples /// ``` diff --git a/src/ldlt.rs b/src/ldlt.rs index 414140f..f4c19cd 100644 --- a/src/ldlt.rs +++ b/src/ldlt.rs @@ -2,9 +2,12 @@ //! LDLT factorization and solves. //! -//! This module provides a stack-allocated LDLT factorization (`A = L D Lᵀ`) intended for -//! symmetric positive definite (SPD) and positive semi-definite (PSD) matrices (e.g. Gram -//! matrices) without pivoting. +//! This module provides a stack-allocated LDLT factorization (`A = L D Lᵀ`) +//! without pivoting. Successful factors require an exactly symmetric input and +//! every computed diagonal pivot to be positive and above the caller's +//! tolerance. Computed zero and tolerance-small positive pivots are diagnosed +//! rather than returned in a usable factor. See `REFERENCES.md` \[4-6, 11-12\] for +//! Cholesky/LDLT background and pivoted symmetric-indefinite alternatives. //! //! # Preconditions //! The input matrix must be **symmetric**. This is a correctness contract, not a hint: @@ -16,17 +19,22 @@ use core::hint::cold_path; +use crate::matrix::SymmetricMatrix; use crate::scaled_product::{RangeCheckedProduct, ScaledProduct, range_checked_product}; use crate::vector::Vector; -use crate::{ArithmeticOperation, FactorizationKind, LaError, SymmetricMatrix, Tolerance}; +use crate::{ArithmeticOperation, FactorizationKind, LaError, Tolerance}; -/// LDLT factorization (`A = L D Lᵀ`) for symmetric positive (semi)definite matrices. +/// LDLT factorization (`A = L D Lᵀ`) for exactly symmetric positive-definite matrices. /// /// `Ldlt<0>` represents the empty factorization. Its determinant is the empty /// product `1.0`, and solving against [`Vector<0>`] returns [`Vector<0>`]. /// /// This factorization is **not** a general-purpose symmetric-indefinite LDLT (no pivoting). -/// It assumes the input matrix is symmetric and (numerically) SPD/PSD. +/// It assumes the input matrix is exactly symmetric and numerically positive +/// definite under the caller's absolute pivot tolerance. Computed zero and +/// tolerance-small positive pivots return [`LaError::Singular`]. Because pivots +/// are computed in binary64, success is not an exact proof that the stored +/// matrix is positive definite. /// /// # Preconditions /// The source matrix passed to [`Matrix::ldlt`](crate::Matrix::ldlt) must be @@ -235,7 +243,8 @@ impl Ldlt { /// Determinant of the original matrix. /// - /// For SPD/PSD matrices, this is the product of the diagonal terms of `D`. + /// For a successfully constructed factorization, this is the product of + /// the diagonal terms of `D`. /// /// # Examples /// ``` diff --git a/src/lib.rs b/src/lib.rs index b5dbe11..ba5ab44 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -371,24 +371,6 @@ pub use matrix::{DeterminantWithErrorBound, Matrix}; pub use tolerance::{DEFAULT_SINGULAR_TOL, Tolerance}; pub use vector::Vector; -/// A finite [`Matrix`] proven exactly symmetric for LDLT factorization. -/// -/// Mirrored entries have equal numeric values; IEEE-754 signed zeros may have -/// different bit patterns because `+0.0 == -0.0`. -#[must_use] -#[derive(Clone, Copy, Debug, PartialEq)] -struct SymmetricMatrix { - matrix: Matrix, -} - -impl SymmetricMatrix { - /// Consume the wrapper and return the underlying matrix. - #[inline] - const fn into_matrix(self) -> Matrix { - self.matrix - } -} - /// Fallibly dispatch a runtime dimension to a concrete stack-allocated matrix. /// /// The macro creates a zero matrix with type `Matrix` for the selected diff --git a/src/lu.rs b/src/lu.rs index f18f3dd..8e799a1 100644 --- a/src/lu.rs +++ b/src/lu.rs @@ -1,6 +1,11 @@ #![forbid(unsafe_code)] //! LU decomposition and solves. +//! +//! The implementation computes `P A = L U` with partial pivoting. Partial +//! pivoting is a practical finite-precision strategy rather than an +//! unconditional accuracy guarantee; see `REFERENCES.md` \[1-3, 11-12\] for +//! stability analysis and standard algorithmic background. use core::hint::cold_path; @@ -13,6 +18,8 @@ use crate::{ArithmeticOperation, FactorizationKind, LaError, Tolerance}; /// /// `Lu<0>` represents the empty factorization. Its determinant is the empty /// product `1.0`, and solving against [`Vector<0>`] returns [`Vector<0>`]. +/// Numerical solves and determinants remain subject to binary64 rounding and +/// matrix conditioning; this type does not provide a certified error bound. #[must_use] #[derive(Clone, Copy, Debug, PartialEq)] pub struct Lu { diff --git a/src/matrix.rs b/src/matrix.rs index c775b80..a6665ac 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -6,18 +6,16 @@ use core::hint::cold_path; use crate::ldlt::Ldlt; use crate::lu::Lu; -use crate::{ - ArithmeticOperation, ERR_COEFF_2, ERR_COEFF_3, ERR_COEFF_4, LaError, SymmetricMatrix, Tolerance, -}; +use crate::{ArithmeticOperation, ERR_COEFF_2, ERR_COEFF_3, ERR_COEFF_4, LaError, Tolerance}; /// A closed-form determinant and its certified absolute error bound. /// /// Values of this type are produced by /// [`Matrix::det_direct_with_errbound`]. The paired result guarantees that the -/// determinant and bound came from one traversal of the same rounded -/// arithmetic tree. The guarantee is unavailable when gradual underflow could -/// invalidate the relative-error analysis or when the matrix dimension exceeds -/// the closed-form D ≤ 4 scope. +/// determinant was evaluated once and that its matching bound was computed for +/// the same matrix in one call. The guarantee is unavailable when gradual +/// underflow could invalidate the relative-error analysis or when the matrix +/// dimension exceeds the closed-form D ≤ 4 scope. #[must_use] #[non_exhaustive] #[derive(Clone, Copy, Debug, PartialEq)] @@ -77,6 +75,16 @@ pub struct Matrix { rows: [[f64; D]; D], } +/// A finite [`Matrix`] proven exactly symmetric for LDLT factorization. +/// +/// Mirrored entries have equal numeric values; IEEE-754 signed zeros may have +/// different bit patterns because `+0.0 == -0.0`. +#[must_use] +#[derive(Clone, Copy, Debug, PartialEq)] +pub(crate) struct SymmetricMatrix { + matrix: Matrix, +} + /// Rounded arithmetic result together with proof that gradual underflow could /// not have changed that operation's result. /// @@ -202,6 +210,12 @@ impl<'a, const D: usize> Det4SharedMinorInput<'a, D> { } impl SymmetricMatrix { + /// Consume the wrapper and return the underlying matrix. + #[inline] + pub(crate) const fn into_matrix(self) -> Matrix { + self.matrix + } + /// Construct a symmetric matrix proof without checking the invariant. /// /// This constructor is only for paths that have already validated exact @@ -677,6 +691,8 @@ impl Matrix { /// `D = 0` follows the empty-matrix convention: factorization succeeds, /// [`Lu::det`](crate::Lu::det) returns `1.0`, and solving a length-zero /// right-hand side returns a length-zero [`Vector`](crate::Vector). + /// Partial pivoting is a practical finite-precision strategy, not a + /// certified accuracy guarantee; see `REFERENCES.md` \[1-3, 11-12\]. /// /// # Examples /// ``` @@ -731,8 +747,13 @@ impl Matrix { /// [`Ldlt::det`](crate::Ldlt::det) returns `1.0`, and solving a length-zero /// right-hand side returns a length-zero [`Vector`](crate::Vector). /// - /// This is intended for symmetric positive definite (SPD) and positive semi-definite (PSD) - /// matrices such as Gram matrices. + /// This is intended for exactly symmetric positive-definite matrices such + /// as nonsingular Gram matrices. Computed zero and tolerance-small positive + /// pivots are diagnosed as singular rather than returned in a usable + /// factorization. Because pivots are computed in binary64, success is not + /// an exact proof that the stored matrix is positive definite. + /// See `REFERENCES.md` \[4-6, 11-12\] for Cholesky/LDLT background and the + /// pivoted symmetric-indefinite alternative. /// /// # Symmetry validation /// The input matrix `self` must be exactly symmetric: every mirrored pair @@ -1143,8 +1164,8 @@ impl Matrix { /// closed-form bound is available. /// /// This is the preferred API when both values are needed: it evaluates the - /// determinant arithmetic tree once, so the approximation and bound cannot - /// accidentally come from separate traversals. + /// determinant arithmetic tree once, then computes the matching bound for + /// the same matrix within that call. /// /// # Examples /// ``` diff --git a/tests/exact_bench_config.rs b/tests/exact_bench_config.rs index 7d7498e..d20a706 100644 --- a/tests/exact_bench_config.rs +++ b/tests/exact_bench_config.rs @@ -40,6 +40,7 @@ fn validate_baseline_and_random_corpus() { /// Report whether `det_sign_exact` can certify this fixture through its direct filter. /// /// Both a missing bound and a direct-path error select the exact fallback. +#[cfg(not(la_stack_v0_4_3_api))] fn fast_filter_is_conclusive(input: ExactInput) -> bool { match input.matrix.det_direct_with_errbound() { Ok(Some(estimate)) => estimate.determinant().abs() > estimate.absolute_error_bound(), @@ -47,6 +48,15 @@ fn fast_filter_is_conclusive(input: ExactInput) -> bool { } } +/// Use the separate v0.4.3 determinant and bound APIs for the same filter test. +#[cfg(la_stack_v0_4_3_api)] +fn fast_filter_is_conclusive(input: ExactInput) -> bool { + match (input.matrix.det_direct(), input.matrix.det_errbound()) { + (Ok(Some(determinant)), Ok(Some(bound))) => determinant.abs() > bound, + (Ok(None) | Err(_), _) | (_, Ok(None) | Err(_)) => false, + } +} + #[test] fn i16_range_rejects_unordered_bounds() { let Err(err) = I16Range::try_new(5, 4) else { diff --git a/tests/vs_linalg_inputs.rs b/tests/vs_linalg_inputs.rs index d245bd8..f484583 100644 --- a/tests/vs_linalg_inputs.rs +++ b/tests/vs_linalg_inputs.rs @@ -10,6 +10,8 @@ use faer::perm::PermRef; use faer::{Mat, Side}; use nalgebra::{Const, DimMin, SMatrix, SVector}; +#[cfg(feature = "exact")] +use la_stack::ExactF64Conversion; use la_stack::{DEFAULT_SINGULAR_TOL, Matrix, Vector}; #[path = "../benches/common/vs_linalg.rs"] @@ -283,11 +285,24 @@ fn stress_inputs_exercise_pivoting_conditioning_and_scaled_products() { let pivoting_rows = make_pivoting_matrix_rows::<8>(); assert!(pivoting_rows[1][0].abs() > pivoting_rows[0][0].abs()); - let pivoting_lu = Matrix::<8>::try_from_rows(pivoting_rows) - .unwrap() - .lu(zero_tolerance) - .unwrap(); - assert!(pivoting_lu.det().unwrap().is_finite()); + let pivoting = Matrix::<8>::try_from_rows(pivoting_rows).unwrap(); + let pivoting_lu = pivoting.lu(zero_tolerance).unwrap(); + let pivoting_det = pivoting_lu.det().unwrap(); + assert!(pivoting_det.is_finite()); + + #[cfg(feature = "exact")] + { + // The fixture is the shared baseline matrix with one row swap, so its + // exact determinant is the negation of the baseline determinant. This + // Bareiss-backed oracle is independent of the timed floating LU path. + let baseline = Matrix::<8>::try_from_rows(make_matrix_rows::<8>()).unwrap(); + let expected = -baseline.det_exact().unwrap(); + assert_close( + "pivoting LU determinant against exact row-swap oracle", + pivoting_det, + expected.to_rounded_f64().unwrap(), + ); + } let ill_conditioned = Matrix::<8>::try_from_rows(make_ill_conditioned_matrix_rows()).unwrap(); let expected_ill_conditioned_det = f64::from_bits((1023_u64 - 448) << 52); diff --git a/typos.toml b/typos.toml index 6c98b70..4d0d833 100644 --- a/typos.toml +++ b/typos.toml @@ -31,7 +31,7 @@ Schreiber = "Schreiber" Tikhomirov = "Tikhomirov" Trefethen = "Trefethen" # Foreign-language terms in reference titles -methode = "methode" +Sur = "Sur" # Math / crate terms faer = "faer" frhs = "frhs" From e2573434c541a8115750b3a965882f7b484b6be0 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Sun, 12 Jul 2026 10:01:06 -0700 Subject: [PATCH 2/2] docs(science): clarify determinant and LDLT failure contracts - State the determinant sign condition without ambiguous Markdown delimiters. - Distinguish singular zero pivots from coupled zero-pivot LDLT rejections. - Align citation and contributor guidance with the positive-definite domain. --- AGENTS.md | 3 ++- CITATION.cff | 3 ++- docs/mathematical_basis.md | 2 +- scripts/tests/test_bench_compare.py | 39 ++++++++++++++++++++++++++--- src/ldlt.rs | 10 +++++--- 5 files changed, 46 insertions(+), 11 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index f700bd7..74f8399 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -430,7 +430,8 @@ When creating or updating issues: - `src/vector.rs`: `Vector` (`[f64; D]`) - `src/matrix.rs`: `Matrix` (`[[f64; D]; D]`) + helpers (`get`, `try_get`, `set`, `inf_norm`, `det`, `det_direct`) - `src/lu.rs`: `Lu` factorization with partial pivoting (`solve`, `det`) - - `src/ldlt.rs`: `Ldlt` factorization without pivoting for symmetric SPD/PSD matrices (`solve`, `det`) + - `src/ldlt.rs`: `Ldlt` factorization without pivoting for exactly + symmetric positive-definite matrices (`solve`, `det`) - `src/exact.rs`: exact arithmetic behind `features = ["exact"]`: - Determinants: `det_exact()`, strict `det_exact_f64()`, rounded `det_exact_rounded_f64()`, and `det_sign_exact()` via a scaled `BigInt` diff --git a/CITATION.cff b/CITATION.cff index 9fa1f58..da83df8 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -23,6 +23,7 @@ keywords: - "robust-predicates" abstract: >- la-stack is a Rust library providing fast, stack-allocated linear algebra primitives - and factorization routines (LU with partial pivoting and LDLT for symmetric SPD/PSD matrices) + and factorization routines (LU with partial pivoting and no-pivot LDLT for exactly + symmetric positive-definite matrices) for fixed small dimensions using const generics. license: BSD-3-Clause diff --git a/docs/mathematical_basis.md b/docs/mathematical_basis.md index da24132..4e77fe3 100644 --- a/docs/mathematical_basis.md +++ b/docs/mathematical_basis.md @@ -198,7 +198,7 @@ required by `Matrix::ldlt`. | General floating solve | `lu(tol)` then `solve` | Approximate; absolute pivot policy | | Positive-definite floating solve | `ldlt(tol)` then `solve` | Exact symmetry; computed pivots must exceed tolerance; success is not a certificate | | Floating determinant, any `D` | `det` | No certified bound; zero is not exact singularity | -| `D ≤ 4` error-bounded determinant/sign test | `det_direct_with_errbound` | Sign is certified only when `|estimate| > bound`; otherwise inconclusive | +| `D ≤ 4` error-bounded determinant/sign test | `det_direct_with_errbound` | Sign is certified when estimate magnitude exceeds bound; otherwise inconclusive | | Exact determinant sign | `det_sign_exact` | Exact for stored binary64 entries | | Exact determinant value or solve | `det_exact`, `solve_exact` | Exact for represented inputs | | Binary64 output from an exact result | Strict or rounded conversions | Strict conversion forbids rounding | diff --git a/scripts/tests/test_bench_compare.py b/scripts/tests/test_bench_compare.py index e5a56f1..498eb84 100644 --- a/scripts/tests/test_bench_compare.py +++ b/scripts/tests/test_bench_compare.py @@ -456,13 +456,23 @@ def test_collect_comparisons_records_each_missing_side_in_registry_order(tmp_pat assert (first_three[2].missing_current, first_three[2].missing_baseline) == (True, True) -def test_v043_adapter_requires_bound_only_baselines_but_allows_missing_paired_rows(tmp_path: Path) -> None: +def _write_v043_exact_comparison_samples( + criterion_dir: Path, + *, + missing_current: frozenset[tuple[str, str]] = frozenset(), +) -> None: unavailable = bench_compare._V0_4_3_UNAVAILABLE_BASELINE_ROWS for group, benches in bench_compare.EXACT_GROUPS.items(): for bench in benches: - _write_estimates(tmp_path / group / bench / "new" / "estimates.json", "median", 10.0) - if (group, bench) not in unavailable: - _write_estimates(tmp_path / group / bench / "v0.4.3" / "estimates.json", "median", 20.0) + row = (group, bench) + if row not in missing_current: + _write_estimates(criterion_dir / group / bench / "new" / "estimates.json", "median", 10.0) + if row not in unavailable: + _write_estimates(criterion_dir / group / bench / "v0.4.3" / "estimates.json", "median", 20.0) + + +def test_v043_adapter_requires_bound_only_baselines_but_allows_missing_paired_rows(tmp_path: Path) -> None: + _write_v043_exact_comparison_samples(tmp_path) collection = bench_compare._collect_comparisons( tmp_path, @@ -478,6 +488,27 @@ def test_v043_adapter_requires_bound_only_baselines_but_allows_missing_paired_ro assert all((f"exact_d{dimension}", "det_direct_with_errbound") not in comparison_rows for dimension in (2, 3, 4)) +def test_v043_adapter_reports_missing_current_for_unavailable_paired_row(tmp_path: Path) -> None: + target = ("exact_d2", "det_direct_with_errbound") + _write_v043_exact_comparison_samples(tmp_path, missing_current=frozenset({target})) + + collection = bench_compare._collect_comparisons( + tmp_path, + "v0.4.3", + "median", + suite="exact", + policy=bench_compare.ComparisonPolicy(baseline_api_compatibility="la_stack_v0_4_3_api"), + ) + + assert len(collection.gaps) == 1 + gap = collection.gaps[0] + assert gap.suite == "exact" + assert gap.group == target[0] + assert gap.bench == target[1] + assert gap.missing_current + assert not gap.missing_baseline + + def test_collect_comparisons_reports_wholly_absent_selected_suite(tmp_path: Path) -> None: group = tmp_path / "exact_d2" _write_estimates(group / "det" / "new" / "estimates.json", "median", 10.0) diff --git a/src/ldlt.rs b/src/ldlt.rs index f4c19cd..3fb9f04 100644 --- a/src/ldlt.rs +++ b/src/ldlt.rs @@ -31,10 +31,12 @@ use crate::{ArithmeticOperation, FactorizationKind, LaError, Tolerance}; /// /// This factorization is **not** a general-purpose symmetric-indefinite LDLT (no pivoting). /// It assumes the input matrix is exactly symmetric and numerically positive -/// definite under the caller's absolute pivot tolerance. Computed zero and -/// tolerance-small positive pivots return [`LaError::Singular`]. Because pivots -/// are computed in binary64, success is not an exact proof that the stored -/// matrix is positive definite. +/// definite under the caller's absolute pivot tolerance. An uncoupled computed +/// zero or a tolerance-small positive pivot returns [`LaError::Singular`]; a +/// computed zero with non-zero remaining coupling returns +/// [`LaError::NotPositiveSemidefinite`]. Because pivots are computed in +/// binary64, success is not an exact proof that the stored matrix is positive +/// definite. /// /// # Preconditions /// The source matrix passed to [`Matrix::ldlt`](crate::Matrix::ldlt) must be