Optimize the hyperparameters of a sparse Gaussian process

Estimates kernel parameters, optionally the noise, and optionally the inducing inputs of a sparse model (fit_sparse_gp()) by maximizing the VFE bound or the FITC approximate marginal likelihood, with the optimizer of optimize_gp(): L-BFGS-B on unconstrained coordinates, deterministic multiple starts, and analytical gradients.

Usage

optimize_sparse_gp(
  x,
  y,
  kernel,
  noise_variance = 1e-06,
  mean = zero_mean(),
  n_inducing = NULL,
  inducing_points = NULL,
  selection = c("farthest", "quantile", "variance"),
  method = c("vfe", "fitc"),
  optimize_inducing = FALSE,
  optimize_noise = TRUE,
  n_starts = 3L,
  start_spread = 1,
  lower = 1e-08,
  upper = 1e+08,
  control = list(maxit = 500),
  fitc_diagonal_floor = 1e-10,
  initial_jitter = 0,
  fallback_jitter = 1e-10,
  jitter_multiplier = 10,
  max_attempts = 8L,
  symmetry_tolerance = sqrt(.Machine$double.eps),
  gradient = c("analytical", "numerical")
)

Arguments

x: Numeric vector or matrix of training inputs.
y: Numeric response vector.
kernel: Initial kernel specification.
noise_variance: Noise variance: a scalar starting value, or one value per observation forming a known pattern whose scale is estimated, as in optimize_gp(). With optimize_noise = FALSE it is fixed.
mean: Gaussian-process mean specification with fixed coefficients.
n_inducing: Optional number of inducing points selected from x.
inducing_points: Optional explicit inducing locations. If supplied, n_inducing is ignored.
selection: Inducing-point selection method used when explicit points are not supplied: "farthest", "quantile", or "variance" (greedy conditional-variance selection under kernel; see select_inducing_points()).
method: The approximation: "vfe" (the default here) or "fitc".
optimize_inducing: If TRUE, optimize the inducing inputs jointly with the hyperparameters, starting from the selected or supplied points. The kernel must be differentiable in its inputs (kernel_input_gradient()).
optimize_noise, n_starts, start_spread, lower, upper, control, gradient: As in optimize_gp().
fitc_diagonal_floor: Strictly positive numerical floor applied to the FITC diagonal term before inversion.
initial_jitter: Non-negative jitter tried first in the inducing-system Cholesky factorizations, relative to the scale of each factorized matrix (the mean of its diagonal).
fallback_jitter: Positive relative jitter tried after the first failed factorization.
jitter_multiplier: Multiplicative jitter escalation factor.
max_attempts: Maximum Cholesky attempts.
symmetry_tolerance: Relative tolerance for checking that covariance matrices are symmetric.

Value

A fitted gaussianprocesses_sparse_model with optimizer diagnostics in $optimization, which has the fields of optimize_gp().

Details

VFE is the default because its objective is a lower bound on the exact log marginal likelihood: maximizing it fits the hyperparameters of the exact model as well as the inducing points allow, and pays for poor coverage through the trace term. FITC's objective is not a bound; its diagonal correction can absorb noise, and maximizing it tends to underestimate the noise variance (Bauer, van der Wilk, and Rasmussen, 2016). The research note on sparse approximations compares the two.

The gradient is analytical for both approximations, with the inducing points fixed: with $C = Q_{ff} + \Lambda$ and $M = C^{-1} r r^\top C^{-1} - C^{-1}$, $$\frac{\partial}{\partial \theta_j} \log N(r \mid 0, C) = \frac12 \mathrm{tr}\left(M \left(\frac{\partial Q_{ff}}{\partial \theta_j} + \frac{\partial \Lambda}{\partial \theta_j}\right)\right),$$ where $\partial Q_{ff}$ follows from the derivatives of $K_{uf}$ and $K_{uu}$, plus the derivative of the VFE trace term. It is computed in $O(n m^2)$ time without forming an $n \times n$ matrix, processing the observations in blocks.

Optimizing the inducing inputs

With optimize_inducing = TRUE the $m \times D$ inducing inputs $Z$ are optimized as well, on the real line, as parameters named inducing.<j>.<d>. Their gradient follows from the same terms: $Z_{jd}$ changes only row $j$ of $K_{uf}$ and row and column $j$ of $K_{uu}$, $$\frac{\partial \mathcal F}{\partial Z_{jd}} = \sum_i G^{uf}_{ji} \frac{\partial k(z_j, x_i)}{\partial z_{jd}} + 2 \sum_l G^{uu}_{jl} \frac{\partial k(z_j, z_l)}{\partial z_{jd}},$$ with the kernel input derivatives of kernel_input_gradient() and the sensitivities $G^{uf}$ and $G^{uu}$ of the objective to $K_{uf}$ and $K_{uu}$, at the same $O(n m^2)$ cost.

Optimizing $Z$ is recommended with VFE: the bound can only move towards the exact marginal likelihood, so better inducing inputs make the approximation better. With FITC it is prone to overfitting: the inducing inputs can move or collapse so that FITC's diagonal correction absorbs the noise, and the noise variance is then underestimated (Bauer, van der Wilk, and Rasmussen, 2016). Inducing points that become nearly identical under the kernel are reported in $inducing_collapse with a warning of class gaussianprocesses_inducing_collapse_warning; $K_{uu}$ is factorized with the package's jitter policy, so the fit stays defined. Greedy variance selection (selection = "variance") gives a good starting point.

Stability

Experimental: this interface may change in a minor release, and every change is listed in NEWS. See gaussianprocesses-package for the policy.

Examples

set.seed(1)
x <- seq(-3, 3, length.out = 300)
y <- sin(2 * x) + rnorm(300, sd = 0.2)

model <- optimize_sparse_gp(x, y, rbf_kernel(), noise_variance = 0.1,
  n_inducing = 20, n_starts = 1)
model
#> Sparse Gaussian-process regression model (VFE)
#>   observations: 300
#>   inducing points: 20
#>   inducing selection: farthest
#>   trace term tr(K - Q) / noise: 7.207e-06
#>   variational lower bound on the log marginal likelihood: 41.51277
#>   optimization: converged (L-BFGS-B, analytical gradient)
model$noise_variance
#> [1] 0.03733125