Allow sorting SVD according to singular values

src/linalg/decomposition.rs

@@ -74,7 +74,31 @@ impl<T: ComplexField, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     }
 
     /// Computes the Singular Value Decomposition using implicit shift.
+    /// The singular values are guaranteed to be sorted in descending order.
+    /// If this order is not required consider using `svd_unordered`.
     pub fn svd(self, compute_u: bool, compute_v: bool) -> SVD<T, R, C>
+    where
+        R: DimMin<C>,
+        DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
+        DefaultAllocator: Allocator<T, R, C>
+            + Allocator<T, C>
+            + Allocator<T, R>
+            + Allocator<T, DimDiff<DimMinimum<R, C>, U1>>
+            + Allocator<T, DimMinimum<R, C>, C>
+            + Allocator<T, R, DimMinimum<R, C>>
+            + Allocator<T, DimMinimum<R, C>>
+            + Allocator<T::RealField, DimMinimum<R, C>>
+            + Allocator<T::RealField, DimDiff<DimMinimum<R, C>, U1>>
+            + Allocator<(usize, usize), DimMinimum<R, C>>
+            + Allocator<(T::RealField, usize), DimMinimum<R, C>>,
+    {
+        SVD::new(self.into_owned(), compute_u, compute_v)
+    }
+
+    /// Computes the Singular Value Decomposition using implicit shift.
+    /// The singular values are not guaranteed to be sorted in any particular order.
+    /// If a descending order is required, consider using `svd` instead.
+    pub fn svd_unordered(self, compute_u: bool, compute_v: bool) -> SVD<T, R, C>
     where
         R: DimMin<C>,
         DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
@@ -88,10 +112,12 @@ impl<T: ComplexField, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
             + Allocator<T::RealField, DimMinimum<R, C>>
             + Allocator<T::RealField, DimDiff<DimMinimum<R, C>, U1>>,
     {
-        SVD::new(self.into_owned(), compute_u, compute_v)
+        SVD::new_unordered(self.into_owned(), compute_u, compute_v)
     }
 
     /// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift.
+    /// The singular values are guaranteed to be sorted in descending order.
+    /// If this order is not required consider using `try_svd_unordered`.
     ///
     /// # Arguments
     ///
@@ -119,10 +145,47 @@ impl<T: ComplexField, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
             + Allocator<T, R, DimMinimum<R, C>>
             + Allocator<T, DimMinimum<R, C>>
             + Allocator<T::RealField, DimMinimum<R, C>>
-            + Allocator<T::RealField, DimDiff<DimMinimum<R, C>, U1>>,
+            + Allocator<T::RealField, DimDiff<DimMinimum<R, C>, U1>>
+            + Allocator<(usize, usize), DimMinimum<R, C>>
+            + Allocator<(T::RealField, usize), DimMinimum<R, C>>,
     {
         SVD::try_new(self.into_owned(), compute_u, compute_v, eps, max_niter)
     }
+
+    /// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift.
+    /// The singular values are not guaranteed to be sorted in any particular order.
+    /// If a descending order is required, consider using `try_svd` instead.
+    ///
+    /// # Arguments
+    ///
+    /// * `compute_u` − set this to `true` to enable the computation of left-singular vectors.
+    /// * `compute_v` − set this to `true` to enable the computation of right-singular vectors.
+    /// * `eps`       − tolerance used to determine when a value converged to 0.
+    /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this
+    /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
+    /// continues indefinitely until convergence.
+    pub fn try_svd_unordered(
+        self,
+        compute_u: bool,
+        compute_v: bool,
+        eps: T::RealField,
+        max_niter: usize,
+    ) -> Option<SVD<T, R, C>>
+    where
+        R: DimMin<C>,
+        DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
+        DefaultAllocator: Allocator<T, R, C>
+            + Allocator<T, C>
+            + Allocator<T, R>
+            + Allocator<T, DimDiff<DimMinimum<R, C>, U1>>
+            + Allocator<T, DimMinimum<R, C>, C>
+            + Allocator<T, R, DimMinimum<R, C>>
+            + Allocator<T, DimMinimum<R, C>>
+            + Allocator<T::RealField, DimMinimum<R, C>>
+            + Allocator<T::RealField, DimDiff<DimMinimum<R, C>, U1>>,
+    {
+        SVD::try_new_unordered(self.into_owned(), compute_u, compute_v, eps, max_niter)
+    }
 }
 
 /// # Square matrix decomposition

src/linalg/svd.rs

@@ -9,6 +9,7 @@ use crate::base::{DefaultAllocator, Matrix, Matrix2x3, OMatrix, OVector, Vector2
 use crate::constraint::{SameNumberOfRows, ShapeConstraint};
 use crate::dimension::{Dim, DimDiff, DimMin, DimMinimum, DimSub, U1};
 use crate::storage::Storage;
+use crate::RawStorage;
 use simba::scalar::{ComplexField, RealField};
 
 use crate::linalg::givens::GivensRotation;
@@ -79,8 +80,10 @@ where
         + Allocator<T::RealField, DimDiff<DimMinimum<R, C>, U1>>,
 {
     /// Computes the Singular Value Decomposition of `matrix` using implicit shift.
-    pub fn new(matrix: OMatrix<T, R, C>, compute_u: bool, compute_v: bool) -> Self {
+    /// The singular values are not guaranteed to be sorted in any particular order.
-        Self::try_new(
+    /// If a descending order is required, consider using `new` instead.
+    pub fn new_unordered(matrix: OMatrix<T, R, C>, compute_u: bool, compute_v: bool) -> Self {
+        Self::try_new_unordered(
             matrix,
             compute_u,
             compute_v,
@@ -91,6 +94,8 @@ where
     }
 
     /// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift.
+    /// The singular values are not guaranteed to be sorted in any particular order.
+    /// If a descending order is required, consider using `try_new` instead.
     ///
     /// # Arguments
     ///
@@ -100,7 +105,7 @@ where
     /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this
     /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
     /// continues indefinitely until convergence.
-    pub fn try_new(
+    pub fn try_new_unordered(
         mut matrix: OMatrix<T, R, C>,
         compute_u: bool,
         compute_v: bool,
@@ -612,6 +617,114 @@ where
     }
 }
 
+impl<T: ComplexField, R: DimMin<C>, C: Dim> SVD<T, R, C>
+where
+    DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
+    DefaultAllocator: Allocator<T, R, C>
+        + Allocator<T, C>
+        + Allocator<T, R>
+        + Allocator<T, DimDiff<DimMinimum<R, C>, U1>>
+        + Allocator<T, DimMinimum<R, C>, C>
+        + Allocator<T, R, DimMinimum<R, C>>
+        + Allocator<T, DimMinimum<R, C>>
+        + Allocator<T::RealField, DimMinimum<R, C>>
+        + Allocator<T::RealField, DimDiff<DimMinimum<R, C>, U1>>
+        + Allocator<(usize, usize), DimMinimum<R, C>> // for sorted singular values
+        + Allocator<(T::RealField, usize), DimMinimum<R, C>>, // for sorted singular values
+{
+    /// Computes the Singular Value Decomposition of `matrix` using implicit shift.
+    /// The singular values are guaranteed to be sorted in descending order.
+    /// If this order is not required consider using `new_unordered`.
+    pub fn new(matrix: OMatrix<T, R, C>, compute_u: bool, compute_v: bool) -> Self {
+        let mut svd = Self::new_unordered(matrix, compute_u, compute_v);
+        svd.sort_by_singular_values();
+        svd
+    }
+
+    /// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift.
+    /// The singular values are guaranteed to be sorted in descending order.
+    /// If this order is not required consider using `try_new_unordered`.
+    ///
+    /// # Arguments
+    ///
+    /// * `compute_u` − set this to `true` to enable the computation of left-singular vectors.
+    /// * `compute_v` − set this to `true` to enable the computation of right-singular vectors.
+    /// * `eps`       − tolerance used to determine when a value converged to 0.
+    /// * `max_niter` − maximum total number of iterations performed by the algorithm. If this
+    /// number of iteration is exceeded, `None` is returned. If `niter == 0`, then the algorithm
+    /// continues indefinitely until convergence.
+    pub fn try_new(
+        matrix: OMatrix<T, R, C>,
+        compute_u: bool,
+        compute_v: bool,
+        eps: T::RealField,
+        max_niter: usize,
+    ) -> Option<Self> {
+        Self::try_new_unordered(matrix, compute_u, compute_v, eps, max_niter).map(|mut svd| {
+            svd.sort_by_singular_values();
+            svd
+        })
+    }
+
+    /// Sort the estimated components of the SVD by its singular values in descending order.
+    /// Such an ordering is often implicitly required when the decompositions are used for estimation or fitting purposes.
+    /// Using this function is only required if `new_unordered` or `try_new_unorderd` were used and the specific sorting is required afterward.
+    pub fn sort_by_singular_values(&mut self) {
+        const VALUE_PROCESSED: usize = usize::MAX;
+
+        // Collect the singular values with their original index, ...
+        let mut singular_values = self.singular_values.map_with_location(|r, _, e| (e, r));
+        assert_ne!(
+            singular_values.data.shape().0.value(),
+            VALUE_PROCESSED,
+            "Too many singular values"
+        );
+
+        // ... sort the singular values, ...
+        singular_values
+            .as_mut_slice()
+            .sort_unstable_by(|(a, _), (b, _)| b.partial_cmp(a).expect("Singular value was NaN"));
+
+        // ... and store them.
+        self.singular_values
+            .zip_apply(&singular_values, |value, (new_value, _)| {
+                value.clone_from(&new_value)
+            });
+
+        // Calculate required permutations given the sorted indices.
+        // We need to identify all circles to calculate the required swaps.
+        let mut permutations =
+            crate::PermutationSequence::identity_generic(singular_values.data.shape().0);
+
+        for i in 0..singular_values.len() {
+            let mut index_1 = i;
+            let mut index_2 = singular_values[i].1;
+
+            // Check whether the value was already visited ...
+            while index_2 != VALUE_PROCESSED // ... or a "double swap" must be avoided.
+                && singular_values[index_2].1 != VALUE_PROCESSED
+            {
+                // Add the permutation ...
+                permutations.append_permutation(index_1, index_2);
+                // ... and mark the value as visited.
+                singular_values[index_1].1 = VALUE_PROCESSED;
+
+                index_1 = index_2;
+                index_2 = singular_values[index_1].1;
+            }
+        }
+
+        // Permute the optional components
+        if let Some(u) = self.u.as_mut() {
+            permutations.permute_columns(u);
+        }
+
+        if let Some(v_t) = self.v_t.as_mut() {
+            permutations.permute_rows(v_t);
+        }
+    }
+}
+
 impl<T: ComplexField, R: DimMin<C>, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S>
 where
     DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
@@ -626,9 +739,11 @@ where
         + Allocator<T::RealField, DimDiff<DimMinimum<R, C>, U1>>,
 {
     /// Computes the singular values of this matrix.
+    /// The singular values are not guaranteed to be sorted in any particular order.
+    /// If a descending order is required, consider using `singular_values` instead.
     #[must_use]
-    pub fn singular_values(&self) -> OVector<T::RealField, DimMinimum<R, C>> {
+    pub fn singular_values_unordered(&self) -> OVector<T::RealField, DimMinimum<R, C>> {
-        SVD::new(self.clone_owned(), false, false).singular_values
+        SVD::new_unordered(self.clone_owned(), false, false).singular_values
     }
 
     /// Computes the rank of this matrix.
@@ -636,7 +751,7 @@ where
     /// All singular values below `eps` are considered equal to 0.
     #[must_use]
     pub fn rank(&self, eps: T::RealField) -> usize {
-        let svd = SVD::new(self.clone_owned(), false, false);
+        let svd = SVD::new_unordered(self.clone_owned(), false, false);
         svd.rank(eps)
     }
 
@@ -647,7 +762,31 @@ where
     where
         DefaultAllocator: Allocator<T, C, R>,
     {
-        SVD::new(self.clone_owned(), true, true).pseudo_inverse(eps)
+        SVD::new_unordered(self.clone_owned(), true, true).pseudo_inverse(eps)
+    }
+}
+
+impl<T: ComplexField, R: DimMin<C>, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S>
+where
+    DimMinimum<R, C>: DimSub<U1>,
+    DefaultAllocator: Allocator<T, R, C>
+        + Allocator<T, C>
+        + Allocator<T, R>
+        + Allocator<T, DimDiff<DimMinimum<R, C>, U1>>
+        + Allocator<T, DimMinimum<R, C>, C>
+        + Allocator<T, R, DimMinimum<R, C>>
+        + Allocator<T, DimMinimum<R, C>>
+        + Allocator<T::RealField, DimMinimum<R, C>>
+        + Allocator<T::RealField, DimDiff<DimMinimum<R, C>, U1>>
+        + Allocator<(usize, usize), DimMinimum<R, C>>
+        + Allocator<(T::RealField, usize), DimMinimum<R, C>>,
+{
+    /// Computes the singular values of this matrix.
+    /// The singular values are guaranteed to be sorted in descending order.
+    /// If this order is not required consider using `singular_values_unordered`.
+    #[must_use]
+    pub fn singular_values(&self) -> OVector<T::RealField, DimMinimum<R, C>> {
+        SVD::new(self.clone_owned(), false, false).singular_values
     }
 }
 

tests/linalg/svd.rs

@@ -326,6 +326,13 @@ fn svd_fail() {
         0.07311092531259344,   0.5579247949052946,  0.14518764691585773,  0.03502980663114896,   0.7991329455957719,   0.4929930019965745,
         0.12293810556077789,   0.6617084679545999,   0.9002240700227326, 0.027153062135304884,   0.3630189466989524,  0.18207502727558866,
         0.843196731466686,  0.08951878746549924,   0.7533450877576973, 0.009558876499740077,   0.9429679490873482,   0.9355764454129878);
+    
+    // Check unordered ...
+    let svd = m.clone().svd_unordered(true, true);
+    let recomp = svd.recompose().unwrap();
+    assert_relative_eq!(m, recomp, epsilon = 1.0e-5);
+
+    // ... and ordered SVD.
     let svd = m.clone().svd(true, true);
     let recomp = svd.recompose().unwrap();
     assert_relative_eq!(m, recomp, epsilon = 1.0e-5);
@@ -344,3 +351,45 @@ fn svd_err() {
         svd.clone().pseudo_inverse(-1.0)
     );
 }
+
+#[test]
+#[rustfmt::skip]
+fn svd_sorted() {
+    let reference = nalgebra::matrix![
+        1.0, 2.0, 3.0, 4.0;
+        5.0, 6.0, 7.0, 8.0;
+        9.0, 10.0, 11.0, 12.0
+    ];
+
+    let mut svd = nalgebra::SVD {
+        singular_values: nalgebra::matrix![1.72261225; 2.54368356e+01; 5.14037515e-16],
+        u: Some(nalgebra::matrix![
+            -0.88915331, -0.20673589, 0.40824829;
+            -0.25438183, -0.51828874, -0.81649658;
+            0.38038964, -0.82984158, 0.40824829
+        ]),
+        v_t: Some(nalgebra::matrix![
+            0.73286619,  0.28984978, -0.15316664, -0.59618305;
+            -0.40361757, -0.46474413, -0.52587069, -0.58699725;
+            0.44527162, -0.83143156,  0.32704826,  0.05911168
+        ]),
+    };
+
+    assert_relative_eq!(
+        svd.recompose().expect("valid SVD"),
+        reference,
+        epsilon = 1.0e-5
+    );
+
+    svd.sort_by_singular_values();
+
+    // Ensure successful sorting
+    assert_relative_eq!(svd.singular_values.x, 2.54368356e+01, epsilon = 1.0e-5);
+
+    // Ensure that the sorted components represent the same decomposition
+    assert_relative_eq!(
+        svd.recompose().expect("valid SVD"),
+        reference,
+        epsilon = 1.0e-5
+    );
+}