la_stack/
gram.rs

1#![forbid(unsafe_code)]
2
3//! Fixed-size Gram construction.
4
5use crate::{LaError, Matrix, Vector};
6
7/// Construct a stack-backed [`Matrix<M>`] of pairwise vector dot products.
8///
9/// A Gram matrix records pairwise inner products: diagonal entries are squared
10/// vector lengths, and off-diagonal entries encode their relative angles.
11/// The input contains `M` finite-by-construction [`Vector<N>`] values.
12/// If `V` has these vectors as rows, the mathematical Gram matrix is `G = V Vᵀ`,
13/// with `G[i,j] = vectors[i] · vectors[j]`. In exact real arithmetic and for
14/// `M ≤ N`, its determinant is the squared volume of the spanned
15/// parallelotope. For edges from one simplex vertex, the simplex volume is
16/// `sqrt(det(G)) / M!`; this also handles facets embedded in higher dimensions.
17/// See `REFERENCES.md` \[16\] for the Gram determinant and volume interpretation.
18///
19/// Each upper-triangle dot product is computed once using [`Vector::dot`]'s
20/// left-to-right fused multiply-add reduction and copied to the other triangle,
21/// giving bit-for-bit symmetry. No absolute rounding-error bound is provided.
22/// Rounding and underflow can destroy positive semidefiniteness or rank; this
23/// operation proves neither positive definiteness nor affine independence.
24/// [`Matrix::ldlt`] retains its symmetry and positive-definiteness preconditions.
25/// Forming a Gram matrix squares the spectral condition number of an exact input
26/// with full row rank. See the floating-point discussion in `REFERENCES.md` \[9-11\].
27///
28/// `M` and `N` are independent, with no dimension cap (including dimensions
29/// through 8); storage is `O(M²)` and work is `O(M²(N + 1))`, including output
30/// initialization when `N = 0`. `M = 0` returns an empty
31/// matrix; `N = 0` returns an all-zero matrix. No optional feature is required.
32///
33/// # Errors
34/// Returns [`LaError::NonFinite`] if a dot-product accumulator overflows, even
35/// when the exact result would be finite after cancellation. The error preserves
36/// [`Vector::dot`]'s [`Computation`](crate::NonFiniteOrigin::Computation) origin,
37/// [`VectorDotProduct`](crate::ArithmeticOperation::VectorDotProduct) operation,
38/// and first failing reduction [`Step`](crate::NonFiniteLocation::Step).
39/// The step indexes the vector coordinate, not the output matrix cell.
40/// Pairs are visited in upper-triangle row order.
41///
42/// # Examples
43/// ```
44/// use la_stack::prelude::*;
45/// # fn main() -> Result<(), LaError> {
46/// let vectors = [Vector::try_new([1.0, 0.0, 0.0])?,
47///                Vector::try_new([0.0, 2.0, 0.0])?];
48/// let gram = gram_matrix(&vectors)?;
49/// assert_eq!(gram.as_rows(), &[[1.0, 0.0], [0.0, 4.0]]);
50/// assert_eq!(gram.ldlt(Tolerance::try_new(0.0)?)?.det()?, 4.0);
51/// # Ok(())
52/// # }
53/// ```
54pub const fn gram_matrix<const M: usize, const N: usize>(
55    vectors: &[Vector<N>; M],
56) -> Result<Matrix<M>, LaError> {
57    let mut rows = [[0.0; M]; M];
58    let mut i = 0;
59    while i < M {
60        let mut j = i;
61        while j < M {
62            let value = match vectors[i].dot(&vectors[j]) {
63                Ok(value) => value,
64                Err(error) => return Err(error),
65            };
66            rows[i][j] = value;
67            rows[j][i] = value;
68            j += 1;
69        }
70        i += 1;
71    }
72    Matrix::try_from_rows(rows)
73}
la_stack/gram.rs

la_stack/
gram.rs
