Merge branch 'dev' into master-public
# Conflicts: # src/linalg/svd.rs
This commit is contained in:
commit
82106caa9e
|
@ -484,33 +484,37 @@ where
|
||||||
|
|
||||||
/// Rebuild the original matrix.
|
/// Rebuild the original matrix.
|
||||||
///
|
///
|
||||||
/// This is useful if some of the singular values have been manually modified. Panics if the
|
/// This is useful if some of the singular values have been manually modified.
|
||||||
/// right- and left- singular vectors have not been computed at construction-time.
|
/// Returns `Err` if the right- and left- singular vectors have not been
|
||||||
pub fn recompose(self) -> MatrixMN<N, R, C> {
|
/// computed at construction-time.
|
||||||
let mut u = self.u.expect("SVD recomposition: U has not been computed.");
|
pub fn recompose(self) -> Result<MatrixMN<N, R, C>, &'static str> {
|
||||||
let v_t = self
|
match (self.u, self.v_t) {
|
||||||
.v_t
|
(Some(mut u), Some(v_t)) => {
|
||||||
.expect("SVD recomposition: V^t has not been computed.");
|
|
||||||
|
|
||||||
for i in 0..self.singular_values.len() {
|
for i in 0..self.singular_values.len() {
|
||||||
let val = self.singular_values[i];
|
let val = self.singular_values[i];
|
||||||
u.column_mut(i).mul_assign(val);
|
u.column_mut(i).mul_assign(val);
|
||||||
}
|
}
|
||||||
|
Ok(u * v_t)
|
||||||
u * v_t
|
}
|
||||||
|
(None, None) => Err("SVD recomposition: U and V^t have not been computed."),
|
||||||
|
(None, _) => Err("SVD recomposition: U has not been computed."),
|
||||||
|
(_, None) => Err("SVD recomposition: V^t has not been computed.")
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Computes the pseudo-inverse of the decomposed matrix.
|
/// Computes the pseudo-inverse of the decomposed matrix.
|
||||||
///
|
///
|
||||||
/// Any singular value smaller than `eps` is assumed to be zero.
|
/// Any singular value smaller than `eps` is assumed to be zero.
|
||||||
/// Panics if the right- and left- singular vectors have not been computed at
|
/// Returns `Err` if the right- and left- singular vectors have not
|
||||||
/// construction-time.
|
/// been computed at construction-time.
|
||||||
pub fn pseudo_inverse(mut self, eps: N) -> MatrixMN<N, C, R>
|
pub fn pseudo_inverse(mut self, eps: N) -> Result<MatrixMN<N, C, R>, &'static str>
|
||||||
where DefaultAllocator: Allocator<N, C, R> {
|
where
|
||||||
assert!(
|
DefaultAllocator: Allocator<N, C, R>,
|
||||||
eps >= N::zero(),
|
{
|
||||||
"SVD pseudo inverse: the epsilon must be non-negative."
|
if eps < N::zero() {
|
||||||
);
|
Err("SVD pseudo inverse: the epsilon must be non-negative.")
|
||||||
|
}
|
||||||
|
else {
|
||||||
for i in 0..self.singular_values.len() {
|
for i in 0..self.singular_values.len() {
|
||||||
let val = self.singular_values[i];
|
let val = self.singular_values[i];
|
||||||
|
|
||||||
|
@ -521,37 +525,31 @@ where
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
self.recompose().transpose()
|
self.recompose().map(|m| m.transpose())
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Solves the system `self * x = b` where `self` is the decomposed matrix and `x` the unknown.
|
/// Solves the system `self * x = b` where `self` is the decomposed matrix and `x` the unknown.
|
||||||
///
|
///
|
||||||
/// Any singular value smaller than `eps` is assumed to be zero.
|
/// Any singular value smaller than `eps` is assumed to be zero.
|
||||||
/// Returns `None` if the singular vectors `U` and `V` have not been computed.
|
/// Returns `Err` if the singular vectors `U` and `V` have not been computed.
|
||||||
// FIXME: make this more generic wrt the storage types and the dimensions for `b`.
|
// FIXME: make this more generic wrt the storage types and the dimensions for `b`.
|
||||||
pub fn solve<R2: Dim, C2: Dim, S2>(
|
pub fn solve<R2: Dim, C2: Dim, S2>(
|
||||||
&self,
|
&self,
|
||||||
b: &Matrix<N, R2, C2, S2>,
|
b: &Matrix<N, R2, C2, S2>,
|
||||||
eps: N,
|
eps: N,
|
||||||
) -> MatrixMN<N, C, C2>
|
) -> Result<MatrixMN<N, C, C2>, &'static str>
|
||||||
where
|
where
|
||||||
S2: Storage<N, R2, C2>,
|
S2: Storage<N, R2, C2>,
|
||||||
DefaultAllocator: Allocator<N, C, C2> + Allocator<N, DimMinimum<R, C>, C2>,
|
DefaultAllocator: Allocator<N, C, C2> + Allocator<N, DimMinimum<R, C>, C2>,
|
||||||
ShapeConstraint: SameNumberOfRows<R, R2>,
|
ShapeConstraint: SameNumberOfRows<R, R2>,
|
||||||
{
|
{
|
||||||
assert!(
|
if eps < N::zero() {
|
||||||
eps >= N::zero(),
|
Err("SVD solve: the epsilon must be non-negative.")
|
||||||
"SVD solve: the epsilon must be non-negative."
|
}
|
||||||
);
|
else {
|
||||||
let u = self
|
match (&self.u, &self.v_t) {
|
||||||
.u
|
(Some(u), Some(v_t)) => {
|
||||||
.as_ref()
|
|
||||||
.expect("SVD solve: U has not been computed.");
|
|
||||||
let v_t = self
|
|
||||||
.v_t
|
|
||||||
.as_ref()
|
|
||||||
.expect("SVD solve: V^t has not been computed.");
|
|
||||||
|
|
||||||
let mut ut_b = u.tr_mul(b);
|
let mut ut_b = u.tr_mul(b);
|
||||||
|
|
||||||
for j in 0..ut_b.ncols() {
|
for j in 0..ut_b.ncols() {
|
||||||
|
@ -567,7 +565,13 @@ where
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
v_t.tr_mul(&ut_b)
|
Ok(v_t.tr_mul(&ut_b))
|
||||||
|
}
|
||||||
|
(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.")
|
||||||
|
}
|
||||||
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -624,8 +628,10 @@ where
|
||||||
/// Computes the pseudo-inverse of this matrix.
|
/// Computes the pseudo-inverse of this matrix.
|
||||||
///
|
///
|
||||||
/// All singular values below `eps` are considered equal to 0.
|
/// All singular values below `eps` are considered equal to 0.
|
||||||
pub fn pseudo_inverse(self, eps: N) -> MatrixMN<N, C, R>
|
pub fn pseudo_inverse(self, eps: N) -> Result<MatrixMN<N, C, R>, &'static str>
|
||||||
where DefaultAllocator: Allocator<N, C, R> {
|
where
|
||||||
|
DefaultAllocator: Allocator<N, C, R>,
|
||||||
|
{
|
||||||
SVD::new(self.clone_owned(), true, true).pseudo_inverse(eps)
|
SVD::new(self.clone_owned(), true, true).pseudo_inverse(eps)
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
|
@ -10,7 +10,7 @@ mod quickcheck_tests {
|
||||||
fn svd(m: DMatrix<f64>) -> bool {
|
fn svd(m: DMatrix<f64>) -> bool {
|
||||||
if m.len() > 0 {
|
if m.len() > 0 {
|
||||||
let svd = m.clone().svd(true, true);
|
let svd = m.clone().svd(true, true);
|
||||||
let recomp_m = svd.clone().recompose();
|
let recomp_m = svd.clone().recompose().unwrap();
|
||||||
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
|
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
|
||||||
let ds = DMatrix::from_diagonal(&s);
|
let ds = DMatrix::from_diagonal(&s);
|
||||||
|
|
||||||
|
@ -91,7 +91,7 @@ mod quickcheck_tests {
|
||||||
fn svd_pseudo_inverse(m: DMatrix<f64>) -> bool {
|
fn svd_pseudo_inverse(m: DMatrix<f64>) -> bool {
|
||||||
if m.len() > 0 {
|
if m.len() > 0 {
|
||||||
let svd = m.clone().svd(true, true);
|
let svd = m.clone().svd(true, true);
|
||||||
let pinv = svd.pseudo_inverse(1.0e-10);
|
let pinv = svd.pseudo_inverse(1.0e-10).unwrap();
|
||||||
|
|
||||||
if m.nrows() > m.ncols() {
|
if m.nrows() > m.ncols() {
|
||||||
println!("{}", &pinv * &m);
|
println!("{}", &pinv * &m);
|
||||||
|
@ -118,10 +118,10 @@ mod quickcheck_tests {
|
||||||
let b1 = DVector::new_random(n);
|
let b1 = DVector::new_random(n);
|
||||||
let b2 = DMatrix::new_random(n, nb);
|
let b2 = DMatrix::new_random(n, nb);
|
||||||
|
|
||||||
let sol1 = svd.solve(&b1, 1.0e-7);
|
let sol1 = svd.solve(&b1, 1.0e-7).unwrap();
|
||||||
let sol2 = svd.solve(&b2, 1.0e-7);
|
let sol2 = svd.solve(&b2, 1.0e-7).unwrap();
|
||||||
|
|
||||||
let recomp = svd.recompose();
|
let recomp = svd.recompose().unwrap();
|
||||||
if !relative_eq!(m, recomp, epsilon = 1.0e-6) {
|
if !relative_eq!(m, recomp, epsilon = 1.0e-6) {
|
||||||
println!("{}{}", m, recomp);
|
println!("{}{}", m, recomp);
|
||||||
}
|
}
|
||||||
|
@ -263,22 +263,22 @@ fn svd_singular_horizontal() {
|
||||||
fn svd_zeros() {
|
fn svd_zeros() {
|
||||||
let m = DMatrix::from_element(10, 10, 0.0);
|
let m = DMatrix::from_element(10, 10, 0.0);
|
||||||
let svd = m.clone().svd(true, true);
|
let svd = m.clone().svd(true, true);
|
||||||
assert_eq!(m, svd.recompose());
|
assert_eq!(Ok(m), svd.recompose());
|
||||||
}
|
}
|
||||||
|
|
||||||
#[test]
|
#[test]
|
||||||
fn svd_identity() {
|
fn svd_identity() {
|
||||||
let m = DMatrix::<f64>::identity(10, 10);
|
let m = DMatrix::<f64>::identity(10, 10);
|
||||||
let svd = m.clone().svd(true, true);
|
let svd = m.clone().svd(true, true);
|
||||||
assert_eq!(m, svd.recompose());
|
assert_eq!(Ok(m), svd.recompose());
|
||||||
|
|
||||||
let m = DMatrix::<f64>::identity(10, 15);
|
let m = DMatrix::<f64>::identity(10, 15);
|
||||||
let svd = m.clone().svd(true, true);
|
let svd = m.clone().svd(true, true);
|
||||||
assert_eq!(m, svd.recompose());
|
assert_eq!(Ok(m), svd.recompose());
|
||||||
|
|
||||||
let m = DMatrix::<f64>::identity(15, 10);
|
let m = DMatrix::<f64>::identity(15, 10);
|
||||||
let svd = m.clone().svd(true, true);
|
let svd = m.clone().svd(true, true);
|
||||||
assert_eq!(m, svd.recompose());
|
assert_eq!(Ok(m), svd.recompose());
|
||||||
}
|
}
|
||||||
|
|
||||||
#[test]
|
#[test]
|
||||||
|
@ -295,7 +295,7 @@ fn svd_with_delimited_subproblem() {
|
||||||
m[(8,8)] = 16.0; m[(3,9)] = 17.0;
|
m[(8,8)] = 16.0; m[(3,9)] = 17.0;
|
||||||
m[(9,9)] = 18.0;
|
m[(9,9)] = 18.0;
|
||||||
let svd = m.clone().svd(true, true);
|
let svd = m.clone().svd(true, true);
|
||||||
assert!(relative_eq!(m, svd.recompose(), epsilon = 1.0e-7));
|
assert!(relative_eq!(m, svd.recompose().unwrap(), epsilon = 1.0e-7));
|
||||||
|
|
||||||
// Rectangular versions.
|
// Rectangular versions.
|
||||||
let mut m = DMatrix::<f64>::from_element(15, 10, 0.0);
|
let mut m = DMatrix::<f64>::from_element(15, 10, 0.0);
|
||||||
|
@ -310,10 +310,10 @@ fn svd_with_delimited_subproblem() {
|
||||||
m[(8,8)] = 16.0; m[(3,9)] = 17.0;
|
m[(8,8)] = 16.0; m[(3,9)] = 17.0;
|
||||||
m[(9,9)] = 18.0;
|
m[(9,9)] = 18.0;
|
||||||
let svd = m.clone().svd(true, true);
|
let svd = m.clone().svd(true, true);
|
||||||
assert!(relative_eq!(m, svd.recompose(), epsilon = 1.0e-7));
|
assert!(relative_eq!(m, svd.recompose().unwrap(), epsilon = 1.0e-7));
|
||||||
|
|
||||||
let svd = m.transpose().svd(true, true);
|
let svd = m.transpose().svd(true, true);
|
||||||
assert!(relative_eq!(m.transpose(), svd.recompose(), epsilon = 1.0e-7));
|
assert!(relative_eq!(m.transpose(), svd.recompose().unwrap(), epsilon = 1.0e-7));
|
||||||
}
|
}
|
||||||
|
|
||||||
#[test]
|
#[test]
|
||||||
|
@ -329,7 +329,15 @@ fn svd_fail() {
|
||||||
println!("Singular values: {}", svd.singular_values);
|
println!("Singular values: {}", svd.singular_values);
|
||||||
println!("u: {:.5}", svd.u.unwrap());
|
println!("u: {:.5}", svd.u.unwrap());
|
||||||
println!("v: {:.5}", svd.v_t.unwrap());
|
println!("v: {:.5}", svd.v_t.unwrap());
|
||||||
let recomp = svd.recompose();
|
let recomp = svd.recompose().unwrap();
|
||||||
println!("{:.5}{:.5}", m, recomp);
|
println!("{:.5}{:.5}", m, recomp);
|
||||||
assert!(relative_eq!(m, recomp, epsilon = 1.0e-5));
|
assert!(relative_eq!(m, recomp, epsilon = 1.0e-5));
|
||||||
}
|
}
|
||||||
|
|
||||||
|
#[test]
|
||||||
|
fn svd_err() {
|
||||||
|
let m = DMatrix::from_element(10, 10, 0.0);
|
||||||
|
let svd = m.clone().svd(false, false);
|
||||||
|
assert_eq!(Err("SVD recomposition: U and V^t have not been computed."), svd.clone().recompose());
|
||||||
|
assert_eq!(Err("SVD pseudo inverse: the epsilon must be non-negative."), svd.clone().pseudo_inverse(-1.0));
|
||||||
|
}
|
||||||
|
|
Loading…
Reference in New Issue