From 45f2fc4f9297d3665f780f2516a59c23cd1f8888 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Crozet=20S=C3=A9bastien?= Date: Fri, 13 Nov 2020 17:26:47 +0100 Subject: [PATCH] Move all matrix decomposition methods under a single impl. --- src/base/matrix.rs | 15 ++ src/linalg/bidiagonal.rs | 17 +- src/linalg/cholesky.rs | 15 +- src/linalg/decomposition.rs | 232 ++++++++++++++++++++++++++++ src/linalg/full_piv_lu.rs | 12 -- src/linalg/hessenberg.rs | 12 +- src/linalg/lu.rs | 10 -- src/linalg/mod.rs | 1 + src/linalg/qr.rs | 10 -- src/linalg/schur.rs | 20 --- src/linalg/svd.rs | 25 --- src/linalg/symmetric_eigen.rs | 26 ---- src/linalg/symmetric_tridiagonal.rs | 14 +- 13 files changed, 252 insertions(+), 157 deletions(-) create mode 100644 src/linalg/decomposition.rs diff --git a/src/base/matrix.rs b/src/base/matrix.rs index ca7fa77b..2ba91826 100644 --- a/src/base/matrix.rs +++ b/src/base/matrix.rs @@ -75,6 +75,21 @@ pub type MatrixCross = /// Note that mixing `Dynamic` with type-level unsigned integers is allowed. Actually, a /// dynamically-sized column vector should be represented as a `Matrix` (given /// some concrete types for `N` and a compatible data storage type `S`). +/// +/// # Documentation by feature +/// Because `Matrix` is the most generic types that groups all matrix and vectors of **nalgebra** +/// this documentation page contains every single matrix/vector-related method. In order to make +/// browsing this page simpler, the next subsections contain direct links to groups of methods +/// related to a specific topic. +/// +/// #### Matrix decomposition +/// - [Rectangular matrix decomposition](#rectangular-matrix-decomposition). +/// - [Square matrix decomposition](#square-matrix-decomposition). +/// +/// #### Matrix slicing +/// - [Slicing](#slicing) +/// - [Mutable slicing](#mutable-slicing) +/// - [Range-based slicing](#range-based-slicing), [mutable range-based slicing](#mutable-range-based-slicing). #[repr(C)] #[derive(Clone, Copy)] pub struct Matrix { diff --git a/src/linalg/bidiagonal.rs b/src/linalg/bidiagonal.rs index c19a0662..3a789978 100644 --- a/src/linalg/bidiagonal.rs +++ b/src/linalg/bidiagonal.rs @@ -40,7 +40,7 @@ where + Allocator> + Allocator, U1>>, { - // FIXME: perhaps we should pack the axises into different vectors so that axises for `v_t` are + // FIXME: perhaps we should pack the axes into different vectors so that axes for `v_t` are // contiguous. This prevents some useless copies. uv: MatrixMN, /// The diagonal elements of the decomposed matrix. @@ -359,18 +359,3 @@ where // // res self.q_determinant() // // } // } - -impl, C: Dim, S: Storage> Matrix -where - DimMinimum: DimSub, - DefaultAllocator: Allocator - + Allocator - + Allocator - + Allocator> - + Allocator, U1>>, -{ - /// Computes the bidiagonalization using householder reflections. - pub fn bidiagonalize(self) -> Bidiagonal { - Bidiagonal::new(self.into_owned()) - } -} diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 17187fdf..bd2f9281 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -6,7 +6,7 @@ use simba::scalar::ComplexField; use simba::simd::SimdComplexField; use crate::allocator::Allocator; -use crate::base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, SquareMatrix, Vector}; +use crate::base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, Vector}; use crate::constraint::{SameNumberOfRows, ShapeConstraint}; use crate::dimension::{Dim, DimAdd, DimDiff, DimSub, DimSum, U1}; use crate::storage::{Storage, StorageMut}; @@ -363,16 +363,3 @@ where } } } - -impl> SquareMatrix -where - DefaultAllocator: Allocator, -{ - /// Attempts to compute the Cholesky decomposition of this matrix. - /// - /// Returns `None` if the input matrix is not definite-positive. The input matrix is assumed - /// to be symmetric and only the lower-triangular part is read. - pub fn cholesky(self) -> Option> { - Cholesky::new(self.into_owned()) - } -} diff --git a/src/linalg/decomposition.rs b/src/linalg/decomposition.rs new file mode 100644 index 00000000..67cc4c6a --- /dev/null +++ b/src/linalg/decomposition.rs @@ -0,0 +1,232 @@ +use crate::storage::Storage; +use crate::{ + Allocator, Bidiagonal, Cholesky, ComplexField, DefaultAllocator, Dim, DimDiff, DimMin, + DimMinimum, DimSub, FullPivLU, Hessenberg, Matrix, Schur, SymmetricEigen, SymmetricTridiagonal, + LU, QR, SVD, U1, +}; + +/// # Rectangular matrix decomposition +/// +/// This section contains the methods for computing some common decompositions of rectangular +/// matrices with real or complex components. The following are currently supported: +/// +/// | Decomposition | Factors | Details | +/// | -------------------------|---------------------|--------------| +/// | QR | `Q * R` | `Q` is an unitary matrix, and `R` is upper-triangular. | +/// | LU with partial pivoting | `P⁻¹ * L * U` | `L` is lower-triangular with a diagonal filled with `1` and `U` is upper-triangular. `P` is a permutation matrix. | +/// | LU with full pivoting | `P⁻¹ * L * U ~ Q⁻¹` | `L` is lower-triangular with a diagonal filled with `1` and `U` is upper-triangular. `P` and `Q` are permutation matrices. | +/// | SVD | `U * Σ * Vᵀ` | `U` and `V` are two orthogonal matrices and `Σ` is a diagonal matrix containing the singular values. | +impl> Matrix { + /// Computes the bidiagonalization using householder reflections. + pub fn bidiagonalize(self) -> Bidiagonal + where + R: DimMin, + DimMinimum: DimSub, + DefaultAllocator: Allocator + + Allocator + + Allocator + + Allocator> + + Allocator, U1>>, + { + Bidiagonal::new(self.into_owned()) + } + + /// Computes the LU decomposition with full pivoting of `matrix`. + /// + /// This effectively computes `P, L, U, Q` such that `P * matrix * Q = LU`. + pub fn full_piv_lu(self) -> FullPivLU + where + R: DimMin, + DefaultAllocator: Allocator + Allocator<(usize, usize), DimMinimum>, + { + FullPivLU::new(self.into_owned()) + } + + /// Computes the LU decomposition with partial (row) pivoting of `matrix`. + pub fn lu(self) -> LU + where + R: DimMin, + DefaultAllocator: Allocator + Allocator<(usize, usize), DimMinimum>, + { + LU::new(self.into_owned()) + } + + /// Computes the QR decomposition of this matrix. + pub fn qr(self) -> QR + where + R: DimMin, + DefaultAllocator: Allocator + Allocator + Allocator>, + { + QR::new(self.into_owned()) + } + + /// Computes the Singular Value Decomposition using implicit shift. + pub fn svd(self, compute_u: bool, compute_v: bool) -> SVD + where + R: DimMin, + DimMinimum: DimSub, // for Bidiagonal. + DefaultAllocator: Allocator + + Allocator + + Allocator + + Allocator, U1>> + + Allocator, C> + + Allocator> + + Allocator> + + Allocator> + + Allocator, U1>>, + { + SVD::new(self.into_owned(), compute_u, compute_v) + } + + /// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift. + /// + /// # Arguments + /// + /// * `compute_u` − set this to `true` to enable the computation of left-singular vectors. + /// * `compute_v` − set this to `true` to enable the computation of right-singular vectors. + /// * `eps` − tolerance used to determine when a value converged to 0. + /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this + /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm + /// continues indefinitely until convergence. + pub fn try_svd( + self, + compute_u: bool, + compute_v: bool, + eps: N::RealField, + max_niter: usize, + ) -> Option> + where + R: DimMin, + DimMinimum: DimSub, // for Bidiagonal. + DefaultAllocator: Allocator + + Allocator + + Allocator + + Allocator, U1>> + + Allocator, C> + + Allocator> + + Allocator> + + Allocator> + + Allocator, U1>>, + { + SVD::try_new(self.into_owned(), compute_u, compute_v, eps, max_niter) + } +} + +/// # Square matrix decomposition +/// +/// This section contains the methods for computing some common decompositions of square +/// matrices with real or complex components. The following are currently supported: +/// +/// | Decomposition | Factors | Details | +/// | -------------------------|---------------------------|--------------| +/// | Hessenberg | `Q * H * Qᵀ` | `Q` is a unitary matrix and `H` an upper-Hessenberg matrix. | +/// | Cholesky | `L * Lᵀ` | `L` is a lower-triangular 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. | +/// | Symmetric tridiagonalization | `Q ~ T ~ Qᵀ` | `Q` is an unitary matrix, and `T` is a tridiagonal matrix. | +impl> Matrix { + /// Attempts to compute the Cholesky decomposition of this matrix. + /// + /// Returns `None` if the input matrix is not definite-positive. The input matrix is assumed + /// to be symmetric and only the lower-triangular part is read. + pub fn cholesky(self) -> Option> + where + DefaultAllocator: Allocator, + { + Cholesky::new(self.into_owned()) + } + + /// Computes the Hessenberg decomposition of this matrix using householder reflections. + pub fn hessenberg(self) -> Hessenberg + where + D: DimSub, + DefaultAllocator: Allocator + Allocator + Allocator>, + { + Hessenberg::new(self.into_owned()) + } + + /// Computes the Schur decomposition of a square matrix. + pub fn schur(self) -> Schur + where + D: DimSub, // For Hessenberg. + DefaultAllocator: Allocator> + + Allocator> + + Allocator + + Allocator, + { + Schur::new(self.into_owned()) + } + + /// Attempts to compute the Schur decomposition of a square matrix. + /// + /// If only eigenvalues are needed, it is more efficient to call the matrix method + /// `.eigenvalues()` instead. + /// + /// # Arguments + /// + /// * `eps` − tolerance used to determine when a value converged to 0. + /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this + /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm + /// continues indefinitely until convergence. + pub fn try_schur(self, eps: N::RealField, max_niter: usize) -> Option> + where + D: DimSub, // For Hessenberg. + DefaultAllocator: Allocator> + + Allocator> + + Allocator + + Allocator, + { + Schur::try_new(self.into_owned(), eps, max_niter) + } + + /// Computes the eigendecomposition of this symmetric matrix. + /// + /// Only the lower-triangular part (including the diagonal) of `m` is read. + pub fn symmetric_eigen(self) -> SymmetricEigen + where + D: DimSub, + DefaultAllocator: Allocator + + Allocator> + + Allocator + + Allocator>, + { + SymmetricEigen::new(self.into_owned()) + } + + /// Computes the eigendecomposition of the given symmetric matrix with user-specified + /// convergence parameters. + /// + /// Only the lower-triangular part (including the diagonal) of `m` is read. + /// + /// # Arguments + /// + /// * `eps` − tolerance used to determine when a value converged to 0. + /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this + /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm + /// continues indefinitely until convergence. + pub fn try_symmetric_eigen( + self, + eps: N::RealField, + max_niter: usize, + ) -> Option> + where + D: DimSub, + DefaultAllocator: Allocator + + Allocator> + + Allocator + + Allocator>, + { + SymmetricEigen::try_new(self.into_owned(), eps, max_niter) + } + + /// Computes the tridiagonalization of this symmetric matrix. + /// + /// Only the lower-triangular part (including the diagonal) of `m` is read. + pub fn symmetric_tridiagonalize(self) -> SymmetricTridiagonal + where + D: DimSub, + DefaultAllocator: Allocator + Allocator>, + { + SymmetricTridiagonal::new(self.into_owned()) + } +} diff --git a/src/linalg/full_piv_lu.rs b/src/linalg/full_piv_lu.rs index 146b2cb7..fb8ff670 100644 --- a/src/linalg/full_piv_lu.rs +++ b/src/linalg/full_piv_lu.rs @@ -257,15 +257,3 @@ where } } } - -impl, C: Dim, S: Storage> Matrix -where - DefaultAllocator: Allocator + Allocator<(usize, usize), DimMinimum>, -{ - /// Computes the LU decomposition with full pivoting of `matrix`. - /// - /// This effectively computes `P, L, U, Q` such that `P * matrix * Q = LU`. - pub fn full_piv_lu(self) -> FullPivLU { - FullPivLU::new(self.into_owned()) - } -} diff --git a/src/linalg/hessenberg.rs b/src/linalg/hessenberg.rs index 57a24097..1e795b41 100644 --- a/src/linalg/hessenberg.rs +++ b/src/linalg/hessenberg.rs @@ -2,7 +2,7 @@ use serde::{Deserialize, Serialize}; use crate::allocator::Allocator; -use crate::base::{DefaultAllocator, MatrixMN, MatrixN, SquareMatrix, VectorN}; +use crate::base::{DefaultAllocator, MatrixMN, MatrixN, VectorN}; use crate::dimension::{DimDiff, DimSub, U1}; use crate::storage::Storage; use simba::scalar::ComplexField; @@ -131,13 +131,3 @@ where &self.hess } } - -impl, S: Storage> SquareMatrix -where - DefaultAllocator: Allocator + Allocator + Allocator>, -{ - /// Computes the Hessenberg decomposition of this matrix using householder reflections. - pub fn hessenberg(self) -> Hessenberg { - Hessenberg::new(self.into_owned()) - } -} diff --git a/src/linalg/lu.rs b/src/linalg/lu.rs index 6aa84d62..9e417f10 100644 --- a/src/linalg/lu.rs +++ b/src/linalg/lu.rs @@ -379,13 +379,3 @@ pub fn gauss_step_swap( .axpy(-pivot_row[k].inlined_clone(), &coeffs, N::one()); } } - -impl, C: Dim, S: Storage> Matrix -where - DefaultAllocator: Allocator + Allocator<(usize, usize), DimMinimum>, -{ - /// Computes the LU decomposition with partial (row) pivoting of `matrix`. - pub fn lu(self) -> LU { - LU::new(self.into_owned()) - } -} diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs index 4a1edbcf..df265498 100644 --- a/src/linalg/mod.rs +++ b/src/linalg/mod.rs @@ -8,6 +8,7 @@ mod determinant; // FIXME: this should not be needed. However, the exp uses // explicit float operations on `f32` and `f64`. We need to // get rid of these to allow exp to be used on a no-std context. +mod decomposition; #[cfg(feature = "std")] mod exp; mod full_piv_lu; diff --git a/src/linalg/qr.rs b/src/linalg/qr.rs index 48d7a896..d3cf3275 100644 --- a/src/linalg/qr.rs +++ b/src/linalg/qr.rs @@ -288,13 +288,3 @@ where // res self.q_determinant() // } } - -impl, C: Dim, S: Storage> Matrix -where - DefaultAllocator: Allocator + Allocator + Allocator>, -{ - /// Computes the QR decomposition of this matrix. - pub fn qr(self) -> QR { - QR::new(self.into_owned()) - } -} diff --git a/src/linalg/schur.rs b/src/linalg/schur.rs index 8d8be72d..1292bad7 100644 --- a/src/linalg/schur.rs +++ b/src/linalg/schur.rs @@ -496,26 +496,6 @@ where + Allocator + Allocator, { - /// Computes the Schur decomposition of a square matrix. - pub fn schur(self) -> Schur { - Schur::new(self.into_owned()) - } - - /// Attempts to compute the Schur decomposition of a square matrix. - /// - /// If only eigenvalues are needed, it is more efficient to call the matrix method - /// `.eigenvalues()` instead. - /// - /// # Arguments - /// - /// * `eps` − tolerance used to determine when a value converged to 0. - /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this - /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm - /// continues indefinitely until convergence. - pub fn try_schur(self, eps: N::RealField, max_niter: usize) -> Option> { - Schur::try_new(self.into_owned(), eps, max_niter) - } - /// Computes the eigenvalues of this matrix. pub fn eigenvalues(&self) -> Option> { assert!( diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs index 94ec4d96..93cba072 100644 --- a/src/linalg/svd.rs +++ b/src/linalg/svd.rs @@ -616,31 +616,6 @@ where + Allocator> + Allocator, U1>>, { - /// Computes the Singular Value Decomposition using implicit shift. - pub fn svd(self, compute_u: bool, compute_v: bool) -> SVD { - SVD::new(self.into_owned(), compute_u, compute_v) - } - - /// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift. - /// - /// # Arguments - /// - /// * `compute_u` − set this to `true` to enable the computation of left-singular vectors. - /// * `compute_v` − set this to `true` to enable the computation of right-singular vectors. - /// * `eps` − tolerance used to determine when a value converged to 0. - /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this - /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm - /// continues indefinitely until convergence. - pub fn try_svd( - self, - compute_u: bool, - compute_v: bool, - eps: N::RealField, - max_niter: usize, - ) -> Option> { - SVD::try_new(self.into_owned(), compute_u, compute_v, eps, max_niter) - } - /// Computes the singular values of this matrix. pub fn singular_values(&self) -> VectorN> { SVD::new(self.clone_owned(), false, false).singular_values diff --git a/src/linalg/symmetric_eigen.rs b/src/linalg/symmetric_eigen.rs index 2d62816e..64d3515a 100644 --- a/src/linalg/symmetric_eigen.rs +++ b/src/linalg/symmetric_eigen.rs @@ -308,32 +308,6 @@ where + Allocator + Allocator>, { - /// Computes the eigendecomposition of this symmetric matrix. - /// - /// Only the lower-triangular part (including the diagonal) of `m` is read. - pub fn symmetric_eigen(self) -> SymmetricEigen { - SymmetricEigen::new(self.into_owned()) - } - - /// Computes the eigendecomposition of the given symmetric matrix with user-specified - /// convergence parameters. - /// - /// Only the lower-triangular part (including the diagonal) of `m` is read. - /// - /// # Arguments - /// - /// * `eps` − tolerance used to determine when a value converged to 0. - /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this - /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm - /// continues indefinitely until convergence. - pub fn try_symmetric_eigen( - self, - eps: N::RealField, - max_niter: usize, - ) -> Option> { - SymmetricEigen::try_new(self.into_owned(), eps, max_niter) - } - /// Computes the eigenvalues of this symmetric matrix. /// /// Only the lower-triangular part of the matrix is read. diff --git a/src/linalg/symmetric_tridiagonal.rs b/src/linalg/symmetric_tridiagonal.rs index f4dc13ba..e8d9fb5d 100644 --- a/src/linalg/symmetric_tridiagonal.rs +++ b/src/linalg/symmetric_tridiagonal.rs @@ -2,7 +2,7 @@ use serde::{Deserialize, Serialize}; use crate::allocator::Allocator; -use crate::base::{DefaultAllocator, MatrixMN, MatrixN, SquareMatrix, VectorN}; +use crate::base::{DefaultAllocator, MatrixMN, MatrixN, VectorN}; use crate::dimension::{DimDiff, DimSub, U1}; use crate::storage::Storage; use simba::scalar::ComplexField; @@ -162,15 +162,3 @@ where &q * self.tri * q.adjoint() } } - -impl, S: Storage> SquareMatrix -where - DefaultAllocator: Allocator + Allocator>, -{ - /// Computes the tridiagonalization of this symmetric matrix. - /// - /// Only the lower-triangular part (including the diagonal) of `m` is read. - pub fn symmetric_tridiagonalize(self) -> SymmetricTridiagonal { - SymmetricTridiagonal::new(self.into_owned()) - } -}