nalgebra/linalg/
svd.rs

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

nalgebra/linalg/
svd.rs
