First compiling commit for take-2 of polar-decomposition:

Code inspired by this thread: https://github.com/dimforge/nalgebra/pull/656
Main person behind this is LucasCampos
This commit is contained in:
metric-space 2021-12-22 00:12:27 -05:00
parent 507ead2f10
commit 6ac6e7f75e
2 changed files with 40 additions and 0 deletions

View File

@ -643,6 +643,35 @@ where
}
}
impl<T: ComplexField, R: DimMin<C>, C: Dim> SVD<T, R, C>
where
DefaultAllocator: Allocator<T, DimMinimum<R, C>, C>
+ Allocator<T, R, DimMinimum<R, C>>
+ Allocator<T::RealField, DimMinimum<R, C>>,
{
/// converts SVD results to a polar form
pub fn to_polar(&self) -> Result<(OMatrix<T, R, R>, OMatrix<T, R, C>), &'static str>
where DefaultAllocator: Allocator<T, R, C> //result
+ Allocator<T, DimMinimum<R, C>, R> // adjoint
+ Allocator<T, DimMinimum<R, C>> // mapped vals
+ Allocator<T, R, R> // square matrix & result
+ Allocator<T, DimMinimum<R, C>, DimMinimum<R, C>> // ?
,
{
match (&self.u, &self.v_t) {
(Some(u), Some(v_t)) => Ok((
u * OMatrix::from_diagonal(&self.singular_values.map(|e| T::from_real(e)))
* u.adjoint(),
u * v_t,
)),
(None, None) => Err("SVD solve: U and V^t have not been computed."),
(None, _) => Err("SVD solve: U has not been computed."),
(_, None) => Err("SVD solve: V^t has not been computed."),
}
}
}
impl<T: ComplexField, R: DimMin<C>, C: Dim> SVD<T, R, C>
where
DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.

View File

@ -441,3 +441,14 @@ fn svd_sorted() {
epsilon = 1.0e-5
);
}
#[test]
fn svd_polar_decomposition() {
let m = DMatrix::<f64>::new_random(4, 4);
let svd = m.clone().svd(true, true);
let (p,u) = svd.to_polar().unwrap();
assert_relative_eq!(m, p*u, epsilon = 1.0e-5);
}