From 6831027573ff9eebcb1bf4e0385c25b1718ba25f Mon Sep 17 00:00:00 2001 From: Andrew Cassidy Date: Tue, 9 May 2023 22:27:05 -0700 Subject: [PATCH] Document and refine `decompose` module --- src/decompose.rs | 269 ++++++++++++++++++++++++++++++++++----------- tests/decompose.rs | 3 +- 2 files changed, 209 insertions(+), 63 deletions(-) diff --git a/src/decompose.rs b/src/decompose.rs index 2cf12c0..dc12979 100644 --- a/src/decompose.rs +++ b/src/decompose.rs @@ -1,79 +1,93 @@ use crate::util::checked_inv; -use crate::Matrix; -use crate::Vector; +use crate::{Matrix, Vector}; use num_traits::real::Real; +use num_traits::One; use std::iter::{Product, Sum}; +use std::ops::{Mul, Neg, Not}; +/// The parity of an [LU decomposition](LUDecomposition). In other words, how many times the +/// source matrix has to have rows swapped before the decomposition takes place #[derive(Copy, Clone, Debug, PartialEq)] -pub struct LUDecomposition { - pub lu: Matrix, - pub idx: Vector, - pub parity: T, +pub enum Parity { + Even, + Odd, } -impl LUDecomposition +impl Mul for Parity where - T: Real + Default + Sum + Product, + T: Neg + One, { - #[must_use] - pub fn decompose(m: &Matrix) -> Option { - // Implementation from Numerical Recipes §2.3 - let mut lu = m.clone(); - let mut idx: Vector = (0..N).collect(); - let mut parity = T::one(); + type Output = T; - let mut vv: Vector = m - .rows() - .map(|row| { - let m = row.elements().cloned().reduce(|acc, x| acc.max(x.abs()))?; - match m < T::epsilon() { - true => None, - false => Some(T::one() / m), - } - }) - .collect::>()?; // get the inverse maxabs value in each row + fn mul(self, rhs: T) -> Self::Output { + rhs * match self { + Parity::Even => T::one(), + Parity::Odd => -T::one(), + } + } +} - for k in 0..N { - // search for the pivot element and its index - let (ipivot, _) = (lu.col(k) * vv) - .abs() - .elements() - .enumerate() - .skip(k) // below the diagonal - .reduce(|(imax, xmax), (i, x)| match x > xmax { - // Is the figure of merit for the pivot better than the best so far? - true => (i, x), - false => (imax, xmax), - })?; +impl Not for Parity { + type Output = Parity; - // do we need to interchange rows? - if k != ipivot { - lu.pivot_row(ipivot, k); // yes, we do - idx.pivot_row(ipivot, k); - parity = -parity; // change parity of d - vv[ipivot] = vv[k] //interchange scale factor - } + fn not(self) -> Self::Output { + match self { + Parity::Even => Parity::Odd, + Parity::Odd => Parity::Even, + } + } +} - let pivot = lu[(k, k)]; - if pivot.abs() < T::epsilon() { - // if the pivot is zero, the matrix is singular - return None; - }; +/// The result of the [LU decomposition](LUDecomposable::lu) of a matrix. +/// +/// This struct provides a convenient way to reuse one LU decomposition to solve multiple +/// matrix equations. You likely do not need to worry about its contents. +/// +/// See [LU decomposition](https://en.wikipedia.org/wiki/LU_decomposition) +/// on wikipedia for more information +#[derive(Copy, Clone, Debug, PartialEq)] +pub struct LUDecomposition { + /// The $L$ and $U$ matrices combined into one + /// + /// for example if + /// + /// $ U = [[u_{11}, u_{12}, cdots, u_{1n} ], + /// [0, u_{22}, cdots, u_{2n} ], + /// [vdots, vdots, ddots, vdots ], + /// [0, 0, cdots, u_{mn} ]] $ + /// and + /// $ L = [[1, 0, cdots, 0 ], + /// [l_{21}, 1, cdots, 0 ], + /// [vdots, vdots, ddots, vdots ], + /// [l_{m1}, l_{m2}, cdots, 1 ]] $, + /// then + /// $ LU = [[u_{11}, u_{12}, cdots, u_{1n} ], + /// [l_{21}, u_{22}, cdots, u_{2n} ], + /// [vdots, vdots, ddots, vdots ], + /// [l_{m1}, l_{m2}, cdots, u_{mn} ]] $ + /// + /// note that the diagonals of the $L$ matrix are always 1, so no information is lost + pub lu: Matrix, - for i in (k + 1)..N { - // divide by the pivot element - let dpivot = lu[(i, k)] / pivot; - lu[(i, k)] = dpivot; - for j in (k + 1)..N { - // reduce remaining submatrix - lu[(i, j)] = lu[(i, j)] - (dpivot * lu[(k, j)]); - } - } - } + /// The indices of the permutation matrix $P$, such that $PxxA$ = $LxxU$ + /// + /// The permutation matrix rearranges the rows of the original matrix in order to produce + /// the LU decomposition. This makes calculation simpler, but makes the result + /// (known as an LUP decomposition) no longer unique + pub idx: Vector, - return Some(Self { lu, idx, parity }); - } + /// The parity of the decomposition. + pub parity: Parity, +} +impl LUDecomposition +where + T: Real + Default + Sum + Product, +{ + /// Solve for $x$ in $M xx x = b$, where $M$ is the original matrix this is a decomposition of. + /// + /// This is equivalent to [`LUDecomposable::solve`] while allowing the LU decomposition + /// to be reused #[must_use] pub fn solve(&self, b: &Matrix) -> Matrix { let b_permuted = b.permute_rows(&self.idx); @@ -81,7 +95,7 @@ where Matrix::from_cols(b_permuted.cols().map(|mut x| { // Implementation from Numerical Recipes §2.3 // When ii is set to a positive value, - // it will become the index of the first nonvanishing element of b + // it will become the index of the first non-vanishing element of b let mut ii = 0usize; for i in 0..N { // forward substitution using L @@ -107,14 +121,25 @@ where })) } + /// Calculate the determinant $|M|$ of the matrix $M$. + /// If the matrix is singular, the determinant is 0. + /// + /// This is equivalent to [`LUDecomposable::det`] while allowing the LU decomposition + /// to be reused pub fn det(&self) -> T { self.parity * self.lu.diagonals().product() } + /// Calculate the inverse of the original matrix, such that $MxxM^{-1} = I$ + /// + /// This is equivalent to [`Matrix::inverse`] while allowing the LU decomposition to be reused + #[must_use] pub fn inverse(&self) -> Matrix { return self.solve(&Matrix::::identity()); } + /// Separate the $L$ and $U$ sides of the $LU$ matrix. + /// See [the `lu` field](LUDecomposition::lu) for more information pub fn separate(&self) -> (Matrix, Matrix) { let mut l = Matrix::::identity(); let mut u = self.lu; // lu @@ -131,19 +156,83 @@ where } } +/// A Matrix that can be decomposed into an upper and lower diagonal matrix, +/// known as an [LU Decomposition](LUDecomposition). +/// +/// See [LU decomposition](https://en.wikipedia.org/wiki/LU_decomposition) +/// on wikipedia for more information pub trait LUDecomposable where T: Copy + Default + Real + Product + Sum, { + /// return this matrix's [`LUDecomposition`], or [`None`] if the matrix is singular. + /// This can be used to solve for multiple results + /// + /// ``` + /// # use vector_victor::decompose::LUDecomposable; + /// # use vector_victor::{Matrix, Vector}; + /// let m = Matrix::new([[1.0,3.0],[2.0,4.0]]); + /// let lu = m.lu().expect("Cannot decompose a signular matrix"); + /// + /// let b = Vector::vec([7.0,10.0]); + /// assert_eq!(lu.solve(&b), Vector::vec([1.0,2.0])); + /// + /// let c = Vector::vec([10.0, 14.0]); + /// assert_eq!(lu.solve(&c), Vector::vec([1.0,3.0])); + /// + /// ``` #[must_use] fn lu(&self) -> Option>; + /// Calculate the inverse of the matrix, such that $MxxM^{-1} = I$, or [`None`] if the matrix is singular. + /// + /// ``` + /// # use vector_victor::decompose::LUDecomposable; + /// # use vector_victor::Matrix; + /// let m = Matrix::new([[1.0,3.0],[2.0,4.0]]); + /// let mi = m.inverse().expect("Cannot invert a singular matrix"); + /// + /// assert_eq!(mi, Matrix::new([[-2.0, 1.5],[1.0, -0.5]]), "unexpected inverse matrix"); + /// + /// // multiplying a matrix by its inverse yields the identity matrix + /// assert_eq!(m.mmul(&mi), Matrix::identity()) + /// ``` #[must_use] fn inverse(&self) -> Option>; + /// Calculate the determinant $|M|$ of the matrix $M$. + /// If the matrix is singular, the determinant is 0 #[must_use] fn det(&self) -> T; + /// Solve for $x$ in $M xx x = b$, where $M$ is the original matrix this is a decomposition of. + /// + /// ``` + /// # use vector_victor::decompose::LUDecomposable; + /// # use vector_victor::{Matrix, Vector}; + /// + /// let m = Matrix::new([[1.0,3.0],[2.0,4.0]]); + /// let b = Vector::vec([7.0,10.0]); + /// let x = m.solve(&b).expect("Cannot solve a singular matrix"); + /// + /// assert_eq!(x, Vector::vec([1.0,2.0]), "x = [1,2]"); + /// assert_eq!(m.mmul(&x), b, "Mx = b"); + /// ``` + /// + /// $x$ does not need to be a column-vector, it can also be a 2D matrix. For example, + /// the following is another way to calculate the [`inverse`] by solving for the identity matrix $I$. + /// + /// ``` + /// # use vector_victor::decompose::LUDecomposable; + /// # use vector_victor::{Matrix, Vector}; + /// + /// let m = Matrix::new([[1.0,3.0],[2.0,4.0]]); + /// let i = Matrix::::identity(); + /// let mi = m.solve(&i).expect("Cannot solve a singular matrix"); + /// + /// assert_eq!(mi, Matrix::new([[-2.0, 1.5],[1.0, -0.5]])); + /// assert_eq!(m.mmul(&mi), i, "M x M^-1 = I"); + /// ``` #[must_use] fn solve(&self, b: &Matrix) -> Option>; } @@ -153,7 +242,63 @@ where T: Copy + Default + Real + Sum + Product, { fn lu(&self) -> Option> { - LUDecomposition::decompose(self) + // Implementation from Numerical Recipes §2.3 + let mut lu = self.clone(); + let mut idx: Vector = (0..N).collect(); + let mut parity = Parity::Even; + + let mut vv: Vector = self + .rows() + .map(|row| { + let m = row.elements().cloned().reduce(|acc, x| acc.max(x.abs()))?; + checked_inv(m) + }) + .collect::>()?; // get the inverse max abs value in each row + + // for each column in the matrix... + for k in 0..N { + // search for the pivot element and its index + let (ipivot, _) = (lu.col(k) * vv) + .abs() + .elements() + .enumerate() + .skip(k) // below the diagonal + .reduce(|(imax, xmax), (i, x)| match x > xmax { + // Is the figure of merit for the pivot better than the best so far? + true => (i, x), + false => (imax, xmax), + })?; + + // do we need to interchange rows? + if k != ipivot { + lu.pivot_row(ipivot, k); // yes, we do + idx.pivot_row(ipivot, k); + parity = !parity; // swap parity + vv[ipivot] = vv[k] // interchange scale factor + } + + // select our pivot, which is now on the diagonal + let pivot = lu[(k, k)]; + if pivot.abs() < T::epsilon() { + // if the pivot is zero, the matrix is singular + return None; + }; + + // for each element in the column k below the diagonal... + // this is called outer product Gaussian elimination + for i in (k + 1)..N { + // divide by the pivot element + lu[(i, k)] = lu[(i, k)] / pivot; + + // for each element in the column k below the diagonal... + for j in (k + 1)..N { + // reduce remaining submatrix + lu[(i, j)] = lu[(i, j)] - (lu[(i, k)] * lu[(k, j)]); + } + } + } + + return Some(LUDecomposition { lu, idx, parity }); } fn inverse(&self) -> Option> { diff --git a/tests/decompose.rs b/tests/decompose.rs index e159e50..5daa705 100644 --- a/tests/decompose.rs +++ b/tests/decompose.rs @@ -7,6 +7,7 @@ use num_traits::real::Real; use num_traits::Zero; use std::fmt::Debug; use std::iter::{Product, Sum}; +use vector_victor::decompose::Parity::Even; use vector_victor::decompose::{LUDecomposable, LUDecomposition}; use vector_victor::{Matrix, Vector}; @@ -25,7 +26,7 @@ fn test_lu_identity