diff --git a/src/pyrecest/filters/candidate_mixture.py b/src/pyrecest/filters/candidate_mixture.py index a5f2f5d6a..2ed8021a9 100644 --- a/src/pyrecest/filters/candidate_mixture.py +++ b/src/pyrecest/filters/candidate_mixture.py @@ -99,6 +99,17 @@ def __post_init__(self) -> None: covariances = covariances.copy() if not np.all(np.isfinite(covariances)): raise ValueError("covariances must contain only finite values") + covariance_transposes = np.swapaxes(covariances, -1, -2) + covariance_scale = np.maximum( + np.maximum(np.abs(covariances), np.abs(covariance_transposes)), 1.0 + ) + with np.errstate(over="ignore", under="ignore", invalid="ignore"): + relative_asymmetry = np.abs( + covariances / covariance_scale + - covariance_transposes / covariance_scale + ) + if np.any(relative_asymmetry > 1e-12): + raise ValueError("covariances must be symmetric") covariances = np.stack([_symmetrize(value) for value in covariances]) try: cholesky = np.linalg.cholesky(covariances) diff --git a/tests/filters/test_candidate_mixture_covariance_symmetry.py b/tests/filters/test_candidate_mixture_covariance_symmetry.py new file mode 100644 index 000000000..d3f217590 --- /dev/null +++ b/tests/filters/test_candidate_mixture_covariance_symmetry.py @@ -0,0 +1,48 @@ +import unittest + +import numpy as np +from pyrecest.filters.candidate_mixture import GaussianMixtureMeasurementFactor + + +class GaussianMixtureMeasurementCovarianceSymmetryTest(unittest.TestCase): + def test_rejects_nonsymmetric_shared_covariance(self): + nonsymmetric = np.array([[2.0, 1.0], [0.0, 2.0]]) + + with self.assertRaisesRegex(ValueError, "covariances must be symmetric"): + GaussianMixtureMeasurementFactor( + means=np.zeros((2, 2)), + covariances=nonsymmetric, + ) + + def test_rejects_nonsymmetric_component_covariance(self): + covariances = np.array( + [ + np.eye(2), + [[2.0, 1.0], [0.0, 2.0]], + ] + ) + + with self.assertRaisesRegex(ValueError, "covariances must be symmetric"): + GaussianMixtureMeasurementFactor( + means=np.zeros((2, 2)), + covariances=covariances, + ) + + def test_tolerates_roundoff_scale_asymmetry(self): + covariance = np.array([[2.0, 1.0 + 1e-13], [1.0, 2.0]]) + + factor = GaussianMixtureMeasurementFactor( + means=np.zeros((1, 2)), + covariances=covariance, + ) + + np.testing.assert_allclose( + factor.covariances[0], + 0.5 * (covariance + covariance.T), + rtol=0.0, + atol=0.0, + ) + + +if __name__ == "__main__": + unittest.main()