docs: strengthen LDLT symmetry precondition and add is_symmetric API

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<D>`):
- `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 <oz-agent@warp.dev>

Author: acgetchell

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,

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<const D: usize> Ldlt<D> {
 
 #[cfg(debug_assertions)]
 fn debug_assert_symmetric<const D: usize>(a: &Matrix<D>) {
-    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);
+    }
 }

src/matrix.rs

@@ -152,6 +152,94 @@ impl<const D: usize> Matrix<D> {
         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<const D: usize> Matrix<D> {
     /// 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<const D: usize> Matrix<D> {
     /// 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<Ldlt<D>, 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 [<is_symmetric_true_for_identity_ $d d>]() {
+                    let m = Matrix::<$d>::identity();
+                    assert!(m.is_symmetric(1e-12));
+                    assert_eq!(m.first_asymmetry(1e-12), None);
+                }
+
+                #[test]
+                fn [<is_symmetric_true_for_zero_ $d d>]() {
+                    let m = Matrix::<$d>::zero();
+                    assert!(m.is_symmetric(1e-12));
+                    assert_eq!(m.first_asymmetry(1e-12), None);
+                }
+
+                #[test]
+                fn [<is_symmetric_true_for_constructed_symmetric_ $d d>]() {
+                    // 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 [<is_symmetric_false_for_asymmetric_offdiagonal_ $d d>]() {
+                    // 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 [<is_symmetric_false_for_nan_offdiagonal_ $d d>]() {
+                    // 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);
+    }
 }