Merge pull request #1236 from vasilNnikolov/fix_bug_1218

Fix bug 1218
This commit is contained in:
Sébastien Crozet 2023-04-30 15:27:44 +02:00 committed by GitHub
commit 1a271ac541
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 30 additions and 6 deletions

View File

@ -4,6 +4,10 @@ documented here.
This project adheres to [Semantic Versioning](https://semver.org/).
## Unreleased
### Fixed
- Fixed severe catastrophic cancellation issue in variance calculation.
## [0.32.2] (07 March 2023)

View File

@ -111,6 +111,7 @@ serde_json = "1.0"
rand_xorshift = "0.3"
rand_isaac = "0.3"
criterion = { version = "0.4", features = ["html_reports"] }
nalgebra = { path = ".", features = ["debug", "compare", "rand", "macros"]}
# For matrix comparison macro
matrixcompare = "0.3.0"

View File

@ -335,12 +335,12 @@ impl<T: Scalar, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> {
if self.is_empty() {
T::zero()
} else {
let val = self.iter().cloned().fold((T::zero(), T::zero()), |a, b| {
(a.0 + b.clone() * b.clone(), a.1 + b)
});
let denom = T::one() / crate::convert::<_, T>(self.len() as f64);
let vd = val.1 * denom.clone();
val.0 * denom - vd.clone() * vd
let n_elements: T = crate::convert(self.len() as f64);
let mean = self.mean();
self.iter().cloned().fold(T::zero(), |acc, x| {
acc + (x.clone() - mean.clone()) * (x.clone() - mean.clone())
}) / n_elements
}
}

View File

@ -11,6 +11,7 @@ mod reshape;
#[cfg(feature = "rkyv-serialize-no-std")]
mod rkyv;
mod serde;
mod variance;
#[cfg(feature = "compare")]
mod matrixcompare;

18
tests/core/variance.rs Normal file
View File

@ -0,0 +1,18 @@
use nalgebra::DVector;
#[test]
fn test_variance_catastrophic_cancellation() {
let long_repeating_vector = DVector::repeat(10_000, 100000000.0);
assert_eq!(long_repeating_vector.variance(), 0.0);
let short_vec = DVector::from_vec(vec![1., 2., 3.]);
assert_eq!(short_vec.variance(), 2.0 / 3.0);
let short_vec =
DVector::<f64>::from_vec(vec![1.0e8 + 4.0, 1.0e8 + 7.0, 1.0e8 + 13.0, 1.0e8 + 16.0]);
assert_eq!(short_vec.variance(), 22.5);
let short_vec =
DVector::<f64>::from_vec(vec![1.0e9 + 4.0, 1.0e9 + 7.0, 1.0e9 + 13.0, 1.0e9 + 16.0]);
assert_eq!(short_vec.variance(), 22.5);
}