diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index bbd233eb..12674e4c 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -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(&mut self, x: &Vector, sigma: N::RealField) where S2: Storage, DefaultAllocator: Allocator, ShapeConstraint: SameNumberOfRows, { - // 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::(); - 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::() { - // 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( 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(chol : &mut Matrix, x: &Vector, 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(chol : &mut Matrix, x: &Vector, sigma: N::RealField) where N: ComplexField, D: Dim,