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

Start fixing SVD.

This commit is contained in:
sebcrozet 2019-03-18 11:23:19 +01:00
parent 010c009cff
commit e4748c69ce
13 changed files with 248 additions and 206 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/geometry/reflection.rs
View File

@ -61,10 +61,10 @@ impl<N: Complex, D: Dim, S: Storage<N, D>> Reflection<N, D, S> {
}
}
/// Applies the reflection to the rows of `rhs`.
/// Applies the reflection to the rows of `lhs`.
pub fn reflect_rows<R2: Dim, C2: Dim, S2, S3>(
&self,
rhs: &mut Matrix<N, R2, C2, S2>,
lhs: &mut Matrix<N, R2, C2, S2>,
work: &mut Vector<N, R2, S3>,
) where
S2: StorageMut<N, R2, C2>,
@ -72,13 +72,13 @@ impl<N: Complex, D: Dim, S: Storage<N, D>> Reflection<N, D, S> {
ShapeConstraint: DimEq<C2, D> + AreMultipliable<R2, C2, D, U1>,
DefaultAllocator: Allocator<N, D>
{
rhs.mul_to(&self.axis, work);
lhs.mul_to(&self.axis, work);
if !self.bias.is_zero() {
work.add_scalar_mut(-self.bias);
}
let m_two: N = ::convert(-2.0f64);
rhs.ger(m_two, &work, &self.axis.conjugate(), N::one());
lhs.ger(m_two, &work, &self.axis.conjugate(), N::one());
}
}

20
src/linalg/bidiagonal.rs
View File

@ -49,9 +49,9 @@ where
// contiguous. This prevents some useless copies.
uv: MatrixMN<N, R, C>,
/// The diagonal elements of the decomposed matrix.
pub diagonal: VectorN<N, DimMinimum<R, C>>,
diagonal: VectorN<N, DimMinimum<R, C>>,
/// The off-diagonal elements of the decomposed matrix.
pub off_diagonal: VectorN<N, DimDiff<DimMinimum<R, C>, U1>>,
off_diagonal: VectorN<N, DimDiff<DimMinimum<R, C>, U1>>,
upper_diagonal: bool,
}
@ -143,9 +143,9 @@ where
Bidiagonal {
uv: matrix,
diagonal: diagonal,
off_diagonal: off_diagonal,
upper_diagonal: upper_diagonal,
diagonal,
off_diagonal,
upper_diagonal,
}
}
@ -231,7 +231,7 @@ where
res
}
/// Computes the orthogonal matrix `V` of this `U * D * V` decomposition.
/// Computes the orthogonal matrix `V_t` of this `U * D * V_t` decomposition.
pub fn v_t(&self) -> MatrixMN<N, DimMinimum<R, C>, C>
where DefaultAllocator: Allocator<N, DimMinimum<R, C>, C> {
let (nrows, ncols) = self.uv.data.shape();
@ -258,13 +258,13 @@ where
}
/// The diagonal part of this decomposed matrix.
pub fn diagonal(&self) -> &VectorN<N, DimMinimum<R, C>> {
&self.diagonal
pub fn diagonal(&self) -> VectorN<N, DimMinimum<R, C>> {
self.diagonal.map(|e| N::from_real(e.modulus()))
}
/// The off-diagonal part of this decomposed matrix.
pub fn off_diagonal(&self) -> &VectorN<N, DimDiff<DimMinimum<R, C>, U1>> {
&self.off_diagonal
pub fn off_diagonal(&self) -> VectorN<N, DimDiff<DimMinimum<R, C>, U1>> {
self.off_diagonal.map(|e| N::from_real(e.modulus()))
}
#[doc(hidden)]

83
src/linalg/givens.rs
View File

@ -1,7 +1,7 @@
//! Construction of givens rotations.
use alga::general::{Complex, Real};
use num::Zero;
use num::{Zero, One};
use num_complex::Complex as NumComplex;
use base::dimension::{Dim, U2};
@ -12,60 +12,37 @@ use base::{Vector, Matrix};
use geometry::UnitComplex;
/// A Givens rotation.
#[derive(Debug)]
#[derive(Debug, Clone, Copy)]
pub struct GivensRotation<N: Complex> {
// FIXME: c should be a `N::Real`.
c: N,
c: N::Real,
s: N
}
// XXX: remove this
/// Computes the rotation `R` required such that the `y` component of `R * v` is zero.
///
/// Returns `None` if no rotation is needed (i.e. if `v.y == 0`). Otherwise, this returns the norm
/// of `v` and the rotation `r` such that `R * v = [ |v|, 0.0 ]^t` where `|v|` is the norm of `v`.
pub fn cancel_y<N: Real, S: Storage<N, U2>>(v: &Vector<N, U2, S>) -> Option<(UnitComplex<N>, N)> {
if !v[1].is_zero() {
let c = NumComplex::new(v[0], -v[1]);
Some(UnitComplex::from_complex_and_get(c))
} else {
None
}
}
// XXX: remove this
/// Computes the rotation `R` required such that the `x` component of `R * v` is zero.
///
/// Returns `None` if no rotation is needed (i.e. if `v.x == 0`). Otherwise, this returns the norm
/// of `v` and the rotation `r` such that `R * v = [ 0.0, |v| ]^t` where `|v|` is the norm of `v`.
pub fn cancel_x<N: Real, S: Storage<N, U2>>(v: &Vector<N, U2, S>) -> Option<(UnitComplex<N>, N)> {
if !v[0].is_zero() {
let c = NumComplex::new(v[1], v[0]);
Some(UnitComplex::from_complex_and_get(c))
} else {
None
}
}
// Matrix = UnitComplex * Matrix
impl<N: Complex> GivensRotation<N> {
/// The Givents rotation that does nothing.
pub fn identity() -> Self {
Self {
c: N::Real::one(),
s: N::zero()
}
}
/// Initializes a Givens rotation from its non-normalized cosine an sine components.
pub fn new(c: N, s: N) -> Self {
let res = Self::try_new(c, s, N::Real::zero()).unwrap();
println!("The rot: {:?}", res);
res
pub fn new(c: N, s: N) -> (Self, N) {
Self::try_new(c, s, N::Real::zero()).unwrap()
}
/// Initializes a Givens rotation form its non-normalized cosine an sine components.
pub fn try_new(c: N, s: N, eps: N::Real) -> Option<Self> {
pub fn try_new(c: N, s: N, eps: N::Real) -> Option<(Self, N)> {
let (mod0, sign0) = c.to_exp();
let denom = (mod0 * mod0 + s.modulus_squared()).sqrt();
if denom > eps {
let c = N::from_real(mod0 / denom);
let s = s / sign0.scale(denom);
Some(Self { c, s })
let norm = sign0.scale(denom);
let c = mod0 / denom;
let s = s / norm;
Some((Self { c, s }, norm))
} else {
None
}
@ -79,8 +56,8 @@ impl<N: Complex> GivensRotation<N> {
if !v[1].is_zero() {
let (mod0, sign0) = v[0].to_exp();
let denom = (mod0 * mod0 + v[1].modulus_squared()).sqrt();
let c = N::from_real(mod0 / denom);
let s = (sign0 * v[1].conjugate()).unscale(-denom);
let c = mod0 / denom;
let s = -v[1] / sign0.scale(denom);
let r = sign0.scale(denom);
Some((Self { c, s }, r))
} else {
@ -94,11 +71,11 @@ impl<N: Complex> GivensRotation<N> {
/// of `v` and the rotation `r` such that `R * v = [ 0.0, |v| ]^t` where `|v|` is the norm of `v`.
pub fn cancel_x<S: Storage<N, U2>>(v: &Vector<N, U2, S>) -> Option<(Self, N)> {
if !v[0].is_zero() {
let (mod0, sign0) = v[0].to_exp();
let denom = (mod0 * mod0 + v[1].modulus_squared()).sqrt();
let c = N::from_real(mod0 / denom);
let s = (sign0 * v[1].conjugate()).unscale(denom);
let r = sign0.scale(denom);
let (mod1, sign1) = v[1].to_exp();
let denom = (mod1 * mod1 + v[0].modulus_squared()).sqrt();
let c = mod1 / denom;
let s = (v[0].conjugate() * sign1).unscale(denom);
let r = sign1.scale(denom);
Some((Self { c, s }, r))
} else {
None
@ -106,7 +83,7 @@ impl<N: Complex> GivensRotation<N> {
}
/// The cos part of this roration.
pub fn c(&self) -> N {
pub fn c(&self) -> N::Real {
self.c
}
@ -140,8 +117,8 @@ impl<N: Complex> GivensRotation<N> {
let a = *rhs.get_unchecked((0, j));
let b = *rhs.get_unchecked((1, j));
*rhs.get_unchecked_mut((0, j)) = c * a + -s.conjugate() * b;
*rhs.get_unchecked_mut((1, j)) = s * a + c * b;
*rhs.get_unchecked_mut((0, j)) = a.scale(c) - s.conjugate() * b;
*rhs.get_unchecked_mut((1, j)) = s * a + b.scale(c);
}
}
}
@ -167,8 +144,8 @@ impl<N: Complex> GivensRotation<N> {
let a = *lhs.get_unchecked((j, 0));
let b = *lhs.get_unchecked((j, 1));
*lhs.get_unchecked_mut((j, 0)) = c * a + s * b;
*lhs.get_unchecked_mut((j, 1)) = -s.conjugate() * a + c * b;
*lhs.get_unchecked_mut((j, 0)) = a.scale(c) + s * b;
*lhs.get_unchecked_mut((j, 1)) = -s.conjugate() * a + b.scale(c);
}
}
}

15
src/linalg/householder.rs
View File

@ -23,20 +23,20 @@ pub fn reflection_axis_mut<N: Complex, D: Dim, S: StorageMut<N, D>>(
let reflection_norm = reflection_sq_norm.sqrt();
let factor;
let scaled_norm;
let signed_norm;
unsafe {
let (modulus, exp) = column.vget_unchecked(0).to_exp();
scaled_norm = exp.scale(reflection_norm);
let (modulus, sign) = column.vget_unchecked(0).to_exp();
signed_norm = sign.scale(reflection_norm);
factor = (reflection_sq_norm + modulus * reflection_norm) * ::convert(2.0);
*column.vget_unchecked_mut(0) += scaled_norm;
*column.vget_unchecked_mut(0) += signed_norm;
};
if !factor.is_zero() {
column.unscale_mut(factor.sqrt());
(-scaled_norm, true)
(-signed_norm, true)
} else {
(-scaled_norm, false)
(-signed_norm, false)
}
}
@ -67,7 +67,7 @@ pub fn clear_column_unchecked<N: Complex, R: Dim, C: Dim>(
}
}
/// Uses an hoseholder reflection to zero out the `irow`-th row, ending before the `shift + 1`-th
/// Uses an householder reflection to zero out the `irow`-th row, ending before the `shift + 1`-th
/// superdiagonal element.
#[doc(hidden)]
pub fn clear_row_unchecked<N: Complex, R: Dim, C: Dim>(
@ -85,6 +85,7 @@ pub fn clear_row_unchecked<N: Complex, R: Dim, C: Dim>(
axis.tr_copy_from(&top.columns_range(irow + shift..));
let (reflection_norm, not_zero) = reflection_axis_mut(&mut axis);
axis.conjugate_mut(); // So that reflect_rows actually cancels the first row.
*diag_elt = reflection_norm;
if not_zero {

2
src/linalg/mod.rs
View File

@ -30,6 +30,6 @@ pub use self::lu::*;
pub use self::permutation_sequence::*;
pub use self::qr::*;
pub use self::schur::*;
pub use self::svd::*;
//pub use self::svd::*;
pub use self::symmetric_eigen::*;
pub use self::symmetric_tridiagonal::*;

7
src/linalg/schur.rs
View File

@ -417,10 +417,11 @@ where
if compute_q {
// XXX: we have to build the matrix manually because
// rot.to_rotation_matrix().unwrap() causes an ICE.
let c = N::from_real(rot.c());
q = Some(MatrixN::from_column_slice_generic(
dim,
dim,
&[rot.c(), rot.s(), -rot.s().conjugate(), rot.c()],
&[c, rot.s(), -rot.s().conjugate(), c],
));
}
}
@ -479,9 +480,9 @@ fn compute_2x2_basis<N: Complex, S: Storage<N, U2, U2>>(
// This is necessary for numerical stability of the normalization of the complex
// number.
if x1.modulus() > x2.modulus() {
Some(GivensRotation::new(x1, h10))
Some(GivensRotation::new(x1, h10).0)
} else {
Some(GivensRotation::new(x2, h10))
Some(GivensRotation::new(x2, h10).0)
}
} else {
None

224
src/linalg/svd.rs
View File

@ -1,20 +1,22 @@
#[cfg(feature = "serde-serialize")]
use serde::{Deserialize, Serialize};
use num_complex::Complex;
use num_complex::Complex as NumComplex;
use num::Zero;
use std::ops::MulAssign;
use approx::AbsDiffEq;
use alga::general::Real;
use alga::general::Complex;
use allocator::Allocator;
use base::{DefaultAllocator, Matrix, Matrix2x3, MatrixMN, Vector2, VectorN};
use constraint::{SameNumberOfRows, ShapeConstraint};
use dimension::{Dim, DimDiff, DimMin, DimMinimum, DimSub, U1, U2};
use storage::Storage;
use geometry::UnitComplex;
use linalg::givens;
use linalg::symmetric_eigen;
use linalg::Bidiagonal;
use linalg::givens::GivensRotation;
/// Singular Value Decomposition of a general matrix.
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
@ -43,7 +45,7 @@ use linalg::Bidiagonal;
))
)]
#[derive(Clone, Debug)]
pub struct SVD<N: Real, R: DimMin<C>, C: Dim>
pub struct SVD<N: Complex, R: DimMin<C>, C: Dim>
where DefaultAllocator: Allocator<N, DimMinimum<R, C>, C>
+ Allocator<N, R, DimMinimum<R, C>>
+ Allocator<N, DimMinimum<R, C>>
@ -52,11 +54,13 @@ where DefaultAllocator: Allocator<N, DimMinimum<R, C>, C>
pub u: Option<MatrixMN<N, R, DimMinimum<R, C>>>,
/// The right-singular vectors `V^t` of this SVD.
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.
pub singular_values: VectorN<N, DimMinimum<R, C>>,
}
impl<N: Real, 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
DefaultAllocator: Allocator<N, DimMinimum<R, C>, C>
+ Allocator<N, R, DimMinimum<R, C>>
@ -66,7 +70,7 @@ where
VectorN<N, DimMinimum<R, C>>: Copy,
{}
impl<N: Real, R: DimMin<C>, C: Dim> SVD<N, R, C>
impl<N: Complex, R: DimMin<C>, C: Dim> SVD<N, R, C>
where
DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
DefaultAllocator: Allocator<N, R, C>
@ -79,7 +83,7 @@ where
{
/// 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 {
Self::try_new(matrix, compute_u, compute_v, N::default_epsilon(), 0).unwrap()
Self::try_new(matrix, compute_u, compute_v, N::Real::default_epsilon(), 0).unwrap()
}
/// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift.
@ -96,7 +100,7 @@ where
mut matrix: MatrixMN<N, R, C>,
compute_u: bool,
compute_v: bool,
eps: N,
eps: N::Real,
max_niter: usize,
) -> Option<Self>
{
@ -108,18 +112,20 @@ where
let min_nrows_ncols = nrows.min(ncols);
let dim = min_nrows_ncols.value();
let m_amax = matrix.amax();
let m_amax = matrix.camax();
if !m_amax.is_zero() {
matrix /= m_amax;
matrix.unscale_mut(m_amax);
}
let mut b = Bidiagonal::new(matrix);
let mut u = if compute_u { Some(b.u()) } else { None };
let mut v_t = if compute_v { Some(b.v_t()) } else { None };
let mut diagonal = b.diagonal();
let mut off_diagonal = b.off_diagonal();
let mut niter = 0;
let (mut start, mut end) = Self::delimit_subproblem(&mut b, &mut u, &mut v_t, dim - 1, eps);
let (mut start, mut end) = Self::delimit_subproblem(&mut diagonal, &mut off_diagonal, &mut u, &mut v_t, b.is_upper_diagonal(), dim - 1, eps);
while end != start {
let subdim = end - start + 1;
@ -131,19 +137,19 @@ where
let mut vec;
{
let dm = b.diagonal[m];
let dn = b.diagonal[n];
let fm = b.off_diagonal[m];
let dm = diagonal[m];
let dn = diagonal[n];
let fm = off_diagonal[m];
let tmm = dm * dm + b.off_diagonal[m - 1] * b.off_diagonal[m - 1];
let tmm = dm * dm + off_diagonal[m - 1] * off_diagonal[m - 1];
let tmn = dm * fm;
let tnn = dn * dn + fm * fm;
let shift = symmetric_eigen::wilkinson_shift(tmm, tnn, tmn);
vec = Vector2::new(
b.diagonal[start] * b.diagonal[start] - shift,
b.diagonal[start] * b.off_diagonal[start],
diagonal[start] * diagonal[start] - shift,
diagonal[start] * off_diagonal[start],
);
}
@ -151,31 +157,31 @@ where
let m12 = if k == n - 1 {
N::zero()
} else {
b.off_diagonal[k + 1]
off_diagonal[k + 1]
};
let mut subm = Matrix2x3::new(
b.diagonal[k],
b.off_diagonal[k],
diagonal[k],
off_diagonal[k],
N::zero(),
N::zero(),
b.diagonal[k + 1],
diagonal[k + 1],
m12,
);
if let Some((rot1, norm1)) = givens::cancel_y(&vec) {
rot1.conjugate()
if let Some((rot1, norm1)) = GivensRotation::cancel_y(&vec) {
rot1.inverse()
.rotate_rows(&mut subm.fixed_columns_mut::<U2>(0));
if k > start {
// This is not the first iteration.
b.off_diagonal[k - 1] = norm1;
off_diagonal[k - 1] = norm1;
}
let v = Vector2::new(subm[(0, 0)], subm[(1, 0)]);
// FIXME: does the case `v.y == 0` ever happen?
let (rot2, norm2) =
givens::cancel_y(&v).unwrap_or((UnitComplex::identity(), subm[(0, 0)]));
GivensRotation::cancel_y(&v).unwrap_or((GivensRotation::identity(), subm[(0, 0)]));
rot2.rotate(&mut subm.fixed_columns_mut::<U2>(1));
subm[(0, 0)] = norm2;
@ -197,12 +203,12 @@ where
}
}
b.diagonal[k + 0] = subm[(0, 0)];
b.diagonal[k + 1] = subm[(1, 1)];
b.off_diagonal[k + 0] = subm[(0, 1)];
diagonal[k + 0] = subm[(0, 0)];
diagonal[k + 1] = subm[(1, 1)];
off_diagonal[k + 0] = subm[(0, 1)];
if k != n - 1 {
b.off_diagonal[k + 1] = subm[(1, 2)];
off_diagonal[k + 1] = subm[(1, 2)];
}
vec.x = subm[(0, 1)];
@ -214,16 +220,16 @@ where
} else if subdim == 2 {
// Solve the remaining 2x2 subproblem.
let (u2, s, v2) = Self::compute_2x2_uptrig_svd(
b.diagonal[start],
b.off_diagonal[start],
b.diagonal[start + 1],
diagonal[start],
off_diagonal[start],
diagonal[start + 1],
compute_u && b.is_upper_diagonal() || compute_v && !b.is_upper_diagonal(),
compute_v && b.is_upper_diagonal() || compute_u && !b.is_upper_diagonal(),
);
b.diagonal[start + 0] = s[0];
b.diagonal[start + 1] = s[1];
b.off_diagonal[start] = N::zero();
diagonal[start + 0] = s[0];
diagonal[start + 1] = s[1];
off_diagonal[start] = N::zero();
if let Some(ref mut u) = u {
let rot = if b.is_upper_diagonal() {
@ -247,7 +253,7 @@ where
}
// Re-delimit the subproblem in case some decoupling occurred.
let sub = Self::delimit_subproblem(&mut b, &mut u, &mut v_t, end, eps);
let sub = Self::delimit_subproblem(&mut diagonal, &mut off_diagonal, &mut u, &mut v_t, b.is_upper_diagonal(), end, eps);
start = sub.0;
end = sub.1;
@ -257,24 +263,25 @@ where
}
}
b.diagonal *= m_amax;
diagonal.scale_mut(m_amax);
// Ensure all singular value are non-negative.
for i in 0..dim {
let sval = b.diagonal[i];
if sval < N::zero() {
b.diagonal[i] = -sval;
let sval = diagonal[i];
let (modulus, sign) = sval.to_exp();
if modulus != N::Real::zero() {
diagonal[i] = N::from_real(modulus);
if let Some(ref mut u) = u {
u.column_mut(i).neg_mut();
u.column_mut(i).mul_assign(sign);
}
}
}
Some(Self {
u: u,
v_t: v_t,
singular_values: b.diagonal,
u,
v_t,
singular_values: diagonal,
})
}
@ -287,10 +294,10 @@ where
m22: N,
compute_u: bool,
compute_v: bool,
) -> (Option<UnitComplex<N>>, Vector2<N>, Option<UnitComplex<N>>)
) -> (Option<GivensRotation<N>>, Vector2<N>, Option<GivensRotation<N>>)
{
let two: N = ::convert(2.0f64);
let half: N = ::convert(0.5f64);
let two: N::Real = ::convert(2.0f64);
let half: N::Real = ::convert(0.5f64);
let denom = (m11 + m22).hypot(m12) + (m11 - m22).hypot(m12);
@ -298,24 +305,29 @@ where
// 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 v1 = two * m11 * m22 / denom;
let v2 = half * denom;
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 = Vector2::new(m11 * m12, v1 * v1 - m11 * m11).normalize();
let (csv, sgn_v) = GivensRotation::new(m11 * m12, v1 * v1 - m11 * m11);
v1 *= sgn_v;
v2 *= sgn_v;
if compute_v {
v_t = Some(UnitComplex::new_unchecked(Complex::new(csv.x, csv.y)));
v_t = Some(csv);
}
if compute_u {
let cu = (m11 * csv.x + m12 * csv.y) / v1;
let su = (m22 * csv.y) / v1;
let cu = (m11.scale(csv.c()) + m12 * csv.s()) / v1;
let su = (m22 * csv.s()) / v1;
let (csu, sgn_u) = GivensRotation::new(cu, su);
u = Some(UnitComplex::new_unchecked(Complex::new(cu, su)));
v1 *= sgn_u;
v2 *= sgn_u;
u = Some(csu);
}
}
@ -338,11 +350,13 @@ where
*/
fn delimit_subproblem(
b: &mut Bidiagonal<N, R, C>,
diagonal: &mut VectorN<N, DimMinimum<R, C>>,
off_diagonal: &mut VectorN<N, DimDiff<DimMinimum<R, C>, U1>>,
u: &mut Option<MatrixMN<N, R, DimMinimum<R, C>>>,
v_t: &mut Option<MatrixMN<N, DimMinimum<R, C>, C>>,
is_upper_diagonal: bool,
end: usize,
eps: N,
eps: N::Real,
) -> (usize, usize)
{
let mut n = end;
@ -350,20 +364,20 @@ where
while n > 0 {
let m = n - 1;
if b.off_diagonal[m].is_zero()
|| b.off_diagonal[m].abs() <= eps * (b.diagonal[n].abs() + b.diagonal[m].abs())
if off_diagonal[m].is_zero()
|| off_diagonal[m].modulus() <= eps * (diagonal[n].modulus() + diagonal[m].modulus())
{
b.off_diagonal[m] = N::zero();
} else if b.diagonal[m].abs() <= eps {
b.diagonal[m] = N::zero();
Self::cancel_horizontal_off_diagonal_elt(b, u, v_t, m, m + 1);
off_diagonal[m] = N::zero();
} else if diagonal[m].modulus() <= eps {
diagonal[m] = N::zero();
Self::cancel_horizontal_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m, m + 1);
if m != 0 {
Self::cancel_vertical_off_diagonal_elt(b, u, v_t, m - 1);
Self::cancel_vertical_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m - 1);
}
} else if b.diagonal[n].abs() <= eps {
b.diagonal[n] = N::zero();
Self::cancel_vertical_off_diagonal_elt(b, u, v_t, m);
} else if diagonal[n].modulus() <= eps {
diagonal[n] = N::zero();
Self::cancel_vertical_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m);
} else {
break;
}
@ -379,18 +393,18 @@ where
while new_start > 0 {
let m = new_start - 1;
if b.off_diagonal[m].abs() <= eps * (b.diagonal[new_start].abs() + b.diagonal[m].abs())
if off_diagonal[m].modulus() <= eps * (diagonal[new_start].modulus() + diagonal[m].modulus())
{
b.off_diagonal[m] = N::zero();
off_diagonal[m] = N::zero();
break;
}
// FIXME: write a test that enters this case.
else if b.diagonal[m].abs() <= eps {
b.diagonal[m] = N::zero();
Self::cancel_horizontal_off_diagonal_elt(b, u, v_t, m, n);
else if diagonal[m].modulus() <= eps {
diagonal[m] = N::zero();
Self::cancel_horizontal_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m, n);
if m != 0 {
Self::cancel_vertical_off_diagonal_elt(b, u, v_t, m - 1);
Self::cancel_vertical_off_diagonal_elt(diagonal, off_diagonal, u, v_t, is_upper_diagonal, m - 1);
}
break;
}
@ -403,21 +417,23 @@ where
// Cancels the i-th off-diagonal element using givens rotations.
fn cancel_horizontal_off_diagonal_elt(
b: &mut Bidiagonal<N, R, C>,
diagonal: &mut VectorN<N, DimMinimum<R, C>>,
off_diagonal: &mut VectorN<N, DimDiff<DimMinimum<R, C>, U1>>,
u: &mut Option<MatrixMN<N, R, DimMinimum<R, C>>>,
v_t: &mut Option<MatrixMN<N, DimMinimum<R, C>, C>>,
is_upper_diagonal: bool,
i: usize,
end: usize,
)
{
let mut v = Vector2::new(b.off_diagonal[i], b.diagonal[i + 1]);
b.off_diagonal[i] = N::zero();
let mut v = Vector2::new(off_diagonal[i], diagonal[i + 1]);
off_diagonal[i] = N::zero();
for k in i..end {
if let Some((rot, norm)) = givens::cancel_x(&v) {
b.diagonal[k + 1] = norm;
if let Some((rot, norm)) = GivensRotation::cancel_x(&v) {
diagonal[k + 1] = norm;
if b.is_upper_diagonal() {
if is_upper_diagonal {
if let Some(ref mut u) = *u {
rot.inverse()
.rotate_rows(&mut u.fixed_columns_with_step_mut::<U2>(i, k - i));
@ -427,9 +443,9 @@ where
}
if k + 1 != end {
v.x = -rot.sin_angle() * b.off_diagonal[k + 1];
v.y = b.diagonal[k + 2];
b.off_diagonal[k + 1] *= rot.cos_angle();
v.x = -rot.s() * off_diagonal[k + 1];
v.y = diagonal[k + 2];
off_diagonal[k + 1] = off_diagonal[k + 1].scale(rot.c());
}
} else {
break;
@ -439,20 +455,22 @@ where
// Cancels the i-th off-diagonal element using givens rotations.
fn cancel_vertical_off_diagonal_elt(
b: &mut Bidiagonal<N, R, C>,
diagonal: &mut VectorN<N, DimMinimum<R, C>>,
off_diagonal: &mut VectorN<N, DimDiff<DimMinimum<R, C>, U1>>,
u: &mut Option<MatrixMN<N, R, DimMinimum<R, C>>>,
v_t: &mut Option<MatrixMN<N, DimMinimum<R, C>, C>>,
is_upper_diagonal: bool,
i: usize,
)
{
let mut v = Vector2::new(b.diagonal[i], b.off_diagonal[i]);
b.off_diagonal[i] = N::zero();
let mut v = Vector2::new(diagonal[i], off_diagonal[i]);
off_diagonal[i] = N::zero();
for k in (0..i + 1).rev() {
if let Some((rot, norm)) = givens::cancel_y(&v) {
b.diagonal[k] = norm;
if let Some((rot, norm)) = GivensRotation::cancel_y(&v) {
diagonal[k] = norm;
if b.is_upper_diagonal() {
if is_upper_diagonal {
if let Some(ref mut v_t) = *v_t {
rot.rotate(&mut v_t.fixed_rows_with_step_mut::<U2>(k, i - k));
}
@ -462,9 +480,9 @@ where
}
if k > 0 {
v.x = b.diagonal[k - 1];
v.y = rot.sin_angle() * b.off_diagonal[k - 1];
b.off_diagonal[k - 1] *= rot.cos_angle();
v.x = diagonal[k - 1];
v.y = rot.s() * off_diagonal[k - 1];
off_diagonal[k - 1] = off_diagonal[k - 1].scale(rot.c());
}
} else {
break;
@ -474,12 +492,12 @@ where
/// Computes the rank of the decomposed matrix, i.e., the number of singular values greater
/// than `eps`.
pub fn rank(&self, eps: N) -> usize {
pub fn rank(&self, eps: N::Real) -> usize {
assert!(
eps >= N::zero(),
eps >= N::Real::zero(),
"SVD rank: the epsilon must be non-negative."
);
self.singular_values.iter().filter(|e| **e > eps).count()
self.singular_values.iter().filter(|e| e.asum() > eps).count()
}
/// Rebuild the original matrix.
@ -507,18 +525,18 @@ where
/// Any singular value smaller than `eps` is assumed to be zero.
/// Returns `Err` if the right- and left- singular vectors have not
/// been computed at construction-time.
pub fn pseudo_inverse(mut self, eps: N) -> Result<MatrixMN<N, C, R>, &'static str>
pub fn pseudo_inverse(mut self, eps: N::Real) -> Result<MatrixMN<N, C, R>, &'static str>
where
DefaultAllocator: Allocator<N, C, R>,
{
if eps < N::zero() {
if eps < N::Real::zero() {
Err("SVD pseudo inverse: the epsilon must be non-negative.")
}
else {
for i in 0..self.singular_values.len() {
let val = self.singular_values[i];
if val > eps {
if val.asum() > eps {
self.singular_values[i] = N::one() / val;
} else {
self.singular_values[i] = N::zero();
@ -537,14 +555,14 @@ where
pub fn solve<R2: Dim, C2: Dim, S2>(
&self,
b: &Matrix<N, R2, C2, S2>,
eps: N,
eps: N::Real,
) -> Result<MatrixMN<N, C, C2>, &'static str>
where
S2: Storage<N, R2, C2>,
DefaultAllocator: Allocator<N, C, C2> + Allocator<N, DimMinimum<R, C>, C2>,
ShapeConstraint: SameNumberOfRows<R, R2>,
{
if eps < N::zero() {
if eps < N::Real::zero() {
Err("SVD solve: the epsilon must be non-negative.")
}
else {
@ -557,7 +575,7 @@ where
for i in 0..self.singular_values.len() {
let val = self.singular_values[i];
if val > eps {
if val.asum() > eps {
col[i] /= val;
} else {
col[i] = N::zero();
@ -575,7 +593,7 @@ where
}
}
impl<N: Real, R: DimMin<C>, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S>
impl<N: Complex, R: DimMin<C>, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S>
where
DimMinimum<R, C>: DimSub<U1>, // for Bidiagonal.
DefaultAllocator: Allocator<N, R, C>
@ -605,7 +623,7 @@ where
self,
compute_u: bool,
compute_v: bool,
eps: N,
eps: N::Real,
max_niter: usize,
) -> Option<SVD<N, R, C>>
{
@ -620,7 +638,7 @@ where
/// Computes the rank of this matrix.
///
/// All singular values below `eps` are considered equal to 0.
pub fn rank(&self, eps: N) -> usize {
pub fn rank(&self, eps: N::Real) -> usize {
let svd = SVD::new(self.clone_owned(), false, false);
svd.rank(eps)
}
@ -628,7 +646,7 @@ where
/// Computes the pseudo-inverse of this matrix.
///
/// All singular values below `eps` are considered equal to 0.
pub fn pseudo_inverse(self, eps: N) -> Result<MatrixMN<N, C, R>, &'static str>
pub fn pseudo_inverse(self, eps: N::Real) -> Result<MatrixMN<N, C, R>, &'static str>
where
DefaultAllocator: Allocator<N, C, R>,
{

41
src/linalg/symmetric_eigen.rs
View File

@ -1,7 +1,7 @@
#[cfg(feature = "serde-serialize")]
use serde::{Deserialize, Serialize};
use num::Zero;
use num::{Zero, One};
use num_complex::Complex as NumComplex;
use approx::AbsDiffEq;
use std::ops::MulAssign;
@ -119,6 +119,8 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
q = Some(res.0);
diag = res.1;
off_diag = res.2;
println!("Tridiagonalization q: {:.5?}", q);
} else {
let res = SymmetricTridiagonal::new(m).unpack_tridiagonal();
q = None;
@ -150,6 +152,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
let j = i + 1;
if let Some((rot, norm)) = GivensRotation::cancel_y(&v) {
println!("Canceling: {:.5?} with norm: {:.5?}", rot, norm);
if i > start {
// Not the first iteration.
off_diag[i - 1] = norm;
@ -160,19 +163,33 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
let mij = off_diag[i];
let cc = rot.c() * rot.c();
let ss = rot.s() * rot.s().conjugate();
let cs = rot.c() * rot.s();
let ss = rot.s().modulus_squared(); // rot.s() * rot.s().conjugate()
let cs = rot.s().scale(rot.c());
let b = cs * mij.conjugate() + cs.conjugate() * mij;
// b = cs * mij.conjugate() + cs.conjugate() * mij
let b = N::from_real((cs * mij.conjugate()).real() * ::convert(2.0));
diag[i] = (cc * mii + ss * mjj) - b;
diag[j] = (ss * mii + cc * mjj) + b;
off_diag[i] = cs * (mii - mjj) + mij * cc - mij.conjugate() * rot.s() * rot.s();
diag[i] = (mii.scale(cc) + mjj.scale(ss)) - b;
diag[j] = (mii.scale(ss) + mjj.scale(cc)) + b;
off_diag[i] = cs * (mii - mjj) + mij.scale(cc) - mij.conjugate() * rot.s() * rot.s();
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 {
v.x = off_diag[i];
v.y = -rot.s() * off_diag[i + 1];
off_diag[i + 1] *= rot.c();
off_diag[i + 1] = off_diag[i + 1].scale(rot.c());
}
if let Some(ref mut q) = q {
@ -197,12 +214,12 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
diag[start + 0] = eigvals[0];
diag[start + 1] = eigvals[1];
println!("Eigvals: {:?}", eigvals);
println!("m: {}", m);
println!("Curr q: {:?}", q);
println!("Eigvals: {:.5?}", eigvals);
println!("m: {:.5}", m);
println!("Curr q: {:.5?}", 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) {
rot.rotate_rows(&mut q.fixed_columns_mut::<U2>(start));
}
}

1
src/linalg/symmetric_tridiagonal.rs
View File

@ -83,6 +83,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>
m.ger_symm(-N::one(), &p, &axis.conjugate(), N::one());
m.ger_symm(-N::one(), &axis, &p.conjugate(), N::one());
m.ger_symm(dot * ::convert(2.0), &axis, &axis.conjugate(), N::one());
println!("The m: {}", m);
}
}

1
tests/linalg/bidiagonal.rs
View File

@ -41,6 +41,7 @@ quickcheck! {
relative_eq!(m, &u * d * &v_t, epsilon = 1.0e-7)
}
fn bidiagonal_static_square(m: Matrix4<RandComplex<f64>>) -> bool {
let m = m.map(|e| e.0);
let bidiagonal = m.bidiagonalize();

13
tests/linalg/eigen.rs
View File

@ -9,10 +9,9 @@ mod quickcheck_tests {
use std::cmp;
quickcheck! {
/*
fn symmetric_eigen(n: usize) -> bool {
let n = cmp::max(1, cmp::min(n, 10));
let m = DMatrix::<RandComplex<f64>>::new_random(n, n).map(|e| e.0);
let m = DMatrix::<RandComplex<f64>>::new_random(n, n).map(|e| e.0).hermitian_part();
let eig = m.clone().symmetric_eigen();
let recomp = eig.recompose();
@ -23,7 +22,7 @@ mod quickcheck_tests {
fn symmetric_eigen_singular(n: usize) -> bool {
let n = cmp::max(1, cmp::min(n, 10));
let mut m = DMatrix::<RandComplex<f64>>::new_random(n, n).map(|e| e.0);
let mut m = DMatrix::<RandComplex<f64>>::new_random(n, n).map(|e| e.0).hermitian_part();
m.row_mut(n / 2).fill(na::zero());
m.column_mut(n / 2).fill(na::zero());
let eig = m.clone().symmetric_eigen();
@ -35,7 +34,7 @@ mod quickcheck_tests {
}
fn symmetric_eigen_static_square_4x4(m: Matrix4<RandComplex<f64>>) -> bool {
let m = m.map(|e| e.0);
let m = m.map(|e| e.0).hermitian_part();
let eig = m.symmetric_eigen();
let recomp = eig.recompose();
@ -43,7 +42,6 @@ mod quickcheck_tests {
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
}
*/
fn symmetric_eigen_static_square_3x3(m: Matrix3<RandComplex<f64>>) -> bool {
let m = m.map(|e| e.0).hermitian_part();
@ -57,7 +55,6 @@ mod quickcheck_tests {
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
}
/*
fn symmetric_eigen_static_square_2x2(m: Matrix2<RandComplex<f64>>) -> bool {
let m = m.map(|e| e.0).hermitian_part();
let eig = m.symmetric_eigen();
@ -68,11 +65,9 @@ mod quickcheck_tests {
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5)
}
*/
}
}
/*
// Test proposed on the issue #176 of rulinalg.
#[test]
fn symmetric_eigen_singular_24x24() {
@ -115,7 +110,7 @@ fn symmetric_eigen_singular_24x24() {
recomp.lower_triangle(),
epsilon = 1.0e-5
));
}*/
}
// #[cfg(feature = "arbitrary")]
// quickcheck! {

21
tests/linalg/svd.rs
View File

@ -7,8 +7,11 @@ mod quickcheck_tests {
DMatrix, DVector, Matrix2, Matrix2x5, Matrix3, Matrix3x5, Matrix4, Matrix5x2, Matrix5x3,
};
use std::cmp;
use core::helper::{RandScalar, RandComplex};
quickcheck! {
/*
fn svd(m: DMatrix<f64>) -> bool {
if m.len() > 0 {
let svd = m.clone().svd(true, true);
@ -68,6 +71,7 @@ mod quickcheck_tests {
relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5)
}
fn svd_static_square(m: Matrix4<f64>) -> bool {
let svd = m.svd(true, true);
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
@ -78,18 +82,27 @@ mod quickcheck_tests {
u.is_orthogonal(1.0e-5) &&
v_t.is_orthogonal(1.0e-5)
}
*/
fn svd_static_square_2x2(m: Matrix2<f64>) -> bool {
fn svd_static_square_2x2(m: Matrix2<RandComplex<f64>>) -> bool {
let m = m.map(|e| e.0);
let svd = m.svd(true, true);
let (u, s, v_t) = (svd.u.unwrap(), svd.singular_values, svd.v_t.unwrap());
let ds = Matrix2::from_diagonal(&s);
s.iter().all(|e| *e >= 0.0) &&
println!("u, s, v_t: {}{}{}", u, s, v_t);
println!("m: {}", m);
println!("recomp: {}", u * ds * 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) &&
relative_eq!(m, u * ds * v_t, epsilon = 1.0e-5) &&
u.is_orthogonal(1.0e-5) &&
v_t.is_orthogonal(1.0e-5)
}
/*
fn svd_pseudo_inverse(m: DMatrix<f64>) -> bool {
if m.len() > 0 {
let svd = m.clone().svd(true, true);
@ -140,9 +153,11 @@ mod quickcheck_tests {
true
}
*/
}
}
/*
// Test proposed on the issue #176 of rulinalg.
#[test]
fn svd_singular() {
@ -342,3 +357,5 @@ fn svd_err() {
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));
}
*/

14
tests/linalg/tridiagonal.rs
View File

@ -15,6 +15,20 @@ quickcheck! {
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7)
}
fn symm_tridiagonal_singular(n: usize) -> bool {
let n = cmp::max(1, cmp::min(n, 4));
let mut m = DMatrix::<RandComplex<f64>>::new_random(n, n).map(|e| e.0).hermitian_part();
m.row_mut(n / 2).fill(na::zero());
m.column_mut(n / 2).fill(na::zero());
let tri = m.clone().symmetric_tridiagonalize();
println!("Tri: {:?}", tri);
let recomp = tri.recompose();
println!("Recomp: {:?}", recomp);
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7)
}
fn symm_tridiagonal_static_square(m: Matrix4<RandComplex<f64>>) -> bool {
let m = m.map(|e| e.0).hermitian_part();
let tri = m.symmetric_tridiagonalize();