Merge pull request #1018 from Christopher22/sorted_svd
Allow descending sorting of estimated SVD
This commit is contained in:
commit
640ab4b12d
|
@ -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
|
||||
|
|
|
@ -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 {
|
||||
Self::try_new(
|
||||
/// The singular values are not guaranteed to be sorted in any particular order.
|
||||
/// 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>> {
|
||||
SVD::new(self.clone_owned(), false, false).singular_values
|
||||
pub fn singular_values_unordered(&self) -> OVector<T::RealField, DimMinimum<R, C>> {
|
||||
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
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
@ -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
|
||||
);
|
||||
}
|
||||
|
|
Loading…
Reference in New Issue