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

Fix Schur decomposition.

Browse Source
This commit is contained in:
sebcrozet 2019-03-12 13:15:02 +01:00
parent 77f048b6b9
commit 010c009cff
11 changed files with 181 additions and 84 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

79
src/base/blas.rs
View File

@ -684,7 +684,7 @@ where
assert!(
a.is_square(),
"Syetric gemv: the input matrix must be square."
"Symmetric gemv: the input matrix must be square."
);
assert!(
dim2 == dim3 && dim1 == dim2,
@ -715,6 +715,83 @@ where
}
}
/// Computes `self = alpha * a * x + beta * self`, where `a` is a **symmetric** matrix, `x` a
/// vector, and `alpha, beta` two scalars.
///
/// If `beta` is zero, `self` is never read. If `self` is read, only its lower-triangular part
/// (including the diagonal) is actually read.
///
/// # Examples:
///
/// ```
/// # use nalgebra::{Matrix2, Vector2};
/// let mat = Matrix2::new(1.0, 2.0,
/// 2.0, 4.0);
/// let mut vec1 = Vector2::new(1.0, 2.0);
/// let vec2 = Vector2::new(0.1, 0.2);
/// vec1.gemv_symm(10.0, &mat, &vec2, 5.0);
/// assert_eq!(vec1, Vector2::new(10.0, 20.0));
///
///
/// // The matrix upper-triangular elements can be garbage because it is never
/// // read by this method. Therefore, it is not necessary for the caller to
/// // fill the matrix struct upper-triangle.
/// let mat = Matrix2::new(1.0, 9999999.9999999,
/// 2.0, 4.0);
/// let mut vec1 = Vector2::new(1.0, 2.0);
/// vec1.gemv_symm(10.0, &mat, &vec2, 5.0);
/// assert_eq!(vec1, Vector2::new(10.0, 20.0));
/// ```
#[inline]
pub fn cgemv_symm<D2: Dim, D3: Dim, SB, SC>(
&mut self,
alpha: N,
a: &SquareMatrix<N, D2, SB>,
x: &Vector<N, D3, SC>,
beta: N,
) where
N: Complex,
SB: Storage<N, D2, D2>,
SC: Storage<N, D3>,
ShapeConstraint: DimEq<D, D2> + AreMultipliable<D2, D2, D3, U1>,
{
let dim1 = self.nrows();
let dim2 = a.nrows();
let dim3 = x.nrows();
assert!(
a.is_square(),
"Symmetric cgemv: the input matrix must be square."
);
assert!(
dim2 == dim3 && dim1 == dim2,
"Symmetric cgemv: dimensions mismatch."
);
if dim2 == 0 {
return;
}
// FIXME: avoid bound checks.
let col2 = a.column(0);
let val = unsafe { *x.vget_unchecked(0) };
self.axpy(alpha * val, &col2, beta);
self[0] += alpha * a.slice_range(1.., 0).cdot(&x.rows_range(1..));
for j in 1..dim2 {
let col2 = a.column(j);
let dot = col2.rows_range(j..).cdot(&x.rows_range(j..));
let val;
unsafe {
val = *x.vget_unchecked(j);
*self.vget_unchecked_mut(j) += alpha * dot;
}
self.rows_range_mut(j + 1..)
.axpy(alpha * val, &col2.rows_range(j + 1..), N::one());
}
}
/// Computes `self = alpha * a.transpose() * x + beta * self`, where `a` is a matrix, `x` a vector, and
/// `alpha, beta` two scalars.
///

7
src/base/matrix.rs
View File

@ -997,7 +997,12 @@ impl<N: Complex, D: Dim, S: StorageMut<N, D, D>> Matrix<N, D, D, S> {
let dim = self.shape().0;
for i in 1..dim {
for i in 0..dim {
{
let diag = unsafe { self.get_unchecked_mut((i, i)) };
*diag = diag.conjugate();
}
for j in 0..i {
unsafe {
let ref_ij = self.get_unchecked_mut((i, j)) as *mut N;

4
src/geometry/reflection.rs
View File

@ -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.conjugate(), work);
rhs.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, N::one());
rhs.ger(m_two, &work, &self.axis.conjugate(), N::one());
}
}

31
src/linalg/givens.rs
View File

@ -1,6 +1,7 @@
//! Construction of givens rotations.
use alga::general::{Complex, Real};
use num::Zero;
use num_complex::Complex as NumComplex;
use base::dimension::{Dim, U2};
@ -11,7 +12,9 @@ use base::{Vector, Matrix};
use geometry::UnitComplex;
/// A Givens rotation.
#[derive(Debug)]
pub struct GivensRotation<N: Complex> {
// FIXME: c should be a `N::Real`.
c: N,
s: N
}
@ -47,24 +50,22 @@ pub fn cancel_x<N: Real, S: Storage<N, U2>>(v: &Vector<N, U2, S>) -> Option<(Uni
// Matrix = UnitComplex * Matrix
impl<N: Complex> GivensRotation<N> {
/// Initializes a Givens rotation form 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 {
let denom = (c.modulus_squared() + s.modulus_squared()).sqrt();
Self {
c: c.unscale(denom),
s: s.unscale(denom)
}
let res = Self::try_new(c, s, N::Real::zero()).unwrap();
println!("The rot: {:?}", res);
res
}
/// 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> {
let denom = (c.modulus_squared() + s.modulus_squared()).sqrt();
let (mod0, sign0) = c.to_exp();
let denom = (mod0 * mod0 + s.modulus_squared()).sqrt();
if denom > eps {
Some(Self {
c: c.unscale(denom),
s: s.unscale(denom)
})
let c = N::from_real(mod0 / denom);
let s = s / sign0.scale(denom);
Some(Self { c, s })
} else {
None
}
@ -116,7 +117,7 @@ impl<N: Complex> GivensRotation<N> {
/// The inverse of this givens rotation.
pub fn inverse(&self) -> Self {
Self { c: self.c, s: -self.s.conjugate() }
Self { c: self.c, s: -self.s }
}
/// Performs the multiplication `rhs = self * rhs` in-place.
@ -139,8 +140,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.conjugate() * b;
*rhs.get_unchecked_mut((0, j)) = c * a + -s.conjugate() * b;
*rhs.get_unchecked_mut((1, j)) = s * a + c * b;
}
}
}
@ -167,7 +168,7 @@ impl<N: Complex> GivensRotation<N> {
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.conjugate() * b;
*lhs.get_unchecked_mut((j, 1)) = -s.conjugate() * a + c * b;
}
}
}

21
src/linalg/schur.rs
View File

@ -96,6 +96,7 @@ where
let dim = m.data.shape().0;
// Specialization would make this easier.
if dim.value() == 0 {
let vecs = Some(MatrixN::from_element_generic(dim, dim, N::zero()));
let vals = MatrixN::from_element_generic(dim, dim, N::zero());
@ -107,9 +108,7 @@ where
} else {
return Some((None, m));
}
}
// Specialization would make this easier.
else if dim.value() == 2 {
} else if dim.value() == 2 {
return decompose_2x2(m, compute_q);
}
@ -421,7 +420,7 @@ where
q = Some(MatrixN::from_column_slice_generic(
dim,
dim,
&[rot.c(), rot.s(), -rot.s().conjugate(), rot.c().conjugate()],
&[rot.c(), rot.s(), -rot.s().conjugate(), rot.c()],
));
}
}
@ -444,9 +443,9 @@ fn compute_2x2_eigvals<N: Complex, S: Storage<N, U2, U2>>(
let h01 = m[(0, 1)];
let h11 = m[(1, 1)];
// NOTE: this discriminant computation is mor stable than the
// NOTE: this discriminant computation is more stable than the
// one based on the trace and determinant: 0.25 * tra * tra - det
// because et ensures positiveness for symmetric matrices.
// because it ensures positiveness for symmetric matrices.
let val = (h00 - h11) * ::convert(0.5);
let discr = h10 * h01 + val * val;
@ -471,16 +470,18 @@ fn compute_2x2_basis<N: Complex, S: Storage<N, U2, U2>>(
}
if let Some((eigval1, eigval2)) = compute_2x2_eigvals(m) {
let x1 = m[(1, 1)] - eigval1;
let x2 = m[(1, 1)] - eigval2;
let x1 = eigval1 - m[(1, 1)];
let x2 = eigval2 - m[(1, 1)];
println!("eigval1: {}, eigval2: {}, h10: {}", eigval1, eigval2, h10);
// NOTE: Choose the one that yields a larger x component.
// 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))
} else {
Some(GivensRotation::new(x2, -h10))
Some(GivensRotation::new(x2, h10))
}
} else {
None

17
src/linalg/symmetric_eigen.rs
View File

@ -43,6 +43,7 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
/// The eigenvectors of the decomposed matrix.
pub eigenvectors: MatrixN<N, D>,
// FIXME: this should be a VectorN<N::Real, D>
/// The unsorted eigenvalues of the decomposed matrix.
pub eigenvalues: VectorN<N, D>,
}
@ -159,14 +160,14 @@ 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();
let ss = rot.s() * rot.s().conjugate();
let cs = rot.c() * rot.s();
let b = cs * ::convert(2.0) * mij;
let b = cs * mij.conjugate() + cs.conjugate() * mij;
diag[i] = (cc * mii + ss * mjj) - b;
diag[j] = (ss * mii + cc * mjj) + b;
off_diag[i] = cs * (mii - mjj) + mij * (cc - ss);
off_diag[i] = cs * (mii - mjj) + mij * cc - mij.conjugate() * rot.s() * rot.s();
if i != n - 1 {
v.x = off_diag[i];
@ -187,10 +188,8 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, D>
}
} else if subdim == 2 {
let m = Matrix2::new(
diag[start],
off_diag[start],
off_diag[start],
diag[start + 1],
diag[start], off_diag[start].conjugate(),
off_diag[start], diag[start + 1],
);
let eigvals = m.eigenvalues().unwrap();
let basis = Vector2::new(eigvals.x - diag[start + 1], off_diag[start]);
@ -198,6 +197,10 @@ 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);
if let Some(ref mut q) = q {
if let Some(rot) = GivensRotation::try_new(basis.x, basis.y, eps) {
rot.rotate_rows(&mut q.fixed_columns_mut::<U2>(start));

8
src/linalg/symmetric_tridiagonal.rs
View File

@ -53,6 +53,8 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>
pub fn new(mut m: MatrixN<N, D>) -> Self {
let dim = m.data.shape().0;
println!("Input m: {}", m.index((0.., 0..)));
assert!(
m.is_square(),
"Unable to compute the symmetric tridiagonal decomposition of a non-square matrix."
@ -75,11 +77,11 @@ where DefaultAllocator: Allocator<N, D, D> + Allocator<N, DimDiff<D, U1>>
if not_zero {
let mut p = p.rows_range_mut(i..);
p.gemv_symm(::convert(2.0), &m, &axis.conjugate(), N::zero());
let dot = axis.dot(&p);
p.cgemv_symm(::convert(2.0), &m, &axis, N::zero());
let dot = axis.cdot(&p);
// p.axpy(-dot, &axis.conjugate(), N::one());
m.ger_symm(-N::one(), &p, &axis.conjugate(), N::one());
m.ger_symm(-N::one(), &axis, &p, N::one());
m.ger_symm(-N::one(), &axis, &p.conjugate(), N::one());
m.ger_symm(dot * ::convert(2.0), &axis, &axis.conjugate(), N::one());
}
}

14
tests/linalg/eigen.rs
View File

@ -9,6 +9,7 @@ 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);
@ -42,29 +43,36 @@ 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);
let m = m.map(|e| e.0).hermitian_part();
let eig = m.symmetric_eigen();
let recomp = eig.recompose();
println!("Eigenvectors: {}", eig.eigenvectors);
println!("Eigenvalues: {}", eig.eigenvalues);
println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
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);
let m = m.map(|e| e.0).hermitian_part();
let eig = m.symmetric_eigen();
let recomp = eig.recompose();
println!("Eigenvectors: {}", eig.eigenvectors);
println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
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() {
@ -107,7 +115,7 @@ fn symmetric_eigen_singular_24x24() {
recomp.lower_triangle(),
epsilon = 1.0e-5
));
}
}*/
// #[cfg(feature = "arbitrary")]
// quickcheck! {

7
tests/linalg/hessenberg.rs
View File

@ -4,6 +4,7 @@ use na::{DMatrix, Matrix2, Matrix4};
use core::helper::{RandScalar, RandComplex};
use std::cmp;
#[test]
fn hessenberg_simple() {
let m = Matrix2::new(1.0, 0.0, 1.0, 3.0);
@ -19,20 +20,20 @@ quickcheck! {
let hess = m.clone().hessenberg();
let (p, h) = hess.unpack();
relative_eq!(m, &p * h * p.transpose(), epsilon = 1.0e-7)
relative_eq!(m, &p * h * p.conjugate_transpose(), epsilon = 1.0e-7)
}
fn hessenberg_static_mat2(m: Matrix2<RandComplex<f64>>) -> bool {
let m = m.map(|e| e.0);
let hess = m.hessenberg();
let (p, h) = hess.unpack();
relative_eq!(m, p * h * p.transpose(), epsilon = 1.0e-7)
relative_eq!(m, p * h * p.conjugate_transpose(), epsilon = 1.0e-7)
}
fn hessenberg_static(m: Matrix4<RandComplex<f64>>) -> bool {
let m = m.map(|e| e.0);
let hess = m.hessenberg();
let (p, h) = hess.unpack();
relative_eq!(m, p * h * p.transpose(), epsilon = 1.0e-7)
relative_eq!(m, p * h * p.conjugate_transpose(), epsilon = 1.0e-7)
}
}

34
tests/linalg/real_schur.rs
View File

@ -2,7 +2,6 @@
use na::{DMatrix, Matrix3, Matrix4};
#[test]
fn schur_simpl_mat3() {
let m = Matrix3::new(-2.0, -4.0, 2.0,
@ -19,48 +18,53 @@ fn schur_simpl_mat3() {
mod quickcheck_tests {
use std::cmp;
use na::{DMatrix, Matrix2, Matrix3, Matrix4};
use core::helper::{RandScalar, RandComplex};
quickcheck! {
fn schur(n: usize) -> bool {
let n = cmp::max(1, cmp::min(n, 10));
let m = DMatrix::<f64>::new_random(n, n);
let m = DMatrix::<RandComplex<f64>>::new_random(n, n).map(|e| e.0);
let (vecs, vals) = m.clone().real_schur().unpack();
if !relative_eq!(&vecs * &vals * vecs.transpose(), m, epsilon = 1.0e-7) {
println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose());
if !relative_eq!(&vecs * &vals * vecs.conjugate_transpose(), m, epsilon = 1.0e-7) {
println!("{:.5}{:.5}", m, &vecs * &vals * vecs.conjugate_transpose());
}
relative_eq!(&vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7)
relative_eq!(&vecs * vals * vecs.conjugate_transpose(), m, epsilon = 1.0e-7)
}
fn schur_static_mat2(m: Matrix2<f64>) -> bool {
fn schur_static_mat2(m: Matrix2<RandComplex<f64>>) -> bool {
let m = m.map(|e| e.0);
let (vecs, vals) = m.clone().real_schur().unpack();
let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7);
let ok = relative_eq!(vecs * vals * vecs.conjugate_transpose(), m, epsilon = 1.0e-7);
if !ok {
println!("{:.5}{:.5}", vecs, vals);
println!("Reconstruction:{}{}", m, &vecs * &vals * vecs.transpose());
println!("Vecs: {:.5} Vals: {:.5}", vecs, vals);
println!("Reconstruction:{}{}", m, &vecs * &vals * vecs.conjugate_transpose());
}
ok
}
fn schur_static_mat3(m: Matrix3<f64>) -> bool {
fn schur_static_mat3(m: Matrix3<RandComplex<f64>>) -> bool {
let m = m.map(|e| e.0);
let (vecs, vals) = m.clone().real_schur().unpack();
let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7);
let ok = relative_eq!(vecs * vals * vecs.conjugate_transpose(), m, epsilon = 1.0e-7);
if !ok {
println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose());
println!("Vecs: {:.5} Vals: {:.5}", vecs, vals);
println!("{:.5}{:.5}", m, &vecs * &vals * vecs.conjugate_transpose());
}
ok
}
fn schur_static_mat4(m: Matrix4<f64>) -> bool {
fn schur_static_mat4(m: Matrix4<RandComplex<f64>>) -> bool {
let m = m.map(|e| e.0);
let (vecs, vals) = m.clone().real_schur().unpack();
let ok = relative_eq!(vecs * vals * vecs.transpose(), m, epsilon = 1.0e-7);
let ok = relative_eq!(vecs * vals * vecs.conjugate_transpose(), m, epsilon = 1.0e-7);
if !ok {
println!("{:.5}{:.5}", m, &vecs * &vals * vecs.transpose());
println!("{:.5}{:.5}", m, &vecs * &vals * vecs.conjugate_transpose());
}
ok
}

43
tests/linalg/tridiagonal.rs
View File

@ -6,33 +6,28 @@ use na::{DMatrix, Matrix2, Matrix4};
use core::helper::{RandScalar, RandComplex};
quickcheck! {
// fn symm_tridiagonal(n: usize) -> bool {
// let n = cmp::max(1, cmp::min(n, 50));
// let m = DMatrix::<RandComplex<f64>>::new_random(n, n).map(|e| e.0).hermitian_part();
// let tri = m.clone().symmetric_tridiagonalize();
// let recomp = tri.recompose();
//
// println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
//
// 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();
println!("Internal tri: {}{}", tri.internal_tri(), tri.off_diagonal());
fn symm_tridiagonal(n: usize) -> bool {
let n = cmp::max(1, cmp::min(n, 50));
let m = DMatrix::<RandComplex<f64>>::new_random(n, n).map(|e| e.0).hermitian_part();
let tri = m.clone().symmetric_tridiagonalize();
let recomp = tri.recompose();
println!("{}{}", m.lower_triangle(), recomp.lower_triangle());
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7)
}
// fn symm_tridiagonal_static_square_2x2(m: Matrix2<RandComplex<f64>>) -> bool {
// let m = m.map(|e| e.0).hermitian_part();
// let tri = m.symmetric_tridiagonalize();
// let recomp = tri.recompose();
//
// 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();
let recomp = tri.recompose();
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7)
}
fn symm_tridiagonal_static_square_2x2(m: Matrix2<RandComplex<f64>>) -> bool {
let m = m.map(|e| e.0).hermitian_part();
let tri = m.symmetric_tridiagonalize();
let recomp = tri.recompose();
relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-7)
}
}