From 010c009cfffd5ba748e02b32c803bad35ce8fc47 Mon Sep 17 00:00:00 2001 From: sebcrozet Date: Tue, 12 Mar 2019 13:15:02 +0100 Subject: [PATCH] Fix Schur decomposition. --- src/base/blas.rs | 79 ++++++++++++++++++++++++++++- src/base/matrix.rs | 7 ++- src/geometry/reflection.rs | 4 +- src/linalg/givens.rs | 31 +++++------ src/linalg/schur.rs | 21 ++++---- src/linalg/symmetric_eigen.rs | 17 ++++--- src/linalg/symmetric_tridiagonal.rs | 8 +-- tests/linalg/eigen.rs | 14 +++-- tests/linalg/hessenberg.rs | 7 +-- tests/linalg/real_schur.rs | 34 +++++++------ tests/linalg/tridiagonal.rs | 43 +++++++--------- 11 files changed, 181 insertions(+), 84 deletions(-) diff --git a/src/base/blas.rs b/src/base/blas.rs index 6c0c029d..5eddf2b5 100644 --- a/src/base/blas.rs +++ b/src/base/blas.rs @@ -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( + &mut self, + alpha: N, + a: &SquareMatrix, + x: &Vector, + beta: N, + ) where + N: Complex, + SB: Storage, + SC: Storage, + ShapeConstraint: DimEq + AreMultipliable, + { + 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. /// diff --git a/src/base/matrix.rs b/src/base/matrix.rs index 3fa276f0..e7cf0026 100644 --- a/src/base/matrix.rs +++ b/src/base/matrix.rs @@ -997,7 +997,12 @@ impl> Matrix { 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; diff --git a/src/geometry/reflection.rs b/src/geometry/reflection.rs index ac63b7da..fca0f548 100644 --- a/src/geometry/reflection.rs +++ b/src/geometry/reflection.rs @@ -72,13 +72,13 @@ impl> Reflection { ShapeConstraint: DimEq + AreMultipliable, DefaultAllocator: Allocator { - 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()); } } diff --git a/src/linalg/givens.rs b/src/linalg/givens.rs index 7e6a881e..472f3bb3 100644 --- a/src/linalg/givens.rs +++ b/src/linalg/givens.rs @@ -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 { + // FIXME: c should be a `N::Real`. c: N, s: N } @@ -47,24 +50,22 @@ pub fn cancel_x>(v: &Vector) -> Option<(Uni // Matrix = UnitComplex * Matrix impl GivensRotation { - /// 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 { - 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 GivensRotation { /// 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 GivensRotation { 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 GivensRotation { 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; } } } diff --git a/src/linalg/schur.rs b/src/linalg/schur.rs index fa7e1736..2eca9843 100644 --- a/src/linalg/schur.rs +++ b/src/linalg/schur.rs @@ -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>( 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>( } 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 diff --git a/src/linalg/symmetric_eigen.rs b/src/linalg/symmetric_eigen.rs index 8d16232f..6084c1b4 100644 --- a/src/linalg/symmetric_eigen.rs +++ b/src/linalg/symmetric_eigen.rs @@ -43,6 +43,7 @@ where DefaultAllocator: Allocator + Allocator /// The eigenvectors of the decomposed matrix. pub eigenvectors: MatrixN, + // FIXME: this should be a VectorN /// The unsorted eigenvalues of the decomposed matrix. pub eigenvalues: VectorN, } @@ -159,14 +160,14 @@ where DefaultAllocator: Allocator + Allocator 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 + Allocator } } 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 + Allocator 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::(start)); diff --git a/src/linalg/symmetric_tridiagonal.rs b/src/linalg/symmetric_tridiagonal.rs index 276dc6fb..9b54c596 100644 --- a/src/linalg/symmetric_tridiagonal.rs +++ b/src/linalg/symmetric_tridiagonal.rs @@ -53,6 +53,8 @@ where DefaultAllocator: Allocator + Allocator> pub fn new(mut m: MatrixN) -> 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 + Allocator> 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()); } } diff --git a/tests/linalg/eigen.rs b/tests/linalg/eigen.rs index c44bb098..28a2b0d1 100644 --- a/tests/linalg/eigen.rs +++ b/tests/linalg/eigen.rs @@ -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::>::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>) -> 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>) -> 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! { diff --git a/tests/linalg/hessenberg.rs b/tests/linalg/hessenberg.rs index d9245110..72c2baab 100644 --- a/tests/linalg/hessenberg.rs +++ b/tests/linalg/hessenberg.rs @@ -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>) -> 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>) -> 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) } } diff --git a/tests/linalg/real_schur.rs b/tests/linalg/real_schur.rs index 554bbdef..54ad9ce4 100644 --- a/tests/linalg/real_schur.rs +++ b/tests/linalg/real_schur.rs @@ -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::::new_random(n, n); + let m = DMatrix::>::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) -> bool { + fn schur_static_mat2(m: Matrix2>) -> 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) -> bool { + fn schur_static_mat3(m: Matrix3>) -> 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) -> bool { + fn schur_static_mat4(m: Matrix4>) -> 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 } diff --git a/tests/linalg/tridiagonal.rs b/tests/linalg/tridiagonal.rs index 132edba5..5ecf32ff 100644 --- a/tests/linalg/tridiagonal.rs +++ b/tests/linalg/tridiagonal.rs @@ -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::>::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>) -> 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::>::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>) -> 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>) -> 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>) -> 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) + } }