Optimize Gaussian-process hyperparameters — optimize

Estimates kernel parameters and, optionally, observation-noise variance by maximizing the exact log marginal likelihood. Optimization is performed in log-parameter space with deterministic multiple starts.

Usage

optimize_gp(
  x,
  y,
  kernel,
  noise_variance = 1e-06,
  mean = zero_mean(),
  optimize_noise = TRUE,
  n_starts = 3L,
  start_spread = 1,
  lower = 1e-08,
  upper = 1e+08,
  control = list(maxit = 500),
  initial_jitter = 0,
  fallback_jitter = 1e-10,
  jitter_multiplier = 10,
  max_attempts = 8L,
  symmetry_tolerance = sqrt(.Machine$double.eps),
  gradient = c("analytical", "numerical"),
  derivative = NULL,
  noise_groups = NULL
)

Arguments

x: Numeric vector or matrix of training inputs.
y: Numeric response vector.
kernel: Initial Gaussian-process kernel specification.
noise_variance: Non-negative observation-noise variance: a scalar, or one value per observation. With optimize_noise = TRUE, a scalar is the starting value of a homoscedastic noise variance, and a vector is a known relative noise pattern $\Sigma_0$ whose overall scale $s$ is estimated, $\Sigma = s \Sigma_0$, starting from $s = 1$. With optimize_noise = FALSE the noise is fixed at the given value or values.
mean: Gaussian-process mean specification. Estimated or marginalized coefficients are not optimizer parameters: for each value of the hyperparameters they are computed in closed form, and the objective is the profile, restricted, or marginal likelihood (mean_coefficients).
optimize_noise: If TRUE, optimize the noise variance, or the scale of a noise pattern, jointly with the kernel parameters.
n_starts: Number of deterministic optimization starts.
start_spread: Log-space displacement used to construct additional deterministic starts around the supplied parameters.
lower: Strictly positive lower bound applied to all optimized parameters on the natural scale.
upper: Strictly positive upper bound applied to all optimized parameters on the natural scale.
control: List passed to stats::optim(). The default allows 500 L-BFGS-B iterations per start: simple kernels converge in a few dozen, but composites with periodic components can need a few hundred. Each iteration costs about one GP fit with analytical gradients.
initial_jitter: Non-negative jitter tried first, relative to the covariance scale (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 per objective evaluation.
symmetry_tolerance: Relative tolerance for checking that the covariance matrix is symmetric.
gradient: How L-BFGS-B obtains the gradient of the objective. "analytical" (the default) uses log_marginal_likelihood_gradient(), which costs one GP fit per optimizer step. "numerical" lets stats::optim() use finite differences, which costs two fits per parameter per step.
derivative: Optional observation types for derivative observations, as in fit_gp().
noise_groups: Optional group of each observation, such as the output of a multi-output model (gp_stack_outputs()): with optimize_noise = TRUE, one noise variance is estimated per group, as noise_variance[g] for the gth of the sorted unique groups. noise_variance then gives their starting values: one value for every group, or one per group in that order. Experimental.

Value

A fitted gaussianprocesses_model with optimized kernel and noise parameters. Optimization diagnostics are stored in $optimization.

Details

Kernel parameters are identified through the stable paths returned by kernel_parameters(kernel, flatten = TRUE). Every parameter is optimized on an unconstrained coordinate determined by its constraint. Positive parameters, including the noise variance, are optimized on the log scale and stay positive. Real-valued parameters, such as the location of a changepoint_kernel(), are optimized on the identity scale. $optimization$coordinates records the coordinate of each parameter. The deterministic starts do not depend on gradient.

A changepoint location is bounded by the training range of its input column, and the starts after the first place it at evenly spaced quantiles of that column, because the likelihood is often multimodal in the location. A spectral-mixture frequency is bounded by 0 and the Nyquist frequency of its column, and later starts displace it by one frequency-resolution bin (spectral_mixture_kernel()).

The objective is the likelihood that log_marginal_likelihood() reports for the mean: $optimization$likelihood records which, and $optimization$mean_coefficients the coefficients at the optimum.

The returned diagnostics retain every start, convergence code, objective value, final parameters, and the gradient method. Evaluation counts are reported per start in start_results and summed over all starts:

function_evaluations: Evaluations of the negative log marginal likelihood, including those stats::optim() makes to approximate gradients by finite differences.
gradient_evaluations: Analytical gradient evaluations, or the number of finite-difference gradient approximations.
model_fits: GP fits, the dominant cost. With analytical gradients the value and gradient at the same point share one fit.

counts is the unmodified stats::optim() result for the best start; with numerical gradients it excludes the finite-difference evaluations.

$optimization$optima groups the starts by the optimum they reached: starts whose log marginal likelihoods differ by at most $10^{-3}$ and whose coordinates differ by at most $10^{-2}$ share one. It reports, for each distinct optimum, its starts, log marginal likelihood, and gap to the best. Several optima with large gaps mean that the starts matter; several with nearly equal likelihoods but different parameters suggest a ridge of poorly identified parameters (gp_hyperparameter_uncertainty()). summary() reports them.

Stability

Stable: from version 1.0.0 this interface changes incompatibly only in a major release, after a deprecation period. Results and options that concern an experimental model class, kernel, or argument follow that interface's tier. See gaussianprocesses-package for the policy. The noise_groups argument (one noise variance per group) is experimental.

Examples

x <- seq(-2, 2, length.out = 20)
y <- sin(2 * x) + 0.1 * cos(9 * x)

model <- optimize_gp(x, y, kernel = rbf_kernel(), noise_variance = 0.05, n_starts = 1)
model
#> Exact Gaussian-process regression model
#>   observations: 20
#>   input dimensions: 1
#>   mean: ZeroMean()
#>   noise structure: homoscedastic
#>   noise variance: 0.00737229
#>   numerical jitter: 0
#>   optimization: converged (L-BFGS-B, analytical gradient)
#>   evaluations: 20 objective, 20 gradient, 20 fits
#>   log marginal likelihood: 4.349764
#>   kernel:
#>     RBF(variance=1.47029, length_scale=1.04508)
#> 

kernel_parameters(model$kernel, flatten = TRUE)
#>     variance length_scale 
#>     1.470291     1.045076 
model$noise_variance
#> [1] 0.007372292