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

needs faster matrix initialization

This commit is contained in:
Nestor Demeure 2019-11-03 20:00:15 +01:00
parent 46687b7cde
commit 72834186d8
1 changed files with 24 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

80
src/linalg/cholesky.rs
View File

@ -149,48 +149,17 @@ where
/// 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()`.
#[inline]
pub fn rank_one_update<R2: Dim, S2>(&mut self, x: &Vector<N, R2, S2>, sigma: N::RealField)
where
S2: Storage<N, R2, U1>,
DefaultAllocator: Allocator<N, R2, U1>,
ShapeConstraint: SameNumberOfRows<R2, D>,
{
// 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,
self.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 { *self.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 { *self.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 = self.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),
);
}
}
rank_one_update(&mut self.chol, x, sigma)
}
/// Updates the decomposition such that we get the decomposition of a matrix with the given column `c` 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.
/// Since the matrix is square, an identical row will be added in the `j`th row.
pub fn insert_column<R2, S2>(
self,
@ -206,37 +175,32 @@ where
{
// for an explanation of the formulas, see https://en.wikipedia.org/wiki/Cholesky_decomposition#Updating_the_decomposition
let n = col.nrows();
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.");
// TODO what is the fastest way to produce the new matrix ?
let mut chol= self.chol.clone().insert_column(j, N::zero()).insert_row(j, N::zero());
// update the jth row
let top_left_corner = chol.slice_range(..j, ..j);
let colj_minus = col.rows_range(..j);
let rowj = top_left_corner.solve_lower_triangular(&colj_minus).unwrap().adjoint(); // TODO both the row and its adjoint seem to be usefull
chol.slice_range_mut(j, ..j).copy_from(&rowj);
// TODO
//println!("dotc:{} norm2:{}", rowj.dotc(&rowj), rowj.norm_squared());
let top_left_corner = self.chol.slice_range(..j, ..j);
let col_jminus = col.rows_range(..j);
let new_rowj_adjoint = top_left_corner.solve_lower_triangular(&col_jminus).expect("Cholesky::insert_column : Unable to solve lower triangular system!");
new_rowj_adjoint.adjoint_to(&mut chol.slice_range_mut(j, ..j));
// update the center element
let center_element = N::sqrt(col[j] - rowj.dotc(&rowj) );
let center_element = N::sqrt(col[j] - N::from_real(new_rowj_adjoint.norm_squared()));
chol[(j,j)] = center_element;
// update the jth column
let colj_plus = col.rows_range(j+1..);
let bottom_left_corner = chol.slice_range(j+1.., ..j);
let colj = (colj_plus - bottom_left_corner*rowj.adjoint()) / center_element; // TODO that can probably be done with a single optimized operation
chol.slice_range_mut(j+1.., j).copy_from(&colj);
let bottom_left_corner = self.chol.slice_range(j.., ..j);
// new_colj = (col_jplus - bottom_left_corner * new_rowj.adjoint()) / center_element;
let mut new_colj = col.rows_range(j+1..).clone_owned();
new_colj.gemm(-N::one() / center_element, &bottom_left_corner, &new_rowj_adjoint, N::one() / center_element );
chol.slice_range_mut(j+1.., j).copy_from(&new_colj);
// update the bottom right corner
let mut bottom_right_corner = chol.slice_range_mut(j+1.., j+1..);
rank_one_update_helper(&mut bottom_right_corner, &colj, -N::real(N::one()));
rank_one_update(&mut bottom_right_corner, &new_colj, -N::real(N::one()));
Cholesky { chol }
}
@ -254,13 +218,14 @@ where
let n = self.chol.nrows();
assert!(n > 0, "The matrix needs at least one column.");
assert!(j < n, "j needs to be within the bound of the matrix.");
// TODO what is the fastest way to produce the new matrix ?
let mut chol= self.chol.clone().remove_column(j).remove_row(j);
// updates the bottom right corner
let mut corner = chol.slice_range_mut(j.., j..);
let colj = self.chol.slice_range(j+1.., j);
rank_one_update_helper(&mut corner, &colj, N::real(N::one()));
let mut bottom_right_corner = chol.slice_range_mut(j.., j..);
let old_colj = self.chol.slice_range(j+1.., j);
rank_one_update(&mut bottom_right_corner, &old_colj, N::real(N::one()));
Cholesky { chol }
}
@ -281,7 +246,10 @@ where
/// 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()`.
fn rank_one_update_helper<N, D, S, Rx, Sx>(chol : &mut Matrix<N, D, D, S>, x: &Vector<N, Rx, Sx>, sigma: N::RealField)
///
/// 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,