Skip to content

Commit

Permalink
First compiling commit for take-2 of polar-decomposition:
Browse files Browse the repository at this point in the history
Code inspired by this thread: dimforge#656
Main person behind this is LucasCampos
  • Loading branch information
metric-space committed Dec 22, 2021
1 parent 507ead2 commit 6ac6e7f
Show file tree
Hide file tree
Showing 2 changed files with 40 additions and 0 deletions.
29 changes: 29 additions & 0 deletions src/linalg/svd.rs
Expand Up @@ -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.
Expand Down
11 changes: 11 additions & 0 deletions tests/linalg/svd.rs
Expand Up @@ -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);

}

0 comments on commit 6ac6e7f

Please sign in to comment.