nalgebra/tests/linalg/cholesky.rs

175 lines
7.6 KiB
Rust
Raw Permalink Normal View History

2021-03-01 00:52:14 +08:00
#![cfg(all(feature = "proptest-support", feature = "debug"))]
#[test]
// #[rustfmt::skip]
fn cholesky_with_substitute() {
// Make a tiny covariance matrix with a small covariance value.
let m = na::Matrix2::new(1.0, f64::NAN, 1.0, 1e-32);
// Show that the cholesky fails for our matrix. We then try again with a substitute.
assert!(na::Cholesky::new(m).is_none());
// ...and show that we get some result this time around.
assert!(na::Cholesky::new_with_substitute(m, 1e-8).is_some());
}
macro_rules! gen_tests(
($module: ident, $scalar: ty) => {
mod $module {
use na::debug::RandomSDP;
2023-01-14 23:22:27 +08:00
use na::dimension::{Const, Dyn};
use na::{DMatrix, DVector, Matrix4x3, Vector4};
use rand::random;
use simba::scalar::ComplexField;
#[allow(unused_imports)]
2019-03-23 21:29:07 +08:00
use crate::core::helper::{RandScalar, RandComplex};
2021-03-01 00:52:14 +08:00
use crate::proptest::*;
use proptest::{prop_assert, proptest};
proptest! {
#[test]
fn cholesky(n in PROPTEST_MATRIX_DIM) {
2023-01-14 23:22:27 +08:00
let m = RandomSDP::new(Dyn(n), || random::<$scalar>().0).unwrap();
let l = m.clone().cholesky().unwrap().unpack();
2021-03-01 00:52:14 +08:00
prop_assert!(relative_eq!(m, &l * l.adjoint(), epsilon = 1.0e-7));
}
2021-03-01 00:52:14 +08:00
#[test]
fn cholesky_static(_n in PROPTEST_MATRIX_DIM) {
let m = RandomSDP::new(Const::<4>, || random::<$scalar>().0).unwrap();
let chol = m.cholesky().unwrap();
let l = chol.unpack();
2021-03-01 00:52:14 +08:00
prop_assert!(relative_eq!(m, &l * l.adjoint(), epsilon = 1.0e-7));
}
2021-03-01 00:52:14 +08:00
#[test]
fn cholesky_solve(n in PROPTEST_MATRIX_DIM, nb in PROPTEST_MATRIX_DIM) {
2023-01-14 23:22:27 +08:00
let m = RandomSDP::new(Dyn(n), || random::<$scalar>().0).unwrap();
let chol = m.clone().cholesky().unwrap();
let b1 = DVector::<$scalar>::new_random(n).map(|e| e.0);
let b2 = DMatrix::<$scalar>::new_random(n, nb).map(|e| e.0);
let sol1 = chol.solve(&b1);
let sol2 = chol.solve(&b2);
2021-03-01 00:52:14 +08:00
prop_assert!(relative_eq!(&m * &sol1, b1, epsilon = 1.0e-7));
prop_assert!(relative_eq!(&m * &sol2, b2, epsilon = 1.0e-7));
}
2021-03-01 00:52:14 +08:00
#[test]
fn cholesky_solve_static(_n in PROPTEST_MATRIX_DIM) {
let m = RandomSDP::new(Const::<4>, || random::<$scalar>().0).unwrap();
let chol = m.clone().cholesky().unwrap();
let b1 = Vector4::<$scalar>::new_random().map(|e| e.0);
let b2 = Matrix4x3::<$scalar>::new_random().map(|e| e.0);
let sol1 = chol.solve(&b1);
let sol2 = chol.solve(&b2);
2021-03-01 00:52:14 +08:00
prop_assert!(relative_eq!(m * sol1, b1, epsilon = 1.0e-7));
prop_assert!(relative_eq!(m * sol2, b2, epsilon = 1.0e-7));
}
2021-03-01 00:52:14 +08:00
#[test]
fn cholesky_inverse(n in PROPTEST_MATRIX_DIM) {
2023-01-14 23:22:27 +08:00
let m = RandomSDP::new(Dyn(n), || random::<$scalar>().0).unwrap();
let m1 = m.clone().cholesky().unwrap().inverse();
let id1 = &m * &m1;
let id2 = &m1 * &m;
2021-03-01 00:52:14 +08:00
prop_assert!(id1.is_identity(1.0e-7) && id2.is_identity(1.0e-7));
}
2021-03-01 00:52:14 +08:00
#[test]
fn cholesky_inverse_static(_n in PROPTEST_MATRIX_DIM) {
let m = RandomSDP::new(Const::<4>, || random::<$scalar>().0).unwrap();
let m1 = m.clone().cholesky().unwrap().inverse();
let id1 = &m * &m1;
let id2 = &m1 * &m;
2021-03-01 00:52:14 +08:00
prop_assert!(id1.is_identity(1.0e-7) && id2.is_identity(1.0e-7));
}
2019-11-02 22:11:14 +08:00
2021-04-08 10:43:27 +08:00
#[test]
fn cholesky_determinant(n in PROPTEST_MATRIX_DIM) {
2023-01-14 23:22:27 +08:00
let m = RandomSDP::new(Dyn(n), || random::<$scalar>().0).unwrap();
2021-04-08 10:43:27 +08:00
let lu_det = m.clone().lu().determinant();
assert_relative_eq!(lu_det.imaginary(), 0., epsilon = 1.0e-7);
2021-04-08 10:43:27 +08:00
let chol_det = m.cholesky().unwrap().determinant();
prop_assert!(relative_eq!(lu_det.real(), chol_det, epsilon = 1.0e-7));
2021-04-08 10:43:27 +08:00
}
#[test]
fn cholesky_determinant_static(_n in PROPTEST_MATRIX_DIM) {
let m = RandomSDP::new(Const::<4>, || random::<$scalar>().0).unwrap();
2021-04-08 10:43:27 +08:00
let lu_det = m.clone().lu().determinant();
assert_relative_eq!(lu_det.imaginary(), 0., epsilon = 1.0e-7);
2021-04-08 10:43:27 +08:00
let chol_det = m.cholesky().unwrap().determinant();
prop_assert!(relative_eq!(lu_det.real(), chol_det, epsilon = 1.0e-7));
2021-04-08 10:43:27 +08:00
}
2021-03-01 00:52:14 +08:00
#[test]
fn cholesky_rank_one_update(_n in PROPTEST_MATRIX_DIM) {
let mut m = RandomSDP::new(Const::<4>, || random::<$scalar>().0).unwrap();
2019-11-03 01:27:01 +08:00
let x = Vector4::<$scalar>::new_random().map(|e| e.0);
2019-11-02 23:45:30 +08:00
2019-11-03 02:04:07 +08:00
// this is dirty but $scalar is not a scalar type (its a Rand) in this file
2019-11-02 23:45:30 +08:00
let zero = random::<$scalar>().0 * 0.;
let one = zero + 1.;
let sigma = random::<f64>(); // needs to be a real
let sigma_scalar = zero + sigma;
2019-11-02 22:56:59 +08:00
2019-11-03 02:04:07 +08:00
// updates cholesky decomposition and reconstructs m updated
2019-11-02 22:56:59 +08:00
let mut chol = m.clone().cholesky().unwrap();
chol.rank_one_update(&x, sigma);
let m_chol_updated = chol.l() * chol.l().adjoint();
2019-11-02 22:11:14 +08:00
// updates m manually
2019-11-03 01:27:01 +08:00
m.gerc(sigma_scalar, &x, &x, one); // m += sigma * x * x.adjoint()
2019-11-02 22:11:14 +08:00
2021-03-01 00:52:14 +08:00
prop_assert!(relative_eq!(m, m_chol_updated, epsilon = 1.0e-7));
2019-11-02 22:11:14 +08:00
}
2019-11-03 21:33:35 +08:00
2021-03-01 00:52:14 +08:00
#[test]
fn cholesky_insert_column(n in PROPTEST_MATRIX_DIM) {
2019-11-04 04:24:44 +08:00
let n = n.max(1).min(10);
2019-11-03 22:17:20 +08:00
let j = random::<usize>() % n;
2023-01-14 23:22:27 +08:00
let m_updated = RandomSDP::new(Dyn(n), || random::<$scalar>().0).unwrap();
2019-11-03 22:17:20 +08:00
// 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
2019-11-17 20:10:50 +08:00
let chol = m.clone().cholesky().unwrap().insert_column(j, col);
2019-11-03 22:17:20 +08:00
let m_chol_updated = chol.l() * chol.l().adjoint();
2021-03-01 00:52:14 +08:00
prop_assert!(relative_eq!(m_updated, m_chol_updated, epsilon = 1.0e-7));
2019-11-03 22:17:20 +08:00
}
2021-03-01 00:52:14 +08:00
#[test]
fn cholesky_remove_column(n in PROPTEST_MATRIX_DIM) {
2019-11-04 04:24:44 +08:00
let n = n.max(1).min(10);
2019-11-03 21:33:35 +08:00
let j = random::<usize>() % n;
2023-01-14 23:22:27 +08:00
let m = RandomSDP::new(Dyn(n), || random::<$scalar>().0).unwrap();
2019-11-03 21:33:35 +08:00
// 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);
2021-03-01 00:52:14 +08:00
prop_assert!(relative_eq!(m_updated, m_chol_updated, epsilon = 1.0e-7));
2019-11-03 21:33:35 +08:00
}
}
}
}
);
gen_tests!(complex, RandComplex<f64>);
gen_tests!(f64, RandScalar<f64>);