nalgebra/linalg/
svd.rs

Help
1#[cfg(feature = "serde-serialize-no-std")]
2use serde::{Deserialize, Serialize};
3use std::any::TypeId;
4
5use approx::AbsDiffEq;
6use num::{One, Zero};
7
8use crate::allocator::Allocator;
9use crate::base::{DefaultAllocator, Matrix, Matrix2x3, OMatrix, OVector, Vector2};
10use crate::constraint::{SameNumberOfRows, ShapeConstraint};
11use crate::dimension::{Dim, DimDiff, DimMin, DimMinimum, DimSub, U1};
12use crate::storage::Storage;
13use crate::{Matrix2, Matrix3, RawStorage, U2, U3};
14use simba::scalar::{ComplexField, RealField};
15
16use crate::linalg::givens::GivensRotation;
17use crate::linalg::symmetric_eigen;
18use crate::linalg::Bidiagonal;
19
20/// Singular Value Decomposition of a general matrix.
21#[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))]
22#[cfg_attr(
23    feature = "serde-serialize-no-std",
24    serde(bound(
25        serialize = "DefaultAllocator: Allocator<T::RealField, DimMinimum<R, C>>    +
26                           Allocator<T, DimMinimum<R, C>, C> +
27                           Allocator<T, R, DimMinimum<R, C>>,
28         OMatrix<T, R, DimMinimum<R, C>>: Serialize,
29         OMatrix<T, DimMinimum<R, C>, C>: Serialize,
30         OVector<T::RealField, DimMinimum<R, C>>: Serialize"
31    ))
32)]
33#[cfg_attr(
34    feature = "serde-serialize-no-std",
35    serde(bound(
36        deserialize = "DefaultAllocator: Allocator<T::RealField, DimMinimum<R, C>>    +
37                           Allocator<T, DimMinimum<R, C>, C> +
38                           Allocator<T, R, DimMinimum<R, C>>,
39         OMatrix<T, R, DimMinimum<R, C>>: Deserialize<'de>,
40         OMatrix<T, DimMinimum<R, C>, C>: Deserialize<'de>,
41         OVector<T::RealField, DimMinimum<R, C>>: Deserialize<'de>"
42    ))
43)]
44#[derive(Clone, Debug)]
45pub struct SVD<T: ComplexField, R: DimMin<C>, C: Dim>
46where
47    DefaultAllocator: Allocator<T, DimMinimum<R, C>, C>
48        + Allocator<T, R, DimMinimum<R, C>>
49        + Allocator<T::RealField, DimMinimum<R, C>>,
50{
51    /// The left-singular vectors `U` of this SVD.
52    pub u: Option<OMatrix<T, R, DimMinimum<R, C>>>,
53    /// The right-singular vectors `V^t` of this SVD.
54    pub v_t: Option<OMatrix<T, DimMinimum<R, C>, C>>,
55    /// The singular values of this SVD.
56    pub singular_values: OVector<T::RealField, DimMinimum<R, C>>,
57}
58
59impl<T: ComplexField, R: DimMin<C>, C: Dim> Copy for SVD<T, R, C>
60where
61    DefaultAllocator: Allocator<T, DimMinimum<R, C>, C>
62        + Allocator<T, R, DimMinimum<R, C>>
63        + Allocator<T::RealField, DimMinimum<R, C>>,
64    OMatrix<T, R, DimMinimum<R, C>>: Copy,
65    OMatrix<T, DimMinimum<R, C>, C>: Copy,
66    OVector<T::RealField, DimMinimum<R, C>>: Copy,
67{
68}
69
70impl<T: ComplexField, R: DimMin<C>, C: Dim> SVD<T, R, C>
71where
72    DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
73    DefaultAllocator: Allocator<T, R, C>
74        + Allocator<T, C>
75        + Allocator<T, R>
76        + Allocator<T, DimDiff<DimMinimum<R, C>, U1>>
77        + Allocator<T, DimMinimum<R, C>, C>
78        + Allocator<T, R, DimMinimum<R, C>>
79        + Allocator<T, DimMinimum<R, C>>
80        + Allocator<T::RealField, DimMinimum<R, C>>
81        + Allocator<T::RealField, DimDiff<DimMinimum<R, C>, U1>>,
82{
83    fn use_special_always_ordered_svd2() -> bool {
84        TypeId::of::<OMatrix<T, R, C>>() == TypeId::of::<Matrix2<T::RealField>>()
85            && TypeId::of::<Self>() == TypeId::of::<SVD<T::RealField, U2, U2>>()
86    }
87
88    fn use_special_always_ordered_svd3() -> bool {
89        TypeId::of::<OMatrix<T, R, C>>() == TypeId::of::<Matrix3<T::RealField>>()
90            && TypeId::of::<Self>() == TypeId::of::<SVD<T::RealField, U3, U3>>()
91    }
92
93    /// Computes the Singular Value Decomposition of `matrix` using implicit shift.
94    /// The singular values are not guaranteed to be sorted in any particular order.
95    /// If a descending order is required, consider using `new` instead.
96    pub fn new_unordered(matrix: OMatrix<T, R, C>, compute_u: bool, compute_v: bool) -> Self {
97        Self::try_new_unordered(
98            matrix,
99            compute_u,
100            compute_v,
101            T::RealField::default_epsilon(),
102            0,
103        )
104        .unwrap()
105    }
106
107    /// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift.
108    /// The singular values are not guaranteed to be sorted in any particular order.
109    /// If a descending order is required, consider using `try_new` instead.
110    ///
111    /// # Arguments
112    ///
113    /// * `compute_u` − set this to `true` to enable the computation of left-singular vectors.
114    /// * `compute_v` − set this to `true` to enable the computation of right-singular vectors.
115    /// * `eps`       − tolerance used to determine when a value converged to 0.
116    /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this
117    /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
118    /// continues indefinitely until convergence.
119    pub fn try_new_unordered(
120        mut matrix: OMatrix<T, R, C>,
121        compute_u: bool,
122        compute_v: bool,
123        eps: T::RealField,
124        max_niter: usize,
125    ) -> Option<Self> {
126        assert!(
127            !matrix.is_empty(),
128            "Cannot compute the SVD of an empty matrix."
129        );
130        let (nrows, ncols) = matrix.shape_generic();
131        let min_nrows_ncols = nrows.min(ncols);
132
133        if Self::use_special_always_ordered_svd2() {
134            // SAFETY: the reference transmutes are OK since we checked that the types match exactly.
135            let matrix: &Matrix2<T::RealField> = unsafe { std::mem::transmute(&matrix) };
136            let result = super::svd2::svd_ordered2(matrix, compute_u, compute_v);
137            let typed_result: &Self = unsafe { std::mem::transmute(&result) };
138            return Some(typed_result.clone());
139        } else if Self::use_special_always_ordered_svd3() {
140            // SAFETY: the reference transmutes are OK since we checked that the types match exactly.
141            let matrix: &Matrix3<T::RealField> = unsafe { std::mem::transmute(&matrix) };
142            let result = super::svd3::svd_ordered3(matrix, compute_u, compute_v, eps, max_niter);
143            let typed_result: &Self = unsafe { std::mem::transmute(&result) };
144            return Some(typed_result.clone());
145        }
146
147        let dim = min_nrows_ncols.value();
148
149        let m_amax = matrix.camax();
150
151        if !m_amax.is_zero() {
152            matrix.unscale_mut(m_amax.clone());
153        }
154
155        let bi_matrix = Bidiagonal::new(matrix);
156        let mut u = if compute_u { Some(bi_matrix.u()) } else { None };
157        let mut v_t = if compute_v {
158            Some(bi_matrix.v_t())
159        } else {
160            None
161        };
162        let mut diagonal = bi_matrix.diagonal();
163        let mut off_diagonal = bi_matrix.off_diagonal();
164
165        let mut niter = 0;
166        let (mut start, mut end) = Self::delimit_subproblem(
167            &mut diagonal,
168            &mut off_diagonal,
169            &mut u,
170            &mut v_t,
171            bi_matrix.is_upper_diagonal(),
172            dim - 1,
173            eps.clone(),
174        );
175
176        while end != start {
177            let subdim = end - start + 1;
178
179            // Solve the subproblem.
180            #[allow(clippy::comparison_chain)]
181            if subdim > 2 {
182                let m = end - 1;
183                let n = end;
184
185                let mut vec;
186                {
187                    let dm = diagonal[m].clone();
188                    let dn = diagonal[n].clone();
189                    let fm = off_diagonal[m].clone();
190
191                    let tmm = dm.clone() * dm.clone()
192                        + off_diagonal[m - 1].clone() * off_diagonal[m - 1].clone();
193                    let tmn = dm * fm.clone();
194                    let tnn = dn.clone() * dn + fm.clone() * fm;
195
196                    let shift = symmetric_eigen::wilkinson_shift(tmm, tnn, tmn);
197
198                    vec = Vector2::new(
199                        diagonal[start].clone() * diagonal[start].clone() - shift,
200                        diagonal[start].clone() * off_diagonal[start].clone(),
201                    );
202                }
203
204                for k in start..n {
205                    let m12 = if k == n - 1 {
206                        T::RealField::zero()
207                    } else {
208                        off_diagonal[k + 1].clone()
209                    };
210
211                    let mut subm = Matrix2x3::new(
212                        diagonal[k].clone(),
213                        off_diagonal[k].clone(),
214                        T::RealField::zero(),
215                        T::RealField::zero(),
216                        diagonal[k + 1].clone(),
217                        m12,
218                    );
219
220                    if let Some((rot1, norm1)) = GivensRotation::cancel_y(&vec) {
221                        rot1.inverse()
222                            .rotate_rows(&mut subm.fixed_columns_mut::<2>(0));
223                        let rot1 = GivensRotation::new_unchecked(rot1.c(), T::from_real(rot1.s()));
224
225                        if k > start {
226                            // This is not the first iteration.
227                            off_diagonal[k - 1] = norm1;
228                        }
229
230                        let v = Vector2::new(subm[(0, 0)].clone(), subm[(1, 0)].clone());
231                        // TODO: does the case `v.y == 0` ever happen?
232                        let (rot2, norm2) = GivensRotation::cancel_y(&v)
233                            .unwrap_or((GivensRotation::identity(), subm[(0, 0)].clone()));
234
235                        rot2.rotate(&mut subm.fixed_columns_mut::<2>(1));
236                        let rot2 = GivensRotation::new_unchecked(rot2.c(), T::from_real(rot2.s()));
237
238                        subm[(0, 0)] = norm2;
239
240                        if let Some(ref mut v_t) = v_t {
241                            if bi_matrix.is_upper_diagonal() {
242                                rot1.rotate(&mut v_t.fixed_rows_mut::<2>(k));
243                            } else {
244                                rot2.rotate(&mut v_t.fixed_rows_mut::<2>(k));
245                            }
246                        }
247
248                        if let Some(ref mut u) = u {
249                            if bi_matrix.is_upper_diagonal() {
250                                rot2.inverse().rotate_rows(&mut u.fixed_columns_mut::<2>(k));
251                            } else {
252                                rot1.inverse().rotate_rows(&mut u.fixed_columns_mut::<2>(k));
253                            }
254                        }
255
256                        diagonal[k] = subm[(0, 0)].clone();
257                        diagonal[k + 1] = subm[(1, 1)].clone();
258                        off_diagonal[k] = subm[(0, 1)].clone();
259
260                        if k != n - 1 {
261                            off_diagonal[k + 1] = subm[(1, 2)].clone();
262                        }
263
264                        vec.x = subm[(0, 1)].clone();
265                        vec.y = subm[(0, 2)].clone();
266                    } else {
267                        break;
268                    }
269                }
270            } else if subdim == 2 {
271                // Solve the remaining 2x2 subproblem.
272                let (u2, s, v2) = compute_2x2_uptrig_svd(
273                    diagonal[start].clone(),
274                    off_diagonal[start].clone(),
275                    diagonal[start + 1].clone(),
276                    compute_u && bi_matrix.is_upper_diagonal()
277                        || compute_v && !bi_matrix.is_upper_diagonal(),
278                    compute_v && bi_matrix.is_upper_diagonal()
279                        || compute_u && !bi_matrix.is_upper_diagonal(),
280                );
281                let u2 = u2.map(|u2| GivensRotation::new_unchecked(u2.c(), T::from_real(u2.s())));
282                let v2 = v2.map(|v2| GivensRotation::new_unchecked(v2.c(), T::from_real(v2.s())));
283
284                diagonal[start] = s[0].clone();
285                diagonal[start + 1] = s[1].clone();
286                off_diagonal[start] = T::RealField::zero();
287
288                if let Some(ref mut u) = u {
289                    let rot = if bi_matrix.is_upper_diagonal() {
290                        u2.clone().unwrap()
291                    } else {
292                        v2.clone().unwrap()
293                    };
294                    rot.rotate_rows(&mut u.fixed_columns_mut::<2>(start));
295                }
296
297                if let Some(ref mut v_t) = v_t {
298                    let rot = if bi_matrix.is_upper_diagonal() {
299                        v2.unwrap()
300                    } else {
301                        u2.unwrap()
302                    };
303                    rot.inverse().rotate(&mut v_t.fixed_rows_mut::<2>(start));
304                }
305
306                end -= 1;
307            }
308
309            // Re-delimit the subproblem in case some decoupling occurred.
310            let sub = Self::delimit_subproblem(
311                &mut diagonal,
312                &mut off_diagonal,
313                &mut u,
314                &mut v_t,
315                bi_matrix.is_upper_diagonal(),
316                end,
317                eps.clone(),
318            );
319            start = sub.0;
320            end = sub.1;
321
322            niter += 1;
323            if niter == max_niter {
324                return None;
325            }
326        }
327
328        diagonal *= m_amax;
329
330        // Ensure all singular value are non-negative.
331        for i in 0..dim {
332            let sval = diagonal[i].clone();
333
334            if sval < T::RealField::zero() {
335                diagonal[i] = -sval;
336
337                if let Some(ref mut u) = u {
338                    u.column_mut(i).neg_mut();
339                }
340            }
341        }
342
343        Some(Self {
344            u,
345            v_t,
346            singular_values: diagonal,
347        })
348    }
349
350    /*
351    fn display_bidiag(b: &Bidiagonal<T, R, C>, begin: usize, end: usize) {
352        for i in begin .. end {
353            for k in begin .. i {
354                print!("    ");
355            }
356            println!("{}  {}", b.diagonal[i], b.off_diagonal[i]);
357        }
358        for k in begin .. end {
359            print!("    ");
360        }
361        println!("{}", b.diagonal[end]);
362    }
363    */
364
365    fn delimit_subproblem(
366        diagonal: &mut OVector<T::RealField, DimMinimum<R, C>>,
367        off_diagonal: &mut OVector<T::RealField, DimDiff<DimMinimum<R, C>, U1>>,
368        u: &mut Option<OMatrix<T, R, DimMinimum<R, C>>>,
369        v_t: &mut Option<OMatrix<T, DimMinimum<R, C>, C>>,
370        is_upper_diagonal: bool,
371        end: usize,
372        eps: T::RealField,
373    ) -> (usize, usize) {
374        let mut n = end;
375
376        while n > 0 {
377            let m = n - 1;
378
379            if off_diagonal[m].is_zero()
380                || off_diagonal[m].clone().norm1()
381                    <= eps.clone() * (diagonal[n].clone().norm1() + diagonal[m].clone().norm1())
382            {
383                off_diagonal[m] = T::RealField::zero();
384            } else if diagonal[m].clone().norm1() <= eps {
385                diagonal[m] = T::RealField::zero();
386                Self::cancel_horizontal_off_diagonal_elt(
387                    diagonal,
388                    off_diagonal,
389                    u,
390                    v_t,
391                    is_upper_diagonal,
392                    m,
393                    m + 1,
394                );
395
396                if m != 0 {
397                    Self::cancel_vertical_off_diagonal_elt(
398                        diagonal,
399                        off_diagonal,
400                        u,
401                        v_t,
402                        is_upper_diagonal,
403                        m - 1,
404                    );
405                }
406            } else if diagonal[n].clone().norm1() <= eps {
407                diagonal[n] = T::RealField::zero();
408                Self::cancel_vertical_off_diagonal_elt(
409                    diagonal,
410                    off_diagonal,
411                    u,
412                    v_t,
413                    is_upper_diagonal,
414                    m,
415                );
416            } else {
417                break;
418            }
419
420            n -= 1;
421        }
422
423        if n == 0 {
424            return (0, 0);
425        }
426
427        let mut new_start = n - 1;
428        while new_start > 0 {
429            let m = new_start - 1;
430
431            if off_diagonal[m].clone().norm1()
432                <= eps.clone() * (diagonal[new_start].clone().norm1() + diagonal[m].clone().norm1())
433            {
434                off_diagonal[m] = T::RealField::zero();
435                break;
436            }
437            // TODO: write a test that enters this case.
438            else if diagonal[m].clone().norm1() <= eps {
439                diagonal[m] = T::RealField::zero();
440                Self::cancel_horizontal_off_diagonal_elt(
441                    diagonal,
442                    off_diagonal,
443                    u,
444                    v_t,
445                    is_upper_diagonal,
446                    m,
447                    n,
448                );
449
450                if m != 0 {
451                    Self::cancel_vertical_off_diagonal_elt(
452                        diagonal,
453                        off_diagonal,
454                        u,
455                        v_t,
456                        is_upper_diagonal,
457                        m - 1,
458                    );
459                }
460                break;
461            }
462
463            new_start -= 1;
464        }
465
466        (new_start, n)
467    }
468
469    // Cancels the i-th off-diagonal element using givens rotations.
470    fn cancel_horizontal_off_diagonal_elt(
471        diagonal: &mut OVector<T::RealField, DimMinimum<R, C>>,
472        off_diagonal: &mut OVector<T::RealField, DimDiff<DimMinimum<R, C>, U1>>,
473        u: &mut Option<OMatrix<T, R, DimMinimum<R, C>>>,
474        v_t: &mut Option<OMatrix<T, DimMinimum<R, C>, C>>,
475        is_upper_diagonal: bool,
476        i: usize,
477        end: usize,
478    ) {
479        let mut v = Vector2::new(off_diagonal[i].clone(), diagonal[i + 1].clone());
480        off_diagonal[i] = T::RealField::zero();
481
482        for k in i..end {
483            if let Some((rot, norm)) = GivensRotation::cancel_x(&v) {
484                let rot = GivensRotation::new_unchecked(rot.c(), T::from_real(rot.s()));
485                diagonal[k + 1] = norm;
486
487                if is_upper_diagonal {
488                    if let Some(ref mut u) = *u {
489                        rot.inverse()
490                            .rotate_rows(&mut u.fixed_columns_with_step_mut::<2>(i, k - i));
491                    }
492                } else if let Some(ref mut v_t) = *v_t {
493                    rot.rotate(&mut v_t.fixed_rows_with_step_mut::<2>(i, k - i));
494                }
495
496                if k + 1 != end {
497                    v.x = -rot.s().real() * off_diagonal[k + 1].clone();
498                    v.y = diagonal[k + 2].clone();
499                    off_diagonal[k + 1] *= rot.c();
500                }
501            } else {
502                break;
503            }
504        }
505    }
506
507    // Cancels the i-th off-diagonal element using givens rotations.
508    fn cancel_vertical_off_diagonal_elt(
509        diagonal: &mut OVector<T::RealField, DimMinimum<R, C>>,
510        off_diagonal: &mut OVector<T::RealField, DimDiff<DimMinimum<R, C>, U1>>,
511        u: &mut Option<OMatrix<T, R, DimMinimum<R, C>>>,
512        v_t: &mut Option<OMatrix<T, DimMinimum<R, C>, C>>,
513        is_upper_diagonal: bool,
514        i: usize,
515    ) {
516        let mut v = Vector2::new(diagonal[i].clone(), off_diagonal[i].clone());
517        off_diagonal[i] = T::RealField::zero();
518
519        for k in (0..i + 1).rev() {
520            if let Some((rot, norm)) = GivensRotation::cancel_y(&v) {
521                let rot = GivensRotation::new_unchecked(rot.c(), T::from_real(rot.s()));
522                diagonal[k] = norm;
523
524                if is_upper_diagonal {
525                    if let Some(ref mut v_t) = *v_t {
526                        rot.rotate(&mut v_t.fixed_rows_with_step_mut::<2>(k, i - k));
527                    }
528                } else if let Some(ref mut u) = *u {
529                    rot.inverse()
530                        .rotate_rows(&mut u.fixed_columns_with_step_mut::<2>(k, i - k));
531                }
532
533                if k > 0 {
534                    v.x = diagonal[k - 1].clone();
535                    v.y = rot.s().real() * off_diagonal[k - 1].clone();
536                    off_diagonal[k - 1] *= rot.c();
537                }
538            } else {
539                break;
540            }
541        }
542    }
543
544    /// Computes the rank of the decomposed matrix, i.e., the number of singular values greater
545    /// than `eps`.
546    #[must_use]
547    pub fn rank(&self, eps: T::RealField) -> usize {
548        assert!(
549            eps >= T::RealField::zero(),
550            "SVD rank: the epsilon must be non-negative."
551        );
552        self.singular_values.iter().filter(|e| **e > eps).count()
553    }
554
555    /// Rebuild the original matrix.
556    ///
557    /// This is useful if some of the singular values have been manually modified.
558    /// Returns `Err` if the right- and left- singular vectors have not been
559    /// computed at construction-time.
560    pub fn recompose(self) -> Result<OMatrix<T, R, C>, &'static str> {
561        match (self.u, self.v_t) {
562            (Some(mut u), Some(v_t)) => {
563                for i in 0..self.singular_values.len() {
564                    let val = self.singular_values[i].clone();
565                    u.column_mut(i).scale_mut(val);
566                }
567                Ok(u * v_t)
568            }
569            (None, None) => Err("SVD recomposition: U and V^t have not been computed."),
570            (None, _) => Err("SVD recomposition: U has not been computed."),
571            (_, None) => Err("SVD recomposition: V^t has not been computed."),
572        }
573    }
574
575    /// Computes the pseudo-inverse of the decomposed matrix.
576    ///
577    /// Any singular value smaller than `eps` is assumed to be zero.
578    /// Returns `Err` if the right- and left- singular vectors have not
579    /// been computed at construction-time.
580    pub fn pseudo_inverse(mut self, eps: T::RealField) -> Result<OMatrix<T, C, R>, &'static str>
581    where
582        DefaultAllocator: Allocator<T, C, R>,
583    {
584        if eps < T::RealField::zero() {
585            Err("SVD pseudo inverse: the epsilon must be non-negative.")
586        } else {
587            for i in 0..self.singular_values.len() {
588                let val = self.singular_values[i].clone();
589
590                if val > eps {
591                    self.singular_values[i] = T::RealField::one() / val;
592                } else {
593                    self.singular_values[i] = T::RealField::zero();
594                }
595            }
596
597            self.recompose().map(|m| m.adjoint())
598        }
599    }
600
601    /// Solves the system `self * x = b` where `self` is the decomposed matrix and `x` the unknown.
602    ///
603    /// Any singular value smaller than `eps` is assumed to be zero.
604    /// Returns `Err` if the singular vectors `U` and `V` have not been computed.
605    // TODO: make this more generic wrt the storage types and the dimensions for `b`.
606    pub fn solve<R2: Dim, C2: Dim, S2>(
607        &self,
608        b: &Matrix<T, R2, C2, S2>,
609        eps: T::RealField,
610    ) -> Result<OMatrix<T, C, C2>, &'static str>
611    where
612        S2: Storage<T, R2, C2>,
613        DefaultAllocator: Allocator<T, C, C2> + Allocator<T, DimMinimum<R, C>, C2>,
614        ShapeConstraint: SameNumberOfRows<R, R2>,
615    {
616        if eps < T::RealField::zero() {
617            Err("SVD solve: the epsilon must be non-negative.")
618        } else {
619            match (&self.u, &self.v_t) {
620                (Some(u), Some(v_t)) => {
621                    let mut ut_b = u.ad_mul(b);
622
623                    for j in 0..ut_b.ncols() {
624                        let mut col = ut_b.column_mut(j);
625
626                        for i in 0..self.singular_values.len() {
627                            let val = self.singular_values[i].clone();
628                            if val > eps {
629                                col[i] = col[i].clone().unscale(val);
630                            } else {
631                                col[i] = T::zero();
632                            }
633                        }
634                    }
635
636                    Ok(v_t.ad_mul(&ut_b))
637                }
638                (None, None) => Err("SVD solve: U and V^t have not been computed."),
639                (None, _) => Err("SVD solve: U has not been computed."),
640                (_, None) => Err("SVD solve: V^t has not been computed."),
641            }
642        }
643    }
644
645    /// converts SVD results to Polar decomposition form of the original Matrix: `A = P' * U`.
646    ///
647    /// The polar decomposition used here is Left Polar Decomposition (or Reverse Polar Decomposition)
648    /// Returns None if the singular vectors of the SVD haven't been calculated
649    pub fn to_polar(&self) -> Option<(OMatrix<T, R, R>, OMatrix<T, R, C>)>
650    where
651        DefaultAllocator: Allocator<T, R, C> //result
652            + Allocator<T, DimMinimum<R, C>, R> // adjoint
653            + Allocator<T, DimMinimum<R, C>> // mapped vals
654            + Allocator<T, R, R> // result
655            + Allocator<T, DimMinimum<R, C>, DimMinimum<R, C>>, // square matrix
656    {
657        match (&self.u, &self.v_t) {
658            (Some(u), Some(v_t)) => Some((
659                u * OMatrix::from_diagonal(&self.singular_values.map(|e| T::from_real(e)))
660                    * u.adjoint(),
661                u * v_t,
662            )),
663            _ => None,
664        }
665    }
666}
667
668impl<T: ComplexField, R: DimMin<C>, C: Dim> SVD<T, R, C>
669where
670    DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
671    DefaultAllocator: Allocator<T, R, C>
672        + Allocator<T, C>
673        + Allocator<T, R>
674        + Allocator<T, DimDiff<DimMinimum<R, C>, U1>>
675        + Allocator<T, DimMinimum<R, C>, C>
676        + Allocator<T, R, DimMinimum<R, C>>
677        + Allocator<T, DimMinimum<R, C>>
678        + Allocator<T::RealField, DimMinimum<R, C>>
679        + Allocator<T::RealField, DimDiff<DimMinimum<R, C>, U1>>
680        + Allocator<(usize, usize), DimMinimum<R, C>> // for sorted singular values
681        + Allocator<(T::RealField, usize), DimMinimum<R, C>>, // for sorted singular values
682{
683    /// Computes the Singular Value Decomposition of `matrix` using implicit shift.
684    /// The singular values are guaranteed to be sorted in descending order.
685    /// If this order is not required consider using `new_unordered`.
686    pub fn new(matrix: OMatrix<T, R, C>, compute_u: bool, compute_v: bool) -> Self {
687        let mut svd = Self::new_unordered(matrix, compute_u, compute_v);
688
689        if !Self::use_special_always_ordered_svd3() && !Self::use_special_always_ordered_svd2() {
690            svd.sort_by_singular_values();
691        }
692
693        svd
694    }
695
696    /// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift.
697    /// The singular values are guaranteed to be sorted in descending order.
698    /// If this order is not required consider using `try_new_unordered`.
699    ///
700    /// # Arguments
701    ///
702    /// * `compute_u` − set this to `true` to enable the computation of left-singular vectors.
703    /// * `compute_v` − set this to `true` to enable the computation of right-singular vectors.
704    /// * `eps`       − tolerance used to determine when a value converged to 0.
705    /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this
706    /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
707    /// continues indefinitely until convergence.
708    pub fn try_new(
709        matrix: OMatrix<T, R, C>,
710        compute_u: bool,
711        compute_v: bool,
712        eps: T::RealField,
713        max_niter: usize,
714    ) -> Option<Self> {
715        Self::try_new_unordered(matrix, compute_u, compute_v, eps, max_niter).map(|mut svd| {
716            if !Self::use_special_always_ordered_svd3() && !Self::use_special_always_ordered_svd2()
717            {
718                svd.sort_by_singular_values();
719            }
720
721            svd
722        })
723    }
724
725    /// Sort the estimated components of the SVD by its singular values in descending order.
726    /// Such an ordering is often implicitly required when the decompositions are used for estimation or fitting purposes.
727    /// Using this function is only required if `new_unordered` or `try_new_unorderd` were used and the specific sorting is required afterward.
728    pub fn sort_by_singular_values(&mut self) {
729        const VALUE_PROCESSED: usize = usize::MAX;
730
731        // Collect the singular values with their original index, ...
732        let mut singular_values = self.singular_values.map_with_location(|r, _, e| (e, r));
733        assert_ne!(
734            singular_values.data.shape().0.value(),
735            VALUE_PROCESSED,
736            "Too many singular values"
737        );
738
739        // ... sort the singular values, ...
740        singular_values
741            .as_mut_slice()
742            .sort_unstable_by(|(a, _), (b, _)| b.partial_cmp(a).expect("Singular value was NaN"));
743
744        // ... and store them.
745        self.singular_values
746            .zip_apply(&singular_values, |value, (new_value, _)| {
747                value.clone_from(&new_value)
748            });
749
750        // Calculate required permutations given the sorted indices.
751        // We need to identify all circles to calculate the required swaps.
752        let mut permutations =
753            crate::PermutationSequence::identity_generic(singular_values.data.shape().0);
754
755        for i in 0..singular_values.len() {
756            let mut index_1 = i;
757            let mut index_2 = singular_values[i].1;
758
759            // Check whether the value was already visited ...
760            while index_2 != VALUE_PROCESSED // ... or a "double swap" must be avoided.
761                && singular_values[index_2].1 != VALUE_PROCESSED
762            {
763                // Add the permutation ...
764                permutations.append_permutation(index_1, index_2);
765                // ... and mark the value as visited.
766                singular_values[index_1].1 = VALUE_PROCESSED;
767
768                index_1 = index_2;
769                index_2 = singular_values[index_1].1;
770            }
771        }
772
773        // Permute the optional components
774        if let Some(u) = self.u.as_mut() {
775            permutations.permute_columns(u);
776        }
777
778        if let Some(v_t) = self.v_t.as_mut() {
779            permutations.permute_rows(v_t);
780        }
781    }
782}
783
784impl<T: ComplexField, R: DimMin<C>, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S>
785where
786    DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
787    DefaultAllocator: Allocator<T, R, C>
788        + Allocator<T, C>
789        + Allocator<T, R>
790        + Allocator<T, DimDiff<DimMinimum<R, C>, U1>>
791        + Allocator<T, DimMinimum<R, C>, C>
792        + Allocator<T, R, DimMinimum<R, C>>
793        + Allocator<T, DimMinimum<R, C>>
794        + Allocator<T::RealField, DimMinimum<R, C>>
795        + Allocator<T::RealField, DimDiff<DimMinimum<R, C>, U1>>,
796{
797    /// Computes the singular values of this matrix.
798    /// The singular values are not guaranteed to be sorted in any particular order.
799    /// If a descending order is required, consider using `singular_values` instead.
800    #[must_use]
801    pub fn singular_values_unordered(&self) -> OVector<T::RealField, DimMinimum<R, C>> {
802        SVD::new_unordered(self.clone_owned(), false, false).singular_values
803    }
804
805    /// Computes the rank of this matrix.
806    ///
807    /// All singular values below `eps` are considered equal to 0.
808    #[must_use]
809    pub fn rank(&self, eps: T::RealField) -> usize {
810        let svd = SVD::new_unordered(self.clone_owned(), false, false);
811        svd.rank(eps)
812    }
813
814    /// Computes the pseudo-inverse of this matrix.
815    ///
816    /// All singular values below `eps` are considered equal to 0.
817    pub fn pseudo_inverse(self, eps: T::RealField) -> Result<OMatrix<T, C, R>, &'static str>
818    where
819        DefaultAllocator: Allocator<T, C, R>,
820    {
821        SVD::new_unordered(self.clone_owned(), true, true).pseudo_inverse(eps)
822    }
823}
824
825impl<T: ComplexField, R: DimMin<C>, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S>
826where
827    DimMinimum<R, C>: DimSub<U1>,
828    DefaultAllocator: Allocator<T, R, C>
829        + Allocator<T, C>
830        + Allocator<T, R>
831        + Allocator<T, DimDiff<DimMinimum<R, C>, U1>>
832        + Allocator<T, DimMinimum<R, C>, C>
833        + Allocator<T, R, DimMinimum<R, C>>
834        + Allocator<T, DimMinimum<R, C>>
835        + Allocator<T::RealField, DimMinimum<R, C>>
836        + Allocator<T::RealField, DimDiff<DimMinimum<R, C>, U1>>
837        + Allocator<(usize, usize), DimMinimum<R, C>>
838        + Allocator<(T::RealField, usize), DimMinimum<R, C>>,
839{
840    /// Computes the singular values of this matrix.
841    /// The singular values are guaranteed to be sorted in descending order.
842    /// If this order is not required consider using `singular_values_unordered`.
843    #[must_use]
844    pub fn singular_values(&self) -> OVector<T::RealField, DimMinimum<R, C>> {
845        SVD::new(self.clone_owned(), false, false).singular_values
846    }
847}
848
849// Explicit formulae inspired from the paper "Computing the Singular Values of 2-by-2 Complex
850// Matrices", Sanzheng Qiao and Xiaohong Wang.
851// http://www.cas.mcmaster.ca/sqrl/papers/sqrl5.pdf
852fn compute_2x2_uptrig_svd<T: RealField>(
853    m11: T,
854    m12: T,
855    m22: T,
856    compute_u: bool,
857    compute_v: bool,
858) -> (
859    Option<GivensRotation<T>>,
860    Vector2<T>,
861    Option<GivensRotation<T>>,
862) {
863    let two: T::RealField = crate::convert(2.0f64);
864    let half: T::RealField = crate::convert(0.5f64);
865
866    let denom = (m11.clone() + m22.clone()).hypot(m12.clone())
867        + (m11.clone() - m22.clone()).hypot(m12.clone());
868
869    // NOTE: v1 is the singular value that is the closest to m22.
870    // This prevents cancellation issues when constructing the vector `csv` below. If we chose
871    // otherwise, we would have v1 ~= m11 when m12 is small. This would cause catastrophic
872    // cancellation on `v1 * v1 - m11 * m11` below.
873    let mut v1 = m11.clone() * m22.clone() * two / denom.clone();
874    let mut v2 = half * denom;
875
876    let mut u = None;
877    let mut v_t = None;
878
879    if compute_u || compute_v {
880        let (csv, sgn_v) = GivensRotation::new(
881            m11.clone() * m12.clone(),
882            v1.clone() * v1.clone() - m11.clone() * m11.clone(),
883        );
884        v1 *= sgn_v.clone();
885        v2 *= sgn_v;
886
887        if compute_v {
888            v_t = Some(csv.clone());
889        }
890
891        if compute_u {
892            let cu = (m11.scale(csv.c()) + m12 * csv.s()) / v1.clone();
893            let su = (m22 * csv.s()) / v1.clone();
894            let (csu, sgn_u) = GivensRotation::new(cu, su);
895
896            v1 *= sgn_u.clone();
897            v2 *= sgn_u;
898            u = Some(csu);
899        }
900    }
901
902    (u, Vector2::new(v1, v2), v_t)
903}
nalgebra/linalg/svd.rs

nalgebra/linalg/
svd.rs
