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