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.
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.
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] 2The statistical model still has zero observation noise even though the factorization required numerical stabilization.
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] 1.281905e-09Earlier 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 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().
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.915504e-14
#>
#> $condition_number
#> [1] 5.220558e+13
#>
#> $condition_tolerance
#> [1] 1.490116e-08
#>
#> $potentially_ill_conditioned
#> [1] TRUE
#>
#> $min_cholesky_diagonal
#> [1] 1.412941e-06
#>
#> $max_cholesky_diagonal
#> [1] 1
#>
#> $observation_types
#> type observations covariance_scale reciprocal_condition_number
#> 1 value 30 1 1.915504e-14The diagnostic output separates:
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.
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.291551e-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.00 0
#> 2 near_singular FITC 0.02 0
#> retained_model_bytes n_training n_inducing
#> 1 10672 20 20
#> 2 14128 20 6
suite$accuracy
#> scenario model latent_truth_rmse latent_truth_mae
#> 1 near_singular exact 7.663762e-06 5.759123e-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.004145559 0.009472644
#> 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-05Those 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.