nalgebra/linalg/decomposition.rs
1use crate::storage::Storage;
2use crate::{
3 Allocator, Bidiagonal, Cholesky, ColPivQR, ComplexField, DefaultAllocator, Dim, DimDiff,
4 DimMin, DimMinimum, DimSub, FullPivLU, Hessenberg, LBLT, LU, Matrix, OMatrix, QR, RealField,
5 SVD, Schur, SymmetricEigen, SymmetricTridiagonal, U1, UDU,
6};
7
8/// # Rectangular matrix decomposition
9///
10/// This section contains the methods for computing some common decompositions of rectangular
11/// matrices with real or complex components. The following are currently supported:
12///
13/// | Decomposition | Factors | Details |
14/// | -------------------------|---------------------|--------------|
15/// | QR | `Q * R` | `Q` is an unitary matrix, and `R` is upper-triangular. |
16/// | QR with column pivoting | `Q * R * P⁻¹` | `Q` is an unitary matrix, and `R` is upper-triangular. `P` is a permutation matrix. |
17/// | 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. |
18/// | 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. |
19/// | SVD | `U * Σ * Vᵀ` | `U` and `V` are two orthogonal matrices and `Σ` is a diagonal matrix containing the singular values. |
20/// | Polar (Left Polar) | `P' * U` | `U` is semi-unitary/unitary and `P'` is a positive semi-definite Hermitian Matrix
21impl<T: ComplexField, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
22 /// Computes the bidiagonalization using householder reflections.
23 pub fn bidiagonalize(self) -> Bidiagonal<T, R, C>
24 where
25 R: DimMin<C>,
26 DimMinimum<R, C>: DimSub<U1>,
27 DefaultAllocator: Allocator<R, C>
28 + Allocator<C>
29 + Allocator<R>
30 + Allocator<DimMinimum<R, C>>
31 + Allocator<DimDiff<DimMinimum<R, C>, U1>>,
32 {
33 Bidiagonal::new(self.into_owned())
34 }
35
36 /// Computes the LU decomposition with full pivoting of `matrix`.
37 ///
38 /// This effectively computes `P, L, U, Q` such that `P * matrix * Q = LU`.
39 pub fn full_piv_lu(self) -> FullPivLU<T, R, C>
40 where
41 R: DimMin<C>,
42 DefaultAllocator: Allocator<R, C> + Allocator<DimMinimum<R, C>>,
43 {
44 FullPivLU::new(self.into_owned())
45 }
46
47 /// Computes the LU decomposition with partial (row) pivoting of `matrix`.
48 pub fn lu(self) -> LU<T, R, C>
49 where
50 R: DimMin<C>,
51 DefaultAllocator: Allocator<R, C> + Allocator<DimMinimum<R, C>>,
52 {
53 LU::new(self.into_owned())
54 }
55
56 /// Computes the QR decomposition of this matrix.
57 pub fn qr(self) -> QR<T, R, C>
58 where
59 R: DimMin<C>,
60 DefaultAllocator: Allocator<R, C> + Allocator<R> + Allocator<DimMinimum<R, C>>,
61 {
62 QR::new(self.into_owned())
63 }
64
65 /// Computes the QR decomposition (with column pivoting) of this matrix.
66 pub fn col_piv_qr(self) -> ColPivQR<T, R, C>
67 where
68 R: DimMin<C>,
69 DefaultAllocator: Allocator<R, C>
70 + Allocator<R>
71 + Allocator<DimMinimum<R, C>>
72 + Allocator<DimMinimum<R, C>>,
73 {
74 ColPivQR::new(self.into_owned())
75 }
76
77 /// Computes the Singular Value Decomposition using implicit shift.
78 /// The singular values are guaranteed to be sorted in descending order.
79 /// If this order is not required consider using `svd_unordered`.
80 pub fn svd(self, compute_u: bool, compute_v: bool) -> SVD<T, R, C>
81 where
82 R: DimMin<C>,
83 DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
84 DefaultAllocator: Allocator<R, C>
85 + Allocator<C>
86 + Allocator<R>
87 + Allocator<DimDiff<DimMinimum<R, C>, U1>>
88 + Allocator<DimMinimum<R, C>, C>
89 + Allocator<R, DimMinimum<R, C>>
90 + Allocator<DimMinimum<R, C>>,
91 {
92 SVD::new(self.into_owned(), compute_u, compute_v)
93 }
94
95 /// Computes the Singular Value Decomposition using implicit shift.
96 /// The singular values are not guaranteed to be sorted in any particular order.
97 /// If a descending order is required, consider using `svd` instead.
98 pub fn svd_unordered(self, compute_u: bool, compute_v: bool) -> SVD<T, R, C>
99 where
100 R: DimMin<C>,
101 DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
102 DefaultAllocator: Allocator<R, C>
103 + Allocator<C>
104 + Allocator<R>
105 + Allocator<DimDiff<DimMinimum<R, C>, U1>>
106 + Allocator<DimMinimum<R, C>, C>
107 + Allocator<R, DimMinimum<R, C>>
108 + Allocator<DimMinimum<R, C>>,
109 {
110 SVD::new_unordered(self.into_owned(), compute_u, compute_v)
111 }
112
113 /// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift.
114 /// The singular values are guaranteed to be sorted in descending order.
115 /// If this order is not required consider using `try_svd_unordered`.
116 ///
117 /// # Arguments
118 ///
119 /// * `compute_u` − set this to `true` to enable the computation of left-singular vectors.
120 /// * `compute_v` − set this to `true` to enable the computation of right-singular vectors.
121 /// * `eps` − tolerance used to determine when a value converged to 0.
122 /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this
123 /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
124 /// continues indefinitely until convergence.
125 pub fn try_svd(
126 self,
127 compute_u: bool,
128 compute_v: bool,
129 eps: T::RealField,
130 max_niter: usize,
131 ) -> Option<SVD<T, R, C>>
132 where
133 R: DimMin<C>,
134 DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
135 DefaultAllocator: Allocator<R, C>
136 + Allocator<C>
137 + Allocator<R>
138 + Allocator<DimDiff<DimMinimum<R, C>, U1>>
139 + Allocator<DimMinimum<R, C>, C>
140 + Allocator<R, DimMinimum<R, C>>
141 + Allocator<DimMinimum<R, C>>,
142 {
143 SVD::try_new(self.into_owned(), compute_u, compute_v, eps, max_niter)
144 }
145
146 /// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift.
147 /// The singular values are not guaranteed to be sorted in any particular order.
148 /// If a descending order is required, consider using `try_svd` instead.
149 ///
150 /// # Arguments
151 ///
152 /// * `compute_u` − set this to `true` to enable the computation of left-singular vectors.
153 /// * `compute_v` − set this to `true` to enable the computation of right-singular vectors.
154 /// * `eps` − tolerance used to determine when a value converged to 0.
155 /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this
156 /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
157 /// continues indefinitely until convergence.
158 pub fn try_svd_unordered(
159 self,
160 compute_u: bool,
161 compute_v: bool,
162 eps: T::RealField,
163 max_niter: usize,
164 ) -> Option<SVD<T, R, C>>
165 where
166 R: DimMin<C>,
167 DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
168 DefaultAllocator: Allocator<R, C>
169 + Allocator<C>
170 + Allocator<R>
171 + Allocator<DimDiff<DimMinimum<R, C>, U1>>
172 + Allocator<DimMinimum<R, C>, C>
173 + Allocator<R, DimMinimum<R, C>>
174 + Allocator<DimMinimum<R, C>>
175 + Allocator<DimMinimum<R, C>>
176 + Allocator<DimDiff<DimMinimum<R, C>, U1>>,
177 {
178 SVD::try_new_unordered(self.into_owned(), compute_u, compute_v, eps, max_niter)
179 }
180
181 /// Computes the Polar Decomposition of a `matrix` (indirectly uses SVD).
182 pub fn polar(self) -> (OMatrix<T, R, R>, OMatrix<T, R, C>)
183 where
184 R: DimMin<C>,
185 DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
186 DefaultAllocator: Allocator<R, C>
187 + Allocator<DimMinimum<R, C>, R>
188 + Allocator<DimMinimum<R, C>>
189 + Allocator<R, R>
190 + Allocator<DimMinimum<R, C>, DimMinimum<R, C>>
191 + Allocator<C>
192 + Allocator<R>
193 + Allocator<DimDiff<DimMinimum<R, C>, U1>>
194 + Allocator<DimMinimum<R, C>, C>
195 + Allocator<R, DimMinimum<R, C>>
196 + Allocator<DimMinimum<R, C>>
197 + Allocator<DimMinimum<R, C>>
198 + Allocator<DimDiff<DimMinimum<R, C>, U1>>,
199 {
200 SVD::new_unordered(self.into_owned(), true, true)
201 .to_polar()
202 .unwrap()
203 }
204
205 /// Attempts to compute the Polar Decomposition of a `matrix` (indirectly uses SVD).
206 ///
207 /// # Arguments
208 ///
209 /// * `eps` − tolerance used to determine when a value converged to 0 when computing the SVD.
210 /// * `max_niter` − maximum total number of iterations performed by the SVD computation algorithm.
211 pub fn try_polar(
212 self,
213 eps: T::RealField,
214 max_niter: usize,
215 ) -> Option<(OMatrix<T, R, R>, OMatrix<T, R, C>)>
216 where
217 R: DimMin<C>,
218 DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
219 DefaultAllocator: Allocator<R, C>
220 + Allocator<DimMinimum<R, C>, R>
221 + Allocator<DimMinimum<R, C>>
222 + Allocator<R, R>
223 + Allocator<DimMinimum<R, C>, DimMinimum<R, C>>
224 + Allocator<C>
225 + Allocator<R>
226 + Allocator<DimDiff<DimMinimum<R, C>, U1>>
227 + Allocator<DimMinimum<R, C>, C>
228 + Allocator<R, DimMinimum<R, C>>
229 + Allocator<DimMinimum<R, C>>
230 + Allocator<DimMinimum<R, C>>
231 + Allocator<DimDiff<DimMinimum<R, C>, U1>>,
232 {
233 SVD::try_new_unordered(self.into_owned(), true, true, eps, max_niter)
234 .and_then(|svd| svd.to_polar())
235 }
236}
237
238/// # Square matrix decomposition
239///
240/// This section contains the methods for computing some common decompositions of square
241/// matrices with real or complex components. The following are currently supported:
242///
243/// | Decomposition | Factors | Details |
244/// | -------------------------|---------------------------|--------------|
245/// | Hessenberg | `Q * H * Qᵀ` | `Q` is a unitary matrix and `H` an upper-Hessenberg matrix. |
246/// | Cholesky | `L * Lᵀ` | `L` is a lower-triangular matrix. |
247/// | LBLT decomposition | `Pᵀ * L * B * Lᴴ * P` | `L` is unit lower-triangular, `B` is Hermitian block-diagonal, and `P` is a permutation matrix. |
248/// | UDU | `U * D * Uᵀ` | `U` is a upper-triangular matrix, and `D` a diagonal matrix. |
249/// | Schur decomposition | `Q * T * Qᵀ` | `Q` is an unitary matrix and `T` a quasi-upper-triangular matrix. |
250/// | Symmetric eigendecomposition | `Q ~ Λ ~ Qᵀ` | `Q` is an unitary matrix, and `Λ` is a real diagonal matrix. |
251/// | Symmetric tridiagonalization | `Q ~ T ~ Qᵀ` | `Q` is an unitary matrix, and `T` is a tridiagonal matrix. |
252impl<T: ComplexField, D: Dim, S: Storage<T, D, D>> Matrix<T, D, D, S> {
253 /// Attempts to compute the Cholesky decomposition of this matrix.
254 ///
255 /// Returns `None` if the input matrix is not definite-positive. The input matrix is assumed
256 /// to be symmetric and only the lower-triangular part is read.
257 pub fn cholesky(self) -> Option<Cholesky<T, D>>
258 where
259 DefaultAllocator: Allocator<D, D>,
260 {
261 Cholesky::new(self.into_owned())
262 }
263
264 /// Computes the LBLT decomposition of this matrix.
265 /// The input matrix `self` is assumed to be Hermitian (symmetric) and the decomposition will
266 /// only read the lower-triangular part of `self`.
267 pub fn lblt(self) -> LBLT<T, D>
268 where
269 T: Copy,
270 T::RealField: Copy,
271 DefaultAllocator: Allocator<D> + Allocator<D, D>,
272 {
273 LBLT::new(self.into_owned())
274 }
275
276 /// Attempts to compute the UDU decomposition of this matrix.
277 ///
278 /// The input matrix `self` is assumed to be symmetric and this decomposition will only read
279 /// the upper-triangular part of `self`.
280 pub fn udu(self) -> Option<UDU<T, D>>
281 where
282 T: RealField,
283 DefaultAllocator: Allocator<D> + Allocator<D, D>,
284 {
285 UDU::new(self.into_owned())
286 }
287
288 /// Computes the Hessenberg decomposition of this matrix using householder reflections.
289 pub fn hessenberg(self) -> Hessenberg<T, D>
290 where
291 D: DimSub<U1>,
292 DefaultAllocator: Allocator<D, D> + Allocator<D> + Allocator<DimDiff<D, U1>>,
293 {
294 Hessenberg::new(self.into_owned())
295 }
296
297 /// Computes the Schur decomposition of a square matrix.
298 pub fn schur(self) -> Schur<T, D>
299 where
300 D: DimSub<U1>, // For Hessenberg.
301 DefaultAllocator: Allocator<D, DimDiff<D, U1>>
302 + Allocator<DimDiff<D, U1>>
303 + Allocator<D, D>
304 + Allocator<D>,
305 {
306 Schur::new(self.into_owned())
307 }
308
309 /// Attempts to compute the Schur decomposition of a square matrix.
310 ///
311 /// If only eigenvalues are needed, it is more efficient to call the matrix method
312 /// `.eigenvalues()` instead.
313 ///
314 /// # Arguments
315 ///
316 /// * `eps` − tolerance used to determine when a value converged to 0.
317 /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this
318 /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
319 /// continues indefinitely until convergence.
320 pub fn try_schur(self, eps: T::RealField, max_niter: usize) -> Option<Schur<T, D>>
321 where
322 D: DimSub<U1>, // For Hessenberg.
323 DefaultAllocator: Allocator<D, DimDiff<D, U1>>
324 + Allocator<DimDiff<D, U1>>
325 + Allocator<D, D>
326 + Allocator<D>,
327 {
328 Schur::try_new(self.into_owned(), eps, max_niter)
329 }
330
331 /// Computes the eigendecomposition of this symmetric matrix.
332 ///
333 /// Only the lower-triangular part (including the diagonal) of `m` is read.
334 pub fn symmetric_eigen(self) -> SymmetricEigen<T, D>
335 where
336 D: DimSub<U1>,
337 DefaultAllocator:
338 Allocator<D, D> + Allocator<DimDiff<D, U1>> + Allocator<D> + Allocator<DimDiff<D, U1>>,
339 {
340 SymmetricEigen::new(self.into_owned())
341 }
342
343 /// Computes the eigendecomposition of the given symmetric matrix with user-specified
344 /// convergence parameters.
345 ///
346 /// Only the lower-triangular part (including the diagonal) of `m` is read.
347 ///
348 /// # Arguments
349 ///
350 /// * `eps` − tolerance used to determine when a value converged to 0.
351 /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this
352 /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
353 /// continues indefinitely until convergence.
354 pub fn try_symmetric_eigen(
355 self,
356 eps: T::RealField,
357 max_niter: usize,
358 ) -> Option<SymmetricEigen<T, D>>
359 where
360 D: DimSub<U1>,
361 DefaultAllocator:
362 Allocator<D, D> + Allocator<DimDiff<D, U1>> + Allocator<D> + Allocator<DimDiff<D, U1>>,
363 {
364 SymmetricEigen::try_new(self.into_owned(), eps, max_niter)
365 }
366
367 /// Computes the tridiagonalization of this symmetric matrix.
368 ///
369 /// Only the lower-triangular part (including the diagonal) of `m` is read.
370 pub fn symmetric_tridiagonalize(self) -> SymmetricTridiagonal<T, D>
371 where
372 D: DimSub<U1>,
373 DefaultAllocator: Allocator<D, D> + Allocator<DimDiff<D, U1>>,
374 {
375 SymmetricTridiagonal::new(self.into_owned())
376 }
377}