feat: Add Bunch–Kaufman LDL Decomposition (#1591)

Author: awinick

src/base/norm.rs

@@ -39,13 +39,44 @@ pub trait Norm<T: SimdComplexField> {
         ShapeConstraint: SameNumberOfRows<R1, R2> + SameNumberOfColumns<C1, C2>;
 }
 
-/// Euclidean norm.
+/// Euclidean norm of a vector, or Frobenius norm of a matrix.
+///
+/// Computes sqrt(sum |a_ij|^2) over all elements.
+///
+/// <div class="warning">
+/// For matrices, this is the Frobenius norm, not the matrix 2-norm (spectral norm).
+/// </div>
 #[derive(Copy, Clone, Debug)]
 pub struct EuclideanNorm;
-/// Lp norm.
+
+/// Entrywise Lp norm of a matrix or vector.
+///
+/// Computes (sum |a_ij|^p)^(1/p) over all elements.
+///
+/// <div class="warning">
+/// This does not match the standard mathematical definition of the matrix Lp norm.
+/// </div>
 #[derive(Copy, Clone, Debug)]
 pub struct LpNorm(pub i32);
-/// L-infinite norm aka. Chebytchev norm aka. uniform norm aka. suppremum norm.
+
+/// Induced matrix 1-norm (maximum absolute column sum).
+///
+/// Computes max_j sum_i |a_ij|.
+///
+/// For a column vector, this is the L1 norm.
+/// For a row vector, this is the L-infinity norm.
+#[derive(Copy, Clone, Debug)]
+pub struct OneNorm;
+
+/// Entrywise L-infinity norm of a matrix or vector.
+///
+/// Computes max |a_ij| over all elements.
+///
+/// For a vector this is the standard L-infinity norm.
+///
+/// <div class="warning">
+/// For matrices, this is the entrywise maximum, not the induced matrix infinity-norm.
+/// </div>
 #[derive(Copy, Clone, Debug)]
 pub struct UniformNorm;
 
@@ -120,6 +151,45 @@ impl<T: SimdComplexField> Norm<T> for LpNorm {
     }
 }
 
+impl<T: SimdComplexField> Norm<T> for OneNorm {
+    #[inline]
+    fn norm<R, C, S>(&self, m: &Matrix<T, R, C, S>) -> T::SimdRealField
+    where
+        R: Dim,
+        C: Dim,
+        S: Storage<T, R, C>,
+    {
+        m.column_iter()
+            .map(|col| col.fold(T::SimdRealField::zero(), |a, b| a + b.simd_modulus()))
+            .fold(T::SimdRealField::zero(), T::SimdRealField::simd_max)
+    }
+
+    #[inline]
+    fn metric_distance<R1, C1, S1, R2, C2, S2>(
+        &self,
+        m1: &Matrix<T, R1, C1, S1>,
+        m2: &Matrix<T, R2, C2, S2>,
+    ) -> T::SimdRealField
+    where
+        R1: Dim,
+        C1: Dim,
+        S1: Storage<T, R1, C1>,
+        R2: Dim,
+        C2: Dim,
+        S2: Storage<T, R2, C2>,
+        ShapeConstraint: SameNumberOfRows<R1, R2> + SameNumberOfColumns<C1, C2>,
+    {
+        m1.column_iter()
+            .zip(m2.column_iter())
+            .map(|(c1, c2)| {
+                c1.zip_fold(&c2, T::SimdRealField::zero(), |acc, a, b| {
+                    acc + (a - b).simd_modulus()
+                })
+            })
+            .fold(T::SimdRealField::zero(), T::SimdRealField::simd_max)
+    }
+}
+
 impl<T: SimdComplexField> Norm<T> for UniformNorm {
     #[inline]
     fn norm<R, C, S>(&self, m: &Matrix<T, R, C, S>) -> T::SimdRealField
@@ -158,8 +228,11 @@ impl<T: SimdComplexField> Norm<T> for UniformNorm {
 }
 
 /// # Magnitude and norms
+///
+/// Unless otherwise noted, the norm used throughout is the L2 norm for vectors and
+/// the Frobenius norm for matrices.
 impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
-    /// The squared L2 norm of this vector.
+    /// Squared L2 norm of this vector, or squared Frobenius norm of this matrix.
     #[inline]
     #[must_use]
     pub fn norm_squared(&self) -> T::SimdRealField
@@ -176,7 +249,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
         res
     }
 
-    /// The L2 norm of this matrix.
+    /// L2 norm of this vector, or Frobenius norm of this matrix.
     ///
     /// Use `.apply_norm` to apply a custom norm.
     #[inline]
@@ -188,7 +261,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
         self.norm_squared().simd_sqrt()
     }
 
-    /// Compute the distance between `self` and `rhs` using the metric induced by the euclidean norm.
+    /// Distance between `self` and `rhs` using the L2 norm for vectors, or the Frobenius for matrices.
     ///
     /// Use `.apply_metric_distance` to apply a custom norm.
     #[inline]
@@ -306,7 +379,13 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
         self.unscale(self.norm())
     }
 
-    /// The Lp norm of this matrix.
+    /// Entrywise Lp norm of this matrix or vector.
+    ///
+    /// Computes (sum |a_ij|^p)^(1/p) over all elements.
+    ///
+    /// <div class="warning">
+    /// For matrices, this does not match the standard mathematical definition of the matrix Lp norm.
+    /// </div>
     #[inline]
     #[must_use]
     pub fn lp_norm(&self, p: i32) -> T::SimdRealField
@@ -316,6 +395,21 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
         self.apply_norm(&LpNorm(p))
     }
 
+    /// Induced matrix 1-norm (maximum absolute column sum).
+    ///
+    /// Computes max_j sum_i |a_ij|.
+    ///
+    /// For a column vector, this is the L1 norm.
+    /// For a row vector, this is the L-infinity norm.
+    #[inline]
+    #[must_use]
+    pub fn one_norm(&self) -> T::SimdRealField
+    where
+        T: SimdComplexField,
+    {
+        self.apply_norm(&OneNorm)
+    }
+
     /// Attempts to normalize `self`.
     ///
     /// The components of this matrix can be SIMD types.

src/linalg/decomposition.rs

@@ -1,8 +1,8 @@
 use crate::storage::Storage;
 use crate::{
     Allocator, Bidiagonal, Cholesky, ColPivQR, ComplexField, DefaultAllocator, Dim, DimDiff,
-    DimMin, DimMinimum, DimSub, FullPivLU, Hessenberg, LU, Matrix, OMatrix, QR, RealField, SVD,
-    Schur, SymmetricEigen, SymmetricTridiagonal, U1, UDU,
+    DimMin, DimMinimum, DimSub, FullPivLU, Hessenberg, LBLT, LU, Matrix, OMatrix, QR, RealField,
+    SVD, Schur, SymmetricEigen, SymmetricTridiagonal, U1, UDU,
 };
 
 /// # Rectangular matrix decomposition
@@ -244,6 +244,7 @@ impl<T: ComplexField, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
 /// | -------------------------|---------------------------|--------------|
 /// | Hessenberg               | `Q * H * Qᵀ`             | `Q` is a unitary matrix and `H` an upper-Hessenberg matrix. |
 /// | Cholesky                 | `L * Lᵀ`                 | `L` is a lower-triangular matrix. |
+/// | LBLT decomposition       | `Pᵀ * L * B * Lᴴ * P`    | `L` is unit lower-triangular, `B` is Hermitian block-diagonal, and `P` is a permutation matrix. |
 /// | UDU                      | `U * D * Uᵀ`             | `U` is a upper-triangular matrix, and `D` a diagonal matrix. |
 /// | Schur decomposition      | `Q * T * Qᵀ`             | `Q` is an unitary matrix and `T` a quasi-upper-triangular matrix. |
 /// | Symmetric eigendecomposition | `Q ~ Λ ~ Qᵀ`   | `Q` is an unitary matrix, and `Λ` is a real diagonal matrix. |
@@ -260,6 +261,18 @@ impl<T: ComplexField, D: Dim, S: Storage<T, D, D>> Matrix<T, D, D, S> {
         Cholesky::new(self.into_owned())
     }
 
+    /// Computes the LBLT decomposition of this matrix.
+    /// The input matrix `self` is assumed to be Hermitian (symmetric) and the decomposition will
+    /// only read the lower-triangular part of `self`.
+    pub fn lblt(self) -> LBLT<T, D>
+    where
+        T: Copy,
+        T::RealField: Copy,
+        DefaultAllocator: Allocator<D> + Allocator<D, D>,
+    {
+        LBLT::new(self.into_owned())
+    }
+
     /// Attempts to compute the UDU decomposition of this matrix.
     ///
     /// The input matrix `self` is assumed to be symmetric and this decomposition will only read

src/linalg/exp.rs

@@ -156,31 +156,49 @@ where
     fn d4_tight(&mut self) -> T::RealField {
         if self.d4_exact.is_none() {
             self.calc_a4();
-            self.d4_exact = Some(one_norm(self.a4.as_ref().unwrap()).powf(convert(0.25)));
+            self.d4_exact = Some(self.a4.as_ref().unwrap().one_norm().powf(convert(0.25)));
         }
         self.d4_exact.clone().unwrap()
     }
 
     fn d6_tight(&mut self) -> T::RealField {
         if self.d6_exact.is_none() {
             self.calc_a6();
-            self.d6_exact = Some(one_norm(self.a6.as_ref().unwrap()).powf(convert(1.0 / 6.0)));
+            self.d6_exact = Some(
+                self.a6
+                    .as_ref()
+                    .unwrap()
+                    .one_norm()
+                    .powf(convert(1.0 / 6.0)),
+            );
         }
         self.d6_exact.clone().unwrap()
     }
 
     fn d8_tight(&mut self) -> T::RealField {
         if self.d8_exact.is_none() {
             self.calc_a8();
-            self.d8_exact = Some(one_norm(self.a8.as_ref().unwrap()).powf(convert(1.0 / 8.0)));
+            self.d8_exact = Some(
+                self.a8
+                    .as_ref()
+                    .unwrap()
+                    .one_norm()
+                    .powf(convert(1.0 / 8.0)),
+            );
         }
         self.d8_exact.clone().unwrap()
     }
 
     fn d10_tight(&mut self) -> T::RealField {
         if self.d10_exact.is_none() {
             self.calc_a10();
-            self.d10_exact = Some(one_norm(self.a10.as_ref().unwrap()).powf(convert(1.0 / 10.0)));
+            self.d10_exact = Some(
+                self.a10
+                    .as_ref()
+                    .unwrap()
+                    .one_norm()
+                    .powf(convert(1.0 / 10.0)),
+            );
         }
         self.d10_exact.clone().unwrap()
     }
@@ -196,7 +214,7 @@ where
 
         if self.d4_approx.is_none() {
             self.calc_a4();
-            self.d4_approx = Some(one_norm(self.a4.as_ref().unwrap()).powf(convert(0.25)));
+            self.d4_approx = Some(self.a4.as_ref().unwrap().one_norm().powf(convert(0.25)));
         }
 
         self.d4_approx.clone().unwrap()
@@ -213,7 +231,13 @@ where
 
         if self.d6_approx.is_none() {
             self.calc_a6();
-            self.d6_approx = Some(one_norm(self.a6.as_ref().unwrap()).powf(convert(1.0 / 6.0)));
+            self.d6_approx = Some(
+                self.a6
+                    .as_ref()
+                    .unwrap()
+                    .one_norm()
+                    .powf(convert(1.0 / 6.0)),
+            );
         }
 
         self.d6_approx.clone().unwrap()
@@ -230,7 +254,13 @@ where
 
         if self.d8_approx.is_none() {
             self.calc_a8();
-            self.d8_approx = Some(one_norm(self.a8.as_ref().unwrap()).powf(convert(1.0 / 8.0)));
+            self.d8_approx = Some(
+                self.a8
+                    .as_ref()
+                    .unwrap()
+                    .one_norm()
+                    .powf(convert(1.0 / 8.0)),
+            );
         }
 
         self.d8_approx.clone().unwrap()
@@ -247,7 +277,13 @@ where
 
         if self.d10_approx.is_none() {
             self.calc_a10();
-            self.d10_approx = Some(one_norm(self.a10.as_ref().unwrap()).powf(convert(1.0 / 10.0)));
+            self.d10_approx = Some(
+                self.a10
+                    .as_ref()
+                    .unwrap()
+                    .one_norm()
+                    .powf(convert(1.0 / 10.0)),
+            );
         }
 
         self.d10_approx.clone().unwrap()
@@ -430,7 +466,7 @@ where
     let choose_2m_m = factorial(2 * m) / (m_factorial * m_factorial);
 
     let abs_c_recip = choose_2m_m * factorial(2 * m + 1);
-    let alpha = a_abs_onenorm / one_norm(a);
+    let alpha = a_abs_onenorm / a.one_norm();
     let alpha: f64 = try_convert::<_, f64>(alpha).unwrap() / abs_c_recip as f64;
 
     let u = 2_f64.powf(-53.0);
@@ -451,27 +487,6 @@ where
     q.lu().solve(&p).unwrap()
 }
 
-fn one_norm<T, D>(m: &OMatrix<T, D, D>) -> T::RealField
-where
-    T: ComplexField,
-    D: Dim,
-    DefaultAllocator: Allocator<D, D>,
-{
-    let mut max = <T as ComplexField>::RealField::zero();
-
-    for i in 0..m.ncols() {
-        let col = m.column(i);
-        max = max.max(
-            col.iter()
-                .fold(<T as ComplexField>::RealField::zero(), |a, b| {
-                    a + b.clone().abs()
-                }),
-        );
-    }
-
-    max
-}
-
 impl<T: ComplexField, D> OMatrix<T, D, D>
 where
     D: DimMin<D, Output = D>,
@@ -539,15 +554,3 @@ where
         x
     }
 }
-
-#[cfg(test)]
-mod tests {
-    #[test]
-    #[allow(clippy::float_cmp)]
-    fn one_norm() {
-        use crate::Matrix3;
-        let m = Matrix3::new(-3.0, 5.0, 7.0, 2.0, 6.0, 4.0, 0.0, 2.0, 8.0);
-
-        assert_eq!(super::one_norm(&m), 19.0);
-    }
-}