srenblad/nalgebra
1
0
Fork 0
forked from M-Labs/nalgebra
Code Pull Requests Activity

finished cleaning

Browse Source
This commit is contained in:
Nestor Demeure 2019-11-03 21:24:44 +01:00 committed by Sébastien Crozet
parent 3d08a80d8d
commit 59c6a98615
2 changed files with 69 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

127
src/linalg/cholesky.rs
View File

@ -8,7 +8,6 @@ use crate::base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, SquareMatrix, Vec
use crate::constraint::{SameNumberOfRows, ShapeConstraint}; use crate::constraint::{SameNumberOfRows, ShapeConstraint};
use crate::dimension::{Dim, DimAdd, DimSum, DimDiff, DimSub, Dynamic, U1}; use crate::dimension::{Dim, DimAdd, DimSum, DimDiff, DimSub, Dynamic, U1};
use crate::storage::{Storage, StorageMut}; use crate::storage::{Storage, StorageMut};
use crate::base::allocator::Reallocator;
/// The Cholesky decomposition of a symmetric-definite-positive matrix. /// The Cholesky decomposition of a symmetric-definite-positive matrix.
#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
@ -156,7 +155,7 @@ where
DefaultAllocator: Allocator<N, R2, U1>, DefaultAllocator: Allocator<N, R2, U1>,
ShapeConstraint: SameNumberOfRows<R2, D>, ShapeConstraint: SameNumberOfRows<R2, D>,
{ {
rank_one_update(&mut self.chol, x, sigma) Self::xx_rank_one_update(&mut self.chol, x, sigma)
} }
/// Updates the decomposition such that we get the decomposition of a matrix with the given column `col` in the `j`th position. /// Updates the decomposition such that we get the decomposition of a matrix with the given column `col` in the `j`th position.
@ -170,7 +169,7 @@ where
D: DimAdd<U1>, D: DimAdd<U1>,
R2: Dim, R2: Dim,
S2: Storage<N, R2, U1>, S2: Storage<N, R2, U1>,
DefaultAllocator: Reallocator<N, D, D, D, DimSum<D, U1>> + Reallocator<N, D, DimSum<D, U1>, DimSum<D, U1>, DimSum<D, U1>>, DefaultAllocator: Allocator<N, DimSum<D, U1>, DimSum<D, U1>>,
ShapeConstraint: SameNumberOfRows<R2, DimSum<D, U1>>, ShapeConstraint: SameNumberOfRows<R2, DimSum<D, U1>>,
{ {
// for an explanation of the formulas, see https://en.wikipedia.org/wiki/Cholesky_decomposition#Updating_the_decomposition // for an explanation of the formulas, see https://en.wikipedia.org/wiki/Cholesky_decomposition#Updating_the_decomposition
@ -178,8 +177,12 @@ where
assert_eq!(n, self.chol.nrows() + 1, "The new column must have the size of the factored matrix plus one."); assert_eq!(n, self.chol.nrows() + 1, "The new column must have the size of the factored matrix plus one.");
assert!(j < n, "j needs to be within the bound of the new matrix."); assert!(j < n, "j needs to be within the bound of the new matrix.");
// TODO what is the fastest way to produce the new matrix ? // loads the data into a new matrix with an additional jth row/column
let mut chol= self.chol.clone().insert_column(j, N::zero()).insert_row(j, N::zero()); let mut chol = unsafe { Matrix::new_uninitialized_generic(self.chol.data.shape().0.add(U1), self.chol.data.shape().1.add(U1)) };
chol.slice_range_mut(..j, ..j).copy_from(&self.chol.slice_range(..j, ..j));
chol.slice_range_mut(..j, j+1..).copy_from(&self.chol.slice_range(..j, j..));
chol.slice_range_mut(j+1.., ..j).copy_from(&self.chol.slice_range(j.., ..j));
chol.slice_range_mut(j+1.., j+1..).copy_from(&self.chol.slice_range(j.., j..));
// update the jth row // update the jth row
let top_left_corner = self.chol.slice_range(..j, ..j); let top_left_corner = self.chol.slice_range(..j, ..j);
@ -200,7 +203,7 @@ where
// update the bottom right corner // update the bottom right corner
let mut bottom_right_corner = chol.slice_range_mut(j+1.., j+1..); let mut bottom_right_corner = chol.slice_range_mut(j+1.., j+1..);
rank_one_update(&mut bottom_right_corner, &new_colj, -N::real(N::one())); Self::xx_rank_one_update(&mut bottom_right_corner, &new_colj, -N::real(N::one()));
Cholesky { chol } Cholesky { chol }
} }
@ -208,27 +211,80 @@ where
/// Updates the decomposition such that we get the decomposition of the factored matrix with its `j`th column removed. /// Updates the decomposition such that we get the decomposition of the factored matrix with its `j`th column removed.
/// Since the matrix is square, the `j`th row will also be removed. /// Since the matrix is square, the `j`th row will also be removed.
pub fn remove_column( pub fn remove_column(
self, &self,
j: usize, j: usize,
) -> Cholesky<N, DimDiff<D, U1>> ) -> Cholesky<N, DimDiff<D, U1>>
where where
D: DimSub<U1>, D: DimSub<U1>,
DefaultAllocator: Reallocator<N, D, D, D, DimDiff<D, U1>> + Reallocator<N, D, DimDiff<D, U1>, DimDiff<D, U1>, DimDiff<D, U1>>, DefaultAllocator: Allocator<N, DimDiff<D, U1>, DimDiff<D, U1>>
{ {
let n = self.chol.nrows(); let n = self.chol.nrows();
assert!(n > 0, "The matrix needs at least one column."); assert!(n > 0, "The matrix needs at least one column.");
assert!(j < n, "j needs to be within the bound of the matrix."); assert!(j < n, "j needs to be within the bound of the matrix.");
// TODO what is the fastest way to produce the new matrix ? // loads the data into a new matrix except for the jth row/column
let mut chol= self.chol.clone().remove_column(j).remove_row(j); let mut chol = unsafe { Matrix::new_uninitialized_generic(self.chol.data.shape().0.sub(U1), self.chol.data.shape().1.sub(U1)) };
chol.slice_range_mut(..j, ..j).copy_from(&self.chol.slice_range(..j, ..j));
chol.slice_range_mut(..j, j..).copy_from(&self.chol.slice_range(..j, j+1..));
chol.slice_range_mut(j.., ..j).copy_from(&self.chol.slice_range(j+1.., ..j));
chol.slice_range_mut(j.., j..).copy_from(&self.chol.slice_range(j+1.., j+1..));
// updates the bottom right corner // updates the bottom right corner
let mut bottom_right_corner = chol.slice_range_mut(j.., j..); let mut bottom_right_corner = chol.slice_range_mut(j.., j..);
let old_colj = self.chol.slice_range(j+1.., j); let old_colj = self.chol.slice_range(j+1.., j);
rank_one_update(&mut bottom_right_corner, &old_colj, N::real(N::one())); Self::xx_rank_one_update(&mut bottom_right_corner, &old_colj, N::real(N::one()));
Cholesky { chol } Cholesky { chol }
} }
/// Given the Cholesky decomposition of a matrix `M`, a scalar `sigma` and a vector `v`,
/// performs a rank one update such that we end up with the decomposition of `M + sigma * v*v.adjoint()`.
///
/// This helper method is calling for by `rank_one_update` but also `insert_column` and `remove_column`
/// where it is used on a square slice of the decomposition
fn xx_rank_one_update<Dm, Sm, Rx, Sx>(chol : &mut Matrix<N, Dm, Dm, Sm>, x: &Vector<N, Rx, Sx>, sigma: N::RealField)
where
//N: ComplexField,
Dm: Dim,
Rx: Dim,
Sm: StorageMut<N, Dm, Dm>,
Sx: Storage<N, Rx, U1>,
DefaultAllocator: Allocator<N, Rx, U1>,
{
// heavily inspired by Eigen's `llt_rank_update_lower` implementation https://eigen.tuxfamily.org/dox/LLT_8h_source.html
let n = x.nrows();
assert_eq!(
n,
chol.nrows(),
"The input vector must be of the same size as the factorized matrix."
);
let mut x = x.clone_owned();
let mut beta = crate::one::<N::RealField>();
for j in 0..n {
// updates the diagonal
let diag = N::real(unsafe { *chol.get_unchecked((j, j)) });
let diag2 = diag * diag;
let xj = unsafe { *x.get_unchecked(j) };
let sigma_xj2 = sigma * N::modulus_squared(xj);
let gamma = diag2 * beta + sigma_xj2;
let new_diag = (diag2 + sigma_xj2 / beta).sqrt();
unsafe { *chol.get_unchecked_mut((j, j)) = N::from_real(new_diag) };
beta += sigma_xj2 / diag2;
// updates the terms of L
let mut xjplus = x.rows_range_mut(j + 1..);
let mut col_j = chol.slice_range_mut(j + 1.., j);
// temp_jplus -= (wj / N::from_real(diag)) * col_j;
xjplus.axpy(-xj / N::from_real(diag), &col_j, N::one());
if gamma != crate::zero::<N::RealField>() {
// col_j = N::from_real(nljj / diag) * col_j + (N::from_real(nljj * sigma / gamma) * N::conjugate(wj)) * temp_jplus;
col_j.axpy(
N::from_real(new_diag * sigma / gamma) * N::conjugate(xj),
&xjplus,
N::from_real(new_diag / diag),
);
}
}
}
} }
impl<N: ComplexField, D: DimSub<Dynamic>, S: Storage<N, D, D>> SquareMatrix<N, D, S> impl<N: ComplexField, D: DimSub<Dynamic>, S: Storage<N, D, D>> SquareMatrix<N, D, S>
@ -243,52 +299,3 @@ where
Cholesky::new(self.into_owned()) Cholesky::new(self.into_owned())
} }
} }
/// Given the Cholesky decomposition of a matrix `M`, a scalar `sigma` and a vector `v`,
/// performs a rank one update such that we end up with the decomposition of `M + sigma * v*v.adjoint()`.
///
/// This helper method is calling for by `rank_one_update` but also `insert_column` and `remove_column`
/// where it is used on a square slice of the decomposition
fn rank_one_update<N, D, S, Rx, Sx>(chol : &mut Matrix<N, D, D, S>, x: &Vector<N, Rx, Sx>, sigma: N::RealField)
where
N: ComplexField,
D: Dim,
Rx: Dim,
S: StorageMut<N, D, D>,
Sx: Storage<N, Rx, U1>,
DefaultAllocator: Allocator<N, Rx, U1>,
{
// heavily inspired by Eigen's `llt_rank_update_lower` implementation https://eigen.tuxfamily.org/dox/LLT_8h_source.html
let n = x.nrows();
assert_eq!(
n,
chol.nrows(),
"The input vector must be of the same size as the factorized matrix."
);
let mut x = x.clone_owned();
let mut beta = crate::one::<N::RealField>();
for j in 0..n {
// updates the diagonal
let diag = N::real(unsafe { *chol.get_unchecked((j, j)) });
let diag2 = diag * diag;
let xj = unsafe { *x.get_unchecked(j) };
let sigma_xj2 = sigma * N::modulus_squared(xj);
let gamma = diag2 * beta + sigma_xj2;
let new_diag = (diag2 + sigma_xj2 / beta).sqrt();
unsafe { *chol.get_unchecked_mut((j, j)) = N::from_real(new_diag) };
beta += sigma_xj2 / diag2;
// updates the terms of L
let mut xjplus = x.rows_range_mut(j + 1..);
let mut col_j = chol.slice_range_mut(j + 1.., j);
// temp_jplus -= (wj / N::from_real(diag)) * col_j;
xjplus.axpy(-xj / N::from_real(diag), &col_j, N::one());
if gamma != crate::zero::<N::RealField>() {
// col_j = N::from_real(nljj / diag) * col_j + (N::from_real(nljj * sigma / gamma) * N::conjugate(wj)) * temp_jplus;
col_j.axpy(
N::from_real(new_diag * sigma / gamma) * N::conjugate(xj),
&xjplus,
N::from_real(new_diag / diag),
);
}
}
}

8
tests/linalg/cholesky.rs
View File

@ -100,7 +100,7 @@ macro_rules! gen_tests(
} }
fn cholesky_insert_column(n: usize) -> bool { fn cholesky_insert_column(n: usize) -> bool {
let n = n.max(1).min(50); let n = n.max(1).min(10);
let j = random::<usize>() % n; let j = random::<usize>() % n;
let m_updated = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap(); let m_updated = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap();
@ -112,15 +112,11 @@ macro_rules! gen_tests(
let chol = m.clone().cholesky().unwrap().insert_column(j, &col); let chol = m.clone().cholesky().unwrap().insert_column(j, &col);
let m_chol_updated = chol.l() * chol.l().adjoint(); let m_chol_updated = chol.l() * chol.l().adjoint();
println!("n={} j={}", n, j);
println!("chol updated:{}", m_chol_updated);
println!("m updated:{}", m_updated);
relative_eq!(m_updated, m_chol_updated, epsilon = 1.0e-7) relative_eq!(m_updated, m_chol_updated, epsilon = 1.0e-7)
} }
fn cholesky_remove_column(n: usize) -> bool { fn cholesky_remove_column(n: usize) -> bool {
let n = n.max(1).min(5); let n = n.max(1).min(10);
let j = random::<usize>() % n; let j = random::<usize>() % n;
let m = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap(); let m = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap();