From 50435c7de24f043371b703476b99780afc5f2dc1 Mon Sep 17 00:00:00 2001 From: Andrea De Vittori Date: Fri, 21 Nov 2025 12:53:33 -0500 Subject: [PATCH 1/3] Fixed cov(self) to handle small covariances and large means in a vectorized fashion --- pyest/gm/gm.py | 28 ++++++++++++++++++++++------ 1 file changed, 22 insertions(+), 6 deletions(-) diff --git a/pyest/gm/gm.py b/pyest/gm/gm.py index ad97d36..f980aa3 100755 --- a/pyest/gm/gm.py +++ b/pyest/gm/gm.py @@ -871,18 +871,33 @@ def mean(self): return np.sum(self.w[:, np.newaxis] * self.m, axis=0) def cov(self): - """ return covariance of the distribution + """Compute the covariance matrix of the distribution. Returns ------- ndarray - (nx,nx) full covariance matrix of distribution + (nx, nx) full covariance matrix of the distribution """ - wsum = np.sum(self.w) + # Compute the mean of the distribution mean = self.mean() - return np.sum([ - w/wsum*(P + np.outer(m, m)) for w, m, P in self - ], axis=0) - np.outer(mean, mean) + + # Extract weights for each sample (assuming self.w exists) + w = self.w + + # Normalise weights to sum to 1 + w = w / np.sum(w) + + # Extract individual covariances (self.P should be shape (n_samples, nx, nx)) + covs = self.P + + # Compute the difference between each sample and the mean + diffs = self.m - mean # shape (n_samples, nx) + + # Compute the outer product of diffs for each sample + outer_diffs = diffs[:, :, None] * diffs[:, None, :] # shape (n_samples, nx, nx) + + # Weighted sum of covariances and outer products + return np.sum(w[:, None, None] * (covs + outer_diffs), axis=0) def cdf(self, x, allow_singular=False): """ evaluate GM CDF at points x""" @@ -1024,3 +1039,4 @@ def __init__( self.state_idxs = np.array(state_idxs) self.spectral_direction_only = spectral_direction_only self.variance_preserving = variance_preserving + From 8d83b7f06d5a812e32768660f90401604b1495ed Mon Sep 17 00:00:00 2001 From: Andrea De Vittori Date: Fri, 21 Nov 2025 12:59:37 -0500 Subject: [PATCH 2/3] Add unit test to ensure covariance is not altered after splitting --- tests/test_gm.py | 52 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/tests/test_gm.py b/tests/test_gm.py index 32bdfbf..3d721ab 100644 --- a/tests/test_gm.py +++ b/tests/test_gm.py @@ -8,6 +8,7 @@ from pyest.sensors import ConvexPolyhedralFieldOfView from pyest.utils import fail import pytest +from pyest.metrics import max_covariance_ratio def test_gm_equal_weights(): @@ -555,6 +556,57 @@ def test_init_mismatch_fail(): fail(gm.GaussianMixture, ValueError, w, m, P) +def test_cov(): + # Initial covariance matrix P0 (6×6) + P0 = np.array([ + [2.25e-06, 0, 0, 0, 0, 0], + [0, 2.25e-06, 0, 0, 0, 0], + [0, 0, 2.25e-06, 0, 0, 0], + [0, 0, 0, 1.77777778e-14, 0, 0], + [0, 0, 0, 0, 1.77777778e-14, 0], + [0, 0, 0, 0, 0, 1.77777778e-14] + ]) + + # Initial large mean vector m0 + m0 = np.array([ + -3.95911227e+05, -1.80422801e+05, -8.03243134e+04, + 5.48813474e-01, -1.02480997e+00, -5.93612009e-01 + ]) + + # Single-component GMM (weight 1) + weights = np.array([1]) + p0 = gm.GaussianMixture(weights, m0, P0) + + # Options for Gaussian splitting + split_opts = gm.GaussSplitOptions( + L=3, # number of components created per split + lam=1e-3, # scaling parameter + recurse_depth=3, # maximum recursion depth + min_weight=-np.inf + ) + + # No weight threshold + split_tol = -np.inf + + # Arguments passed to recursive_split: splitting method + tolerance + recursive_split_args = (gm.split.id_variance, split_tol) + + # Apply recursive splitting to the initial GMM + p0 = gm.split.recursive_split(p0, split_opts, *recursive_split_args) + + # Cholesky decomposition of the resulting GMM covariance + Schol_gmm = np.atleast_2d(np.linalg.cholesky(p0.cov())) + + # Cholesky decomposition of the original covariance + Schol_P0 = np.atleast_2d(np.linalg.cholesky(P0)) + + # Compute the maximum covariance ratio between the two matrices + max_cov_ratio = abs(max_covariance_ratio(Schol_gmm, Schol_P0)) + + # Test: the covariance should not increase by more than 0.01% + assert max_cov_ratio <= 1.0001, f"max_cov_ratio is too large: {max_cov_ratio}" if __name__ == '__main__': pytest.main([__file__]) + + From 0dc1ad27601b30a6aec703597accdaf09a506291 Mon Sep 17 00:00:00 2001 From: Keith LeGrand <26415081+keithalegrand@users.noreply.github.com> Date: Sat, 22 Nov 2025 21:45:37 -0500 Subject: [PATCH 3/3] clean up comments, add copy to prevent modification by reference --- pyest/gm/gm.py | 24 +++++++++++------------- tests/test_gm.py | 4 ++-- 2 files changed, 13 insertions(+), 15 deletions(-) diff --git a/pyest/gm/gm.py b/pyest/gm/gm.py index f980aa3..dfff8b2 100755 --- a/pyest/gm/gm.py +++ b/pyest/gm/gm.py @@ -6,6 +6,7 @@ import numpy as np import pandas as pd from numpy import pi, sqrt +from copy import copy from numpy.random import rand, randn from scipy.stats._multivariate import _LOG_2PI, _PSD, _squeeze_output from scipy.stats import Covariance, _covariance @@ -880,22 +881,19 @@ def cov(self): """ # Compute the mean of the distribution mean = self.mean() - - # Extract weights for each sample (assuming self.w exists) - w = self.w - - # Normalise weights to sum to 1 + + w = copy(self.w) w = w / np.sum(w) - - # Extract individual covariances (self.P should be shape (n_samples, nx, nx)) + + # Extract individual covariances covs = self.P - - # Compute the difference between each sample and the mean - diffs = self.m - mean # shape (n_samples, nx) - + + # Compute the difference between each mixand mean and the distribution mean + diffs = self.m - mean # shape (nC, nx) + # Compute the outer product of diffs for each sample - outer_diffs = diffs[:, :, None] * diffs[:, None, :] # shape (n_samples, nx, nx) - + outer_diffs = diffs[:, :, None] * diffs[:, None, :] # shape (nC, nx, nx) + # Weighted sum of covariances and outer products return np.sum(w[:, None, None] * (covs + outer_diffs), axis=0) diff --git a/tests/test_gm.py b/tests/test_gm.py index 3d721ab..981c9d3 100644 --- a/tests/test_gm.py +++ b/tests/test_gm.py @@ -572,7 +572,7 @@ def test_cov(): -3.95911227e+05, -1.80422801e+05, -8.03243134e+04, 5.48813474e-01, -1.02480997e+00, -5.93612009e-01 ]) - + # Single-component GMM (weight 1) weights = np.array([1]) p0 = gm.GaussianMixture(weights, m0, P0) @@ -584,7 +584,7 @@ def test_cov(): recurse_depth=3, # maximum recursion depth min_weight=-np.inf ) - + # No weight threshold split_tol = -np.inf