diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 286e3f2..e55d6bb 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -1,9 +1,9 @@ //! # Cholesky Decomposition //! -//! every positive definite matrix \\(A \in R^{n \times n}\\) can be factored as +//! every positive definite matrix \\(A \in R^{n \times n}\\) can be factored as //! //! \\[A = R^TR\\] -//! +//! //! where \\(R\\) is upper triangular matrix with positive diagonal elements //! //! Example: @@ -12,8 +12,8 @@ //! use crate::smartcore::linalg::cholesky::*; //! //! let A = DenseMatrix::from_2d_array(&[ -//! &[25., 15., -5.], -//! &[15., 18., 0.], +//! &[25., 15., -5.], +//! &[15., 18., 0.], //! &[-5., 0., 11.] //! ]); //! @@ -41,14 +41,14 @@ use crate::math::num::RealNumber; /// Results of Cholesky decomposition. pub struct Cholesky> { R: M, - t: PhantomData + t: PhantomData, } impl> Cholesky { pub(crate) fn new(R: M) -> Cholesky { Cholesky { R: R, - t: PhantomData + t: PhantomData, } } @@ -65,10 +65,10 @@ impl> Cholesky { } } R - } - + } + /// Get upper triangular matrix. - pub fn U(&self) -> M { + pub fn U(&self) -> M { let (n, _) = self.R.shape(); let mut R = M::zeros(n, n); @@ -80,20 +80,20 @@ impl> Cholesky { } } R - } + } /// Solves Ax = b - pub(crate) fn solve(&self, mut b: M) -> Result { - + pub(crate) fn solve(&self, mut b: M) -> Result { let (bn, m) = b.shape(); let (rn, _) = self.R.shape(); if bn != rn { - return Err(Failed::because(FailedError::SolutionFailed, &format!( - "Can't solve Ax = b for x. Number of rows in b != number of rows in R." - ))); + return Err(Failed::because( + FailedError::SolutionFailed, + &format!("Can't solve Ax = b for x. Number of rows in b != number of rows in R."), + )); } - + for k in 0..bn { for j in 0..m { for i in 0..k { @@ -102,7 +102,7 @@ impl> Cholesky { b.div_element_mut(k, j, self.R.get(k, k)); } } - + for k in (0..bn).rev() { for j in 0..m { for i in k + 1..bn { @@ -128,11 +128,12 @@ pub trait CholeskyDecomposableMatrix: BaseMatrix { let (m, n) = self.shape(); if m != n { - return Err(Failed::because(FailedError::DecompositionFailed, &format!( - "Can't do Cholesky decomposition on a non-square matrix" - ))); + return Err(Failed::because( + FailedError::DecompositionFailed, + &format!("Can't do Cholesky decomposition on a non-square matrix"), + )); } - + for j in 0..n { let mut d = T::zero(); for k in 0..j { @@ -147,9 +148,10 @@ pub trait CholeskyDecomposableMatrix: BaseMatrix { d = self.get(j, j) - d; if d < T::zero() { - return Err(Failed::because(FailedError::DecompositionFailed, &format!( - "The matrix is not positive definite." - ))); + return Err(Failed::because( + FailedError::DecompositionFailed, + &format!("The matrix is not positive definite."), + )); } self.set(j, j, d.sqrt()); @@ -172,35 +174,33 @@ mod tests { #[test] fn cholesky_decompose() { let a = DenseMatrix::from_2d_array(&[&[25., 15., -5.], &[15., 18., 0.], &[-5., 0., 11.]]); - let l = DenseMatrix::from_2d_array(&[ - &[5.0, 0.0, 0.0], - &[3.0, 3.0, 0.0], - &[-1.0, 1.0, 3.0], - ]); - let u = DenseMatrix::from_2d_array(&[ - &[5.0, 3.0, -1.0], - &[0.0, 3.0, 1.0], - &[0.0, 0.0, 3.0], - ]); + let l = + DenseMatrix::from_2d_array(&[&[5.0, 0.0, 0.0], &[3.0, 3.0, 0.0], &[-1.0, 1.0, 3.0]]); + let u = + DenseMatrix::from_2d_array(&[&[5.0, 3.0, -1.0], &[0.0, 3.0, 1.0], &[0.0, 0.0, 3.0]]); let cholesky = a.cholesky().unwrap(); - - assert!(cholesky.L().abs().approximate_eq(&l.abs(), 1e-4)); - assert!(cholesky.U().abs().approximate_eq(&u.abs(), 1e-4)); - assert!(cholesky.L().matmul(&cholesky.U()).abs().approximate_eq(&a.abs(), 1e-4)); + + assert!(cholesky.L().abs().approximate_eq(&l.abs(), 1e-4)); + assert!(cholesky.U().abs().approximate_eq(&u.abs(), 1e-4)); + assert!(cholesky + .L() + .matmul(&cholesky.U()) + .abs() + .approximate_eq(&a.abs(), 1e-4)); } #[test] fn cholesky_solve_mut() { let a = DenseMatrix::from_2d_array(&[&[25., 15., -5.], &[15., 18., 0.], &[-5., 0., 11.]]); let b = DenseMatrix::from_2d_array(&[&[40., 51., 28.]]); - let expected = DenseMatrix::from_2d_array(&[ - &[1.0, 2.0, 3.0] - ]); - + let expected = DenseMatrix::from_2d_array(&[&[1.0, 2.0, 3.0]]); + let cholesky = a.cholesky().unwrap(); - assert!(cholesky.solve(b.transpose()).unwrap().transpose().approximate_eq(&expected, 1e-4)); - + assert!(cholesky + .solve(b.transpose()) + .unwrap() + .transpose() + .approximate_eq(&expected, 1e-4)); } - } diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs index 61749b0..fb12909 100644 --- a/src/linalg/mod.rs +++ b/src/linalg/mod.rs @@ -33,6 +33,7 @@ //! let u: DenseMatrix = svd.U; //! ``` +pub mod cholesky; /// The matrix is represented in terms of its eigenvalues and eigenvectors. pub mod evd; /// Factors a matrix as the product of a lower triangular matrix and an upper triangular matrix. @@ -49,18 +50,17 @@ pub mod ndarray_bindings; pub mod qr; /// Singular value decomposition. pub mod svd; -pub mod cholesky; use std::fmt::{Debug, Display}; use std::marker::PhantomData; use std::ops::Range; use crate::math::num::RealNumber; +use cholesky::CholeskyDecomposableMatrix; use evd::EVDDecomposableMatrix; use lu::LUDecomposableMatrix; use qr::QRDecomposableMatrix; use svd::SVDDecomposableMatrix; -use cholesky::CholeskyDecomposableMatrix; /// Column or row vector pub trait BaseVector: Clone + Debug { diff --git a/src/linalg/naive/dense_matrix.rs b/src/linalg/naive/dense_matrix.rs index b5ecd90..d3d6353 100644 --- a/src/linalg/naive/dense_matrix.rs +++ b/src/linalg/naive/dense_matrix.rs @@ -8,11 +8,11 @@ use serde::de::{Deserializer, MapAccess, SeqAccess, Visitor}; use serde::ser::{SerializeStruct, Serializer}; use serde::{Deserialize, Serialize}; +use crate::linalg::cholesky::CholeskyDecomposableMatrix; use crate::linalg::evd::EVDDecomposableMatrix; use crate::linalg::lu::LUDecomposableMatrix; use crate::linalg::qr::QRDecomposableMatrix; use crate::linalg::svd::SVDDecomposableMatrix; -use crate::linalg::cholesky::CholeskyDecomposableMatrix; use crate::linalg::Matrix; pub use crate::linalg::{BaseMatrix, BaseVector}; use crate::math::num::RealNumber; diff --git a/src/linalg/nalgebra_bindings.rs b/src/linalg/nalgebra_bindings.rs index a400a67..e0b885b 100644 --- a/src/linalg/nalgebra_bindings.rs +++ b/src/linalg/nalgebra_bindings.rs @@ -42,11 +42,11 @@ use std::ops::{AddAssign, DivAssign, MulAssign, Range, SubAssign}; use nalgebra::{DMatrix, Dynamic, Matrix, MatrixMN, RowDVector, Scalar, VecStorage, U1}; +use crate::linalg::cholesky::CholeskyDecomposableMatrix; use crate::linalg::evd::EVDDecomposableMatrix; use crate::linalg::lu::LUDecomposableMatrix; use crate::linalg::qr::QRDecomposableMatrix; use crate::linalg::svd::SVDDecomposableMatrix; -use crate::linalg::cholesky::CholeskyDecomposableMatrix; use crate::linalg::Matrix as SmartCoreMatrix; use crate::linalg::{BaseMatrix, BaseVector}; use crate::math::num::RealNumber; diff --git a/src/linalg/ndarray_bindings.rs b/src/linalg/ndarray_bindings.rs index 76749a7..00c9745 100644 --- a/src/linalg/ndarray_bindings.rs +++ b/src/linalg/ndarray_bindings.rs @@ -49,11 +49,11 @@ use std::ops::SubAssign; use ndarray::ScalarOperand; use ndarray::{s, stack, Array, ArrayBase, Axis, Ix1, Ix2, OwnedRepr}; +use crate::linalg::cholesky::CholeskyDecomposableMatrix; use crate::linalg::evd::EVDDecomposableMatrix; use crate::linalg::lu::LUDecomposableMatrix; use crate::linalg::qr::QRDecomposableMatrix; use crate::linalg::svd::SVDDecomposableMatrix; -use crate::linalg::cholesky::CholeskyDecomposableMatrix; use crate::linalg::Matrix; use crate::linalg::{BaseMatrix, BaseVector}; use crate::math::num::RealNumber;