QR decomposition depends less on DMat internals

This commit is contained in:
Vincent Barrielle 2014-05-09 22:14:37 +02:00
parent 2fd880a62d
commit 5611307b4d
8 changed files with 116 additions and 63 deletions

View File

@ -119,6 +119,7 @@ extern crate test;
pub mod na; pub mod na;
mod structs; mod structs;
mod traits; mod traits;
mod linalg;
// mod lower_triangular; // mod lower_triangular;
// mod chol; // mod chol;

View File

@ -0,0 +1,58 @@
use std::num::{Zero, Float};
use na::DVec;
use na::DMat;
use traits::operations::{Transpose, ColSlice};
use traits::geometry::Norm;
use std::cmp::min;
/// QR decomposition using Householder reflections
/// # Arguments
/// * `m` matrix to decompose
pub fn decomp_qr<N: Clone + Num + Float>(m: &DMat<N>) -> (DMat<N>, DMat<N>) {
let rows = m.nrows();
let cols = m.ncols();
assert!(rows >= cols);
let mut q : DMat<N> = DMat::new_identity(rows);
let mut r = m.clone();
let subtract_reflection = |vec: DVec<N>| -> DMat<N> {
// FIXME: we don't handle the complex case here
let mut qk : DMat<N> = DMat::new_identity(rows);
let start = rows - vec.at.len();
for j in range(start, rows) {
for i in range(start, rows) {
unsafe {
let vv = vec.at_fast(i-start)*vec.at_fast(j-start);
let qkij = qk.at_fast(i,j);
qk.set_fast(i, j, qkij - vv - vv);
}
}
}
qk
};
let iterations = min(rows-1, cols);
for ite in range(0u, iterations) {
let mut v = r.col_slice(ite, ite, rows);
//let mut v = r.col_slice<DVec<N>>(ite, rows-ite, rows);
let alpha =
if unsafe { v.at_fast(ite) } >= Zero::zero() {
-Norm::norm(&v)
}
else {
Norm::norm(&v)
};
unsafe {
let x = v.at_fast(0);
v.set_fast(0, x - alpha);
}
let _ = v.normalize();
let qk = subtract_reflection(v);
r = qk * r;
q = q * Transpose::transpose_cpy(&qk);
}
(q, r)
}

4
src/linalg/mod.rs Normal file
View File

@ -0,0 +1,4 @@
pub use self::decompositions::decomp_qr;
mod decompositions;

View File

@ -39,12 +39,12 @@ pub use traits::{
Transpose, Transpose,
UniformSphereSample, UniformSphereSample,
AnyVec, AnyVec,
VecExt VecExt,
ColSlice, RowSlice
}; };
pub use structs::{ pub use structs::{
Identity, Identity,
decomp_qr,
DMat, DVec, DMat, DVec,
Iso2, Iso3, Iso4, Iso2, Iso3, Iso4,
Mat1, Mat2, Mat3, Mat4, Mat1, Mat2, Mat3, Mat4,
@ -53,6 +53,10 @@ pub use structs::{
Vec0, Vec1, Vec2, Vec3, Vec4, Vec5, Vec6 Vec0, Vec1, Vec2, Vec3, Vec4, Vec5, Vec6
}; };
pub use linalg::{
decomp_qr
};
/// Traits to work around the language limitations related to operator overloading. /// Traits to work around the language limitations related to operator overloading.
/// ///
/// The trait names are formed by: /// The trait names are formed by:

View File

@ -4,14 +4,12 @@
use rand::Rand; use rand::Rand;
use rand; use rand;
use std::num::{One, Zero, Float}; use std::num::{One, Zero};
use traits::operations::ApproxEq; use traits::operations::ApproxEq;
use std::mem; use std::mem;
use structs::dvec::{DVec, DVecMulRhs}; use structs::dvec::{DVec, DVecMulRhs};
use traits::operations::{Inv, Transpose, Mean, Cov}; use traits::operations::{Inv, Transpose, Mean, Cov, ColSlice, RowSlice};
use traits::structure::Cast; use traits::structure::Cast;
use traits::geometry::Norm;
use std::cmp::min;
use std::fmt::{Show, Formatter, Result}; use std::fmt::{Show, Formatter, Result};
#[doc(hidden)] #[doc(hidden)]
@ -497,61 +495,38 @@ impl<N: Clone + Num + Cast<f32> + DMatDivRhs<N, DMat<N>> + ToStr > Cov<DMat<N>>
} }
} }
impl<N: Clone> ColSlice<DVec<N>> for DMat<N> {
fn col_slice(&self, col_id :uint, row_start: uint, row_end: uint) -> DVec<N> {
assert!(col_id < self.ncols);
assert!(row_start < row_end);
assert!(row_end <= self.nrows);
// we can init from slice thanks to the matrix being column major
let start= self.offset(row_start, col_id);
let stop = self.offset(row_end, col_id);
let slice = DVec::from_vec(
row_end - row_start, self.mij.slice(start, stop));
slice
}
}
/// QR decomposition using Householder reflections impl<N: Clone> RowSlice<DVec<N>> for DMat<N> {
/// # Arguments fn row_slice(&self, row_id :uint, col_start: uint, col_end: uint) -> DVec<N> {
/// * `m` matrix to decompose assert!(row_id < self.nrows);
pub fn decomp_qr<N: Clone + Num + Float + Show>(m: &DMat<N>) -> (DMat<N>, DMat<N>) { assert!(col_start < col_end);
let rows = m.nrows(); assert!(col_end <= self.ncols);
let cols = m.ncols(); let mut slice : DVec<N> = unsafe {
assert!(rows >= cols); DVec::new_uninitialized(self.nrows)
let mut q : DMat<N> = DMat::new_identity(rows);
let mut r = m.clone();
let subtract_reflection = |vec: DVec<N>| -> DMat<N> {
// FIXME: we don't handle the complex case here
let mut qk : DMat<N> = DMat::new_identity(rows);
let start = rows - vec.at.len();
for j in range(start, rows) {
for i in range(start, rows) {
unsafe {
let vv = vec.at_fast(i-start)*vec.at_fast(j-start);
let qkij = qk.at_fast(i,j);
qk.set_fast(i, j, qkij - vv - vv);
}
}
}
qk
};
let iterations = min(rows-1, cols);
for ite in range(0u, iterations) {
// we get the ite-th column truncated from its first ite elements,
// this is fast thanks to the matrix being column major
let start= m.offset(ite, ite);
let stop = m.offset(rows, ite);
let mut v = DVec::from_vec(rows - ite, r.mij.slice(start, stop));
let alpha =
if unsafe { v.at_fast(ite) } >= Zero::zero() {
-Norm::norm(&v)
}
else {
Norm::norm(&v)
}; };
let mut slice_idx = 0u;
for col_id in range(col_start, col_end) {
unsafe { unsafe {
let x = v.at_fast(0); slice.set_fast(slice_idx, self.at_fast(row_id, col_id));
v.set_fast(0, x - alpha);
} }
let _ = v.normalize(); slice_idx += 1;
let qk = subtract_reflection(v); }
r = qk * r; slice
q = q * Transpose::transpose_cpy(&qk);
} }
(q, r)
} }
impl<N: ApproxEq<N>> ApproxEq<N> for DMat<N> { impl<N: ApproxEq<N>> ApproxEq<N> for DMat<N> {
#[inline] #[inline]
@ -578,9 +553,9 @@ impl<N: Show + Clone> Show for DMat<N> {
fn fmt(&self, form:&mut Formatter) -> Result { fn fmt(&self, form:&mut Formatter) -> Result {
for i in range(0u, self.nrows()) { for i in range(0u, self.nrows()) {
for j in range(0u, self.ncols()) { for j in range(0u, self.ncols()) {
write!(form.buf, "{} ", self.at(i, j)); let _ = write!(form.buf, "{} ", self.at(i, j));
} }
write!(form.buf, "\n"); let _ = write!(form.buf, "\n");
} }
write!(form.buf, "\n") write!(form.buf, "\n")
} }

View File

@ -1,7 +1,6 @@
//! Data structures and implementations. //! Data structures and implementations.
pub use self::dmat::DMat; pub use self::dmat::DMat;
pub use self::dmat::decomp_qr;
pub use self::dvec::DVec; pub use self::dvec::DVec;
pub use self::vec::{Vec0, Vec1, Vec2, Vec3, Vec4, Vec5, Vec6}; pub use self::vec::{Vec0, Vec1, Vec2, Vec3, Vec4, Vec5, Vec6};
pub use self::mat::{Identity, Mat1, Mat2, Mat3, Mat4, Mat5, Mat6}; pub use self::mat::{Identity, Mat1, Mat2, Mat3, Mat4, Mat5, Mat6};

View File

@ -8,7 +8,7 @@ pub use self::structure::{FloatVec, FloatVecExt, Basis, Cast, Col, Dim, Indexabl
Iterable, IterableMut, Mat, Row, AnyVec, VecExt}; Iterable, IterableMut, Mat, Row, AnyVec, VecExt};
pub use self::operations::{Absolute, ApproxEq, Cov, Inv, LMul, Mean, Outer, PartialOrd, RMul, pub use self::operations::{Absolute, ApproxEq, Cov, Inv, LMul, Mean, Outer, PartialOrd, RMul,
ScalarAdd, ScalarSub, Transpose}; ScalarAdd, ScalarSub, Transpose, ColSlice, RowSlice};
pub use self::operations::{PartialOrdering, PartialLess, PartialEqual, PartialGreater, NotComparable}; pub use self::operations::{PartialOrdering, PartialLess, PartialEqual, PartialGreater, NotComparable};
pub mod geometry; pub mod geometry;

View File

@ -244,6 +244,18 @@ pub trait Mean<N> {
fn mean(&Self) -> N; fn mean(&Self) -> N;
} }
/// Trait for objects that support column slicing
pub trait ColSlice<VecLike> {
/// Returns a view to a slice of a column of a matrix.
fn col_slice(&self, col_id :uint, row_start: uint, row_end: uint) -> VecLike;
}
/// Trait for objects that support column slicing
pub trait RowSlice<VecLike> {
/// Returns a view to a slice of a row of a matrix.
fn row_slice(&self, row_id :uint, col_start: uint, col_end: uint) -> VecLike;
}
// /// Cholesky decomposition. // /// Cholesky decomposition.
// pub trait Chol { // pub trait Chol {
// /// Performs the cholesky decomposition on `self`. The resulting upper-triangular matrix is // /// Performs the cholesky decomposition on `self`. The resulting upper-triangular matrix is