Numerical Stability in Gaussian-Process Computation • gaussianprocesses

Why numerical stability matters

Gaussian-process regression repeatedly solves systems involving

C=K+\Sigma_\varepsilon.

When inputs are duplicated, almost duplicated, or connected by a kernel with a very long length scale, C can become poorly conditioned.

The package treats numerical stabilization as part of the implementation, rather than silently replacing the statistical model.

Cholesky rather than explicit inversion

The core computation is

C = R^\top R.

To compute

\alpha=C^{-1}y,

the implementation solves

R^\top z=y, \qquad R\alpha=z.

No explicit matrix inverse is required.

This is both more stable and usually more efficient.

Jitter is not observation noise

If Cholesky factorization fails, the numerical layer can try

C_\delta = C + \delta I.

The scalar \delta is jitter. It is a numerical device.

Observation noise,

\sigma_\varepsilon^2,

is part of the statistical model.

They are stored separately.

Consider duplicated observations with zero statistical noise:

model <- fit_gp(
  x = c(0, 0),
  y = c(1, 1),
  kernel = rbf_kernel(),
  noise_variance = 0,
  fallback_jitter = 1e-10
)

model$noise_variance
#> [1] 0
model$jitter
#> [1] 1e-10
model$cholesky_attempts
#> [1] 2

The statistical model still has zero observation noise even though the factorization required numerical stabilization.

Jitter is relative to the covariance scale

The jitter is not a fixed number. It is a multiple of the covariance scale, the mean prior variance

s = \frac{1}{n}\operatorname{tr}(C), \qquad \delta = \epsilon\, s,

where the relative jitter \epsilon climbs the ladder 0, then fallback_jitter = 10^{-10}, multiplied by jitter_multiplier = 10 after each failure.

Because \delta is proportional to C, results do not depend on the units of y. Multiplying y by c and every variance parameter by c^2 multiplies C, the jitter, and every posterior variance by c^2, and the posterior mean by c:

# 40 nearly coincident inputs without noise make C numerically singular, so
# every fit needs jitter.
x <- seq(0, 1, length.out = 40)

fit_in_units <- function(c) {
  fit_gp(
    x,
    c * sin(3 * x),
    kernel = rbf_kernel(variance = c^2, length_scale = 0.5),
    noise_variance = 0
  )
}

unit <- fit_in_units(1)
small <- fit_in_units(1e-4)

# The absolute jitter scales with C; the relative jitter is the same.
c(unit = unit$jitter, small = small$jitter)
#>  unit small 
#> 1e-10 1e-18
c(
  unit = gp_numerical_diagnostics(unit)$relative_jitter,
  small = gp_numerical_diagnostics(small)$relative_jitter
)
#>  unit small 
#> 1e-10 1e-10

# Relative difference of the posterior means after undoing the scaling.
x_new <- c(0.25, 0.5, 1.5)
unit_mean <- predict_gp(unit, x_new)$mean
small_mean <- predict_gp(small, x_new)$mean / 1e-4
max(abs(small_mean - unit_mean)) / max(abs(unit_mean))
#> [1] 3.909373e-10

Earlier versions used the same ladder as absolute values. Then the jitter was negligible for a covariance of scale 10^8 and dominant for one of scale 10^{-8}. The study in inst/benchmarks/jitter-scaling.R rescales problems that need jitter by c = 10^{-4} and c = 10^{4}:

Policy	Log marginal likelihood	Posterior mean and variance	Function draws
Absolute jitter	off by 0.4 to 1.7	off by up to 37%	off by up to 110%
Relative jitter	within 7 \times 10^{-8}	within 2 \times 10^{-8}	within 1.2 \times 10^{-6}

The table gives relative differences after undoing the scaling. The remaining differences under relative jitter come from rounding: scaling by c^2 = 10^{\pm 8} is not exact in binary arithmetic, and a nearly singular covariance amplifies the rounding.

The trace term of the likelihood gradient

The gradient of the log marginal likelihood with respect to a log parameter \eta is

\frac12\left(\alpha^\top D \alpha - \operatorname{tr}(C^{-1}D)\right), \qquad D = \frac{\partial C}{\partial \eta}.

The trace needs C^{-1} in some form. log_marginal_likelihood_gradient() computes it once from the Cholesky factor with chol2inv(), and each parameter then costs O(n^2). The alternative uses two triangular solves per parameter and never forms C^{-1}, but each parameter then costs O(n^3).

inst/benchmarks/gradient-numerics.R compares the two. A periodic kernel on a regular grid gives a circulant C whose eigenvalues are sums of Bessel functions, so the exact gradient is known to near machine precision even when C is ill-conditioned:

Condition number of C	`chol2inv()`	Triangular solves
10^4	1.5 \times 10^{-13}	1.5 \times 10^{-13}
10^8	6.0 \times 10^{-11}	6.0 \times 10^{-11}
10^{12}	4.8 \times 10^{-7}	4.8 \times 10^{-7}

These are the largest relative errors of the gradient. The two methods are equally accurate, and both lose accuracy in proportion to the conditioning of C, as any method in double precision must. The solves were 4 to 18 times slower, and the gap grows with the number of parameters, so the package keeps chol2inv().

Near-singular scenario

The benchmark suite contains a reproducible near-singular design:

simulation <- gp_simulation_scenario(
  "near_singular",
  n = 30,
  seed = 1729
)

near_singular_model <- fit_gp(
  simulation$x,
  simulation$observed,
  kernel = simulation$kernel,
  noise_variance = simulation$noise_variance
)

gp_numerical_diagnostics(
  near_singular_model
)
#> $observation_noise_variance
#>  [1] 1e-12 1e-12 1e-12 1e-12 1e-12 1e-12 1e-12 1e-12 1e-12 1e-12 1e-12 1e-12
#> [13] 1e-12 1e-12 1e-12 1e-12 1e-12 1e-12 1e-12 1e-12 1e-12 1e-12 1e-12 1e-12
#> [25] 1e-12 1e-12 1e-12 1e-12 1e-12 1e-12
#> 
#> $numerical_jitter
#> [1] 0
#> 
#> $relative_jitter
#> [1] 0
#> 
#> $covariance_scale
#> [1] 1
#> 
#> $cholesky_attempts
#> [1] 1
#> 
#> $used_numerical_jitter
#> [1] FALSE
#> 
#> $reciprocal_condition_number
#> [1] 1.915273e-14
#> 
#> $condition_number
#> [1] 5.221188e+13
#> 
#> $condition_tolerance
#> [1] 1.490116e-08
#> 
#> $potentially_ill_conditioned
#> [1] TRUE
#> 
#> $min_cholesky_diagonal
#> [1] 1.412981e-06
#> 
#> $max_cholesky_diagonal
#> [1] 1
#> 
#> $observation_types
#>    type observations covariance_scale reciprocal_condition_number
#> 1 value           30                1                1.915273e-14

The diagnostic output separates:

observation-noise variance;
numerical jitter, absolute and relative to the covariance scale;
Cholesky attempts;
reciprocal condition number;
condition-number flag.

Posterior variances

Analytically, posterior variances cannot be negative. Numerically, expressions such as

K_{**}-K_{*f}C^{-1}K_{f*}

can produce tiny negative diagonal values due to floating-point cancellation.

The package distinguishes two cases.

For a value of floating-point scale,

-\varepsilon < 0,

the marginal variance is clipped to zero.

For a materially negative value, the computation stops with a numerical error.

This prevents an instability from being silently converted into an apparently valid uncertainty estimate.

Benchmarking stability separately from speed

suite <- run_gp_benchmark_suite(
  scenarios = "near_singular",
  n = 20,
  n_inducing = 6,
  repeats = 1,
  seed = 42
)

suite$stability
#>        scenario exact_numerical_jitter exact_reciprocal_condition_number
#> 1 near_singular                      0                      4.291552e-14
#>   exact_cholesky_attempts sparse_kuu_jitter sparse_a_jitter
#> 1                       1                 0               0
#>   sparse_fitc_floor_count
#> 1                       6
suite$performance
#>        scenario model fit_elapsed_seconds predict_elapsed_seconds
#> 1 near_singular exact               0.001                   0.001
#> 2 near_singular  FITC               0.004                   0.001
#>   retained_model_bytes n_training n_inducing
#> 1                10680         20         20
#> 2                14136         20          6
suite$accuracy
#>        scenario model latent_truth_rmse latent_truth_mae
#> 1 near_singular exact      7.663259e-06     5.758081e-06
#> 2 near_singular  FITC      4.145902e-03     3.239425e-03
#>   posterior_mean_rmse_vs_exact posterior_mean_max_abs_error_vs_exact
#> 1                  0.000000000                           0.000000000
#> 2                  0.004145557                           0.009472646
#>   latent_variance_mae_vs_exact latent_variance_max_abs_error_vs_exact
#> 1                 0.000000e+00                           0.000000e+00
#> 2                 1.428755e-05                           3.964489e-05

Those are different tables intentionally.

A method can be faster but less accurate. A model can fit the latent truth well while still requiring numerical stabilization. A small reciprocal condition number does not by itself prove statistical misspecification.

Keeping these concepts separate is one of the central design rules of the package.