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