diff --git a/CHANGELOG.md b/CHANGELOG.md index 77300de..4a1a593 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,7 @@ - bugfix DimensionalityEstimator dimensionality initialization - implement 'fixed' gaussian proces type to allow more inducing points than datapoints - implement `copy()` method for `Predictor` class + - fix: uncertainty computation of sparse GP # v1.4.3 diff --git a/mellon/conditional.py b/mellon/conditional.py index 990f330..74ff2bd 100644 --- a/mellon/conditional.py +++ b/mellon/conditional.py @@ -5,7 +5,7 @@ from jax.numpy import diag as diagonal from jax.numpy.linalg import cholesky from jax.scipy.linalg import solve_triangular -from .util import ensure_2d, stabilize, DEFAULT_JITTER, add_variance +from .util import ensure_2d, stabilize, DEFAULT_JITTER, add_variance, add_projected_variance from .base_predictor import Predictor, ExpPredictor, PredictorTime from .decomposition import DEFAULT_SIGMA @@ -278,9 +278,9 @@ def __init__( LLB = stabilize(LLB, jitter) else: logger.debug("Assuming y is not the mean of the GP.") - y_cov_factor = _sigma_to_y_cov_factor(sigma, y_cov_factor, xu.shape[0]) + y_cov_factor = _sigma_to_y_cov_factor(sigma, y_cov_factor, x.shape[0]) sigma = None - LLB = add_variance(LLB, y_cov_factor, jitter=jitter) + LLB = add_projected_variance(LLB, A, y_cov_factor, jitter=jitter) L_B = cholesky(LLB) r = y - mu diff --git a/mellon/util.py b/mellon/util.py index b9b621d..f6ec0a4 100644 --- a/mellon/util.py +++ b/mellon/util.py @@ -315,6 +315,33 @@ def add_variance(K, M=None, jitter=DEFAULT_JITTER): return K +def add_projected_variance(K, A, y_cov_factor, jitter=DEFAULT_JITTER): + """ + Adds the projected observation noise covariance to K and stabilizes it. + + Parameters + ---------- + K : array_like, shape (n_landmarks, n_landmarks) + The initial covariance matrix. + A : array_like, shape (n_landmarks, n_obs) + The projection matrix from observations to inducing points. + y_cov_factor : array_like, shape (n_obs, n_obs) + The observation noise covariance matrix. + jitter : float, optional + A small number to stabilize the covariance matrix. Defaults to 1e-6. + + Returns + ------- + stabilized_K : array_like, shape (n_landmarks, n_landmarks) + The stabilized covariance matrix with added projected variance. + """ + noise = A @ y_cov_factor @ A.T + noise_diag = diagonal(noise) + diff = where(noise_diag < jitter, jitter - noise_diag, 0) + K = K + noise + diagonal(diff) + return K + + def mle(nn_distances, d): R""" Nearest Neighbor distribution maximum likelihood estimate for log density