diff --git a/README.md b/README.md index a5c5a67..d7be54f 100644 --- a/README.md +++ b/README.md @@ -108,6 +108,17 @@ let det = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap().det(); assert!((det - 1.0).abs() <= 1e-12); ``` +> ⚠️ **LDLT precondition:** The input matrix must be **symmetric**. Symmetry is +> verified by a `debug_assert!` in debug builds only; release builds silently accept +> asymmetric inputs and produce a meaningless factorization. Pre-validate with +> [`Matrix::is_symmetric`](https://docs.rs/la-stack/latest/la_stack/struct.Matrix.html#method.is_symmetric) +> (or locate the offending pair with +> [`Matrix::first_asymmetry`](https://docs.rs/la-stack/latest/la_stack/struct.Matrix.html#method.first_asymmetry)) +> when you cannot statically guarantee symmetry, or fall back to `lu()` if your +> matrices may not be symmetric at all. See +> [`Matrix::ldlt`](https://docs.rs/la-stack/latest/la_stack/struct.Matrix.html#method.ldlt) +> for details. + ## ⚡ Compile-time determinants (D ≤ 4) `det_direct()` is a `const fn` providing closed-form determinants for D=0–4, diff --git a/src/ldlt.rs b/src/ldlt.rs index fe022ea..965ab2d 100644 --- a/src/ldlt.rs +++ b/src/ldlt.rs @@ -3,6 +3,16 @@ //! 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. +//! +//! # Preconditions +//! The input matrix must be **symmetric**. This is a correctness contract, not a hint: +//! the factorization algorithm reads only the lower triangle and implicitly assumes the +//! upper triangle mirrors it. Symmetry is verified by a `debug_assert!` in debug builds +//! only; in release builds an asymmetric input will silently produce a meaningless +//! factorization. Callers who cannot statically guarantee symmetry should pre-validate +//! with [`Matrix::is_symmetric`](crate::Matrix::is_symmetric) (or locate the offending +//! pair with [`Matrix::first_asymmetry`](crate::Matrix::first_asymmetry)), or fall back +//! to [`crate::Lu`] if their matrices are not guaranteed to be symmetric at all. use core::hint::cold_path; @@ -15,6 +25,16 @@ use crate::vector::Vector; /// This factorization is **not** a general-purpose symmetric-indefinite LDLT (no pivoting). /// It assumes the input matrix is symmetric and (numerically) SPD/PSD. /// +/// # Preconditions +/// The source matrix passed to [`Matrix::ldlt`](crate::Matrix::ldlt) must be symmetric +/// (`A[i][j] == A[j][i]` within rounding). Asymmetric inputs panic in debug builds via +/// `debug_assert!` and are silently accepted in release builds — producing a +/// mathematically meaningless factorization whose [`Self::det`] and [`Self::solve_vec`] +/// results are wrong without any error being reported. Pre-validate with +/// [`Matrix::is_symmetric`](crate::Matrix::is_symmetric) when the input cannot be +/// statically guaranteed symmetric; see [`Matrix::ldlt`](crate::Matrix::ldlt) for further +/// details and alternatives. +/// /// # Storage /// The factors are stored in a single [`Matrix`]: /// - `D` is stored on the diagonal. @@ -193,17 +213,18 @@ impl Ldlt { #[cfg(debug_assertions)] fn debug_assert_symmetric(a: &Matrix) { - let scale = a.inf_norm().max(1.0); - let eps = 1e-12 * scale; - - for r in 0..D { - for c in (r + 1)..D { - let diff = (a.rows[r][c] - a.rows[c][r]).abs(); - debug_assert!( - diff <= eps, - "matrix must be symmetric (diff={diff}, eps={eps}) at ({r}, {c})" - ); - } + // Delegate to the public predicate so the runtime check and the documented + // contract on `Matrix::ldlt` cannot drift apart. `first_asymmetry` is used + // (rather than `is_symmetric`) so the panic message can name the offending + // pair — which is invaluable for debugging. + if let Some((r, c)) = a.first_asymmetry(1e-12) { + let diff = (a.rows[r][c] - a.rows[c][r]).abs(); + let eps = 1e-12 * a.inf_norm().max(1.0); + debug_assert!( + false, + "matrix must be symmetric (diff={diff}, eps={eps}) at ({r}, {c}); \ + pre-validate with Matrix::is_symmetric before calling ldlt" + ); } } @@ -433,4 +454,74 @@ mod tests { let err = ldlt.solve_vec(b).unwrap_err(); assert_eq!(err, LaError::NonFinite { row: None, col: 1 }); } + + /// Verifies the symmetry precondition documented on [`Matrix::ldlt`] is + /// enforced by `debug_assert_symmetric` in debug builds. The test is + /// gated on `debug_assertions` so `cargo test --release` simply skips it + /// (the assertion is compiled out in release builds — see the + /// Preconditions section of `Matrix::ldlt`). + #[cfg(debug_assertions)] + #[test] + #[should_panic(expected = "matrix must be symmetric")] + fn debug_asymmetric_input_panics() { + // a[0][1] = 2.0 but a[1][0] = -2.0 → clearly asymmetric. + let a = Matrix::<3>::from_rows([[4.0, 2.0, 0.0], [-2.0, 5.0, 1.0], [0.0, 1.0, 3.0]]); + let _ = a.ldlt(DEFAULT_SINGULAR_TOL); + } + + // ----------------------------------------------------------------------- + // Defensive-path coverage for `solve_vec`. + // + // `Ldlt::factor` guarantees that every stored diagonal is finite and + // strictly greater than the recorded `tol`. `solve_vec` still re-checks + // both invariants as a safety net (see the `!diag.is_finite()` and + // `diag <= self.tol` guards in the diagonal solve). Those branches are + // unreachable through the public API, so the only way to exercise them + // is to construct `Ldlt` directly with corrupt factors. The tests below + // document and verify that the safety nets return the documented error + // variants. + // ----------------------------------------------------------------------- + + macro_rules! gen_solve_vec_defensive_tests { + ($d:literal) => { + paste! { + /// `solve_vec` must surface `NonFinite` when a stored + /// diagonal is NaN, even though `factor` cannot produce + /// such a factorization. + #[test] + fn []() { + let mut factors = Matrix::<$d>::identity(); + factors.rows[$d - 1][$d - 1] = f64::NAN; + let ldlt = Ldlt::<$d> { + factors, + tol: DEFAULT_SINGULAR_TOL, + }; + let b = Vector::<$d>::new([1.0; $d]); + let err = ldlt.solve_vec(b).unwrap_err(); + assert_eq!(err, LaError::NonFinite { row: None, col: $d - 1 }); + } + + /// `solve_vec` must surface `Singular` when a stored + /// diagonal is at or below the recorded tolerance, even + /// though `factor` cannot produce such a factorization. + #[test] + fn []() { + let mut factors = Matrix::<$d>::identity(); + factors.rows[$d - 1][$d - 1] = 0.0; + let ldlt = Ldlt::<$d> { + factors, + tol: DEFAULT_SINGULAR_TOL, + }; + let b = Vector::<$d>::new([1.0; $d]); + let err = ldlt.solve_vec(b).unwrap_err(); + assert_eq!(err, LaError::Singular { pivot_col: $d - 1 }); + } + } + }; + } + + gen_solve_vec_defensive_tests!(2); + gen_solve_vec_defensive_tests!(3); + gen_solve_vec_defensive_tests!(4); + gen_solve_vec_defensive_tests!(5); } diff --git a/src/lu.rs b/src/lu.rs index 2771271..62d1817 100644 --- a/src/lu.rs +++ b/src/lu.rs @@ -501,4 +501,50 @@ mod tests { let err = lu.solve_vec(b).unwrap_err(); assert_eq!(err, LaError::NonFinite { row: None, col: 1 }); } + + // ----------------------------------------------------------------------- + // Defensive-path coverage for `solve_vec`. + // + // `Lu::factor` guarantees that every stored U diagonal satisfies + // `|U[i,i]| > tol`. `solve_vec` still re-checks during back-substitution + // as a safety net (see the `diag.abs() <= self.tol` guard). That branch + // is unreachable through the public API, so the only way to exercise it + // is to construct `Lu` directly with a corrupt U. The tests below + // document and verify that the safety net returns `Singular`. + // ----------------------------------------------------------------------- + + macro_rules! gen_solve_vec_defensive_singular_tests { + ($d:literal) => { + paste! { + /// `solve_vec` must surface `Singular` when a stored U + /// diagonal is at or below the recorded tolerance, even + /// though `factor` cannot produce such a factorization. + #[test] + fn []() { + let mut factors = Matrix::<$d>::identity(); + factors.rows[$d - 1][$d - 1] = 0.0; + + let mut piv = [0usize; $d]; + for (i, p) in piv.iter_mut().enumerate() { + *p = i; + } + + let lu = Lu::<$d> { + factors, + piv, + piv_sign: 1.0, + tol: DEFAULT_PIVOT_TOL, + }; + let b = Vector::<$d>::new([0.0; $d]); + let err = lu.solve_vec(b).unwrap_err(); + assert_eq!(err, LaError::Singular { pivot_col: $d - 1 }); + } + } + }; + } + + gen_solve_vec_defensive_singular_tests!(2); + gen_solve_vec_defensive_singular_tests!(3); + gen_solve_vec_defensive_singular_tests!(4); + gen_solve_vec_defensive_singular_tests!(5); } diff --git a/src/matrix.rs b/src/matrix.rs index 6ae4c0b..ee54844 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -152,6 +152,103 @@ impl Matrix { max_row_sum } + /// Returns `true` if the matrix is symmetric within a relative tolerance. + /// + /// Two entries `self[r][c]` and `self[c][r]` are considered equal (for the + /// purposes of symmetry) when + /// `|self[r][c] - self[c][r]| <= rel_tol * max(1.0, self.inf_norm())`. + /// This mirrors the predicate used internally by the debug-build symmetry + /// check inside [`ldlt`](Self::ldlt), so callers can pre-validate matrices + /// that may come from untrusted sources without relying on a debug-only + /// panic. + /// + /// Use [`first_asymmetry`](Self::first_asymmetry) to locate the first + /// offending pair when this returns `false`. + /// + /// # NaN / infinity handling + /// Any non-finite `|self[r][c] - self[c][r]|` (NaN or ±∞) causes this + /// predicate to return `false`. This catches both NaN off-diagonals and + /// asymmetric pairs where one side is infinite and the other is finite + /// (which would otherwise slip through when `inf_norm()` blows `eps` up + /// to `+∞` and makes `diff > eps` trivially false). A matrix whose + /// [`inf_norm`](Self::inf_norm) is `+∞` can still tolerate *finite* + /// asymmetries under an infinite `eps` — callers who need strict equality + /// on large-magnitude finite entries should validate finiteness + /// separately. + /// + /// # Panics + /// In debug builds, panics if `rel_tol` is negative or NaN; in release + /// builds these are silently treated as garbage-in garbage-out, matching + /// the convention of [`lu`](Self::lu) and [`ldlt`](Self::ldlt). + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// let a = Matrix::<2>::from_rows([[4.0, 2.0], [2.0, 3.0]]); + /// assert!(a.is_symmetric(1e-12)); + /// + /// let b = Matrix::<2>::from_rows([[4.0, 2.0], [3.0, 3.0]]); + /// assert!(!b.is_symmetric(1e-12)); + /// ``` + #[inline] + #[must_use] + pub fn is_symmetric(&self, rel_tol: f64) -> bool { + self.first_asymmetry(rel_tol).is_none() + } + + /// Returns the indices `(r, c)` (with `r < c`) of the first off-diagonal + /// pair that violates symmetry, or `None` if the matrix is symmetric + /// within `rel_tol`. + /// + /// Iteration order is row-major over the strict upper triangle, so the + /// returned indices are the lexicographically smallest such pair. The + /// predicate is the same as [`is_symmetric`](Self::is_symmetric): + /// `|self[r][c] - self[c][r]| <= rel_tol * max(1.0, self.inf_norm())`. + /// + /// # Panics + /// In debug builds, panics if `rel_tol` is negative or NaN. + /// + /// # Examples + /// ``` + /// use la_stack::prelude::*; + /// + /// let a = Matrix::<3>::from_rows([ + /// [1.0, 2.0, 0.0], + /// [2.0, 4.0, 5.0], + /// [0.0, 6.0, 9.0], // 6.0 breaks symmetry with a[1][2] = 5.0 + /// ]); + /// assert_eq!(a.first_asymmetry(1e-12), Some((1, 2))); + /// assert_eq!(Matrix::<3>::identity().first_asymmetry(1e-12), None); + /// ``` + #[inline] + #[must_use] + pub fn first_asymmetry(&self, rel_tol: f64) -> Option<(usize, usize)> { + debug_assert!( + rel_tol >= 0.0, + "rel_tol must be non-negative (got {rel_tol})" + ); + let eps = rel_tol * self.inf_norm().max(1.0); + for r in 0..D { + for c in (r + 1)..D { + let diff = (self.rows[r][c] - self.rows[c][r]).abs(); + // Any non-finite `diff` is reported as asymmetric: + // * NaN contaminates one side only, and `diff > eps` would + // silently skip it because ordered comparisons against NaN + // are always `false`. + // * ±∞ arises when exactly one of `self[r][c]` / `self[c][r]` + // is infinite; a naive `diff > eps` misses this when the + // matrix's `inf_norm()` pushes `eps` to `+∞` (because + // `∞ > ∞` is `false`). + if !diff.is_finite() || diff > eps { + cold_path(); + return Some((r, c)); + } + } + } + None + } + /// Compute an LU decomposition with partial pivoting. /// /// # Examples @@ -185,11 +282,31 @@ impl Matrix { /// This is intended for symmetric positive definite (SPD) and positive semi-definite (PSD) /// matrices such as Gram matrices. /// + /// # Preconditions + /// **The input matrix `self` must be symmetric** — that is, `self[i][j] == self[j][i]` + /// (within rounding) for all `i`, `j`. This is a *correctness* precondition, not merely + /// a performance hint. + /// + /// - In **debug builds** a `debug_assert!` verifies symmetry via + /// [`is_symmetric`](Self::is_symmetric) (relative tolerance scaled by the matrix's + /// infinity norm) and panics if it fails. + /// - In **release builds** the check is compiled out for performance. An asymmetric + /// input will be accepted silently and produce a mathematically meaningless + /// factorization — subsequent calls to [`Ldlt::det`] and [`Ldlt::solve_vec`] will + /// return wrong results with no error. + /// + /// Callers who cannot statically guarantee symmetry should pre-validate with + /// [`is_symmetric`](Self::is_symmetric) (or locate the offending pair with + /// [`first_asymmetry`](Self::first_asymmetry)) before calling `ldlt`. If you need a + /// general-purpose factorization that tolerates non-symmetric inputs, use + /// [`lu`](Self::lu) instead. + /// /// # Examples /// ``` /// use la_stack::prelude::*; /// /// # fn main() -> Result<(), LaError> { + /// // Note the symmetric layout: a[0][1] == a[1][0] == 2.0. /// let a = Matrix::<2>::from_rows([[4.0, 2.0], [2.0, 3.0]]); /// let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL)?; /// @@ -210,6 +327,9 @@ impl Matrix { /// is `<= tol` (non-positive or too small). This treats PSD degeneracy (and indefinite inputs) /// as singular/degenerate. /// Returns [`LaError::NonFinite`] if NaN/∞ is detected during factorization. + /// + /// Note that an *asymmetric* input is **not** reported as an error in release builds — + /// see the [Preconditions](#preconditions) section above. #[inline] pub fn ldlt(self, tol: f64) -> Result, LaError> { Ldlt::factor(self, tol) @@ -819,4 +939,121 @@ mod tests { gen_inf_norm_nonfinite_tests!(3); gen_inf_norm_nonfinite_tests!(4); gen_inf_norm_nonfinite_tests!(5); + + // === is_symmetric / first_asymmetry (public LDLT preconditions helpers) === + + macro_rules! gen_is_symmetric_tests { + ($d:literal) => { + paste! { + #[test] + fn []() { + let m = Matrix::<$d>::identity(); + assert!(m.is_symmetric(1e-12)); + assert_eq!(m.first_asymmetry(1e-12), None); + } + + #[test] + fn []() { + let m = Matrix::<$d>::zero(); + assert!(m.is_symmetric(1e-12)); + assert_eq!(m.first_asymmetry(1e-12), None); + } + + #[test] + fn []() { + // Construct A = M + Mᵀ so A is provably symmetric. + let mut m = [[0.0f64; $d]; $d]; + for r in 0..$d { + for c in 0..$d { + #[allow(clippy::cast_precision_loss)] + { + m[r][c] = (r * $d + c) as f64; + } + } + } + let mut sym = [[0.0f64; $d]; $d]; + for r in 0..$d { + for c in 0..$d { + sym[r][c] = m[r][c] + m[c][r]; + } + } + let a = Matrix::<$d>::from_rows(sym); + assert!(a.is_symmetric(1e-12)); + assert_eq!(a.first_asymmetry(1e-12), None); + } + + #[test] + fn []() { + // Perturb a single off-diagonal entry so symmetry fails. + let mut rows = [[0.0f64; $d]; $d]; + for i in 0..$d { + rows[i][i] = 1.0; + } + rows[0][$d - 1] = 1.0; + rows[$d - 1][0] = -1.0; // breaks symmetry + let a = Matrix::<$d>::from_rows(rows); + assert!(!a.is_symmetric(1e-12)); + assert_eq!(a.first_asymmetry(1e-12), Some((0, $d - 1))); + } + + #[test] + fn []() { + // A NaN off-diagonal must be detected as asymmetric. + let mut rows = [[0.0f64; $d]; $d]; + for i in 0..$d { + rows[i][i] = 1.0; + } + rows[0][1] = f64::NAN; + rows[1][0] = f64::NAN; + let a = Matrix::<$d>::from_rows(rows); + assert!(!a.is_symmetric(1e-12)); + // (0, 1) is the first upper-triangular pair involving the NaN. + assert_eq!(a.first_asymmetry(1e-12), Some((0, 1))); + } + } + }; + } + + gen_is_symmetric_tests!(2); + gen_is_symmetric_tests!(3); + gen_is_symmetric_tests!(4); + gen_is_symmetric_tests!(5); + + #[test] + fn is_symmetric_tolerance_scales_with_inf_norm() { + // Off-diagonal entries differ by 1e-6. With inf_norm ≈ 2e6, the + // relative tolerance 1e-12 yields eps ≈ 2e-6, which accepts the gap; + // a stricter tol of 1e-15 rejects it. + let a = Matrix::<2>::from_rows([[1.0e6, 1.0e6 + 1.0e-6], [1.0e6, 1.0e6]]); + assert!(a.is_symmetric(1e-12)); + assert!(!a.is_symmetric(1e-15)); + } + + #[test] + fn first_asymmetry_returns_lexicographically_first_pair() { + // Two asymmetric pairs: (0, 2) and (1, 2). We must get (0, 2) first. + let a = Matrix::<3>::from_rows([[1.0, 0.0, 2.0], [0.0, 1.0, 3.0], [-2.0, -3.0, 1.0]]); + assert_eq!(a.first_asymmetry(1e-12), Some((0, 2))); + } + + /// Regression: a single infinite off-diagonal paired with a finite entry + /// used to slip through as "symmetric" because `inf_norm()` blew `eps` up + /// to `+∞` and `∞ > ∞` evaluates to `false`. After the fix, any + /// non-finite `|a[r][c] - a[c][r]|` is reported as an asymmetry regardless + /// of `eps`. + #[test] + fn first_asymmetry_flags_infinite_offdiagonal_against_finite() { + let a = Matrix::<2>::from_rows([[1.0, f64::INFINITY], [0.0, 1.0]]); + assert_eq!(a.first_asymmetry(1e-12), Some((0, 1))); + assert!(!a.is_symmetric(1e-12)); + } + + #[cfg(debug_assertions)] + #[test] + #[should_panic(expected = "rel_tol must be non-negative")] + fn first_asymmetry_debug_panics_on_negative_tol() { + // Mirrors the `debug_assert!(tol >= 0.0)` convention used by + // `Matrix::lu` / `Matrix::ldlt`. + let _ = Matrix::<2>::identity().first_asymmetry(-1.0); + } }