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<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/// | UDU | `U * D * Uᵀ` | `U` is a upper-triangular matrix, and `D` a diagonal matrix. |
248/// | Schur decomposition | `Q * T * Qᵀ` | `Q` is an unitary matrix and `T` a quasi-upper-triangular matrix. |
249/// | Symmetric eigendecomposition | `Q ~ Λ ~ Qᵀ` | `Q` is an unitary matrix, and `Λ` is a real diagonal matrix. |
250/// | Symmetric tridiagonalization | `Q ~ T ~ Qᵀ` | `Q` is an unitary matrix, and `T` is a tridiagonal matrix. |
251impl<T: ComplexField, D: Dim, S: Storage<T, D, D>> Matrix<T, D, D, S> {
252 /// Attempts to compute the Cholesky decomposition of this matrix.
253 ///
254 /// Returns `None` if the input matrix is not definite-positive. The input matrix is assumed
255 /// to be symmetric and only the lower-triangular part is read.
256 pub fn cholesky(self) -> Option<Cholesky<T, D>>
257 where
258 DefaultAllocator: Allocator<D, D>,
259 {
260 Cholesky::new(self.into_owned())
261 }
262
263 /// Attempts to compute the UDU decomposition of this matrix.
264 ///
265 /// The input matrix `self` is assumed to be symmetric and this decomposition will only read
266 /// the upper-triangular part of `self`.
267 pub fn udu(self) -> Option<UDU<T, D>>
268 where
269 T: RealField,
270 DefaultAllocator: Allocator<D> + Allocator<D, D>,
271 {
272 UDU::new(self.into_owned())
273 }
274
275 /// Computes the Hessenberg decomposition of this matrix using householder reflections.
276 pub fn hessenberg(self) -> Hessenberg<T, D>
277 where
278 D: DimSub<U1>,
279 DefaultAllocator: Allocator<D, D> + Allocator<D> + Allocator<DimDiff<D, U1>>,
280 {
281 Hessenberg::new(self.into_owned())
282 }
283
284 /// Computes the Schur decomposition of a square matrix.
285 pub fn schur(self) -> Schur<T, D>
286 where
287 D: DimSub<U1>, // For Hessenberg.
288 DefaultAllocator: Allocator<D, DimDiff<D, U1>>
289 + Allocator<DimDiff<D, U1>>
290 + Allocator<D, D>
291 + Allocator<D>,
292 {
293 Schur::new(self.into_owned())
294 }
295
296 /// Attempts to compute the Schur decomposition of a square matrix.
297 ///
298 /// If only eigenvalues are needed, it is more efficient to call the matrix method
299 /// `.eigenvalues()` instead.
300 ///
301 /// # Arguments
302 ///
303 /// * `eps` − tolerance used to determine when a value converged to 0.
304 /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this
305 /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
306 /// continues indefinitely until convergence.
307 pub fn try_schur(self, eps: T::RealField, max_niter: usize) -> Option<Schur<T, D>>
308 where
309 D: DimSub<U1>, // For Hessenberg.
310 DefaultAllocator: Allocator<D, DimDiff<D, U1>>
311 + Allocator<DimDiff<D, U1>>
312 + Allocator<D, D>
313 + Allocator<D>,
314 {
315 Schur::try_new(self.into_owned(), eps, max_niter)
316 }
317
318 /// Computes the eigendecomposition of this symmetric matrix.
319 ///
320 /// Only the lower-triangular part (including the diagonal) of `m` is read.
321 pub fn symmetric_eigen(self) -> SymmetricEigen<T, D>
322 where
323 D: DimSub<U1>,
324 DefaultAllocator:
325 Allocator<D, D> + Allocator<DimDiff<D, U1>> + Allocator<D> + Allocator<DimDiff<D, U1>>,
326 {
327 SymmetricEigen::new(self.into_owned())
328 }
329
330 /// Computes the eigendecomposition of the given symmetric matrix with user-specified
331 /// convergence parameters.
332 ///
333 /// Only the lower-triangular part (including the diagonal) of `m` is read.
334 ///
335 /// # Arguments
336 ///
337 /// * `eps` − tolerance used to determine when a value converged to 0.
338 /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this
339 /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
340 /// continues indefinitely until convergence.
341 pub fn try_symmetric_eigen(
342 self,
343 eps: T::RealField,
344 max_niter: usize,
345 ) -> Option<SymmetricEigen<T, D>>
346 where
347 D: DimSub<U1>,
348 DefaultAllocator:
349 Allocator<D, D> + Allocator<DimDiff<D, U1>> + Allocator<D> + Allocator<DimDiff<D, U1>>,
350 {
351 SymmetricEigen::try_new(self.into_owned(), eps, max_niter)
352 }
353
354 /// Computes the tridiagonalization of this symmetric matrix.
355 ///
356 /// Only the lower-triangular part (including the diagonal) of `m` is read.
357 pub fn symmetric_tridiagonalize(self) -> SymmetricTridiagonal<T, D>
358 where
359 D: DimSub<U1>,
360 DefaultAllocator: Allocator<D, D> + Allocator<DimDiff<D, U1>>,
361 {
362 SymmetricTridiagonal::new(self.into_owned())
363 }
364}