diff --git a/src/base/norm.rs b/src/base/norm.rs index 37bdcef85..90c194ad8 100644 --- a/src/base/norm.rs +++ b/src/base/norm.rs @@ -39,13 +39,44 @@ pub trait Norm { ShapeConstraint: SameNumberOfRows + SameNumberOfColumns; } -/// Euclidean norm. +/// Euclidean norm of a vector, or Frobenius norm of a matrix. +/// +/// Computes sqrt(sum |a_ij|^2) over all elements. +/// +///
+/// For matrices, this is the Frobenius norm, not the matrix 2-norm (spectral norm). +///
#[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. +/// +///
+/// This does not match the standard mathematical definition of the matrix Lp norm. +///
#[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. +/// +///
+/// For matrices, this is the entrywise maximum, not the induced matrix infinity-norm. +///
#[derive(Copy, Clone, Debug)] pub struct UniformNorm; @@ -120,6 +151,45 @@ impl Norm for LpNorm { } } +impl Norm for OneNorm { + #[inline] + fn norm(&self, m: &Matrix) -> T::SimdRealField + where + R: Dim, + C: Dim, + S: Storage, + { + 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( + &self, + m1: &Matrix, + m2: &Matrix, + ) -> T::SimdRealField + where + R1: Dim, + C1: Dim, + S1: Storage, + R2: Dim, + C2: Dim, + S2: Storage, + ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, + { + 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 Norm for UniformNorm { #[inline] fn norm(&self, m: &Matrix) -> T::SimdRealField @@ -158,8 +228,11 @@ impl Norm 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> Matrix { - /// 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> Matrix { 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> Matrix { 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> Matrix { 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. + /// + ///
+ /// For matrices, this does not match the standard mathematical definition of the matrix Lp norm. + ///
#[inline] #[must_use] pub fn lp_norm(&self, p: i32) -> T::SimdRealField @@ -316,6 +395,21 @@ impl> Matrix { 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. diff --git a/src/linalg/decomposition.rs b/src/linalg/decomposition.rs index b424b799a..c22ffceae 100644 --- a/src/linalg/decomposition.rs +++ b/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> Matrix { /// | -------------------------|---------------------------|--------------| /// | 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> Matrix { 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 + where + T: Copy, + T::RealField: Copy, + DefaultAllocator: Allocator + Allocator, + { + 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 diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs index be090f643..82b76a2cd 100644 --- a/src/linalg/exp.rs +++ b/src/linalg/exp.rs @@ -156,7 +156,7 @@ 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() } @@ -164,7 +164,13 @@ where 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() } @@ -172,7 +178,13 @@ where 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() } @@ -180,7 +192,13 @@ where 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(m: &OMatrix) -> T::RealField -where - T: ComplexField, - D: Dim, - DefaultAllocator: Allocator, -{ - let mut max = ::RealField::zero(); - - for i in 0..m.ncols() { - let col = m.column(i); - max = max.max( - col.iter() - .fold(::RealField::zero(), |a, b| { - a + b.clone().abs() - }), - ); - } - - max -} - impl OMatrix where D: DimMin, @@ -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); - } -} diff --git a/src/linalg/lblt.rs b/src/linalg/lblt.rs new file mode 100644 index 000000000..52a3ae5e6 --- /dev/null +++ b/src/linalg/lblt.rs @@ -0,0 +1,503 @@ +use crate::{ + ComplexField, DefaultAllocator, Dim, Matrix, OMatrix, OVector, RealField, Storage, U1, + allocator::Allocator, convert, +}; +use num_traits::{One, Zero}; + +#[cfg(feature = "serde-serialize-no-std")] +use serde::{Deserialize, Serialize}; + +/// A pivot represented as `(index, block_size)`, where `block_size` is either 1 or 2. +type Pivot = (usize, usize); + +/// Bunch–Kaufman LBL^H factorization of a Hermitian matrix with symmetric pivoting. +#[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))] +#[cfg_attr( + feature = "serde-serialize-no-std", + serde(bound(serialize = "DefaultAllocator: Allocator, + OMatrix: Serialize, + OVector: Serialize, + Option: Serialize")) +)] +#[cfg_attr( + feature = "serde-serialize-no-std", + serde(bound(deserialize = "DefaultAllocator: Allocator, + OMatrix: Deserialize<'de>, + OVector: Deserialize<'de>, + Option: Deserialize<'de>")) +)] +#[cfg_attr(feature = "defmt", derive(defmt::Format))] +#[derive(Clone, Debug)] +pub struct LBLT +where + DefaultAllocator: Allocator + Allocator, +{ + matrix: OMatrix, + pivots: OVector, + zero_pivot: Option, +} + +impl LBLT +where + T::RealField: Copy, + DefaultAllocator: Allocator + Allocator, +{ + /// Compute the factorization of a complex Hermitian matrix using the Bunch-Kaufman + /// block-diagonal pivoting method: + /// + /// P A P^T = L * B * L^H + /// + /// where P is the permutation induced by the pivot sequence, L is unit lower + /// triangular in the permuted basis, and B is Hermitian block diagonal with + /// 1-by-1 and 2-by-2 diagonal blocks. + /// + /// This implementation follows the partial pivoting (Algorithm A) variant from + /// Bunch & Kaufman (1977), which is also the basis for LAPACK’s `?sytrf/?hetrf` routines. + pub fn new(mut matrix: OMatrix) -> Self { + assert!(matrix.is_square()); + let n = matrix.nrows(); + + let mut pivots = OVector::from_element_generic(matrix.shape_generic().0, U1, (0, 0)); + let mut zero_pivot = None; + + // Bunch–Kaufman pivot threshold: (1 + sqrt(17)) / 8 + let alpha: T::RealField = convert(0.6403882032022076); + + // current pivot position + let mut k = 0; + + while k < n { + let mut block_size = 1; + + // Ensure the diagonal element is real + matrix[(k, k)] = T::from_real(matrix[(k, k)].real()); + let diag_abs = matrix[(k, k)].real().abs(); + + // Row index and magnitude of the largest off-diagonal entry in the active part + // of column k. + let (imax, colmax) = if k + 1 < n { + let mut imax = k + 1; + let mut colmax = matrix[(imax, k)].norm1(); + for i in (k + 2)..n { + let magnitude = matrix[(i, k)].norm1(); + if magnitude > colmax { + imax = i; + colmax = magnitude; + } + } + (imax, colmax) + } else { + // If k is the last column, there is no off-diagonal candidate. + (0, T::RealField::zero()) + }; + + if diag_abs.max(colmax) == T::RealField::zero() { + // Column k is zero: store a 1x1 pivot, and skip all other logic. + if zero_pivot.is_none() { + zero_pivot = Some(k); + } + + pivots[k] = (k, 1); + k += 1; + continue; + } + + let pivot_index: usize; + + if diag_abs < alpha * colmax { + let mut rowmax = T::RealField::zero(); + for j in k..imax { + rowmax = rowmax.max(matrix[(imax, j)].norm1()); + } + for j in (imax + 1)..matrix.nrows() { + rowmax = rowmax.max(matrix[(j, imax)].norm1()); + } + + if diag_abs >= alpha * colmax * (colmax / rowmax) { + // Even though A[k, k] is not diagonally dominant, it is still large enough + // compared with the candidate row/column growth, so keep a 1x1 pivot at k. + pivot_index = k; + } else { + pivot_index = imax; + + if matrix[(imax, imax)].real().abs() < alpha * rowmax { + // The candidate diagonal at imax is also too small relative to its row + // maximum, so use a 2x2 pivot block involving k and k+1. + block_size = 2; + } + } + } else { + // The diagonal dominates column k strongly enough to use A[k, k] as a 1x1 pivot + // without any row/column interchange. + pivot_index = k; + } + + let pivot_target = k + block_size - 1; + + if pivot_index != pivot_target { + // Hermitian two-sided interchange for the chosen pivot. + for i in (pivot_index + 1)..matrix.nrows() { + // Swap entries below both indices + matrix.swap((i, pivot_target), (i, pivot_index)); + } + + for j in (pivot_target + 1)..pivot_index { + // Swap the strip between the two indices. + matrix.swap((j, pivot_target), (pivot_index, j)); + matrix[(j, pivot_target)] = matrix[(j, pivot_target)].conjugate(); + matrix[(pivot_index, j)] = matrix[(pivot_index, j)].conjugate(); + } + + // The cross entry between the swapped indices remains in the same slot. + matrix[(pivot_index, pivot_target)] = + matrix[(pivot_index, pivot_target)].conjugate(); + + // Swap the diagonal entries. + matrix.swap((pivot_target, pivot_target), (pivot_index, pivot_index)); + + if k + 1 == pivot_target { + // For a 2x2 pivot block, move the off-diagonal block entry. + matrix.swap((k + 1, k), (pivot_index, k)); + } + } + + if block_size == 1 { + // 1x1 pivot block D(k) + if k + 1 < n { + let inv_diag = T::RealField::one() / matrix[(k, k)].real(); + + for j in (k + 1)..n { + let jk_conj = matrix[(j, k)].conjugate(); + + for i in j..n { + // Rank-1 Hermitian update of the trailing submatrix. + matrix[(i, j)] = + matrix[(i, j)] + matrix[(i, k)].scale(-inv_diag) * jk_conj; + } + + // Keep the Hermitian diagonal explicitly real. + matrix[(j, j)] = T::from_real(matrix[(j, j)].real()); + } + + // Normalize column k so that it stores the multipliers of L. + for i in (k + 1)..n { + matrix[(i, k)] = matrix[(i, k)].scale(inv_diag); + } + } + + pivots[k] = (pivot_index, 1); + } else { + // 2x2 pivot block D(k:k+1) + if k + 2 < n { + // Form the scaled inverse-coefficient data for the 2x2 Hermitian pivot block. + let d = matrix[(k + 1, k)].abs(); + let d11 = matrix[(k + 1, k + 1)].real() / d; + let d22 = matrix[(k, k)].real() / d; + let d21 = matrix[(k + 1, k)].unscale(d); + let scale = T::RealField::one() / (d * (d11 * d22 - T::RealField::one())); + + for j in (k + 2)..n { + // These are the two transformed entries for row j. Together they represent + // the action of inv(D(k:k+1)) on the stored columns k and k+1. + let work1 = + (matrix[(j, k)].scale(d11) - matrix[(j, k + 1)] * d21).scale(scale); + let work2 = (matrix[(j, k + 1)].scale(d22) + - matrix[(j, k)] * d21.conjugate()) + .scale(scale); + + for i in j..n { + // Rank-2 Hermitian update of the trailing submatrix. + matrix[(i, j)] = matrix[(i, j)] + - matrix[(i, k)] * work1.conjugate() + - matrix[(i, k + 1)] * work2.conjugate(); + } + + matrix[(j, k)] = work1; + matrix[(j, k + 1)] = work2; + + // Keep the Hermitian diagonal explicitly real. + matrix[(j, j)] = T::from_real(matrix[(j, j)].real()); + } + } + + pivots[k] = (pivot_index, 2); + pivots[k + 1] = (pivot_index, 2); + } + + k += block_size; + } + + Self { + matrix, + pivots, + zero_pivot, + } + } + + /// Returns the permutation-aware factor P^T L. + /// + /// This factor can be combined directly with `d()` to reconstruct the original + /// matrix. In general `P^T L` is not lower triangular, even though `L` itself is + /// unit lower triangular in the permuted basis. + pub fn l_permuted(&self) -> OMatrix { + let n = self.matrix.nrows(); + let (nrows, ncols) = self.matrix.shape_generic(); + let mut l_permuted = OMatrix::identity_generic(nrows, ncols); + + let mut k = 0; + while k < n { + let (pivot_index, block_size) = self.pivots[k]; + + if block_size == 1 { + // Right-multiply by the permutation: swap the affected columns. + l_permuted.swap_columns(k, pivot_index); + + // Right-multiply by the unit-lower factor for this 1x1 step. + for row in 0..n { + for i in (k + 1)..n { + l_permuted[(row, k)] = + l_permuted[(row, k)] + l_permuted[(row, i)] * self.matrix[(i, k)]; + } + } + + k += 1; + } else { + // Right-multiply by the permutation: swap the affected columns. + l_permuted.swap_columns(k + 1, pivot_index); + + // Right-multiply by the unit-lower factor for this 2x2 step. + for row in 0..n { + for i in (k + 2)..n { + l_permuted[(row, k)] = + l_permuted[(row, k)] + l_permuted[(row, i)] * self.matrix[(i, k)]; + l_permuted[(row, k + 1)] = l_permuted[(row, k + 1)] + + l_permuted[(row, i)] * self.matrix[(i, k + 1)]; + } + } + + k += 2; + } + } + + l_permuted + } + + /// The block diagonal matrix of this decomposition. + pub fn d(&self) -> OMatrix { + let n = self.matrix.nrows(); + let (nrows, ncols) = self.matrix.shape_generic(); + let mut d = OMatrix::zeros_generic(nrows, ncols); + + let mut k = 0; + while k < n { + d[(k, k)] = self.matrix[(k, k)]; + if self.pivots[k].1 == 2 { + d[(k + 1, k)] = self.matrix[(k + 1, k)]; + d[(k, k + 1)] = self.matrix[(k + 1, k)].conjugate(); + d[(k + 1, k + 1)] = self.matrix[(k + 1, k + 1)]; + k += 1; + } + k += 1; + } + + d + } + + /// Solves the linear system A * x = b using this factorization. + pub fn solve(&self, b: &Matrix) -> Option> + where + S: Storage, + DefaultAllocator: Allocator, + { + let mut result = b.clone_owned(); + + if self.solve_mut(&mut result) { + Some(result) + } else { + None + } + } + + /// Solves the linear system A * x = b in place, overwriting `b` with the solution. + pub fn solve_mut(&self, b: &mut OMatrix) -> bool + where + DefaultAllocator: Allocator, + { + assert_eq!(self.matrix.nrows(), b.nrows()); + + if self.zero_pivot.is_some() { + return false; + } + + let (n, m) = b.shape(); + + // Solve L * y = P^T * b using the stored pivot sequence and multipliers. + let mut k = 0; + while k < n { + let (pivot_index, block_size) = self.pivots[k]; + + if block_size == 1 { + b.swap_rows(k, pivot_index); + + for j in 0..m { + for i in (k + 1)..n { + b[(i, j)] = b[(i, j)] - self.matrix[(i, k)] * b[(k, j)]; + } + } + + k += 1; + } else { + b.swap_rows(k + 1, pivot_index); + + for j in 0..m { + for i in (k + 2)..n { + b[(i, j)] = b[(i, j)] + - self.matrix[(i, k)] * b[(k, j)] + - self.matrix[(i, k + 1)] * b[(k + 1, j)]; + } + } + + k += 2; + } + } + + // Solve D * z = y, handling 1x1 and 2x2 diagonal blocks. + let mut k = 0; + while k < n { + if self.pivots[k].1 == 1 { + for j in 0..m { + b[(k, j)] = b[(k, j)].unscale(self.matrix[(k, k)].real()); + } + k += 1; + } else { + let d11 = self.matrix[(k, k)].real(); + let d22 = self.matrix[(k + 1, k + 1)].real(); + let d21 = self.matrix[(k + 1, k)]; + + let det = d11 * d22 - d21.modulus_squared(); + + for j in 0..m { + let b_k = b[(k, j)]; + let b_k1 = b[(k + 1, j)]; + + b[(k, j)] = (b_k.scale(d22) - b_k1 * d21.conjugate()).unscale(det); + b[(k + 1, j)] = (b_k1.scale(d11) - b_k * d21).unscale(det); + } + k += 2; + } + } + + // Solve L^H * x = z, undoing the pivot sequence in reverse order. + let mut k = n; + while k > 0 { + let k1 = k - 1; + + for j in 0..m { + for i in k..n { + b[(k1, j)] = b[(k1, j)] - self.matrix[(i, k1)].conjugate() * b[(i, j)]; + } + } + + if self.pivots[k1].1 == 1 { + k -= 1; + } else { + let k2 = k - 2; + for j in 0..m { + for i in k..n { + b[(k2, j)] = b[(k2, j)] - self.matrix[(i, k2)].conjugate() * b[(i, j)]; + } + } + k -= 2; + } + + b.swap_rows(k1, self.pivots[k1].0); + } + + true + } + + /// Computes the determinant of the decomposed matrix. + pub fn determinant(&self) -> T::RealField { + let n = self.matrix.nrows(); + let mut determinant = T::RealField::one(); + + let mut k = 0; + while k < n { + if self.pivots[k].1 == 1 { + determinant *= self.matrix[(k, k)].real(); + k += 1; + } else { + determinant *= self.matrix[(k, k)].real() * self.matrix[(k + 1, k + 1)].real() + - self.matrix[(k + 1, k)].modulus_squared(); + k += 2; + } + } + + determinant + } +} + +#[cfg(test)] +mod tests { + use crate::{DMatrix, DVector}; + + use super::*; + + #[test] + fn zero_matrix() { + for n in 1..=5 { + let lblt = DMatrix::::zeros(n, n).lblt(); + assert_eq!(lblt.l_permuted(), DMatrix::identity(n, n)); + assert_eq!(lblt.d(), DMatrix::zeros(n, n)); + assert_eq!(lblt.zero_pivot, Some(0)); + assert_eq!(lblt.pivots, DVector::from_fn(n, |i, _| (i, 1))); + assert!(lblt.determinant().is_zero()); + assert!(lblt.solve(&DVector::from_element(n, 1.0)).is_none()); + } + } + + #[test] + fn identity_matrix() { + for n in 1..=5 { + let identity = DMatrix::::identity(n, n); + let lblt = identity.clone().lblt(); + + assert_eq!(lblt.l_permuted(), identity); + assert_eq!(lblt.d(), identity); + assert_eq!(lblt.zero_pivot, None); + assert_eq!(lblt.pivots, DVector::from_fn(n, |i, _| (i, 1))); + assert!(lblt.determinant().is_one()); + } + } + + #[test] + fn exchange_matrix() { + for n in 1..=15 { + let exchange = DMatrix::from_fn(n, n, |i, j| if i + j + 1 == n { 1.0 } else { 0.0 }); + let lblt = exchange.clone().lblt(); + + let mut expected = Vec::with_capacity(n); + let m = (n + 2) / 4; + for r in 0..m { + let pivot = n - 2 * r - 1; + expected.push((pivot, 2)); + expected.push((pivot, 2)); + } + if !n.is_multiple_of(2) { + expected.push((2 * m, 1)); + } + + for r in m..(n / 2) { + let pivot = 2 * r + n % 2 + 1; + expected.push((pivot, 2)); + expected.push((pivot, 2)); + } + + let l_permuted = lblt.l_permuted(); + let reconstruction = &l_permuted * lblt.d() * l_permuted.adjoint(); + + assert_eq!(exchange, reconstruction); + assert_eq!(lblt.pivots.as_slice(), expected); + } + } +} diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs index 60962d9d1..474674c76 100644 --- a/src/linalg/mod.rs +++ b/src/linalg/mod.rs @@ -17,6 +17,7 @@ pub mod givens; mod hessenberg; pub mod householder; mod inverse; +mod lblt; mod lu; mod permutation_sequence; mod pow; @@ -39,6 +40,7 @@ pub use self::cholesky::*; pub use self::col_piv_qr::*; pub use self::full_piv_lu::*; pub use self::hessenberg::*; +pub use self::lblt::*; pub use self::lu::*; pub use self::permutation_sequence::*; pub use self::qr::*; diff --git a/tests/linalg/lblt.rs b/tests/linalg/lblt.rs new file mode 100644 index 000000000..9fa98f586 --- /dev/null +++ b/tests/linalg/lblt.rs @@ -0,0 +1,121 @@ +use approx::assert_relative_eq; +use na::{ + Complex, ComplexField, DMatrix, DVector, DefaultAllocator, Dim, DimMin, Dyn, LBLT, OMatrix, + OVector, RealField, allocator::Allocator, +}; +use num_traits::FromPrimitive; +use rand::{Rng, SeedableRng, rngs::StdRng}; +use rand_distr::{Distribution, StandardNormal}; + +fn decomposition_error(matrix: &OMatrix, lblt: &LBLT) -> T::RealField +where + T: Copy + ComplexField, + T::RealField: Copy, + N: Dim, + DefaultAllocator: Allocator + Allocator, +{ + let n = T::RealField::from_usize(matrix.nrows()).unwrap(); + let l_permuted = lblt.l_permuted(); + let reconstruction = &l_permuted * lblt.d() * l_permuted.adjoint(); + (reconstruction - matrix).one_norm() / (matrix.one_norm() * n) +} + +/// Samples a Haar-random isometry using a QR factorization of a complex Gaussian matrix. +fn random_isometry(n: N, m: M, rng: &mut impl Rng) -> OMatrix, N, M> +where + T: Copy + RealField, + N: DimMin, + M: Dim, + DefaultAllocator: Allocator + Allocator + Allocator + Allocator, + StandardNormal: Distribution, +{ + let gaussian = OMatrix::from_fn_generic(n, m, |_, _| { + Complex::new(rng.sample(StandardNormal), rng.sample(StandardNormal)) + }); + + let qr = gaussian.qr(); + let mut q = qr.q(); + let r = qr.r(); + + for j in 0..m.value() { + let phase = r[(j, j)] / r[(j, j)].abs(); + for el in q.column_mut(j) { + *el *= phase; + } + } + q +} + +/// Sample a random Hermitian matrix with the given real diagonal spectrum. +fn random_hermitian_from_diag( + diag: &OVector, + rng: &mut impl Rng, +) -> OMatrix, N, N> +where + T: Copy + RealField, + N: DimMin, + DefaultAllocator: Allocator + Allocator, + StandardNormal: Distribution, +{ + let n = diag.shape_generic().0; + let u = random_isometry(n, n, rng); + let d = OMatrix::from_diagonal(&diag.map(|x| Complex::new(x, T::zero()))); + &u * d * u.adjoint() +} + +// This indefinite spectrum reliably exercises diverse 1x1/2x2 pivot patterns across sizes. +#[test] +fn alternating_unit_spectrum() { + let mut rng = StdRng::seed_from_u64(0); + + for n in 2..=10 { + let n_i32 = i32::try_from(n).unwrap(); + let diag = DVector::from_iterator(n, (0..n_i32).map(|i| (-1.0).powi(i))); + + for _ in 0..10 { + let matrix: DMatrix> = random_hermitian_from_diag(&diag, &mut rng); + let lblt = matrix.clone().lblt(); + assert!(decomposition_error(&matrix, &lblt) < 1e-14); + + assert_relative_eq!( + lblt.determinant(), + (-1.0).powi(n_i32 / 2), + max_relative = 1e-12 + ); + + let b = random_isometry(Dyn(n), Dyn(n), &mut rng); + assert!( + (&matrix * &lblt.solve(&b).unwrap() - &b).one_norm() / matrix.one_norm() < 1e-12 + ); + } + } +} + +// An alternating-sign geometric spectrum combines extreme scaling with indefiniteness. +// After Haar randomization this reliably triggers a wide variety of 1×1 and 2×2 pivots. +#[test] +fn alternating_geometric_spectrum() { + let mut rng = StdRng::seed_from_u64(0); + + for n in 2..=10 { + let n_i32 = i32::try_from(n).unwrap(); + let diag = DVector::from_iterator(n, (0..n_i32).map(|i| (-10.0).powi(i))); + + for _ in 0..10 { + let matrix: DMatrix> = random_hermitian_from_diag(&diag, &mut rng); + let lblt = matrix.clone().lblt(); + assert!(decomposition_error(&matrix, &lblt) < 1e-14); + + assert_relative_eq!( + lblt.determinant(), + (-10.0).powi(n_i32 * (n_i32 - 1) / 2), + max_relative = 10.0.powi(-11 + (n_i32 / 2)) + ); + + let b = random_isometry(Dyn(n), Dyn(n), &mut rng); + assert!( + (&matrix * &lblt.solve(&b).unwrap() - &b).one_norm() / matrix.one_norm() < 1e-12 + ); + } + } +} diff --git a/tests/linalg/mod.rs b/tests/linalg/mod.rs index d9bd6cd91..e6330510a 100644 --- a/tests/linalg/mod.rs +++ b/tests/linalg/mod.rs @@ -8,6 +8,7 @@ mod exp; mod full_piv_lu; mod hessenberg; mod inverse; +mod lblt; mod lu; mod pow; mod qr;