-
Notifications
You must be signed in to change notification settings - Fork 103
Open
Description
The exact-equality symmetry validation on MultivariateNormal is not appropriate for floating-point calculations, and is not required for the Cholesky decomposition from nalgebra. Here's a quick script to show an example:
use rand::Rng;
let mut exactly_sym = 0;
let mut almost_sym = 0;
let mut not_sym = 0;
let mut cholesky = 0;
let total = 10_000;
for _ in 0..total {
// Matrix that is nicely behaved, even this doesn't always pass an exact equality check!
let sym =
nalgebra::debug::RandomSDP::new(nalgebra::U2, || thread_rng().r#gen::<f64>()).unwrap();
assert_relative_eq!(
sym.upper_triangle(),
sym.lower_triangle().transpose(),
epsilon = 1e-8
);
// Doesn't really matter how we construct this, we could even just use the identity matrix!
let mut f = nalgebra::Matrix2::<f64>::identity();
f[(0, 1)] = 1.0;
// This should also be symmetric, by mathematical definition
let also_sym = f * sym * f.transpose();
if also_sym.upper_triangle() != also_sym.lower_triangle().transpose() {
if also_sym.upper_triangle().relative_eq(
&also_sym.lower_triangle().transpose(),
1e-10,
1e-10,
) {
almost_sym += 1;
} else {
not_sym += 1;
}
} else {
exactly_sym += 1;
}
match nalgebra::Cholesky::new(also_sym) {
Some(_) => {
cholesky += 1;
}
None => {}
}
}
println!("Almost symmetric: {almost_sym}/{total}");
println!("Exactly symmetric: {exactly_sym}/{total}");
println!("Cholesky ok: {cholesky}/{total}");
assert_eq!(not_sym, 0);Here are the results:
Almost symmetric: 2059/10000 # Notice that these two
Exactly symmetric: 7941/10000 # values sum to 10000
Cholesky ok: 10000/10000It seems to me that it would be more appropriate for MultivariateNormal to either check for approximate equality, or simply pass through nalgebra's assumption that the covariance is (sufficiently) symmetric. Alternatively, a constructor that allows users to pass in a Cholesky struct directly would also suffice.
Metadata
Metadata
Assignees
Labels
No labels