diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 755f5610..0485e9e5 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -4,9 +4,9 @@ use serde::{Deserialize, Serialize}; use alga::general::ComplexField; use crate::allocator::Allocator; -use crate::base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, SquareMatrix}; +use crate::base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, SquareMatrix, Vector}; use crate::constraint::{SameNumberOfRows, ShapeConstraint}; -use crate::dimension::{Dim, DimName, DimAdd, DimSum, DimDiff, DimSub, Dynamic, U1}; +use crate::dimension::{Dim, DimAdd, DimSum, DimDiff, DimSub, Dynamic, U1}; use crate::storage::{Storage, StorageMut}; use crate::base::allocator::Reallocator; @@ -149,7 +149,7 @@ 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()`. - pub fn rank_one_update(&mut self, x: &Matrix, sigma: N::RealField) + pub fn rank_one_update(&mut self, x: &Vector, sigma: N::RealField) where S2: Storage, DefaultAllocator: Allocator, @@ -192,17 +192,19 @@ where /// Updates the decomposition such that we get the decomposition of a matrix with the given column `c` 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( + pub fn insert_column( self, j: usize, - col: &Matrix, + col: &Vector, ) -> Cholesky> where D: DimAdd, - DefaultAllocator: Reallocator> + Reallocator, DimSum, DimSum>, + R2: Dim, S2: Storage, + DefaultAllocator: Reallocator> + Reallocator, DimSum, DimSum>, ShapeConstraint: SameNumberOfRows>, { + // for an explanation of the formulas, see https://en.wikipedia.org/wiki/Cholesky_decomposition#Updating_the_decomposition let n = col.nrows(); assert_eq!( n, @@ -211,7 +213,6 @@ where ); 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 ? - // TODO check for adjoint problems let mut chol= self.chol.clone().insert_column(j, N::zero()).insert_row(j, N::zero()); // update the jth row @@ -225,16 +226,15 @@ where chol[(j,j)] = center_element; // update the jth column - let colj_plus = col.rows_range(j+1..).adjoint(); - let bottom_left_corner = chol.slice_range(j+1, ..j-1); - let colj = (colj_plus - bottom_left_corner*rowj.adjoint()) / center_element; + let colj_plus = col.rows_range(j+1..); + let bottom_left_corner = chol.slice_range(j+1.., ..j-1); + 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); // update 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+1.., j+1..); rank_one_update_helper(&mut bottom_right_corner, &colj, -N::real(N::one())); - // TODO see https://en.wikipedia.org/wiki/Cholesky_decomposition#Updating_the_decomposition Cholesky { chol } } @@ -278,15 +278,14 @@ 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: &Matrix, sigma: N::RealField) +fn rank_one_update_helper(chol : &mut Matrix, x: &Vector, sigma: N::RealField) where N: ComplexField, - D: DimSub, - R2: Dim, + D: Dim, + Rx: Dim, S: StorageMut, - S2: Storage, - DefaultAllocator: Allocator + Allocator, - ShapeConstraint: SameNumberOfRows, + Sx: Storage, + DefaultAllocator: Allocator, { // heavily inspired by Eigen's `llt_rank_update_lower` implementation https://eigen.tuxfamily.org/dox/LLT_8h_source.html let n = x.nrows();