From 2beb09dab29da70e41c6172b9d8d22a143e56487 Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sat, 2 Nov 2019 14:59:07 +0100 Subject: [PATCH 01/23] first version of rank one update --- src/linalg/cholesky.rs | 58 ++++++++++++++++++++++++++++++++---------- 1 file changed, 44 insertions(+), 14 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 0b6e6db5..606434e9 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -6,28 +6,25 @@ use alga::general::ComplexField; use crate::allocator::Allocator; use crate::base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, SquareMatrix}; use crate::constraint::{SameNumberOfRows, ShapeConstraint}; -use crate::dimension::{Dim, DimSub, Dynamic}; +use crate::dimension::{Dim, DimSub, Dynamic, U1}; use crate::storage::{Storage, StorageMut}; /// The Cholesky decomposition of a symmetric-definite-positive matrix. #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] #[cfg_attr( feature = "serde-serialize", - serde(bound( - serialize = "DefaultAllocator: Allocator, - MatrixN: Serialize" - )) + serde(bound(serialize = "DefaultAllocator: Allocator, + MatrixN: Serialize")) )] #[cfg_attr( feature = "serde-serialize", - serde(bound( - deserialize = "DefaultAllocator: Allocator, - MatrixN: Deserialize<'de>" - )) + serde(bound(deserialize = "DefaultAllocator: Allocator, + MatrixN: Deserialize<'de>")) )] #[derive(Clone, Debug)] pub struct Cholesky -where DefaultAllocator: Allocator +where + DefaultAllocator: Allocator, { chol: MatrixN, } @@ -36,10 +33,12 @@ impl Copy for Cholesky where DefaultAllocator: Allocator, MatrixN: Copy, -{} +{ +} impl> Cholesky -where DefaultAllocator: Allocator +where + DefaultAllocator: Allocator, { /// Attempts to compute the Cholesky decomposition of `matrix`. /// @@ -129,7 +128,7 @@ where DefaultAllocator: Allocator /// `x` the unknown. pub fn solve(&self, b: &Matrix) -> MatrixMN where - S2: StorageMut, + S2: Storage, DefaultAllocator: Allocator, ShapeConstraint: SameNumberOfRows, { @@ -146,10 +145,41 @@ where DefaultAllocator: Allocator self.solve_mut(&mut res); res } + + /// Performs a rank one update of the current decomposition. + /// If `M = L * L^T` before the rank one update, then after it we have `L*L^T = M + sigma * v*v^T` where v must be a vector of same dimension. + /// TODO rewrite comment (current version is taken verbatim from eigen) + /// TODO insures that code is correct for complex numbers, eigen uses abs2 and conj + /// https://eigen.tuxfamily.org/dox/LLT_8h_source.html + pub fn rank_one_update(&mut self, x: &Matrix, sigma: N) + where + S2: Storage, + DefaultAllocator: Allocator, + ShapeConstraint: SameNumberOfRows, + { + let n = x.nrows(); + let mut temp = x.clone_owned(); + for k in 0..n { + let lkk = self.chol[(k, k)]; // TODO unsafe { *matrix.get_unchecked((j, j)) } + let xk = temp[k]; + let r = (lkk * lkk + sigma * xk * xk).sqrt(); + let c = r / lkk; + let s = xk / lkk; + self.chol[(k, k)] = r; + // Update the terms of L + if k < n { + for k2 in (k + 1)..n { + self.chol[(k2, k)] = (self.chol[(k2, k)] + sigma * s * temp[k2]) / c; + temp[k2] = c * temp[k2] - s * self.chol[(k2, k)]; + } + } + } + } } impl, S: Storage> SquareMatrix -where DefaultAllocator: Allocator +where + DefaultAllocator: Allocator, { /// Attempts to compute the Cholesky decomposition of this matrix. /// From 1fe4ef956b0d87762e05f1237d690704b3142b31 Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sat, 2 Nov 2019 15:11:14 +0100 Subject: [PATCH 02/23] added test for update --- tests/linalg/cholesky.rs | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs index cefc2630..52b1198c 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -1,6 +1,5 @@ #![cfg(all(feature = "arbitrary", feature = "debug"))] - macro_rules! gen_tests( ($module: ident, $scalar: ty) => { mod $module { @@ -78,6 +77,22 @@ macro_rules! gen_tests( id1.is_identity(1.0e-7) && id2.is_identity(1.0e-7) } + + fn cholesky_rank_one_update(_n: usize) -> bool { + let m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap(); + let x = Vector4::<$scalar>::new_random().map(|e| e.0); + let sigma : $scalar = 1.; + + // updates m manually + let m_updated = m + sigma * x * x.transpose(); + + // updates cholesky deomposition and reconstruct m + let mut chol = m.clone().cholesky().unwrap(); + chol.rank_one_update(x, sigma); + let m_chol_updated = chol.l() * chol.l().transpose(); + + relative_eq!(m_updated, m_chol_updated, epsilon = 1.0e-7) + } } } } From bdccc81874b2b26646428a4ce3785dc511e382fb Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sat, 2 Nov 2019 15:56:59 +0100 Subject: [PATCH 03/23] got test to compile --- src/linalg/cholesky.rs | 2 +- tests/linalg/cholesky.rs | 19 +++++++++++-------- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 606434e9..6a2c9da8 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -151,7 +151,7 @@ where /// TODO rewrite comment (current version is taken verbatim from eigen) /// TODO insures that code is correct for complex numbers, eigen uses abs2 and conj /// https://eigen.tuxfamily.org/dox/LLT_8h_source.html - pub fn rank_one_update(&mut self, x: &Matrix, sigma: N) + pub fn rank_one_update(&mut self, x: &Matrix, sigma: N) where S2: Storage, DefaultAllocator: Allocator, diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs index 52b1198c..80a54e2f 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -79,19 +79,22 @@ macro_rules! gen_tests( } fn cholesky_rank_one_update(_n: usize) -> bool { - let m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap(); + let mut m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap(); let x = Vector4::<$scalar>::new_random().map(|e| e.0); - let sigma : $scalar = 1.; + let sigma = random::<$scalar>().0; // random::<$scalar>().0; + let one = sigma*0. + 1.; // TODO this is dirty but $scalar appears to not be a scalar type + + // updates cholesky decomposition and reconstructs m + let mut chol = m.clone().cholesky().unwrap(); + chol.rank_one_update(&x, sigma); + let m_chol_updated = chol.l() * chol.l().adjoint(); // updates m manually - let m_updated = m + sigma * x * x.transpose(); + m.syger(sigma, &x, &x, one); // m += sigma * x * x.adjoint() - // updates cholesky deomposition and reconstruct m - let mut chol = m.clone().cholesky().unwrap(); - chol.rank_one_update(x, sigma); - let m_chol_updated = chol.l() * chol.l().transpose(); + println!("m : {:?}", m); - relative_eq!(m_updated, m_chol_updated, epsilon = 1.0e-7) + relative_eq!(m, m_chol_updated, epsilon = 1.0e-7) } } } From e49ecdc0a10a9167a544b43aac942b335078b4eb Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sat, 2 Nov 2019 16:36:23 +0100 Subject: [PATCH 04/23] test is now correct --- tests/linalg/cholesky.rs | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs index 80a54e2f..823ec96f 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -79,10 +79,13 @@ macro_rules! gen_tests( } fn cholesky_rank_one_update(_n: usize) -> bool { - let mut m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap(); - let x = Vector4::<$scalar>::new_random().map(|e| e.0); - let sigma = random::<$scalar>().0; // random::<$scalar>().0; - let one = sigma*0. + 1.; // TODO this is dirty but $scalar appears to not be a scalar type + use nalgebra::dimension::U3; + use nalgebra::Vector3; + let mut m = RandomSDP::new(U3, || random::<$scalar>().0).unwrap(); + let x = Vector3::<$scalar>::new_random().map(|e| e.0); + let mut sigma = random::<$scalar>().0; // random::<$scalar>().0; + let one = sigma*0. + 1.; // TODO this is dirty but $scalar appears to not be a scalar type in this file + sigma = one; // TODO placeholder // updates cholesky decomposition and reconstructs m let mut chol = m.clone().cholesky().unwrap(); @@ -90,9 +93,11 @@ macro_rules! gen_tests( let m_chol_updated = chol.l() * chol.l().adjoint(); // updates m manually - m.syger(sigma, &x, &x, one); // m += sigma * x * x.adjoint() + m.ger(sigma, &x, &x, one); // m += sigma * x * x.adjoint() - println!("m : {:?}", m); + println!("sigma : {}", sigma); + println!("m updated : {}", m); + println!("chol : {}", m_chol_updated); relative_eq!(m, m_chol_updated, epsilon = 1.0e-7) } From 80d2efcb5d4ef801bc403a89c6a9c0ec4790b72c Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sat, 2 Nov 2019 16:45:30 +0100 Subject: [PATCH 05/23] added real constraint on sigma --- src/linalg/cholesky.rs | 11 +++++++++-- tests/linalg/cholesky.rs | 11 +++++++---- 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 6a2c9da8..cbbe5ff8 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -8,6 +8,7 @@ use crate::base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, SquareMatrix}; use crate::constraint::{SameNumberOfRows, ShapeConstraint}; use crate::dimension::{Dim, DimSub, Dynamic, U1}; use crate::storage::{Storage, StorageMut}; +use crate::RealField; /// The Cholesky decomposition of a symmetric-definite-positive matrix. #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] @@ -151,12 +152,18 @@ where /// TODO rewrite comment (current version is taken verbatim from eigen) /// TODO insures that code is correct for complex numbers, eigen uses abs2 and conj /// https://eigen.tuxfamily.org/dox/LLT_8h_source.html - pub fn rank_one_update(&mut self, x: &Matrix, sigma: N) - where + /// TODO insure that sigma is a real + pub fn rank_one_update( + &mut self, + x: &Matrix, + sigma: N2, + ) where + N: From, S2: Storage, DefaultAllocator: Allocator, ShapeConstraint: SameNumberOfRows, { + let sigma = ::from(sigma); let n = x.nrows(); let mut temp = x.clone_owned(); for k in 0..n { diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs index 823ec96f..bef5de95 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -83,9 +83,12 @@ macro_rules! gen_tests( use nalgebra::Vector3; let mut m = RandomSDP::new(U3, || random::<$scalar>().0).unwrap(); let x = Vector3::<$scalar>::new_random().map(|e| e.0); - let mut sigma = random::<$scalar>().0; // random::<$scalar>().0; - let one = sigma*0. + 1.; // TODO this is dirty but $scalar appears to not be a scalar type in this file - sigma = one; // TODO placeholder + + // TODO this is dirty but $scalar appears to not be a scalar type in this file + let zero = random::<$scalar>().0 * 0.; + let one = zero + 1.; + let sigma = random::(); // needs to be a real + let sigma_scalar = zero + sigma; // updates cholesky decomposition and reconstructs m let mut chol = m.clone().cholesky().unwrap(); @@ -93,7 +96,7 @@ macro_rules! gen_tests( let m_chol_updated = chol.l() * chol.l().adjoint(); // updates m manually - m.ger(sigma, &x, &x, one); // m += sigma * x * x.adjoint() + m.ger(sigma_scalar, &x, &x, one); // m += sigma * x * x.adjoint() println!("sigma : {}", sigma); println!("m updated : {}", m); From bc7935ac24d7379b6693ae3b70fd522525a827b9 Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sat, 2 Nov 2019 16:49:57 +0100 Subject: [PATCH 06/23] updated comment --- src/linalg/cholesky.rs | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index cbbe5ff8..4063b4ee 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -147,12 +147,10 @@ where res } - /// Performs a rank one update of the current decomposition. - /// If `M = L * L^T` before the rank one update, then after it we have `L*L^T = M + sigma * v*v^T` where v must be a vector of same dimension. - /// TODO rewrite comment (current version is taken verbatim from eigen) + /// 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^*`. /// TODO insures that code is correct for complex numbers, eigen uses abs2 and conj /// https://eigen.tuxfamily.org/dox/LLT_8h_source.html - /// TODO insure that sigma is a real pub fn rank_one_update( &mut self, x: &Matrix, From 3ae88127eee37995cf9cae64b6a4364aea92a2fc Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sat, 2 Nov 2019 18:27:01 +0100 Subject: [PATCH 07/23] rank update passed tests --- src/linalg/cholesky.rs | 43 ++++++++++++++++++++-------------------- tests/linalg/cholesky.rs | 12 +++-------- 2 files changed, 25 insertions(+), 30 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 4063b4ee..c4049504 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -8,7 +8,6 @@ use crate::base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, SquareMatrix}; use crate::constraint::{SameNumberOfRows, ShapeConstraint}; use crate::dimension::{Dim, DimSub, Dynamic, U1}; use crate::storage::{Storage, StorageMut}; -use crate::RealField; /// The Cholesky decomposition of a symmetric-definite-positive matrix. #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] @@ -149,33 +148,35 @@ 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^*`. - /// TODO insures that code is correct for complex numbers, eigen uses abs2 and conj - /// https://eigen.tuxfamily.org/dox/LLT_8h_source.html - pub fn rank_one_update( - &mut self, - x: &Matrix, - sigma: N2, - ) where - N: From, + pub fn rank_one_update(&mut self, x: &Matrix, sigma: N::RealField) + where S2: Storage, DefaultAllocator: Allocator, ShapeConstraint: SameNumberOfRows, { - let sigma = ::from(sigma); + // for a description of the operation, see https://en.wikipedia.org/wiki/Cholesky_decomposition#Updating_the_decomposition + // heavily inspired by Eigen's implementation https://eigen.tuxfamily.org/dox/LLT_8h_source.html + // TODO use unsafe { *matrix.get_unchecked((j, j)) } let n = x.nrows(); let mut temp = x.clone_owned(); - for k in 0..n { - let lkk = self.chol[(k, k)]; // TODO unsafe { *matrix.get_unchecked((j, j)) } - let xk = temp[k]; - let r = (lkk * lkk + sigma * xk * xk).sqrt(); - let c = r / lkk; - let s = xk / lkk; - self.chol[(k, k)] = r; + let mut beta = crate::one::(); + for j in 0..n { + let ljj = N::real(self.chol[(j, j)]); + let dj = ljj * ljj; + let wj = temp[j]; + let swj2 = sigma * N::modulus_squared(wj); + let gamma = dj * beta + swj2; + let nljj = (dj + swj2 / beta).sqrt(); + self.chol[(j, j)] = N::from_real(nljj); + beta += swj2 / dj; // Update the terms of L - if k < n { - for k2 in (k + 1)..n { - self.chol[(k2, k)] = (self.chol[(k2, k)] + sigma * s * temp[k2]) / c; - temp[k2] = c * temp[k2] - s * self.chol[(k2, k)]; + if j < n { + for k in (j + 1)..n { + temp[k] -= (wj / N::from_real(ljj)) * self.chol[(k, j)]; + if gamma != crate::zero::() { + self.chol[(k, j)] = N::from_real(nljj / ljj) * self.chol[(k, j)] + + (N::from_real(nljj * sigma / gamma) * N::conjugate(wj)) * temp[k]; + } } } } diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs index bef5de95..b04ed402 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -79,10 +79,8 @@ macro_rules! gen_tests( } fn cholesky_rank_one_update(_n: usize) -> bool { - use nalgebra::dimension::U3; - use nalgebra::Vector3; - let mut m = RandomSDP::new(U3, || random::<$scalar>().0).unwrap(); - let x = Vector3::<$scalar>::new_random().map(|e| e.0); + let mut m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap(); + let x = Vector4::<$scalar>::new_random().map(|e| e.0); // TODO this is dirty but $scalar appears to not be a scalar type in this file let zero = random::<$scalar>().0 * 0.; @@ -96,11 +94,7 @@ macro_rules! gen_tests( let m_chol_updated = chol.l() * chol.l().adjoint(); // updates m manually - m.ger(sigma_scalar, &x, &x, one); // m += sigma * x * x.adjoint() - - println!("sigma : {}", sigma); - println!("m updated : {}", m); - println!("chol : {}", m_chol_updated); + m.gerc(sigma_scalar, &x, &x, one); // m += sigma * x * x.adjoint() relative_eq!(m, m_chol_updated, epsilon = 1.0e-7) } From 3123da55293b294b15b91964c0c49067bd721a15 Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sat, 2 Nov 2019 19:04:07 +0100 Subject: [PATCH 08/23] code cleaned --- src/linalg/cholesky.rs | 38 +++++++++++++++++++++----------------- tests/linalg/cholesky.rs | 4 ++-- 2 files changed, 23 insertions(+), 19 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index c4049504..d0a9918c 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -147,7 +147,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^*`. + /// 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) where S2: Storage, @@ -156,27 +156,31 @@ where { // for a description of the operation, see https://en.wikipedia.org/wiki/Cholesky_decomposition#Updating_the_decomposition // heavily inspired by Eigen's implementation https://eigen.tuxfamily.org/dox/LLT_8h_source.html - // TODO use unsafe { *matrix.get_unchecked((j, j)) } let n = x.nrows(); - let mut temp = x.clone_owned(); + let mut x = x.clone_owned(); let mut beta = crate::one::(); for j in 0..n { - let ljj = N::real(self.chol[(j, j)]); - let dj = ljj * ljj; - let wj = temp[j]; - let swj2 = sigma * N::modulus_squared(wj); - let gamma = dj * beta + swj2; - let nljj = (dj + swj2 / beta).sqrt(); - self.chol[(j, j)] = N::from_real(nljj); - beta += swj2 / dj; + 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; // Update the terms of L if j < n { - for k in (j + 1)..n { - temp[k] -= (wj / N::from_real(ljj)) * self.chol[(k, j)]; - if gamma != crate::zero::() { - self.chol[(k, j)] = N::from_real(nljj / ljj) * self.chol[(k, j)] - + (N::from_real(nljj * sigma / gamma) * N::conjugate(wj)) * temp[k]; - } + 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), + ); } } } diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs index b04ed402..ea8402a3 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -82,13 +82,13 @@ macro_rules! gen_tests( let mut m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap(); let x = Vector4::<$scalar>::new_random().map(|e| e.0); - // TODO this is dirty but $scalar appears to not be a scalar type in this file + // this is dirty but $scalar is not a scalar type (its a Rand) in this file let zero = random::<$scalar>().0 * 0.; let one = zero + 1.; let sigma = random::(); // needs to be a real let sigma_scalar = zero + sigma; - // updates cholesky decomposition and reconstructs m + // updates cholesky decomposition and reconstructs m updated let mut chol = m.clone().cholesky().unwrap(); chol.rank_one_update(&x, sigma); let m_chol_updated = chol.l() * chol.l().adjoint(); From d4d3f69429a1ae74fb11a6b2469bae24de40c7ab Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sat, 2 Nov 2019 19:05:39 +0100 Subject: [PATCH 09/23] code cleaned --- src/linalg/cholesky.rs | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index d0a9918c..0f453975 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -154,8 +154,7 @@ where DefaultAllocator: Allocator, ShapeConstraint: SameNumberOfRows, { - // for a description of the operation, see https://en.wikipedia.org/wiki/Cholesky_decomposition#Updating_the_decomposition - // heavily inspired by Eigen's implementation https://eigen.tuxfamily.org/dox/LLT_8h_source.html + // heavily inspired by Eigen's `llt_rank_update_lower` implementation https://eigen.tuxfamily.org/dox/LLT_8h_source.html let n = x.nrows(); let mut x = x.clone_owned(); let mut beta = crate::one::(); From 19fce3addca300033fb28a1590f66acdf6baed73 Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sat, 2 Nov 2019 19:28:46 +0100 Subject: [PATCH 10/23] removed useless if --- src/linalg/cholesky.rs | 27 +++++++++++++-------------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 0f453975..25e632c5 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -159,6 +159,7 @@ where 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) }; @@ -167,20 +168,18 @@ where 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; - // Update the terms of L - if j < n { - 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), - ); - } + // 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), + ); } } } From cf2286cb36984024d4f23d6ca96a5cd92f0fe23a Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sun, 3 Nov 2019 09:36:03 +0100 Subject: [PATCH 11/23] added assertion --- src/linalg/cholesky.rs | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 25e632c5..0362f809 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -156,6 +156,11 @@ where { // 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 { From 03730f594b987e01a5809dcf65873f01571b38e0 Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sun, 3 Nov 2019 13:20:56 +0100 Subject: [PATCH 12/23] finally got the correct type for insert column --- src/linalg/cholesky.rs | 31 ++++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 0362f809..7623e247 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -6,8 +6,9 @@ use alga::general::ComplexField; use crate::allocator::Allocator; use crate::base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, SquareMatrix}; use crate::constraint::{SameNumberOfRows, ShapeConstraint}; -use crate::dimension::{Dim, DimSub, Dynamic, U1}; +use crate::dimension::{Dim, DimAdd, DimSum, DimSub, Dynamic, U1}; use crate::storage::{Storage, StorageMut}; +use crate::base::allocator::Reallocator; /// The Cholesky decomposition of a symmetric-definite-positive matrix. #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] @@ -188,6 +189,34 @@ 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( + self, + j: usize, + c: &Matrix, + ) -> Cholesky> + where + D: DimAdd, + DefaultAllocator: Reallocator> + Reallocator, DimSum, DimSum>, + S2: Storage, + ShapeConstraint: SameNumberOfRows>, + { + let n = c.nrows(); + 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 chol= self.chol.insert_column(j, N::zero()).insert_row(j, N::zero()); + + // TODO see https://en.wikipedia.org/wiki/Cholesky_decomposition#Updating_the_decomposition + unimplemented!(); + Cholesky { chol } + } } impl, S: Storage> SquareMatrix From ebbfc84e96d7ea2f1fe27e5fbb620650b21c0647 Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sun, 3 Nov 2019 13:26:18 +0100 Subject: [PATCH 13/23] added template for remove_column --- src/linalg/cholesky.rs | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 7623e247..f63ab826 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -6,7 +6,7 @@ use alga::general::ComplexField; use crate::allocator::Allocator; use crate::base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, SquareMatrix}; use crate::constraint::{SameNumberOfRows, ShapeConstraint}; -use crate::dimension::{Dim, DimAdd, DimSum, DimSub, Dynamic, U1}; +use crate::dimension::{Dim, DimAdd, DimSum, DimDiff, DimSub, Dynamic, U1}; use crate::storage::{Storage, StorageMut}; use crate::base::allocator::Reallocator; @@ -217,6 +217,26 @@ where unimplemented!(); Cholesky { chol } } + + /// 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. + pub fn remove_column( + self, + j: usize, + ) -> Cholesky> + where + D: DimSub, + DefaultAllocator: Reallocator> + Reallocator, DimDiff, DimDiff>, + { + let n = self.chol.nrows(); + 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 chol= self.chol.remove_column(j).remove_row(j); + + // TODO see https://en.wikipedia.org/wiki/Cholesky_decomposition#Updating_the_decomposition + unimplemented!(); + Cholesky { chol } + } } impl, S: Storage> SquareMatrix From fd5cef6609b704c2e267d0b070e69fc948db95d6 Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sun, 3 Nov 2019 14:33:35 +0100 Subject: [PATCH 14/23] remove column is now working --- src/linalg/cholesky.rs | 57 +++++++++++++++++++++++++++++++++++++--- tests/linalg/cholesky.rs | 19 ++++++++++++++ 2 files changed, 72 insertions(+), 4 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index f63ab826..e6a072de 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -211,7 +211,7 @@ 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 ? - let chol= self.chol.insert_column(j, N::zero()).insert_row(j, N::zero()); + let chol= self.chol.clone().insert_column(j, N::zero()).insert_row(j, N::zero()); // TODO see https://en.wikipedia.org/wiki/Cholesky_decomposition#Updating_the_decomposition unimplemented!(); @@ -229,12 +229,16 @@ where DefaultAllocator: Reallocator> + Reallocator, DimDiff, DimDiff>, { 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 chol= self.chol.remove_column(j).remove_row(j); + let mut chol= self.chol.clone().remove_column(j).remove_row(j); + + // updates the 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())); - // TODO see https://en.wikipedia.org/wiki/Cholesky_decomposition#Updating_the_decomposition - unimplemented!(); Cholesky { chol } } } @@ -251,3 +255,48 @@ where 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()`. +fn rank_one_update_helper(chol : &mut Matrix, x: &Matrix, sigma: N::RealField) + where + N: ComplexField, D: DimSub, R2: Dim, + S: StorageMut, + S2: Storage, + DefaultAllocator: Allocator + 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, + 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 { *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::() { + // 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), + ); + } + } +} \ No newline at end of file diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs index ea8402a3..aa411564 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -98,6 +98,25 @@ macro_rules! gen_tests( relative_eq!(m, m_chol_updated, epsilon = 1.0e-7) } + + fn cholesky_remove_column(n: usize) -> bool { + let n = n.max(1).min(5); + let j = random::() % n; + let m = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap(); + + // remove column from cholesky decomposition and rebuild m + let chol = m.clone().cholesky().unwrap().remove_column(j); + let m_chol_updated = chol.l() * chol.l().adjoint(); + + // remove column from m + let m_updated = m.remove_column(j).remove_row(j); + + println!("n={} j={}", n, j); + println!("chol:{}", m_chol_updated); + println!("m up:{}", m_updated); + + relative_eq!(m_updated, m_chol_updated, epsilon = 1.0e-7) + } } } } From 6d1a00218faaf269080f7be306df27c64db82cae Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sun, 3 Nov 2019 15:17:20 +0100 Subject: [PATCH 15/23] found uneeded storagemut --- src/linalg/cholesky.rs | 26 +++++++++++++++++++++----- src/linalg/solve.rs | 12 ++++++------ tests/linalg/cholesky.rs | 24 ++++++++++++++++++++---- 3 files changed, 47 insertions(+), 15 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index e6a072de..45d232f2 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -195,7 +195,7 @@ where pub fn insert_column( self, j: usize, - c: &Matrix, + col: &Matrix, ) -> Cholesky> where D: DimAdd, @@ -203,7 +203,7 @@ where S2: Storage, ShapeConstraint: SameNumberOfRows>, { - let n = c.nrows(); + let n = col.nrows(); assert_eq!( n, self.chol.nrows() + 1, @@ -211,10 +211,26 @@ 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 ? - let chol= self.chol.clone().insert_column(j, N::zero()).insert_row(j, N::zero()); + // TODO check for adjoint problems + let mut chol= self.chol.clone().insert_column(j, N::zero()).insert_row(j, N::zero()); + + // update the top center element S12 + let top_left_corner = chol.slice_range(..j-1, ..j-1); + let colj = col.rows_range(..j-1); // clone_owned needed to get storage mut for b in solve + let new_colj = top_left_corner.ad_solve_lower_triangular(&colj).unwrap(); + chol.slice_range_mut(..j-1, j).copy_from(&new_colj); + + // update the center element S22 + let rowj = chol.slice_range(j, ..j-1); + let center_element = N::sqrt(col[j] + rowj.dot(&rowj.adjoint())); // TODO is there a better way to multiply a vector by its adjoint ? norm_squared ? + chol[(j,j)] = center_element; + + // update the right center element S23 + //chol.slice_range_mut(j+1.., j).copy_from(&new_rowj); + + // update the bottom right corner // TODO see https://en.wikipedia.org/wiki/Cholesky_decomposition#Updating_the_decomposition - unimplemented!(); Cholesky { chol } } @@ -234,7 +250,7 @@ where // 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 corner + // 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())); diff --git a/src/linalg/solve.rs b/src/linalg/solve.rs index f10b1d00..a6b9196f 100644 --- a/src/linalg/solve.rs +++ b/src/linalg/solve.rs @@ -15,7 +15,7 @@ impl> SquareMatrix { b: &Matrix, ) -> Option> where - S2: StorageMut, + S2: Storage, DefaultAllocator: Allocator, ShapeConstraint: SameNumberOfRows, { @@ -35,7 +35,7 @@ impl> SquareMatrix { b: &Matrix, ) -> Option> where - S2: StorageMut, + S2: Storage, DefaultAllocator: Allocator, ShapeConstraint: SameNumberOfRows, { @@ -191,7 +191,7 @@ impl> SquareMatrix { b: &Matrix, ) -> Option> where - S2: StorageMut, + S2: Storage, DefaultAllocator: Allocator, ShapeConstraint: SameNumberOfRows, { @@ -211,7 +211,7 @@ impl> SquareMatrix { b: &Matrix, ) -> Option> where - S2: StorageMut, + S2: Storage, DefaultAllocator: Allocator, ShapeConstraint: SameNumberOfRows, { @@ -273,7 +273,7 @@ impl> SquareMatrix { b: &Matrix, ) -> Option> where - S2: StorageMut, + S2: Storage, DefaultAllocator: Allocator, ShapeConstraint: SameNumberOfRows, { @@ -293,7 +293,7 @@ impl> SquareMatrix { b: &Matrix, ) -> Option> where - S2: StorageMut, + S2: Storage, DefaultAllocator: Allocator, ShapeConstraint: SameNumberOfRows, { diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs index aa411564..e3e5fdc7 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -99,6 +99,26 @@ macro_rules! gen_tests( relative_eq!(m, m_chol_updated, epsilon = 1.0e-7) } + fn cholesky_insert_column(n: usize) -> bool { + let n = n.max(1).min(5); + let j = random::() % n; + let m_updated = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap(); + + // build m and col from m_updated + let col = m_updated.column(j); + let m = m_updated.clone().remove_column(j).remove_row(j); + + // remove column from cholesky decomposition and rebuild m + let chol = m.clone().cholesky().unwrap().insert_column(j, &col); + 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) + } + fn cholesky_remove_column(n: usize) -> bool { let n = n.max(1).min(5); let j = random::() % n; @@ -111,10 +131,6 @@ macro_rules! gen_tests( // remove column from m let m_updated = m.remove_column(j).remove_row(j); - println!("n={} j={}", n, j); - println!("chol:{}", m_chol_updated); - println!("m up:{}", m_updated); - relative_eq!(m_updated, m_chol_updated, epsilon = 1.0e-7) } } From b4c078ca76387261a9980e8e5f834b7a7bcb5c39 Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sun, 3 Nov 2019 15:43:49 +0100 Subject: [PATCH 16/23] insert does not compile yet --- src/linalg/cholesky.rs | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 45d232f2..755f5610 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -6,7 +6,7 @@ use alga::general::ComplexField; use crate::allocator::Allocator; use crate::base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, SquareMatrix}; use crate::constraint::{SameNumberOfRows, ShapeConstraint}; -use crate::dimension::{Dim, DimAdd, DimSum, DimDiff, DimSub, Dynamic, U1}; +use crate::dimension::{Dim, DimName, DimAdd, DimSum, DimDiff, DimSub, Dynamic, U1}; use crate::storage::{Storage, StorageMut}; use crate::base::allocator::Reallocator; @@ -214,21 +214,25 @@ where // TODO check for adjoint problems let mut chol= self.chol.clone().insert_column(j, N::zero()).insert_row(j, N::zero()); - // update the top center element S12 + // update the jth row let top_left_corner = chol.slice_range(..j-1, ..j-1); - let colj = col.rows_range(..j-1); // clone_owned needed to get storage mut for b in solve - let new_colj = top_left_corner.ad_solve_lower_triangular(&colj).unwrap(); - chol.slice_range_mut(..j-1, j).copy_from(&new_colj); + let colj_minus = col.rows_range(..j-1); + 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-1).copy_from(&rowj); - // update the center element S22 - let rowj = chol.slice_range(j, ..j-1); + // update the center element let center_element = N::sqrt(col[j] + rowj.dot(&rowj.adjoint())); // TODO is there a better way to multiply a vector by its adjoint ? norm_squared ? chol[(j,j)] = center_element; - // update the right center element S23 - //chol.slice_range_mut(j+1.., j).copy_from(&new_rowj); + // 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; + 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..); + 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 } @@ -276,7 +280,9 @@ where /// 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) where - N: ComplexField, D: DimSub, R2: Dim, + N: ComplexField, + D: DimSub, + R2: Dim, S: StorageMut, S2: Storage, DefaultAllocator: Allocator + Allocator, From e6467c2a3e5ac794ae709db2d0e1577b011ecae7 Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sun, 3 Nov 2019 18:02:27 +0100 Subject: [PATCH 17/23] insert does compile --- src/linalg/cholesky.rs | 35 +++++++++++++++++------------------ 1 file changed, 17 insertions(+), 18 deletions(-) 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(); From 46687b7cde62079f0e218214d5f2e1dbf1828c83 Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sun, 3 Nov 2019 18:48:04 +0100 Subject: [PATCH 18/23] tests pass, needs cleanup --- src/linalg/cholesky.rs | 13 ++++++++----- tests/linalg/cholesky.rs | 2 +- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 0485e9e5..bbd233eb 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -216,18 +216,21 @@ where 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-1, ..j-1); - let colj_minus = col.rows_range(..j-1); + 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-1).copy_from(&rowj); + chol.slice_range_mut(j, ..j).copy_from(&rowj); + + // TODO + //println!("dotc:{} norm2:{}", rowj.dotc(&rowj), rowj.norm_squared()); // update the center element - let center_element = N::sqrt(col[j] + rowj.dot(&rowj.adjoint())); // TODO is there a better way to multiply a vector by its adjoint ? norm_squared ? + let center_element = N::sqrt(col[j] - rowj.dotc(&rowj) ); 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-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); diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs index e3e5fdc7..e94cd80f 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -100,7 +100,7 @@ macro_rules! gen_tests( } fn cholesky_insert_column(n: usize) -> bool { - let n = n.max(1).min(5); + let n = n.max(1).min(50); let j = random::() % n; let m_updated = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap(); From 72834186d88821268d725b34c5276dce17059bbd Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sun, 3 Nov 2019 20:00:15 +0100 Subject: [PATCH 19/23] needs faster matrix initialization --- src/linalg/cholesky.rs | 80 +++++++++++++----------------------------- 1 file changed, 24 insertions(+), 56 deletions(-) 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, From 7302defe56eb513308195ddb7f7460bce570de78 Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sun, 3 Nov 2019 21:24:44 +0100 Subject: [PATCH 20/23] finished cleaning --- src/linalg/cholesky.rs | 127 +++++++++++++++++++++------------------ tests/linalg/cholesky.rs | 8 +-- 2 files changed, 69 insertions(+), 66 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 12674e4c..262d1225 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -8,7 +8,6 @@ use crate::base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, SquareMatrix, Vec use crate::constraint::{SameNumberOfRows, ShapeConstraint}; use crate::dimension::{Dim, DimAdd, DimSum, DimDiff, DimSub, Dynamic, U1}; use crate::storage::{Storage, StorageMut}; -use crate::base::allocator::Reallocator; /// The Cholesky decomposition of a symmetric-definite-positive matrix. #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] @@ -156,7 +155,7 @@ where DefaultAllocator: Allocator, ShapeConstraint: SameNumberOfRows, { - 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. @@ -170,7 +169,7 @@ where D: DimAdd, R2: Dim, S2: Storage, - DefaultAllocator: Reallocator> + Reallocator, DimSum, DimSum>, + DefaultAllocator: Allocator, DimSum>, ShapeConstraint: SameNumberOfRows>, { // 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!(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()); + // loads the data into a new matrix with an additional jth row/column + 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 let top_left_corner = self.chol.slice_range(..j, ..j); @@ -200,7 +203,7 @@ where // update the bottom right corner 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 } } @@ -208,27 +211,80 @@ where /// 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. pub fn remove_column( - self, + &self, j: usize, ) -> Cholesky> where D: DimSub, - DefaultAllocator: Reallocator> + Reallocator, DimDiff, DimDiff>, + DefaultAllocator: Allocator, DimDiff> { 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); + // loads the data into a new matrix except for the jth row/column + 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 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())); + Self::xx_rank_one_update(&mut bottom_right_corner, &old_colj, N::real(N::one())); 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(chol : &mut Matrix, x: &Vector, sigma: N::RealField) + where + //N: ComplexField, + Dm: Dim, + Rx: Dim, + Sm: StorageMut, + 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(); + 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::(); + 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::() { + // 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, S: Storage> SquareMatrix @@ -243,52 +299,3 @@ where 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(chol : &mut Matrix, x: &Vector, sigma: N::RealField) - where - N: ComplexField, - D: Dim, - Rx: Dim, - S: StorageMut, - 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(); - 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::(); - 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::() { - // 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), - ); - } - } -} \ No newline at end of file diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs index e94cd80f..5f10339d 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -100,7 +100,7 @@ macro_rules! gen_tests( } fn cholesky_insert_column(n: usize) -> bool { - let n = n.max(1).min(50); + let n = n.max(1).min(10); let j = random::() % n; 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 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) } fn cholesky_remove_column(n: usize) -> bool { - let n = n.max(1).min(5); + let n = n.max(1).min(10); let j = random::() % n; let m = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap(); From 246f72ebc08011ee8806c517205b0b87586fac2b Mon Sep 17 00:00:00 2001 From: sebcrozet Date: Sun, 17 Nov 2019 13:10:50 +0100 Subject: [PATCH 21/23] Fix Cholesky for no-std platforms. --- src/linalg/cholesky.rs | 43 ++++++++++++++++++++++------------------ tests/linalg/cholesky.rs | 2 +- 2 files changed, 25 insertions(+), 20 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 262d1225..111f5680 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -1,6 +1,7 @@ #[cfg(feature = "serde-serialize")] use serde::{Deserialize, Serialize}; +use num::One; use alga::general::ComplexField; use crate::allocator::Allocator; @@ -155,23 +156,24 @@ where DefaultAllocator: Allocator, ShapeConstraint: SameNumberOfRows, { - Self::xx_rank_one_update(&mut self.chol, x, sigma) + Self::xx_rank_one_update(&mut self.chol, &mut x.clone_owned(), sigma) } /// 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, + &self, j: usize, - col: &Vector, + col: Vector, ) -> Cholesky> where D: DimAdd, R2: Dim, S2: Storage, - DefaultAllocator: Allocator, DimSum>, + DefaultAllocator: Allocator, DimSum> + Allocator, ShapeConstraint: SameNumberOfRows>, { + let mut col = col.into_owned(); // 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."); @@ -186,24 +188,26 @@ where // update the jth row 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!"); + + let col_j = col[j]; + let (mut new_rowj_adjoint, mut new_colj) = col.rows_range_pair_mut(..j, j + 1..); + assert!(top_left_corner.solve_lower_triangular_mut(&mut new_rowj_adjoint), "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] - N::from_real(new_rowj_adjoint.norm_squared())); - chol[(j,j)] = center_element; + 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 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 ); + 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..); - Self::xx_rank_one_update(&mut bottom_right_corner, &new_colj, -N::real(N::one())); + Self::xx_rank_one_update(&mut bottom_right_corner, &mut new_colj, -N::RealField::one()); Cholesky { chol } } @@ -216,7 +220,7 @@ where ) -> Cholesky> where D: DimSub, - DefaultAllocator: Allocator, DimDiff> + DefaultAllocator: Allocator, DimDiff> + Allocator { let n = self.chol.nrows(); assert!(n > 0, "The matrix needs at least one column."); @@ -231,25 +235,25 @@ where // updates the bottom right corner let mut bottom_right_corner = chol.slice_range_mut(j.., j..); - let old_colj = self.chol.slice_range(j+1.., j); - Self::xx_rank_one_update(&mut bottom_right_corner, &old_colj, N::real(N::one())); + let mut workspace = self.chol.column(j).clone_owned(); + let mut old_colj = workspace.rows_range_mut(j+1..); + Self::xx_rank_one_update(&mut bottom_right_corner, &mut old_colj, N::RealField::one()); 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()`. + /// performs a rank one update such that we end up with the decomposition of `M + sigma * x*x.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(chol : &mut Matrix, x: &Vector, sigma: N::RealField) + fn xx_rank_one_update(chol : &mut Matrix, x: &mut Vector, sigma: N::RealField) where //N: ComplexField, Dm: Dim, Rx: Dim, Sm: StorageMut, - Sx: Storage, - DefaultAllocator: Allocator, + Sx: StorageMut, { // heavily inspired by Eigen's `llt_rank_update_lower` implementation https://eigen.tuxfamily.org/dox/LLT_8h_source.html let n = x.nrows(); @@ -258,8 +262,9 @@ where 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 { *chol.get_unchecked((j, j)) }); diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs index 5f10339d..a89802b2 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -109,7 +109,7 @@ macro_rules! gen_tests( let m = m_updated.clone().remove_column(j).remove_row(j); // remove column from cholesky decomposition and rebuild m - 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(); relative_eq!(m_updated, m_chol_updated, epsilon = 1.0e-7) From 04ebd9254f00107df6f082db7cca20e429c52c48 Mon Sep 17 00:00:00 2001 From: sebcrozet Date: Sun, 17 Nov 2019 13:24:00 +0100 Subject: [PATCH 22/23] Add some missing spaces. --- src/linalg/cholesky.rs | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 111f5680..1a310cbf 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -148,7 +148,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()`. + /// 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 @@ -182,9 +182,9 @@ where // loads the data into a new matrix with an additional jth row/column 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..)); + 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 let top_left_corner = self.chol.slice_range(..j, ..j); @@ -203,10 +203,10 @@ where let bottom_left_corner = self.chol.slice_range(j.., ..j); // new_colj = (col_jplus - bottom_left_corner * new_rowj.adjoint()) / center_element; 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); + 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..); + let mut bottom_right_corner = chol.slice_range_mut(j + 1.., j + 1..); Self::xx_rank_one_update(&mut bottom_right_corner, &mut new_colj, -N::RealField::one()); Cholesky { chol } @@ -229,21 +229,21 @@ where // loads the data into a new matrix except for the jth row/column 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..)); + 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 let mut bottom_right_corner = chol.slice_range_mut(j.., j..); let mut workspace = self.chol.column(j).clone_owned(); - let mut old_colj = workspace.rows_range_mut(j+1..); + let mut old_colj = workspace.rows_range_mut(j + 1..); Self::xx_rank_one_update(&mut bottom_right_corner, &mut old_colj, N::RealField::one()); 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 * x*x.adjoint()`. + /// performs a rank one update such that we end up with the decomposition of `M + sigma * (x * x.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 From e848f76aa7e458f0c8cbf7c8d1d869d322b60eee Mon Sep 17 00:00:00 2001 From: Nestor Demeure Date: Sun, 17 Nov 2019 13:40:19 +0100 Subject: [PATCH 23/23] Update cholesky.rs corrected typo ni doc for `xx_rank_one_update` --- src/linalg/cholesky.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 1a310cbf..67baefb1 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -242,10 +242,10 @@ where Cholesky { chol } } - /// Given the Cholesky decomposition of a matrix `M`, a scalar `sigma` and a vector `v`, + /// Given the Cholesky decomposition of a matrix `M`, a scalar `sigma` and a vector `x`, /// performs a rank one update such that we end up with the decomposition of `M + sigma * (x * x.adjoint())`. /// - /// This helper method is calling for by `rank_one_update` but also `insert_column` and `remove_column` + /// This helper method is called 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(chol : &mut Matrix, x: &mut Vector, sigma: N::RealField) where