M-Labs
/
nalgebra
3
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Decomposition results: return a real vector whenever applicable. 

This includes singular values, eigenvalues of hermitian matrices, tridiagonalization and bidiagonalization diagonal and off-diagonal elements.
This commit is contained in:
sebcrozet 2019-03-19 14:22:59 +01:00
parent 2f0d95bdbb
commit 3edef2f006
7 changed files with 198 additions and 177 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
src/base/matrix.rs
View File

@ -1026,10 +1026,19 @@ impl<N: Complex, D: Dim, S: StorageMut<N, D, D>> Matrix<N, D, D, S> {
} }
impl<N: Scalar, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> { impl<N: Scalar, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
/// Creates a square matrix with its diagonal set to `diag` and all other entries set to 0. /// The diagonal of this matrix.
#[inline] #[inline]
pub fn diagonal(&self) -> VectorN<N, D> pub fn diagonal(&self) -> VectorN<N, D>
where DefaultAllocator: Allocator<N, D> { where DefaultAllocator: Allocator<N, D> {
self.map_diagonal(|e| e)
}
/// Apply the given function to this matrix's diagonal and returns it.
///
/// This is a more efficient version of `self.diagonal().map(f)` since this
/// allocates only once.
pub fn map_diagonal<N2: Scalar>(&self, mut f: impl FnMut(N) -> N2) -> VectorN<N2, D>
where DefaultAllocator: Allocator<N2, D> {
assert!( assert!(
self.is_square(), self.is_square(),
"Unable to get the diagonal of a non-square matrix." "Unable to get the diagonal of a non-square matrix."
@ -1040,7 +1049,7 @@ impl<N: Scalar, D: Dim, S: Storage<N, D, D>> SquareMatrix<N, D, S> {
for i in 0..dim.value() { for i in 0..dim.value() {
unsafe { unsafe {
*res.vget_unchecked_mut(i) = *self.get_unchecked((i, i)); *res.vget_unchecked_mut(i) = f(*self.get_unchecked((i, i)));
} }
} }

10
src/linalg/bidiagonal.rs
View File

@ -272,13 +272,15 @@ where
} }
/// The diagonal part of this decomposed matrix. /// The diagonal part of this decomposed matrix.
pub fn diagonal(&self) -> VectorN<N, DimMinimum<R, C>> { pub fn diagonal(&self) -> VectorN<N::Real, DimMinimum<R, C>>
self.diagonal.map(|e| N::from_real(e.modulus())) where DefaultAllocator: Allocator<N::Real, DimMinimum<R, C>> {
self.diagonal.map(|e| e.modulus())
} }
/// The off-diagonal part of this decomposed matrix. /// The off-diagonal part of this decomposed matrix.
pub fn off_diagonal(&self) -> VectorN<N, DimDiff<DimMinimum<R, C>, U1>> { pub fn off_diagonal(&self) -> VectorN<N::Real, DimDiff<DimMinimum<R, C>, U1>>
self.off_diagonal.map(|e| N::from_real(e.modulus())) where DefaultAllocator: Allocator<N::Real, DimDiff<DimMinimum<R, C>, U1>> {
self.off_diagonal.map(|e| e.modulus())
} }
#[doc(hidden)] #[doc(hidden)]

10
src/linalg/givens.rs
View File

@ -28,6 +28,16 @@ impl<N: Complex> GivensRotation<N> {
} }
} }
/// Initializes a Givens rotation from its components.
///
/// The components are copies as-is. It is not checked whether they describe
/// an actually valid Givens rotation.
pub fn new_unchecked(c: N::Real, s: N) -> Self {
Self {
c, s
}
}
/// Initializes a Givens rotation from its non-normalized cosine an sine components. /// Initializes a Givens rotation from its non-normalized cosine an sine components.
pub fn new(c: N, s: N) -> (Self, N) { pub fn new(c: N, s: N) -> (Self, N) {
Self::try_new(c, s, N::Real::zero()).unwrap() Self::try_new(c, s, N::Real::zero()).unwrap()

212
src/linalg/svd.rs
View File

@ -2,11 +2,11 @@
use serde::{Deserialize, Serialize}; use serde::{Deserialize, Serialize};
use num_complex::Complex as NumComplex; use num_complex::Complex as NumComplex;
use num::Zero; use num::{Zero, One};
use std::ops::MulAssign; use std::ops::MulAssign;
use approx::AbsDiffEq; use approx::AbsDiffEq;
use alga::general::Complex; use alga::general::{Real, Complex};
use allocator::Allocator; use allocator::Allocator;
use base::{DefaultAllocator, Matrix, Matrix2x3, MatrixMN, Vector2, VectorN}; use base::{DefaultAllocator, Matrix, Matrix2x3, MatrixMN, Vector2, VectorN};
use constraint::{SameNumberOfRows, ShapeConstraint}; use constraint::{SameNumberOfRows, ShapeConstraint};
@ -23,51 +23,47 @@ use linalg::givens::GivensRotation;
#[cfg_attr( #[cfg_attr(
feature = "serde-serialize", feature = "serde-serialize",
serde(bound( serde(bound(
serialize = "DefaultAllocator: Allocator<N, R, C> + serialize = "DefaultAllocator: Allocator<N::Real, DimMinimum<R, C>> +
Allocator<N, DimMinimum<R, C>> +
Allocator<N, DimMinimum<R, C>, C> + Allocator<N, DimMinimum<R, C>, C> +
Allocator<N, R, DimMinimum<R, C>>, Allocator<N, R, DimMinimum<R, C>>,
MatrixMN<N, R, DimMinimum<R, C>>: Serialize, MatrixMN<N, R, DimMinimum<R, C>>: Serialize,
MatrixMN<N, DimMinimum<R, C>, C>: Serialize, MatrixMN<N, DimMinimum<R, C>, C>: Serialize,
VectorN<N, DimMinimum<R, C>>: Serialize" VectorN<N::Real, DimMinimum<R, C>>: Serialize"
)) ))
)] )]
#[cfg_attr( #[cfg_attr(
feature = "serde-serialize", feature = "serde-serialize",
serde(bound( serde(bound(
deserialize = "DefaultAllocator: Allocator<N, R, C> + deserialize = "DefaultAllocator: Allocator<N::Real, DimMinimum<R, C>> +
Allocator<N, DimMinimum<R, C>> +
Allocator<N, DimMinimum<R, C>, C> + Allocator<N, DimMinimum<R, C>, C> +
Allocator<N, R, DimMinimum<R, C>>, Allocator<N, R, DimMinimum<R, C>>,
MatrixMN<N, R, DimMinimum<R, C>>: Deserialize<'de>, MatrixMN<N, R, DimMinimum<R, C>>: Deserialize<'de>,
MatrixMN<N, DimMinimum<R, C>, C>: Deserialize<'de>, MatrixMN<N, DimMinimum<R, C>, C>: Deserialize<'de>,
VectorN<N, DimMinimum<R, C>>: Deserialize<'de>" VectorN<N::Real, DimMinimum<R, C>>: Deserialize<'de>"
)) ))
)] )]
#[derive(Clone, Debug)] #[derive(Clone, Debug)]
pub struct SVD<N: Complex, R: DimMin<C>, C: Dim> pub struct SVD<N: Complex, R: DimMin<C>, C: Dim>
where DefaultAllocator: Allocator<N, DimMinimum<R, C>, C> where DefaultAllocator: Allocator<N, DimMinimum<R, C>, C>
+ Allocator<N, R, DimMinimum<R, C>> + Allocator<N, R, DimMinimum<R, C>>
+ Allocator<N, DimMinimum<R, C>> + Allocator<N::Real, DimMinimum<R, C>>
{ {
/// The left-singular vectors `U` of this SVD. /// The left-singular vectors `U` of this SVD.
pub u: Option<MatrixMN<N, R, DimMinimum<R, C>>>, pub u: Option<MatrixMN<N, R, DimMinimum<R, C>>>,
/// The right-singular vectors `V^t` of this SVD. /// The right-singular vectors `V^t` of this SVD.
pub v_t: Option<MatrixMN<N, DimMinimum<R, C>, C>>, pub v_t: Option<MatrixMN<N, DimMinimum<R, C>, C>>,
// FIXME: this should a vector of N::Real
// because singular values are necessarily real.
/// The singular values of this SVD. /// The singular values of this SVD.
pub singular_values: VectorN<N, DimMinimum<R, C>>, pub singular_values: VectorN<N::Real, DimMinimum<R, C>>,
} }
impl<N: Complex, R: DimMin<C>, C: Dim> Copy for SVD<N, R, C> impl<N: Complex, R: DimMin<C>, C: Dim> Copy for SVD<N, R, C>
where where
DefaultAllocator: Allocator<N, DimMinimum<R, C>, C> DefaultAllocator: Allocator<N, DimMinimum<R, C>, C>
+ Allocator<N, R, DimMinimum<R, C>> + Allocator<N, R, DimMinimum<R, C>>
+ Allocator<N, DimMinimum<R, C>>, + Allocator<N::Real, DimMinimum<R, C>>,
MatrixMN<N, R, DimMinimum<R, C>>: Copy, MatrixMN<N, R, DimMinimum<R, C>>: Copy,
MatrixMN<N, DimMinimum<R, C>, C>: Copy, MatrixMN<N, DimMinimum<R, C>, C>: Copy,
VectorN<N, DimMinimum<R, C>>: Copy, VectorN<N::Real, DimMinimum<R, C>>: Copy,
{} {}
impl<N: Complex, R: DimMin<C>, C: Dim> SVD<N, R, C> impl<N: Complex, R: DimMin<C>, C: Dim> SVD<N, R, C>
@ -79,7 +75,9 @@ where
+ Allocator<N, DimDiff<DimMinimum<R, C>, U1>> + Allocator<N, DimDiff<DimMinimum<R, C>, U1>>
+ Allocator<N, DimMinimum<R, C>, C> + Allocator<N, DimMinimum<R, C>, C>
+ Allocator<N, R, DimMinimum<R, C>> + Allocator<N, R, DimMinimum<R, C>>
+ Allocator<N, DimMinimum<R, C>>, + Allocator<N, DimMinimum<R, C>>
+ Allocator<N::Real, DimMinimum<R, C>>
+ Allocator<N::Real, DimDiff<DimMinimum<R, C>, U1>>,
{ {
/// Computes the Singular Value Decomposition of `matrix` using implicit shift. /// Computes the Singular Value Decomposition of `matrix` using implicit shift.
pub fn new(matrix: MatrixMN<N, R, C>, compute_u: bool, compute_v: bool) -> Self { pub fn new(matrix: MatrixMN<N, R, C>, compute_u: bool, compute_v: bool) -> Self {
@ -155,7 +153,7 @@ where
for k in start..n { for k in start..n {
let m12 = if k == n - 1 { let m12 = if k == n - 1 {
N::zero() N::Real::zero()
} else { } else {
off_diagonal[k + 1] off_diagonal[k + 1]
}; };
@ -163,15 +161,15 @@ where
let mut subm = Matrix2x3::new( let mut subm = Matrix2x3::new(
diagonal[k], diagonal[k],
off_diagonal[k], off_diagonal[k],
N::zero(), N::Real::zero(),
N::zero(), N::Real::zero(),
diagonal[k + 1], diagonal[k + 1],
m12, m12,
); );
if let Some((rot1, norm1)) = GivensRotation::cancel_y(&vec) { if let Some((rot1, norm1)) = GivensRotation::cancel_y(&vec) {
rot1.inverse() rot1.inverse().rotate_rows(&mut subm.fixed_columns_mut::<U2>(0));
.rotate_rows(&mut subm.fixed_columns_mut::<U2>(0)); let rot1 = GivensRotation::new_unchecked(rot1.c(), N::from_real(rot1.s()));
if k > start { if k > start {
// This is not the first iteration. // This is not the first iteration.
@ -182,7 +180,10 @@ where
// FIXME: does the case `v.y == 0` ever happen? // FIXME: does the case `v.y == 0` ever happen?
let (rot2, norm2) = let (rot2, norm2) =
GivensRotation::cancel_y(&v).unwrap_or((GivensRotation::identity(), subm[(0, 0)])); GivensRotation::cancel_y(&v).unwrap_or((GivensRotation::identity(), subm[(0, 0)]));
rot2.rotate(&mut subm.fixed_columns_mut::<U2>(1)); rot2.rotate(&mut subm.fixed_columns_mut::<U2>(1));
let rot2 = GivensRotation::new_unchecked(rot2.c(), N::from_real(rot2.s()));
subm[(0, 0)] = norm2; subm[(0, 0)] = norm2;
if let Some(ref mut v_t) = v_t { if let Some(ref mut v_t) = v_t {
@ -219,17 +220,19 @@ where
} }
} else if subdim == 2 { } else if subdim == 2 {
// Solve the remaining 2x2 subproblem. // Solve the remaining 2x2 subproblem.
let (u2, s, v2) = Self::compute_2x2_uptrig_svd( let (u2, s, v2) = compute_2x2_uptrig_svd(
diagonal[start], diagonal[start],
off_diagonal[start], off_diagonal[start],
diagonal[start + 1], diagonal[start + 1],
compute_u && b.is_upper_diagonal() || compute_v && !b.is_upper_diagonal(), compute_u && b.is_upper_diagonal() || compute_v && !b.is_upper_diagonal(),
compute_v && b.is_upper_diagonal() || compute_u && !b.is_upper_diagonal(), compute_v && b.is_upper_diagonal() || compute_u && !b.is_upper_diagonal(),
); );
let u2 = u2.map(|u2| GivensRotation::new_unchecked(u2.c(), N::from_real(u2.s())));
let v2 = v2.map(|v2| GivensRotation::new_unchecked(v2.c(), N::from_real(v2.s())));
diagonal[start + 0] = s[0]; diagonal[start + 0] = s[0];
diagonal[start + 1] = s[1]; diagonal[start + 1] = s[1];
off_diagonal[start] = N::zero(); off_diagonal[start] = N::Real::zero();
if let Some(ref mut u) = u { if let Some(ref mut u) = u {
let rot = if b.is_upper_diagonal() { let rot = if b.is_upper_diagonal() {
@ -263,17 +266,17 @@ where
} }
} }
diagonal.scale_mut(m_amax); diagonal *= m_amax;
// Ensure all singular value are non-negative. // Ensure all singular value are non-negative.
for i in 0..dim { for i in 0..dim {
let sval = diagonal[i]; let sval = diagonal[i];
let (modulus, sign) = sval.to_exp();
if modulus != N::Real::zero() { if sval < N::Real::zero() {
diagonal[i] = N::from_real(modulus); diagonal[i] = -sval;
if let Some(ref mut u) = u { if let Some(ref mut u) = u {
u.column_mut(i).mul_assign(sign); u.column_mut(i).neg_mut();
} }
} }
} }
@ -285,55 +288,6 @@ where
}) })
} }
// Explicit formulaes inspired from the paper "Computing the Singular Values of 2-by-2 Complex
// Matrices", Sanzheng Qiao and Xiaohong Wang.
// http://www.cas.mcmaster.ca/sqrl/papers/sqrl5.pdf
fn compute_2x2_uptrig_svd(
m11: N,
m12: N,
m22: N,
compute_u: bool,
compute_v: bool,
) -> (Option<GivensRotation<N>>, Vector2<N>, Option<GivensRotation<N>>)
{
let two: N::Real = ::convert(2.0f64);
let half: N::Real = ::convert(0.5f64);
let denom = (m11 + m22).hypot(m12) + (m11 - m22).hypot(m12);
// NOTE: v1 is the singular value that is the closest to m22.
// This prevents cancellation issues when constructing the vector `csv` below. If we chose
// otherwise, we would have v1 ~= m11 when m12 is small. This would cause catastrophic
// cancellation on `v1 * v1 - m11 * m11` below.
let mut v1 = (m11 * m22).scale(two / denom);
let mut v2 = N::from_real(half * denom);
let mut u = None;
let mut v_t = None;
if compute_u || compute_v {
let (csv, sgn_v) = GivensRotation::new(m11 * m12, v1 * v1 - m11 * m11);
v1 *= sgn_v;
v2 *= sgn_v;
if compute_v {
v_t = Some(csv);
}
if compute_u {
let cu = (m11.scale(csv.c()) + m12 * csv.s()) / v1;
let su = (m22 * csv.s()) / v1;
let (csu, sgn_u) = GivensRotation::new(cu, su);
v1 *= sgn_u;
v2 *= sgn_u;
u = Some(csu);
}
}
(u, Vector2::new(v1, v2), v_t)
}
/* /*
fn display_bidiag(b: &Bidiagonal<N, R, C>, begin: usize, end: usize) { fn display_bidiag(b: &Bidiagonal<N, R, C>, begin: usize, end: usize) {
for i in begin .. end { for i in begin .. end {
@ -350,8 +304,8 @@ where
*/ */
fn delimit_subproblem( fn delimit_subproblem(
diagonal: &mut VectorN<N, DimMinimum<R, C>>, diagonal: &mut VectorN<N::Real, DimMinimum<R, C>>,
off_diagonal: &mut VectorN<N, DimDiff<DimMinimum<R, C>, U1>>, off_diagonal: &mut VectorN<N::Real, DimDiff<DimMinimum<R, C>, U1>>,
u: &mut Option<MatrixMN<N, R, DimMinimum<R, C>>>, u: &mut Option<MatrixMN<N, R, DimMinimum<R, C>>>,
v_t: &mut Option<MatrixMN<N, DimMinimum<R, C>, C>>, v_t: &mut Option<MatrixMN<N, DimMinimum<R, C>, C>>,
is_upper_diagonal: bool, is_upper_diagonal: bool,
@ -367,16 +321,16 @@ where
if off_diagonal[m].is_zero() if off_diagonal[m].is_zero()
|| off_diagonal[m].modulus() <= eps * (diagonal[n].modulus() + diagonal[m].modulus()) || off_diagonal[m].modulus() <= eps * (diagonal[n].modulus() + diagonal[m].modulus())
{ {
off_diagonal[m] = N::zero(); off_diagonal[m] = N::Real::zero();
} else if diagonal[m].modulus() <= eps { } else if diagonal[m].modulus() <= eps {
diagonal[m] = N::zero(); diagonal[m] = N::Real::zero();
Self::cancel_horizontal_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m, m + 1); Self::cancel_horizontal_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m, m + 1);
if m != 0 { if m != 0 {
Self::cancel_vertical_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m - 1); Self::cancel_vertical_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m - 1);
} }
} else if diagonal[n].modulus() <= eps { } else if diagonal[n].modulus() <= eps {
diagonal[n] = N::zero(); diagonal[n] = N::Real::zero();
Self::cancel_vertical_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m); Self::cancel_vertical_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m);
} else { } else {
break; break;
@ -395,12 +349,12 @@ where
if off_diagonal[m].modulus() <= eps * (diagonal[new_start].modulus() + diagonal[m].modulus()) if off_diagonal[m].modulus() <= eps * (diagonal[new_start].modulus() + diagonal[m].modulus())
{ {
off_diagonal[m] = N::zero(); off_diagonal[m] = N::Real::zero();
break; break;
} }
// FIXME: write a test that enters this case. // FIXME: write a test that enters this case.
else if diagonal[m].modulus() <= eps { else if diagonal[m].modulus() <= eps {
diagonal[m] = N::zero(); diagonal[m] = N::Real::zero();
Self::cancel_horizontal_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m, n); Self::cancel_horizontal_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m, n);
if m != 0 { if m != 0 {
@ -417,8 +371,8 @@ where
// Cancels the i-th off-diagonal element using givens rotations. // Cancels the i-th off-diagonal element using givens rotations.
fn cancel_horizontal_off_diagonal_elt( fn cancel_horizontal_off_diagonal_elt(
diagonal: &mut VectorN<N, DimMinimum<R, C>>, diagonal: &mut VectorN<N::Real, DimMinimum<R, C>>,
off_diagonal: &mut VectorN<N, DimDiff<DimMinimum<R, C>, U1>>, off_diagonal: &mut VectorN<N::Real, DimDiff<DimMinimum<R, C>, U1>>,
u: &mut Option<MatrixMN<N, R, DimMinimum<R, C>>>, u: &mut Option<MatrixMN<N, R, DimMinimum<R, C>>>,
v_t: &mut Option<MatrixMN<N, DimMinimum<R, C>, C>>, v_t: &mut Option<MatrixMN<N, DimMinimum<R, C>, C>>,
is_upper_diagonal: bool, is_upper_diagonal: bool,
@ -427,10 +381,11 @@ where
) )
{ {
let mut v = Vector2::new(off_diagonal[i], diagonal[i + 1]); let mut v = Vector2::new(off_diagonal[i], diagonal[i + 1]);
off_diagonal[i] = N::zero(); off_diagonal[i] = N::Real::zero();
for k in i..end { for k in i..end {
if let Some((rot, norm)) = GivensRotation::cancel_x(&v) { if let Some((rot, norm)) = GivensRotation::cancel_x(&v) {
let rot = GivensRotation::new_unchecked(rot.c(), N::from_real(rot.s()));
diagonal[k + 1] = norm; diagonal[k + 1] = norm;
if is_upper_diagonal { if is_upper_diagonal {
@ -443,9 +398,9 @@ where
} }
if k + 1 != end { if k + 1 != end {
v.x = -rot.s() * off_diagonal[k + 1]; v.x = -rot.s().real() * off_diagonal[k + 1];
v.y = diagonal[k + 2]; v.y = diagonal[k + 2];
off_diagonal[k + 1] = off_diagonal[k + 1].scale(rot.c()); off_diagonal[k + 1] *= rot.c();
} }
} else { } else {
break; break;
@ -455,8 +410,8 @@ where
// Cancels the i-th off-diagonal element using givens rotations. // Cancels the i-th off-diagonal element using givens rotations.
fn cancel_vertical_off_diagonal_elt( fn cancel_vertical_off_diagonal_elt(
diagonal: &mut VectorN<N, DimMinimum<R, C>>, diagonal: &mut VectorN<N::Real, DimMinimum<R, C>>,
off_diagonal: &mut VectorN<N, DimDiff<DimMinimum<R, C>, U1>>, off_diagonal: &mut VectorN<N::Real, DimDiff<DimMinimum<R, C>, U1>>,
u: &mut Option<MatrixMN<N, R, DimMinimum<R, C>>>, u: &mut Option<MatrixMN<N, R, DimMinimum<R, C>>>,
v_t: &mut Option<MatrixMN<N, DimMinimum<R, C>, C>>, v_t: &mut Option<MatrixMN<N, DimMinimum<R, C>, C>>,
is_upper_diagonal: bool, is_upper_diagonal: bool,
@ -464,10 +419,11 @@ where
) )
{ {
let mut v = Vector2::new(diagonal[i], off_diagonal[i]); let mut v = Vector2::new(diagonal[i], off_diagonal[i]);
off_diagonal[i] = N::zero(); off_diagonal[i] = N::Real::zero();
for k in (0..i + 1).rev() { for k in (0..i + 1).rev() {
if let Some((rot, norm)) = GivensRotation::cancel_y(&v) { if let Some((rot, norm)) = GivensRotation::cancel_y(&v) {
let rot = GivensRotation::new_unchecked(rot.c(), N::from_real(rot.s()));
diagonal[k] = norm; diagonal[k] = norm;
if is_upper_diagonal { if is_upper_diagonal {
@ -481,8 +437,8 @@ where
if k > 0 { if k > 0 {
v.x = diagonal[k - 1]; v.x = diagonal[k - 1];
v.y = rot.s() * off_diagonal[k - 1]; v.y = rot.s().real() * off_diagonal[k - 1];
off_diagonal[k - 1] = off_diagonal[k - 1].scale(rot.c()); off_diagonal[k - 1] *= rot.c();
} }
} else { } else {
break; break;
@ -497,7 +453,7 @@ where
eps >= N::Real::zero(), eps >= N::Real::zero(),
"SVD rank: the epsilon must be non-negative." "SVD rank: the epsilon must be non-negative."
); );
self.singular_values.iter().filter(|e| e.asum() > eps).count() self.singular_values.iter().filter(|e| **e > eps).count()
} }
/// Rebuild the original matrix. /// Rebuild the original matrix.
@ -510,7 +466,7 @@ where
(Some(mut u), Some(v_t)) => { (Some(mut u), Some(v_t)) => {
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).scale_mut(val);
} }
Ok(u * v_t) Ok(u * v_t)
} }
@ -536,10 +492,10 @@ where
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];
if val.asum() > eps { if val > eps {
self.singular_values[i] = N::one() / val; self.singular_values[i] = N::Real::one() / val;
} else { } else {
self.singular_values[i] = N::zero(); self.singular_values[i] = N::Real::zero();
} }
} }
@ -575,8 +531,8 @@ where
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];
if val.asum() > eps { if val > eps {
col[i] /= val; col[i] = col[i].unscale(val);
} else { } else {
col[i] = N::zero(); col[i] = N::zero();
} }
@ -602,7 +558,9 @@ where
+ Allocator<N, DimDiff<DimMinimum<R, C>, U1>> + Allocator<N, DimDiff<DimMinimum<R, C>, U1>>
+ Allocator<N, DimMinimum<R, C>, C> + Allocator<N, DimMinimum<R, C>, C>
+ Allocator<N, R, DimMinimum<R, C>> + Allocator<N, R, DimMinimum<R, C>>
+ Allocator<N, DimMinimum<R, C>>, + Allocator<N, DimMinimum<R, C>>
+ Allocator<N::Real, DimMinimum<R, C>>
+ Allocator<N::Real, DimDiff<DimMinimum<R, C>, U1>>,
{ {
/// Computes the Singular Value Decomposition using implicit shift. /// Computes the Singular Value Decomposition using implicit shift.
pub fn svd(self, compute_u: bool, compute_v: bool) -> SVD<N, R, C> { pub fn svd(self, compute_u: bool, compute_v: bool) -> SVD<N, R, C> {
@ -631,7 +589,7 @@ where
} }
/// Computes the singular values of this matrix. /// Computes the singular values of this matrix.
pub fn singular_values(&self) -> VectorN<N, DimMinimum<R, C>> { pub fn singular_values(&self) -> VectorN<N::Real, DimMinimum<R, C>> {
SVD::new(self.clone_owned(), false, false).singular_values SVD::new(self.clone_owned(), false, false).singular_values
} }
@ -653,3 +611,53 @@ where
SVD::new(self.clone_owned(), true, true).pseudo_inverse(eps) SVD::new(self.clone_owned(), true, true).pseudo_inverse(eps)
} }
} }
// Explicit formulae inspired from the paper "Computing the Singular Values of 2-by-2 Complex
// Matrices", Sanzheng Qiao and Xiaohong Wang.
// http://www.cas.mcmaster.ca/sqrl/papers/sqrl5.pdf
fn compute_2x2_uptrig_svd<N: Real>(
m11: N,
m12: N,
m22: N,
compute_u: bool,
compute_v: bool,
) -> (Option<GivensRotation<N>>, Vector2<N>, Option<GivensRotation<N>>)
{
let two: N::Real = ::convert(2.0f64);
let half: N::Real = ::convert(0.5f64);
let denom = (m11 + m22).hypot(m12) + (m11 - m22).hypot(m12);
// NOTE: v1 is the singular value that is the closest to m22.
// This prevents cancellation issues when constructing the vector `csv` below. If we chose
// otherwise, we would have v1 ~= m11 when m12 is small. This would cause catastrophic
// cancellation on `v1 * v1 - m11 * m11` below.
let mut v1 = m11 * m22 * two / denom;
let mut v2 = half * denom;
let mut u = None;
let mut v_t = None;
if compute_u || compute_v {
let (csv, sgn_v) = GivensRotation::new(m11 * m12, v1 * v1 - m11 * m11);
v1 *= sgn_v;
v2 *= sgn_v;
if compute_v {
v_t = Some(csv);
}
if compute_u {
let cu = (m11.scale(csv.c()) + m12 * csv.s()) / v1;
let su = (m22 * csv.s()) / v1;
let (csu, sgn_u) = GivensRotation::new(cu, su);
v1 *= sgn_u;
v2 *= sgn_u;
u = Some(csu);
}
}
(u, Vector2::new(v1, v2), v_t)
}

76
src/linalg/symmetric_eigen.rs
View File

@ -22,8 +22,8 @@ use linalg::SymmetricTridiagonal;
feature = "serde-serialize", feature = "serde-serialize",
serde(bound( serde(bound(
serialize = "DefaultAllocator: Allocator<N, D, D> + serialize = "DefaultAllocator: Allocator<N, D, D> +
Allocator<N, D>, Allocator<N::Real, D>,
VectorN<N, D>: Serialize, VectorN<N::Real, D>: Serialize,
MatrixN<N, D>: Serialize" MatrixN<N, D>: Serialize"
)) ))
)] )]
@ -31,32 +31,31 @@ use linalg::SymmetricTridiagonal;
feature = "serde-serialize", feature = "serde-serialize",
serde(bound( serde(bound(
deserialize = "DefaultAllocator: Allocator<N, D, D> + deserialize = "DefaultAllocator: Allocator<N, D, D> +
Allocator<N, D>, Allocator<N::Real, D>,
VectorN<N, D>: Deserialize<'de>, VectorN<N::Real, D>: Deserialize<'de>,
MatrixN<N, D>: Deserialize<'de>" MatrixN<N, D>: Deserialize<'de>"
)) ))
)] )]
#[derive(Clone, Debug)] #[derive(Clone, Debug)]
pub struct SymmetricEigen<N: Complex, D: Dim> pub struct SymmetricEigen<N: Complex, D: Dim>
where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D> where DefaultAllocator: Allocator<N, D, D> + Allocator<N::Real, D>
{ {
/// The eigenvectors of the decomposed matrix. /// The eigenvectors of the decomposed matrix.
pub eigenvectors: MatrixN<N, D>, pub eigenvectors: MatrixN<N, D>,
// FIXME: this should be a VectorN<N::Real, D>
/// The unsorted eigenvalues of the decomposed matrix. /// The unsorted eigenvalues of the decomposed matrix.
pub eigenvalues: VectorN<N, D>, pub eigenvalues: VectorN<N::Real, D>,
} }
impl<N: Complex, D: Dim> Copy for SymmetricEigen<N, D> impl<N: Complex, D: Dim> Copy for SymmetricEigen<N, D>
where where
DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>, DefaultAllocator: Allocator<N, D, D> + Allocator<N::Real, D>,
MatrixN<N, D>: Copy, MatrixN<N, D>: Copy,
VectorN<N, D>: Copy, VectorN<N::Real, D>: Copy,
{} {}
impl<N: Complex, D: Dim> SymmetricEigen<N, D> impl<N: Complex, D: Dim> SymmetricEigen<N, D>
where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D> where DefaultAllocator: Allocator<N, D, D> + Allocator<N::Real, D>
{ {
/// Computes the eigendecomposition of the given symmetric matrix. /// Computes the eigendecomposition of the given symmetric matrix.
/// ///
@ -64,7 +63,8 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
pub fn new(m: MatrixN<N, D>) -> Self pub fn new(m: MatrixN<N, D>) -> Self
where where
D: DimSub<U1>, D: DimSub<U1>,
DefaultAllocator: Allocator<N, DimDiff<D, U1>>, DefaultAllocator: Allocator<N, DimDiff<D, U1>> + // For tridiagonalization
Allocator<N::Real, DimDiff<D, U1>>,
{ {
Self::try_new(m, N::Real::default_epsilon(), 0).unwrap() Self::try_new(m, N::Real::default_epsilon(), 0).unwrap()
} }
@ -83,7 +83,8 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
pub fn try_new(m: MatrixN<N, D>, eps: N::Real, max_niter: usize) -> Option<Self> pub fn try_new(m: MatrixN<N, D>, eps: N::Real, max_niter: usize) -> Option<Self>
where where
D: DimSub<U1>, D: DimSub<U1>,
DefaultAllocator: Allocator<N, DimDiff<D, U1>>, DefaultAllocator: Allocator<N, DimDiff<D, U1>> + // For tridiagonalization
Allocator<N::Real, DimDiff<D, U1>>,
{ {
Self::do_decompose(m, true, eps, max_niter).map(|(vals, vecs)| SymmetricEigen { Self::do_decompose(m, true, eps, max_niter).map(|(vals, vecs)| SymmetricEigen {
eigenvectors: vecs.unwrap(), eigenvectors: vecs.unwrap(),
@ -96,10 +97,11 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
eigenvectors: bool, eigenvectors: bool,
eps: N::Real, eps: N::Real,
max_niter: usize, max_niter: usize,
) -> Option<(VectorN<N, D>, Option<MatrixN<N, D>>)> ) -> Option<(VectorN<N::Real, D>, Option<MatrixN<N, D>>)>
where where
D: DimSub<U1>, D: DimSub<U1>,
DefaultAllocator: Allocator<N, DimDiff<D, U1>>, DefaultAllocator: Allocator<N, DimDiff<D, U1>> + // For tridiagonalization
Allocator<N::Real, DimDiff<D, U1>>,
{ {
assert!( assert!(
m.is_square(), m.is_square(),
@ -158,41 +160,29 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
off_diag[i - 1] = norm; off_diag[i - 1] = norm;
} }
let mii = diag[i]; let mii = diag[i];
let mjj = diag[j]; let mjj = diag[j];
let mij = off_diag[i]; let mij = off_diag[i];
let cc = rot.c() * rot.c(); let cc = rot.c() * rot.c();
let ss = rot.s().modulus_squared(); // rot.s() * rot.s().conjugate() let ss = rot.s() * rot.s();
let cs = rot.s().scale(rot.c()); let cs = rot.c() * rot.s();
// b = cs * mij.conjugate() + cs.conjugate() * mij let b = cs * ::convert(2.0) * mij;
let b = N::from_real((cs * mij.conjugate()).real() * ::convert(2.0));
diag[i] = (mii.scale(cc) + mjj.scale(ss)) - b; diag[i] = (cc * mii + ss * mjj) - b;
diag[j] = (mii.scale(ss) + mjj.scale(cc)) + b; diag[j] = (ss * mii + cc * mjj) + b;
off_diag[i] = cs * (mii - mjj) + mij.scale(cc) - mij.conjugate() * rot.s() * rot.s(); off_diag[i] = cs * (mii - mjj) + mij * (cc - ss);
let mut mat = Matrix2::new(
mii, mij.conjugate(),
mij, mjj);
println!("The mat before rotate: {:.5}", mat);
println!("The v before rotate: {:.5?}", v);
rot.rotate(&mut mat);
rot.inverse().rotate_rows(&mut mat);
let mut v2 = v.clone();
rot.rotate(&mut v2);
println!("The v: {:.5?}", v2);
println!("The mat: {:.5}", mat);
println!("Its components: {:.5}, {:.5}, {:.5}", diag[i], diag[j], off_diag[i]);
if i != n - 1 { if i != n - 1 {
v.x = off_diag[i]; v.x = off_diag[i];
v.y = -rot.s() * off_diag[i + 1]; v.y = -rot.s() * off_diag[i + 1];
off_diag[i + 1] = off_diag[i + 1].scale(rot.c()); off_diag[i + 1] *= rot.c();
} }
if let Some(ref mut q) = q { if let Some(ref mut q) = q {
let rot = GivensRotation::new_unchecked(rot.c(), N::from_real(rot.s()));
rot.inverse().rotate_rows(&mut q.fixed_columns_mut::<U2>(i)); rot.inverse().rotate_rows(&mut q.fixed_columns_mut::<U2>(i));
} }
} else { } else {
@ -220,6 +210,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
if let Some(ref mut q) = q { if let Some(ref mut q) = q {
if let Some((rot, _)) = GivensRotation::try_new(basis.x, basis.y, eps) { if let Some((rot, _)) = GivensRotation::try_new(basis.x, basis.y, eps) {
let rot = GivensRotation::new_unchecked(rot.c(), N::from_real(rot.s()));
rot.rotate_rows(&mut q.fixed_columns_mut::<U2>(start)); rot.rotate_rows(&mut q.fixed_columns_mut::<U2>(start));
} }
} }
@ -245,14 +236,14 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
} }
fn delimit_subproblem( fn delimit_subproblem(
diag: &VectorN<N, D>, diag: &VectorN<N::Real, D>,
off_diag: &mut VectorN<N, DimDiff<D, U1>>, off_diag: &mut VectorN<N::Real, DimDiff<D, U1>>,
end: usize, end: usize,
eps: N::Real, eps: N::Real,
) -> (usize, usize) ) -> (usize, usize)
where where
D: DimSub<U1>, D: DimSub<U1>,
DefaultAllocator: Allocator<N, DimDiff<D, U1>>, DefaultAllocator: Allocator<N::Real, DimDiff<D, U1>>,
{ {
let mut n = end; let mut n = end;
@ -277,7 +268,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
if off_diag[m].is_zero() if off_diag[m].is_zero()
|| off_diag[m].modulus() <= eps * (diag[new_start].modulus() + diag[m].modulus()) || off_diag[m].modulus() <= eps * (diag[new_start].modulus() + diag[m].modulus())
{ {
off_diag[m] = N::zero(); off_diag[m] = N::Real::zero();
break; break;
} }
@ -294,7 +285,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
let mut u_t = self.eigenvectors.clone(); let mut u_t = self.eigenvectors.clone();
for i in 0..self.eigenvalues.len() { for i in 0..self.eigenvalues.len() {
let val = self.eigenvalues[i]; let val = self.eigenvalues[i];
u_t.column_mut(i).mul_assign(val); u_t.column_mut(i).scale_mut(val);
} }
u_t.conjugate_transpose_mut(); u_t.conjugate_transpose_mut();
&self.eigenvectors * u_t &self.eigenvectors * u_t
@ -324,7 +315,8 @@ pub fn wilkinson_shift<N: Complex>(tmm: N, tnn: N, tmn: N) -> N {
* *
*/ */
impl<N: Complex, D: DimSub<U1>, S: Storage<N, D, D>> SquareMatrix<N, D, S> impl<N: Complex, D: DimSub<U1>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D> + Allocator<N, DimDiff<D, U1>> where DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>> +
Allocator<N::Real, D> + Allocator<N::Real, DimDiff<D, U1>>
{ {
/// Computes the eigendecomposition of this symmetric matrix. /// Computes the eigendecomposition of this symmetric matrix.
/// ///
@ -351,7 +343,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D> + Allocator<N, DimD
/// Computes the eigenvalues of this symmetric matrix. /// Computes the eigenvalues of this symmetric matrix.
/// ///
/// Only the lower-triangular part of the matrix is read. /// Only the lower-triangular part of the matrix is read.
pub fn symmetric_eigenvalues(&self) -> VectorN<N, D> { pub fn symmetric_eigenvalues(&self) -> VectorN<N::Real, D> {
SymmetricEigen::do_decompose(self.clone_owned(), false, N::Real::default_epsilon(), 0) SymmetricEigen::do_decompose(self.clone_owned(), false, N::Real::default_epsilon(), 0)
.unwrap() .unwrap()
.0 .0

26
src/linalg/symmetric_tridiagonal.rs
View File

@ -102,30 +102,30 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>
/// Retrieve the orthogonal transformation, diagonal, and off diagonal elements of this /// Retrieve the orthogonal transformation, diagonal, and off diagonal elements of this
/// decomposition. /// decomposition.
pub fn unpack(self) -> (MatrixN<N, D>, VectorN<N, D>, VectorN<N, DimDiff<D, U1>>) pub fn unpack(self) -> (MatrixN<N, D>, VectorN<N::Real, D>, VectorN<N::Real, DimDiff<D, U1>>)
where DefaultAllocator: Allocator<N, D> { where DefaultAllocator: Allocator<N::Real, D>
+ Allocator<N::Real, DimDiff<D, U1>> {
let diag = self.diagonal(); let diag = self.diagonal();
let q = self.q(); let q = self.q();
(q, diag, self.off_diagonal.apply_into(|e| N::from_real(e.modulus()))) (q, diag, self.off_diagonal.map(N::modulus))
} }
/// Retrieve the diagonal, and off diagonal elements of this decomposition. /// Retrieve the diagonal, and off diagonal elements of this decomposition.
pub fn unpack_tridiagonal(mut self) -> (VectorN<N, D>, VectorN<N, DimDiff<D, U1>>) pub fn unpack_tridiagonal(mut self) -> (VectorN<N::Real, D>, VectorN<N::Real, DimDiff<D, U1>>)
where DefaultAllocator: Allocator<N, D> { where DefaultAllocator: Allocator<N::Real, D>
let diag = self.diagonal(); + Allocator<N::Real, DimDiff<D, U1>> {
(self.diagonal(), self.off_diagonal.map(N::modulus))
(diag, self.off_diagonal.apply_into(|e| N::from_real(e.modulus())))
} }
/// The diagonal components of this decomposition. /// The diagonal components of this decomposition.
pub fn diagonal(&self) -> VectorN<N, D> pub fn diagonal(&self) -> VectorN<N::Real, D>
where DefaultAllocator: Allocator<N, D> { self.tri.diagonal() } where DefaultAllocator: Allocator<N::Real, D> { self.tri.map_diagonal(|e| e.real()) }
/// The off-diagonal components of this decomposition. /// The off-diagonal components of this decomposition.
pub fn off_diagonal(&self) -> VectorN<N, DimDiff<D, U1>> pub fn off_diagonal(&self) -> VectorN<N::Real, DimDiff<D, U1>>
where DefaultAllocator: Allocator<N, D> { where DefaultAllocator: Allocator<N::Real, DimDiff<D, U1>> {
self.off_diagonal.map(|e| N::from_real(e.modulus())) self.off_diagonal.map(N::modulus)
} }
/// Computes the orthogonal matrix `Q` of this decomposition. /// Computes the orthogonal matrix `Q` of this decomposition.

28
tests/linalg/svd.rs
View File

@ -18,11 +18,11 @@ mod quickcheck_tests {
let svd = m.clone().svd(true, true); let svd = m.clone().svd(true, true);
let recomp_m = svd.clone().recompose().unwrap(); 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.map(|e| Complex::from_real(e)));
println!("{}{}", &m, &u * &ds * &v_t); println!("{}{}", &m, &u * &ds * &v_t);
s.iter().all(|e| e.real() >= 0.0) && s.iter().all(|e| *e >= 0.0) &&
relative_eq!(&u * ds * &v_t, recomp_m, epsilon = 1.0e-5) && relative_eq!(&u * ds * &v_t, recomp_m, epsilon = 1.0e-5) &&
relative_eq!(m, recomp_m, epsilon = 1.0e-5) relative_eq!(m, recomp_m, epsilon = 1.0e-5)
} }
@ -35,9 +35,9 @@ mod quickcheck_tests {
let m = m.map(|e| e.0); let m = m.map(|e| e.0);
let svd = m.svd(true, true); let svd = m.svd(true, true);
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 = Matrix3::from_diagonal(&s); let ds = Matrix3::from_diagonal(&s.map(|e| Complex::from_real(e)));
s.iter().all(|e| e.real() >= 0.0) && s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5) && relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5) &&
u.is_orthogonal(1.0e-5) && u.is_orthogonal(1.0e-5) &&
v_t.is_orthogonal(1.0e-5) v_t.is_orthogonal(1.0e-5)
@ -47,9 +47,9 @@ mod quickcheck_tests {
let m = m.map(|e| e.0); let m = m.map(|e| e.0);
let svd = m.svd(true, true); let svd = m.svd(true, true);
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 = Matrix2::from_diagonal(&s); let ds = Matrix2::from_diagonal(&s.map(|e| Complex::from_real(e)));
s.iter().all(|e| e.real() >= 0.0) && s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5) && relative_eq!(m, &u * ds * &v_t, epsilon = 1.0e-5) &&
u.is_orthogonal(1.0e-5) && u.is_orthogonal(1.0e-5) &&
v_t.is_orthogonal(1.0e-5) v_t.is_orthogonal(1.0e-5)
@ -60,9 +60,9 @@ mod quickcheck_tests {
let svd = m.svd(true, true); let svd = m.svd(true, true);
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 = Matrix3::from_diagonal(&s); let ds = Matrix3::from_diagonal(&s.map(|e| Complex::from_real(e)));
s.iter().all(|e| e.real() >= 0.0) && s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)
} }
@ -70,9 +70,9 @@ mod quickcheck_tests {
let m = m.map(|e| e.0); let m = m.map(|e| e.0);
let svd = m.svd(true, true); let svd = m.svd(true, true);
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 = Matrix2::from_diagonal(&s); let ds = Matrix2::from_diagonal(&s.map(|e| Complex::from_real(e)));
s.iter().all(|e| e.real() >= 0.0) && s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)
} }
@ -80,9 +80,9 @@ mod quickcheck_tests {
let m = m.map(|e| e.0); let m = m.map(|e| e.0);
let svd = m.svd(true, true); let svd = m.svd(true, true);
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 = Matrix4::from_diagonal(&s); let ds = Matrix4::from_diagonal(&s.map(|e| Complex::from_real(e)));
s.iter().all(|e| e.real() >= 0.0) && s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) && relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) &&
u.is_orthogonal(1.0e-5) && u.is_orthogonal(1.0e-5) &&
v_t.is_orthogonal(1.0e-5) v_t.is_orthogonal(1.0e-5)
@ -92,14 +92,14 @@ mod quickcheck_tests {
let m = m.map(|e| e.0); let m = m.map(|e| e.0);
let svd = m.svd(true, true); let svd = m.svd(true, true);
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 = Matrix2::from_diagonal(&s); let ds = Matrix2::from_diagonal(&s.map(|e| Complex::from_real(e)));
println!("u, s, v_t: {}{}{}", u, s, v_t); println!("u, s, v_t: {}{}{}", u, s, v_t);
println!("m: {}", m); println!("m: {}", m);
println!("recomp: {}", u * ds * v_t); println!("recomp: {}", u * ds * v_t);
println!("uu_t, vv_t: {}{}", u * u.conjugate_transpose(), v_t.conjugate_transpose() * v_t); println!("uu_t, vv_t: {}{}", u * u.conjugate_transpose(), v_t.conjugate_transpose() * v_t);
s.iter().all(|e| e.re >= 0.0) && s.iter().all(|e| *e >= 0.0) &&
relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) && relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) &&
u.is_orthogonal(1.0e-5) && u.is_orthogonal(1.0e-5) &&
v_t.is_orthogonal(1.0e-5) v_t.is_orthogonal(1.0e-5)