From 1693307c58e30a804e9b79abc783af99be8dc947 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Mon, 20 Apr 2026 12:28:45 -0700 Subject: [PATCH 1/3] docs: strengthen LDLT symmetry precondition and add is_symmetric API MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit LDLT factorization assumes the input matrix is symmetric, but the contract was only implied by one sentence in the module and struct docs. Asymmetric inputs silently produce mathematically meaningless factorizations in release builds (the `debug_assert_symmetric` check is compiled out), and callers had no supported way to validate symmetry up front. Documentation: - `Matrix::ldlt`: add a prominent `# Preconditions` section spelling out the contract, the debug-vs-release split, and pointers to the new `is_symmetric` / `first_asymmetry` predicates and to `lu()` as the fallback for inputs that may not be symmetric at all. - `src/ldlt.rs`: promote the precondition into a module-level and struct-level `# Preconditions` section with back-links. - `README.md`: add a warning call-out under the LDLT example linking the new predicates. New public API (`Matrix`): - `is_symmetric(&self, rel_tol: f64) -> bool` — infallible predicate sharing the `|A[r][c] - A[c][r]| <= rel_tol * max(1, inf_norm)` convention used internally by LDLT. - `first_asymmetry(&self, rel_tol: f64) -> Option<(usize, usize)>` — returns the lexicographically first off-diagonal pair that violates symmetry, or `None` for symmetric matrices. Used by the debug-build check for pinpointed panic messages. - Both `debug_assert!(rel_tol >= 0.0)`, matching `lu(tol)` / `ldlt(tol)`. - NaN off-diagonals are explicitly reported as asymmetric. Refactor: - `debug_assert_symmetric` now delegates to `first_asymmetry`, so the runtime check and the documented contract share one implementation and cannot drift apart. Tests: - Dimension-generic (D=2..=5): identity, zero, `A = M + Mᵀ`, perturbed off-diagonal asymmetric, NaN off-diagonal asymmetric. - Scalar: tolerance scaling with inf_norm, lexicographic-first pair on D=3, debug-only panic on negative `rel_tol`. - `#[cfg(debug_assertions)] #[should_panic]` test for the LDLT debug-build panic still passes against the refactored assertion. No functional change in release builds of existing APIs. Tests: `just ci` — 123 lib, 26 doc, 284 exact-feature, 101 Python, all examples, all linters/validators. Closes #84 Co-Authored-By: Oz --- README.md | 11 +++ src/ldlt.rs | 57 ++++++++++--- src/matrix.rs | 216 ++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 273 insertions(+), 11 deletions(-) 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..9cac458 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,18 @@ 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); + } } diff --git a/src/matrix.rs b/src/matrix.rs index 6ae4c0b..f51605d 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -152,6 +152,94 @@ 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 NaN off-diagonal entry causes this predicate to return `false` + /// (because `NaN <= eps` is `false`). A matrix whose [`inf_norm`](Self::inf_norm) + /// is `+∞` can produce an `eps` of `+∞`, under which any finite asymmetry + /// is tolerated — callers who need strict equality on infinite 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(); + // NaN is reported as asymmetric: a NaN entry contaminates `diff`, + // and `diff > eps` alone would silently skip it because ordered + // comparisons against NaN are always `false`. + if diff.is_nan() || diff > eps { + cold_path(); + return Some((r, c)); + } + } + } + None + } + /// Compute an LU decomposition with partial pivoting. /// /// # Examples @@ -185,11 +273,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 +318,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 +930,109 @@ 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))); + } + + #[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); + } } From 87d426fca1627445b804fd26b62fc7d9d4f0ae48 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Mon, 20 Apr 2026 13:00:45 -0700 Subject: [PATCH 2/3] Added: defensive-path test coverage for LU and LDLT solve_vec Add unit tests to exercise internal safety nets in the LU and LDLT diagonal solve routines. These tests manually construct factors with invalid diagonals (NaN or sub-tolerance) to verify that solve_vec correctly surfaces NonFinite and Singular errors, even though these states are unreachable via the standard factorization APIs. Refs: #84 --- src/ldlt.rs | 56 +++++++++++++++++++++++++++++++++++++++++++++++++++++ src/lu.rs | 46 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 102 insertions(+) diff --git a/src/ldlt.rs b/src/ldlt.rs index 9cac458..965ab2d 100644 --- a/src/ldlt.rs +++ b/src/ldlt.rs @@ -468,4 +468,60 @@ mod tests { 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); } From 1805779dbca49183fbfa95c68ec00984966aa551 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Mon, 20 Apr 2026 13:19:30 -0700 Subject: [PATCH 3/3] Fixed: report infinite vs finite off-diagonal pairs as asymmetric Update Matrix::first_asymmetry to flag any non-finite difference between off-diagonal pairs as an asymmetry. This prevents cases where a single infinite entry paired with a finite entry would incorrectly pass as symmetric because the matrix's infinite norm blew the tolerance up to infinity, making the comparison `diff > eps` return false. Refs: #84 --- src/matrix.rs | 39 ++++++++++++++++++++++++++++++--------- 1 file changed, 30 insertions(+), 9 deletions(-) diff --git a/src/matrix.rs b/src/matrix.rs index f51605d..ee54844 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -166,11 +166,15 @@ impl Matrix { /// offending pair when this returns `false`. /// /// # NaN / infinity handling - /// Any NaN off-diagonal entry causes this predicate to return `false` - /// (because `NaN <= eps` is `false`). A matrix whose [`inf_norm`](Self::inf_norm) - /// is `+∞` can produce an `eps` of `+∞`, under which any finite asymmetry - /// is tolerated — callers who need strict equality on infinite entries - /// should validate finiteness separately. + /// 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 @@ -228,10 +232,15 @@ impl Matrix { for r in 0..D { for c in (r + 1)..D { let diff = (self.rows[r][c] - self.rows[c][r]).abs(); - // NaN is reported as asymmetric: a NaN entry contaminates `diff`, - // and `diff > eps` alone would silently skip it because ordered - // comparisons against NaN are always `false`. - if diff.is_nan() || diff > eps { + // 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)); } @@ -1027,6 +1036,18 @@ mod tests { 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")]