ndarray_linalg/
cholesky.rs

1//! Cholesky decomposition of Hermitian (or real symmetric) positive definite matrices
2//!
3//! See the [Wikipedia page about Cholesky
4//! decomposition](https://en.wikipedia.org/wiki/Cholesky_decomposition) for
5//! more information.
6//!
7//! # Example
8//!
9//! Using the Cholesky decomposition of `A` for various operations, where `A`
10//! is a Hermitian (or real symmetric) positive definite matrix:
11//!
12//! ```
13//! #[macro_use]
14//! extern crate ndarray;
15//! extern crate ndarray_linalg;
16//!
17//! use ndarray::prelude::*;
18//! use ndarray_linalg::cholesky::*;
19//! # fn main() {
20//!
21//! let a: Array2<f64> = array![
22//!     [  4.,  12., -16.],
23//!     [ 12.,  37., -43.],
24//!     [-16., -43.,  98.]
25//! ];
26//!
27//! // Obtain `L`
28//! let lower = a.cholesky(UPLO::Lower).unwrap();
29//! assert!(lower.abs_diff_eq(&array![
30//!     [ 2., 0., 0.],
31//!     [ 6., 1., 0.],
32//!     [-8., 5., 3.]
33//! ], 1e-9));
34//!
35//! // Find the determinant of `A`
36//! let det = a.detc().unwrap();
37//! assert!((det - 36.).abs() < 1e-9);
38//!
39//! // Solve `A * x = b`
40//! let b = array![4., 13., -11.];
41//! let x = a.solvec(&b).unwrap();
42//! assert!(x.abs_diff_eq(&array![-2., 1., 0.], 1e-9));
43//! # }
44//! ```
45
46use ndarray::*;
47use num_traits::Float;
48
49use crate::convert::*;
50use crate::error::*;
51use crate::layout::*;
52use crate::triangular::IntoTriangular;
53use crate::types::*;
54
55pub use lax::UPLO;
56
57/// Cholesky decomposition of Hermitian (or real symmetric) positive definite matrix
58pub struct CholeskyFactorized<S: Data> {
59    /// `L` from the decomposition `A = L * L^H` or `U` from the decomposition
60    /// `A = U^H * U`.
61    pub factor: ArrayBase<S, Ix2>,
62    /// If this is `UPLO::Lower`, then `self.factor` is `L`. If this is
63    /// `UPLO::Upper`, then `self.factor` is `U`.
64    pub uplo: UPLO,
65}
66
67impl<A, S> CholeskyFactorized<S>
68where
69    A: Scalar + Lapack,
70    S: DataMut<Elem = A>,
71{
72    /// Returns `L` from the Cholesky decomposition `A = L * L^H`.
73    ///
74    /// If `self.uplo == UPLO::Lower`, then no computations need to be
75    /// performed; otherwise, the conjugate transpose of `self.factor` is
76    /// calculated.
77    pub fn into_lower(self) -> ArrayBase<S, Ix2> {
78        match self.uplo {
79            UPLO::Lower => self.factor,
80            UPLO::Upper => self.factor.reversed_axes().mapv_into(|elem| elem.conj()),
81        }
82    }
83
84    /// Returns `U` from the Cholesky decomposition `A = U^H * U`.
85    ///
86    /// If `self.uplo == UPLO::Upper`, then no computations need to be
87    /// performed; otherwise, the conjugate transpose of `self.factor` is
88    /// calculated.
89    pub fn into_upper(self) -> ArrayBase<S, Ix2> {
90        match self.uplo {
91            UPLO::Lower => self.factor.reversed_axes().mapv_into(|elem| elem.conj()),
92            UPLO::Upper => self.factor,
93        }
94    }
95}
96
97impl<A, S> DeterminantC for CholeskyFactorized<S>
98where
99    A: Scalar + Lapack,
100    S: Data<Elem = A>,
101{
102    type Output = <A as Scalar>::Real;
103
104    fn detc(&self) -> Self::Output {
105        Float::exp(self.ln_detc())
106    }
107
108    fn ln_detc(&self) -> Self::Output {
109        self.factor
110            .diag()
111            .iter()
112            .map(|elem| Float::ln(elem.square()))
113            .sum::<Self::Output>()
114    }
115}
116
117impl<A, S> DeterminantCInto for CholeskyFactorized<S>
118where
119    A: Scalar + Lapack,
120    S: Data<Elem = A>,
121{
122    type Output = <A as Scalar>::Real;
123
124    fn detc_into(self) -> Self::Output {
125        self.detc()
126    }
127
128    fn ln_detc_into(self) -> Self::Output {
129        self.ln_detc()
130    }
131}
132
133impl<A, S> InverseC for CholeskyFactorized<S>
134where
135    A: Scalar + Lapack,
136    S: Data<Elem = A>,
137{
138    type Output = Array2<A>;
139
140    fn invc(&self) -> Result<Self::Output> {
141        let f = CholeskyFactorized {
142            factor: replicate(&self.factor),
143            uplo: self.uplo,
144        };
145        f.invc_into()
146    }
147}
148
149impl<A, S> InverseCInto for CholeskyFactorized<S>
150where
151    A: Scalar + Lapack,
152    S: DataMut<Elem = A>,
153{
154    type Output = ArrayBase<S, Ix2>;
155
156    fn invc_into(self) -> Result<Self::Output> {
157        let mut a = self.factor;
158        A::inv_cholesky(a.square_layout()?, self.uplo, a.as_allocated_mut()?)?;
159        triangular_fill_hermitian(&mut a, self.uplo);
160        Ok(a)
161    }
162}
163
164impl<A, S> SolveC<A> for CholeskyFactorized<S>
165where
166    A: Scalar + Lapack,
167    S: Data<Elem = A>,
168{
169    fn solvec_inplace<'a>(&self, b: &'a mut ArrayRef<A, Ix1>) -> Result<&'a mut ArrayRef<A, Ix1>> {
170        A::solve_cholesky(
171            self.factor.square_layout()?,
172            self.uplo,
173            self.factor.as_allocated()?,
174            b.as_slice_mut().unwrap(),
175        )?;
176        Ok(b)
177    }
178}
179
180/// Cholesky decomposition of Hermitian (or real symmetric) positive definite matrix reference
181pub trait Cholesky {
182    type Output;
183
184    /// Computes the Cholesky decomposition of the Hermitian (or real
185    /// symmetric) positive definite matrix.
186    ///
187    /// If the argument is `UPLO::Upper`, then computes the decomposition `A =
188    /// U^H * U` using the upper triangular portion of `A` and returns `U`.
189    /// Otherwise, if the argument is `UPLO::Lower`, computes the decomposition
190    /// `A = L * L^H` using the lower triangular portion of `A` and returns
191    /// `L`.
192    fn cholesky(&self, uplo: UPLO) -> Result<Self::Output>;
193}
194
195/// Cholesky decomposition of Hermitian (or real symmetric) positive definite matrix
196pub trait CholeskyInto {
197    type Output;
198    /// Computes the Cholesky decomposition of the Hermitian (or real
199    /// symmetric) positive definite matrix.
200    ///
201    /// If the argument is `UPLO::Upper`, then computes the decomposition `A =
202    /// U^H * U` using the upper triangular portion of `A` and returns `U`.
203    /// Otherwise, if the argument is `UPLO::Lower`, computes the decomposition
204    /// `A = L * L^H` using the lower triangular portion of `A` and returns
205    /// `L`.
206    fn cholesky_into(self, uplo: UPLO) -> Result<Self::Output>;
207}
208
209/// Cholesky decomposition of Hermitian (or real symmetric) positive definite mutable reference of matrix
210pub trait CholeskyInplace {
211    /// Computes the Cholesky decomposition of the Hermitian (or real
212    /// symmetric) positive definite matrix, writing the result (`L` or `U`
213    /// according to the argument) to `self` and returning it.
214    ///
215    /// If the argument is `UPLO::Upper`, then computes the decomposition `A =
216    /// U^H * U` using the upper triangular portion of `A` and writes `U`.
217    /// Otherwise, if the argument is `UPLO::Lower`, computes the decomposition
218    /// `A = L * L^H` using the lower triangular portion of `A` and writes `L`.
219    fn cholesky_inplace(&mut self, uplo: UPLO) -> Result<&mut Self>;
220}
221
222impl<A> Cholesky for ArrayRef<A, Ix2>
223where
224    A: Scalar + Lapack,
225{
226    type Output = Array2<A>;
227
228    fn cholesky(&self, uplo: UPLO) -> Result<Array2<A>> {
229        let a = replicate(self);
230        a.cholesky_into(uplo)
231    }
232}
233
234impl<A, S> CholeskyInto for ArrayBase<S, Ix2>
235where
236    A: Scalar + Lapack,
237    S: DataMut<Elem = A>,
238{
239    type Output = Self;
240
241    fn cholesky_into(mut self, uplo: UPLO) -> Result<Self> {
242        self.cholesky_inplace(uplo)?;
243        Ok(self)
244    }
245}
246
247impl<A> CholeskyInplace for ArrayRef<A, Ix2>
248where
249    A: Scalar + Lapack,
250{
251    fn cholesky_inplace(&mut self, uplo: UPLO) -> Result<&mut Self> {
252        A::cholesky(self.square_layout()?, uplo, self.as_allocated_mut()?)?;
253        Ok(self.into_triangular(uplo))
254    }
255}
256
257/// Cholesky decomposition of Hermitian (or real symmetric) positive definite matrix reference
258pub trait FactorizeC<S: Data> {
259    /// Computes the Cholesky decomposition of the Hermitian (or real
260    /// symmetric) positive definite matrix.
261    ///
262    /// If the argument is `UPLO::Upper`, then computes the decomposition `A =
263    /// U^H * U` using the upper triangular portion of `A` and returns the
264    /// factorization containing `U`. Otherwise, if the argument is
265    /// `UPLO::Lower`, computes the decomposition `A = L * L^H` using the lower
266    /// triangular portion of `A` and returns the factorization containing `L`.
267    fn factorizec(&self, uplo: UPLO) -> Result<CholeskyFactorized<S>>;
268}
269
270/// Cholesky decomposition of Hermitian (or real symmetric) positive definite matrix
271pub trait FactorizeCInto<S: Data> {
272    /// Computes the Cholesky decomposition of the Hermitian (or real
273    /// symmetric) positive definite matrix.
274    ///
275    /// If the argument is `UPLO::Upper`, then computes the decomposition `A =
276    /// U^H * U` using the upper triangular portion of `A` and returns the
277    /// factorization containing `U`. Otherwise, if the argument is
278    /// `UPLO::Lower`, computes the decomposition `A = L * L^H` using the lower
279    /// triangular portion of `A` and returns the factorization containing `L`.
280    fn factorizec_into(self, uplo: UPLO) -> Result<CholeskyFactorized<S>>;
281}
282
283impl<A, S> FactorizeCInto<S> for ArrayBase<S, Ix2>
284where
285    A: Scalar + Lapack,
286    S: DataMut<Elem = A>,
287{
288    fn factorizec_into(self, uplo: UPLO) -> Result<CholeskyFactorized<S>> {
289        Ok(CholeskyFactorized {
290            factor: self.cholesky_into(uplo)?,
291            uplo,
292        })
293    }
294}
295
296impl<A> FactorizeC<OwnedRepr<A>> for ArrayRef<A, Ix2>
297where
298    A: Scalar + Lapack,
299{
300    fn factorizec(&self, uplo: UPLO) -> Result<CholeskyFactorized<OwnedRepr<A>>> {
301        Ok(CholeskyFactorized {
302            factor: self.cholesky(uplo)?,
303            uplo,
304        })
305    }
306}
307
308/// Solve systems of linear equations with Hermitian (or real symmetric)
309/// positive definite coefficient matrices
310pub trait SolveC<A: Scalar> {
311    /// Solves a system of linear equations `A * x = b` with Hermitian (or real
312    /// symmetric) positive definite matrix `A`, where `A` is `self`, `b` is
313    /// the argument, and `x` is the successful result.
314    fn solvec(&self, b: &ArrayRef<A, Ix1>) -> Result<Array1<A>> {
315        let mut b = replicate(b);
316        self.solvec_inplace(&mut b)?;
317        Ok(b)
318    }
319    /// Solves a system of linear equations `A * x = b` with Hermitian (or real
320    /// symmetric) positive definite matrix `A`, where `A` is `self`, `b` is
321    /// the argument, and `x` is the successful result.
322    fn solvec_into<S: DataMut<Elem = A>>(
323        &self,
324        mut b: ArrayBase<S, Ix1>,
325    ) -> Result<ArrayBase<S, Ix1>> {
326        self.solvec_inplace(&mut b)?;
327        Ok(b)
328    }
329    /// Solves a system of linear equations `A * x = b` with Hermitian (or real
330    /// symmetric) positive definite matrix `A`, where `A` is `self`, `b` is
331    /// the argument, and `x` is the successful result. The value of `x` is
332    /// also assigned to the argument.
333    fn solvec_inplace<'a>(&self, b: &'a mut ArrayRef<A, Ix1>) -> Result<&'a mut ArrayRef<A, Ix1>>;
334}
335
336impl<A> SolveC<A> for ArrayRef<A, Ix2>
337where
338    A: Scalar + Lapack,
339{
340    fn solvec_inplace<'a>(&self, b: &'a mut ArrayRef<A, Ix1>) -> Result<&'a mut ArrayRef<A, Ix1>> {
341        self.factorizec(UPLO::Upper)?.solvec_inplace(b)
342    }
343}
344
345/// Inverse of Hermitian (or real symmetric) positive definite matrix ref
346pub trait InverseC {
347    type Output;
348    /// Computes the inverse of the Hermitian (or real symmetric) positive
349    /// definite matrix.
350    fn invc(&self) -> Result<Self::Output>;
351}
352
353/// Inverse of Hermitian (or real symmetric) positive definite matrix
354pub trait InverseCInto {
355    type Output;
356    /// Computes the inverse of the Hermitian (or real symmetric) positive
357    /// definite matrix.
358    fn invc_into(self) -> Result<Self::Output>;
359}
360
361impl<A> InverseC for ArrayRef<A, Ix2>
362where
363    A: Scalar + Lapack,
364{
365    type Output = Array2<A>;
366
367    fn invc(&self) -> Result<Self::Output> {
368        self.factorizec(UPLO::Upper)?.invc_into()
369    }
370}
371
372impl<A, S> InverseCInto for ArrayBase<S, Ix2>
373where
374    A: Scalar + Lapack,
375    S: DataMut<Elem = A>,
376{
377    type Output = Self;
378
379    fn invc_into(self) -> Result<Self::Output> {
380        self.factorizec_into(UPLO::Upper)?.invc_into()
381    }
382}
383
384/// Determinant of Hermitian (or real symmetric) positive definite matrix ref
385pub trait DeterminantC {
386    type Output;
387
388    /// Computes the determinant of the Hermitian (or real symmetric) positive
389    /// definite matrix.
390    fn detc(&self) -> Self::Output;
391
392    /// Computes the natural log of the determinant of the Hermitian (or real
393    /// symmetric) positive definite matrix.
394    ///
395    /// This method is more robust than `.detc()` to very small or very large
396    /// determinants since it returns the natural logarithm of the determinant
397    /// rather than the determinant itself.
398    fn ln_detc(&self) -> Self::Output;
399}
400
401/// Determinant of Hermitian (or real symmetric) positive definite matrix
402pub trait DeterminantCInto {
403    type Output;
404
405    /// Computes the determinant of the Hermitian (or real symmetric) positive
406    /// definite matrix.
407    fn detc_into(self) -> Self::Output;
408
409    /// Computes the natural log of the determinant of the Hermitian (or real
410    /// symmetric) positive definite matrix.
411    ///
412    /// This method is more robust than `.detc_into()` to very small or very
413    /// large determinants since it returns the natural logarithm of the
414    /// determinant rather than the determinant itself.
415    fn ln_detc_into(self) -> Self::Output;
416}
417
418impl<A> DeterminantC for ArrayRef<A, Ix2>
419where
420    A: Scalar + Lapack,
421{
422    type Output = Result<<A as Scalar>::Real>;
423
424    fn detc(&self) -> Self::Output {
425        Ok(Float::exp(self.ln_detc()?))
426    }
427
428    fn ln_detc(&self) -> Self::Output {
429        Ok(self.factorizec(UPLO::Upper)?.ln_detc())
430    }
431}
432
433impl<A, S> DeterminantCInto for ArrayBase<S, Ix2>
434where
435    A: Scalar + Lapack,
436    S: DataMut<Elem = A>,
437{
438    type Output = Result<<A as Scalar>::Real>;
439
440    fn detc_into(self) -> Self::Output {
441        Ok(Float::exp(self.ln_detc_into()?))
442    }
443
444    fn ln_detc_into(self) -> Self::Output {
445        Ok(self.factorizec_into(UPLO::Upper)?.ln_detc_into())
446    }
447}
ndarray_linalg/cholesky.rs

ndarray_linalg/
cholesky.rs
