diff --git a/src/lib.rs b/src/lib.rs index 595da543..e28a4650 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -151,7 +151,8 @@ pub use structs::{ pub use linalg::{ qr, householder_matrix, - cholesky + cholesky, + hessenberg }; mod structs; diff --git a/src/linalg/decompositions.rs b/src/linalg/decompositions.rs index 5cbcc63b..fc996161 100644 --- a/src/linalg/decompositions.rs +++ b/src/linalg/decompositions.rs @@ -1,5 +1,5 @@ use traits::operations::{Transpose, ApproxEq}; -use traits::structure::{ColSlice, Eye, Indexable, Diag, SquareMat, BaseFloat}; +use traits::structure::{ColSlice, Eye, Indexable, Diag, SquareMat, BaseFloat, Cast}; use traits::geometry::Norm; use std::cmp::min; use std::ops::{Mul, Add, Sub}; @@ -71,7 +71,7 @@ pub fn qr(m: &M) -> (M, M) (q, r) } -/// Eigendecomposition of a square matrix using the qr algorithm. +/// Eigendecomposition of a square symmetric matrix using the qr algorithm pub fn eigen_qr(m: &M, eps: &N, niter: usize) -> (M, V) where N: BaseFloat, V: Mul, @@ -79,39 +79,165 @@ pub fn eigen_qr(m: &M, eps: &N, niter: usize) -> (M, V) M: Indexable<(usize, usize), N> + SquareMat + Add + Sub + ColSlice + ApproxEq + Copy { - let mut eigenvectors: M = ::one::(); - let mut eigenvalues = *m; - // let mut shifter: M = Eye::new_identity(rows); + + let (mut eigenvectors, mut eigenvalues) = hessenberg(m); + + // Allocate arrays for Givens rotation components + let mut c = Vec::::with_capacity(::dim::()-1); + let mut s = Vec::::with_capacity(::dim::()-1); + + if ::dim::() == 1 { + return (eigenvectors, eigenvalues.diag()); + } + + unsafe { + c.set_len(::dim::()-1); + s.set_len(::dim::()-1); + } let mut iter = 0; - for _ in 0..niter { - let mut stop = true; + let mut curdim = ::dim::()-1; + + for _ in 0..::dim::() { + + let mut stop = false; + + while !stop && iter < niter { + + let lambda; + + unsafe { + let a = eigenvalues.unsafe_at((curdim-1, curdim-1)); + let b = eigenvalues.unsafe_at((curdim-1, curdim)); + let c = eigenvalues.unsafe_at((curdim, curdim-1)); + let d = eigenvalues.unsafe_at((curdim, curdim)); + + let trace = a + d; + let det = a * d - b * c; + + let constquarter: N = Cast::from(0.25f64); + let consthalf: N = Cast::from(0.5f64); + + let e = (constquarter * trace * trace - det).sqrt(); + + let lambda1 = consthalf * trace + e; + let lambda2 = consthalf * trace - e; + + if (lambda1 - d).abs() < (lambda2 - d).abs() { + lambda = lambda1; + } + else { + lambda = lambda2; + } + + } + + // Shift matrix + for k in 0..curdim+1 { + unsafe { + let a = eigenvalues.unsafe_at((k,k)); + eigenvalues.unsafe_set((k,k), a - lambda); + } + } + + + // Givens rotation from left + for k in 0..curdim { + let x_i = unsafe { eigenvalues.unsafe_at((k,k)) }; + let x_j = unsafe { eigenvalues.unsafe_at((k+1,k)) }; + let r = (x_i*x_i + x_j*x_j).sqrt(); + let ctmp = x_i / r; + let stmp = -x_j / r; + + c[k] = ctmp; + s[k] = stmp; + + for j in k..(curdim+1) { + + unsafe { + let a = eigenvalues.unsafe_at((k,j)); + let b = eigenvalues.unsafe_at((k+1,j)); + + eigenvalues.unsafe_set((k,j), ctmp * a - stmp * b); + eigenvalues.unsafe_set((k+1,j), stmp * a + ctmp * b); + } - for j in 0..::dim::() { - for i in 0..j { - if unsafe { eigenvalues.unsafe_at((i, j)) }.abs() >= *eps { - stop = false; - break; } } - for i in j + 1..::dim::() { - if unsafe { eigenvalues.unsafe_at((i, j)) }.abs() >= *eps { + // Givens rotation from right applied to eigenvalues + for k in 0..curdim { + for i in 0..(k+2) { + + unsafe { + let a = eigenvalues.unsafe_at((i,k)); + let b = eigenvalues.unsafe_at((i,k+1)); + + + eigenvalues.unsafe_set((i,k), c[k] * a - s[k] * b); + eigenvalues.unsafe_set((i,k+1), s[k] * a + c[k] * b); + } + } + + } + + + // Shift back + for k in 0..curdim+1 { + unsafe { + let a = eigenvalues.unsafe_at((k,k)); + eigenvalues.unsafe_set((k,k), a + lambda); + } + } + + + // Givens rotation from right applied to eigenvectors + for k in 0..curdim { + for i in 0..::dim::() { + + unsafe { + let a = eigenvectors.unsafe_at((i,k)); + let b = eigenvectors.unsafe_at((i,k+1)); + + + eigenvectors.unsafe_set((i,k), c[k] * a - s[k] * b); + eigenvectors.unsafe_set((i,k+1), s[k] * a + c[k] * b); + } + } + } + + iter = iter + 1; + + stop = true; + + for j in 0..curdim { + + // Check last row + if unsafe { eigenvalues.unsafe_at((curdim, j)) }.abs() >= *eps { + stop = false; + break; + } + + // Check last column + if unsafe { eigenvalues.unsafe_at((j, curdim)) }.abs() >= *eps { stop = false; break; } } } + if stop { - break; + + if curdim > 1 { + curdim = curdim - 1; + } + else { + break; + } + } - iter = iter + 1; - let (q, r) = qr(&eigenvalues);; - - eigenvalues = r * q; - eigenvectors = eigenvectors * q; } (eigenvectors, eigenvalues.diag()) @@ -164,3 +290,48 @@ pub fn cholesky(m: &M) -> Result return Ok(out); } + +/// Hessenberg +/// Returns the matrix m in Hessenberg form and the corresponding similarity transformation +/// +/// # Arguments +/// * `m` - matrix to transform +/// +/// # Returns +/// * First return value `q` - Similarity matrix p such that q * h * q^T = m +/// * Second return value `h` - Matrix m in Hessenberg form +pub fn hessenberg(m: &M) -> (M, M) + where N: BaseFloat, + V: Indexable + Norm, + M: Copy + Eye + ColSlice + Transpose + Indexable<(usize, usize), N> + + Mul { + + let mut h = m.clone(); + let (rows, cols) = h.shape(); + + let mut q : M = Eye::new_identity(cols); + + if cols <= 2 { + return (q, h); + } + + for ite in 0..(cols-2) { + let mut v = h.col_slice(ite, ite+1, rows); + + let alpha = Norm::norm(&v); + + unsafe { + let x = v.unsafe_at(0); + v.unsafe_set(0, x - alpha); + } + + if !::is_zero(&v.normalize_mut()) { + let p: M = householder_matrix(rows, ite+1, v); + + q = q * p; + h = p * h * p; + } + } + + return (q, h); +} diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs index 80b22827..6c4f53b6 100644 --- a/src/linalg/mod.rs +++ b/src/linalg/mod.rs @@ -1,4 +1,4 @@ -pub use self::decompositions::{qr, eigen_qr, householder_matrix, cholesky}; +pub use self::decompositions::{qr, eigen_qr, householder_matrix, cholesky, hessenberg}; mod decompositions; diff --git a/tests/mat.rs b/tests/mat.rs index 046a0a06..72c14877 100644 --- a/tests/mat.rs +++ b/tests/mat.rs @@ -3,7 +3,7 @@ extern crate rand; use rand::random; use na::{Vec1, Vec3, Mat1, Mat2, Mat3, Mat4, Mat5, Mat6, Rot2, Rot3, Persp3, PerspMat3, Ortho3, - OrthoMat3, DMat, DVec, Row, Col, BaseFloat, Diag, Transpose, RowSlice, ColSlice}; + OrthoMat3, DMat, DVec, Row, Col, BaseFloat, Diag, Transpose, RowSlice, ColSlice, Shape}; macro_rules! test_inv_mat_impl( ($t: ty) => ( @@ -63,28 +63,49 @@ macro_rules! test_cholesky_impl( ); ); -// NOTE: deactivated untile we get a better convergence rate. -// macro_rules! test_eigen_qr_impl( -// ($t: ty) => { -// for _ in (0usize .. 10000) { -// let randmat : $t = random(); -// // Make it symetric so that we can recompose the matrix to test at the end. -// let randmat = na::transpose(&randmat) * randmat; -// -// let (eigenvectors, eigenvalues) = na::eigen_qr(&randmat, &Float::epsilon(), 100); -// -// let diag: $t = Diag::from_diag(&eigenvalues); -// -// let recomp = eigenvectors * diag * na::transpose(&eigenvectors); -// -// println!("eigenvalues: {}", eigenvalues); -// println!(" mat: {}", randmat); -// println!("recomp: {}", recomp); -// -// assert!(na::approx_eq_eps(&randmat, &recomp, &1.0e-2)); -// } -// } -// ) +macro_rules! test_hessenberg_impl( + ($t: ty) => ( + for _ in (0usize .. 10000) { + + let randmat : $t = random(); + + let (q, h) = na::hessenberg(&randmat); + let recomp = q * h * na::transpose(&q); + + let (rows, cols) = h.shape(); + + // Check if `h` has zero entries below the first subdiagonal + if cols > 2 { + for j in 0..(cols-2) { + for i in (j+2)..rows { + assert!(na::approx_eq(&h[(i,j)], &0.0f64)); + } + } + } + + assert!(na::approx_eq(&randmat, &recomp)); + } + ); +); + +macro_rules! test_eigen_qr_impl( + ($t: ty) => { + for _ in (0usize .. 10000) { + let randmat : $t = random(); + // Make it symetric so that we can recompose the matrix to test at the end. + let randmat = na::transpose(&randmat) * randmat; + let (eigenvectors, eigenvalues) = na::eigen_qr(&randmat, &1e-13, 100); + + let diag: $t = Diag::from_diag(&eigenvalues); + let recomp = eigenvectors * diag * na::transpose(&eigenvectors); + println!("eigenvalues: {:?}", eigenvalues); + println!(" mat: {:?}", randmat); + println!("recomp: {:?}", recomp); + + assert!(na::approx_eq_eps(&randmat, &recomp, &1.0e-2)); + } + } +); #[test] fn test_transpose_mat1() { @@ -532,36 +553,35 @@ fn test_qr_mat6() { test_qr_impl!(Mat6); } -// NOTE: deactivated until we get a better convergence rate. -// #[test] -// fn test_eigen_qr_mat1() { -// test_eigen_qr_impl!(Mat1); -// } -// -// #[test] -// fn test_eigen_qr_mat2() { -// test_eigen_qr_impl!(Mat2); -// } -// -// #[test] -// fn test_eigen_qr_mat3() { -// test_eigen_qr_impl!(Mat3); -// } -// -// #[test] -// fn test_eigen_qr_mat4() { -// test_eigen_qr_impl!(Mat4); -// } -// -// #[test] -// fn test_eigen_qr_mat5() { -// test_eigen_qr_impl!(Mat5); -// } -// -// #[test] -// fn test_eigen_qr_mat6() { -// test_eigen_qr_impl!(Mat6); -// } +#[test] +fn test_eigen_qr_mat1() { + test_eigen_qr_impl!(Mat1); +} + +#[test] +fn test_eigen_qr_mat2() { + test_eigen_qr_impl!(Mat2); +} + +#[test] +fn test_eigen_qr_mat3() { + test_eigen_qr_impl!(Mat3); +} + +#[test] +fn test_eigen_qr_mat4() { + test_eigen_qr_impl!(Mat4); +} + +#[test] +fn test_eigen_qr_mat5() { + test_eigen_qr_impl!(Mat5); +} + +#[test] +fn test_eigen_qr_mat6() { + test_eigen_qr_impl!(Mat6); +} #[test] fn test_from_fn() { @@ -733,6 +753,36 @@ fn test_cholesky_mat6() { test_cholesky_impl!(Mat6); } +#[test] +fn test_hessenberg_mat1() { + test_hessenberg_impl!(Mat1); +} + +#[test] +fn test_hessenberg_mat2() { + test_hessenberg_impl!(Mat2); +} + +#[test] +fn test_hessenberg_mat3() { + test_hessenberg_impl!(Mat3); +} + +#[test] +fn test_hessenberg_mat4() { + test_hessenberg_impl!(Mat4); +} + +#[test] +fn test_hessenberg_mat5() { + test_hessenberg_impl!(Mat5); +} + +#[test] +fn test_hessenberg_mat6() { + test_hessenberg_impl!(Mat6); +} + #[test] fn test_transpose_square_mat() { let col_major_mat = &[0, 1, 2, 3,