M-Labs/nalgebra
3
0
Fork 1
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Move all matrix decomposition methods under a single impl.

Browse Source
This commit is contained in:
Crozet Sébastien 2020-11-13 17:26:47 +01:00
parent 4f443b06a9
commit 45f2fc4f92
13 changed files with 252 additions and 157 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
src/base/matrix.rs
View File

@ -75,6 +75,21 @@ pub type MatrixCross<N, R1, C1, R2, C2> =
/// Note that mixing `Dynamic` with type-level unsigned integers is allowed. Actually, a
/// dynamically-sized column vector should be represented as a `Matrix<N, Dynamic, U1, S>` (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<N: Scalar, R: Dim, C: Dim, S> {

17
src/linalg/bidiagonal.rs
View File

@ -40,7 +40,7 @@ where
+ Allocator<N, DimMinimum<R, C>>
+ Allocator<N, DimDiff<DimMinimum<R, C>, 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<N, R, C>,
/// The diagonal elements of the decomposed matrix.
@ -359,18 +359,3 @@ where
// // res self.q_determinant()
// // }
// }
impl<N: ComplexField, R: DimMin<C>, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S>
where
DimMinimum<R, C>: DimSub<U1>,
DefaultAllocator: Allocator<N, R, C>
+ Allocator<N, C>
+ Allocator<N, R>
+ Allocator<N, DimMinimum<R, C>>
+ Allocator<N, DimDiff<DimMinimum<R, C>, U1>>,
{
/// Computes the bidiagonalization using householder reflections.
pub fn bidiagonalize(self) -> Bidiagonal<N, R, C> {
Bidiagonal::new(self.into_owned())
}
}

15
src/linalg/cholesky.rs
View File

@ -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<N: ComplexField, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S>
where
DefaultAllocator: Allocator<N, D, D>,
{
/// 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<N, D>> {
Cholesky::new(self.into_owned())
}
}

232
src/linalg/decomposition.rs Normal file
View File

@ -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<N: ComplexField, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
/// Computes the bidiagonalization using householder reflections.
pub fn bidiagonalize(self) -> Bidiagonal<N, R, C>
where
R: DimMin<C>,
DimMinimum<R, C>: DimSub<U1>,
DefaultAllocator: Allocator<N, R, C>
+ Allocator<N, C>
+ Allocator<N, R>
+ Allocator<N, DimMinimum<R, C>>
+ Allocator<N, DimDiff<DimMinimum<R, C>, 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<N, R, C>
where
R: DimMin<C>,
DefaultAllocator: Allocator<N, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
{
FullPivLU::new(self.into_owned())
}
/// Computes the LU decomposition with partial (row) pivoting of `matrix`.
pub fn lu(self) -> LU<N, R, C>
where
R: DimMin<C>,
DefaultAllocator: Allocator<N, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
{
LU::new(self.into_owned())
}
/// Computes the QR decomposition of this matrix.
pub fn qr(self) -> QR<N, R, C>
where
R: DimMin<C>,
DefaultAllocator: Allocator<N, R, C> + Allocator<N, R> + Allocator<N, DimMinimum<R, C>>,
{
QR::new(self.into_owned())
}
/// Computes the Singular Value Decomposition using implicit shift.
pub fn svd(self, compute_u: bool, compute_v: bool) -> SVD<N, R, C>
where
R: DimMin<C>,
DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
DefaultAllocator: Allocator<N, R, C>
+ Allocator<N, C>
+ Allocator<N, R>
+ Allocator<N, DimDiff<DimMinimum<R, C>, U1>>
+ Allocator<N, DimMinimum<R, C>, C>
+ Allocator<N, R, DimMinimum<R, C>>
+ Allocator<N, DimMinimum<R, C>>
+ Allocator<N::RealField, DimMinimum<R, C>>
+ Allocator<N::RealField, DimDiff<DimMinimum<R, C>, 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<SVD<N, R, C>>
where
R: DimMin<C>,
DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
DefaultAllocator: Allocator<N, R, C>
+ Allocator<N, C>
+ Allocator<N, R>
+ Allocator<N, DimDiff<DimMinimum<R, C>, U1>>
+ Allocator<N, DimMinimum<R, C>, C>
+ Allocator<N, R, DimMinimum<R, C>>
+ Allocator<N, DimMinimum<R, C>>
+ Allocator<N::RealField, DimMinimum<R, C>>
+ Allocator<N::RealField, DimDiff<DimMinimum<R, C>, 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<N: ComplexField, D: Dim, S: Storage<N, D, D>> Matrix<N, D, D, S> {
/// 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<N, D>>
where
DefaultAllocator: Allocator<N, D, D>,
{
Cholesky::new(self.into_owned())
}
/// Computes the Hessenberg decomposition of this matrix using householder reflections.
pub fn hessenberg(self) -> Hessenberg<N, D>
where
D: DimSub<U1>,
DefaultAllocator: Allocator<N, D, D> + Allocator<N, D> + Allocator<N, DimDiff<D, U1>>,
{
Hessenberg::new(self.into_owned())
}
/// Computes the Schur decomposition of a square matrix.
pub fn schur(self) -> Schur<N, D>
where
D: DimSub<U1>, // For Hessenberg.
DefaultAllocator: Allocator<N, D, DimDiff<D, U1>>
+ Allocator<N, DimDiff<D, U1>>
+ Allocator<N, D, D>
+ Allocator<N, D>,
{
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<N, D>>
where
D: DimSub<U1>, // For Hessenberg.
DefaultAllocator: Allocator<N, D, DimDiff<D, U1>>
+ Allocator<N, DimDiff<D, U1>>
+ Allocator<N, D, D>
+ Allocator<N, D>,
{
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<N, D>
where
D: DimSub<U1>,
DefaultAllocator: Allocator<N, D, D>
+ Allocator<N, DimDiff<D, U1>>
+ Allocator<N::RealField, D>
+ Allocator<N::RealField, DimDiff<D, U1>>,
{
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<N, D>>
where
D: DimSub<U1>,
DefaultAllocator: Allocator<N, D, D>
+ Allocator<N, DimDiff<D, U1>>
+ Allocator<N::RealField, D>
+ Allocator<N::RealField, DimDiff<D, U1>>,
{
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<N, D>
where
D: DimSub<U1>,
DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>,
{
SymmetricTridiagonal::new(self.into_owned())
}
}

12
src/linalg/full_piv_lu.rs
View File

@ -257,15 +257,3 @@ where
}
}
}
impl<N: ComplexField, R: DimMin<C>, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S>
where
DefaultAllocator: Allocator<N, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
{
/// 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<N, R, C> {
FullPivLU::new(self.into_owned())
}
}

12
src/linalg/hessenberg.rs
View File

@ -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<N: ComplexField, D: DimSub<U1>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
where
DefaultAllocator: Allocator<N, D, D> + Allocator<N, D> + Allocator<N, DimDiff<D, U1>>,
{
/// Computes the Hessenberg decomposition of this matrix using householder reflections.
pub fn hessenberg(self) -> Hessenberg<N, D> {
Hessenberg::new(self.into_owned())
}
}

10
src/linalg/lu.rs
View File

@ -379,13 +379,3 @@ pub fn gauss_step_swap<N, R: Dim, C: Dim, S>(
.axpy(-pivot_row[k].inlined_clone(), &coeffs, N::one());
}
}
impl<N: ComplexField, R: DimMin<C>, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S>
where
DefaultAllocator: Allocator<N, R, C> + Allocator<(usize, usize), DimMinimum<R, C>>,
{
/// Computes the LU decomposition with partial (row) pivoting of `matrix`.
pub fn lu(self) -> LU<N, R, C> {
LU::new(self.into_owned())
}
}

1
src/linalg/mod.rs
View File

@ -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;

10
src/linalg/qr.rs
View File

@ -288,13 +288,3 @@ where
// res self.q_determinant()
// }
}
impl<N: ComplexField, R: DimMin<C>, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S>
where
DefaultAllocator: Allocator<N, R, C> + Allocator<N, R> + Allocator<N, DimMinimum<R, C>>,
{
/// Computes the QR decomposition of this matrix.
pub fn qr(self) -> QR<N, R, C> {
QR::new(self.into_owned())
}
}

20
src/linalg/schur.rs
View File

@ -496,26 +496,6 @@ where
+ Allocator<N, D, D>
+ Allocator<N, D>,
{
/// Computes the Schur decomposition of a square matrix.
pub fn schur(self) -> Schur<N, D> {
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<N, D>> {
Schur::try_new(self.into_owned(), eps, max_niter)
}
/// Computes the eigenvalues of this matrix.
pub fn eigenvalues(&self) -> Option<VectorN<N, D>> {
assert!(

25
src/linalg/svd.rs
View File

@ -616,31 +616,6 @@ where
+ Allocator<N::RealField, DimMinimum<R, C>>
+ Allocator<N::RealField, DimDiff<DimMinimum<R, C>, U1>>,
{
/// Computes the Singular Value Decomposition using implicit shift.
pub fn svd(self, compute_u: bool, compute_v: bool) -> SVD<N, R, C> {
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<N, R, C>> {
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<N::RealField, DimMinimum<R, C>> {
SVD::new(self.clone_owned(), false, false).singular_values

26
src/linalg/symmetric_eigen.rs
View File

@ -308,32 +308,6 @@ where
+ Allocator<N::RealField, D>
+ Allocator<N::RealField, DimDiff<D, U1>>,
{
/// 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<N, D> {
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<N, D>> {
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.

14
src/linalg/symmetric_tridiagonal.rs
View File

@ -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<N: ComplexField, D: DimSub<U1>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
where
DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>,
{
/// 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<N, D> {
SymmetricTridiagonal::new(self.into_owned())
}
}