Optimize the hyperparameters of a latent Gaussian-process model

Estimates kernel parameters and, optionally, likelihood parameters by maximizing the Laplace approximation to the log marginal likelihood, with the optimizer of optimize_gp(): L-BFGS-B on unconstrained coordinates, deterministic multiple starts, and analytical gradients.

Usage

optimize_latent_gp(
  x,
  y,
  kernel,
  likelihood,
  mean = zero_mean(),
  exposure = NULL,
  optimize_likelihood = TRUE,
  n_starts = 3L,
  start_spread = 1,
  lower = 1e-08,
  upper = 1e+08,
  control = list(maxit = 500),
  method = "laplace",
  max_iterations = 100L,
  convergence_tolerance = 1e-14,
  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. Matrix rows are observations and columns are input dimensions.
y: Responses, one per row of x, as the likelihood accepts them: numbers for Gaussian likelihoods; 0/1, logical values, or a two-level factor whose second level is 1 for Bernoulli likelihoods; and non-negative integer counts for Poisson likelihoods.
kernel: Initial kernel specification.
likelihood: Likelihood specification; its parameters, if any, are the starting values of the likelihood parameters.
mean: Gaussian-process mean specification, with fixed coefficients, coefficients estimated by maximum likelihood, or coefficients with a Gaussian coefficient_prior().
exposure: Optional positive exposure for Poisson likelihoods: one value or one per observation (see poisson_likelihood()).
optimize_likelihood: If TRUE, optimize the likelihood parameters, such as the variance of gaussian_likelihood(), jointly with the kernel parameters. Likelihoods without parameters ignore it.
n_starts, start_spread, lower, upper, control, gradient: As in optimize_gp().
method: The approximation: "laplace".
max_iterations: Maximum number of Newton iterations.
convergence_tolerance: The Newton iterations stop when the Newton decrement is at most this value (see Details).
initial_jitter, fallback_jitter, jitter_multiplier, max_attempts: The jitter policy for the factorization of $B$, relative to its mean diagonal, as in fit_gp().
symmetry_tolerance: Relative tolerance for checking that $B$ is symmetric.

Value

A fitted gaussianprocesses_latent_model with optimized parameters and optimizer diagnostics in $optimization.

Details

The objective is $\log q(y \mid X, \theta)$ of fit_latent_gp(). Its gradient (Rasmussen and Williams, 2006, Algorithm 5.1) has an explicit part, with the mode $\hat f$ held fixed, and an implicit part through the dependence of $\hat f$ on $\theta$: $$\frac{\partial \hat f}{\partial \theta_j} = (I + K W)^{-1} \frac{\partial K}{\partial \theta_j} \nabla \log p(y \mid \hat f),$$ which needs the third derivatives of the log likelihood (evaluate_likelihood()). Likelihood parameters are named likelihood.<name>, for example likelihood.variance.

Mean coefficients estimated by maximum likelihood (estimate_coefficients("ml")) are optimized jointly with the other parameters, on the real line, as mean.<name> (for example mean.intercept), starting from the generalized linear model fitted without the latent process. Their gradient is $H^\top \nabla \log p(y \mid \hat f)$ plus the implicit part through $\partial \hat f / \partial \beta = (I + K W)^{-1} H$.

Each evaluation starts the Newton iterations from the previous mode when that is better than the prior mean, which saves most of the iterations once the optimizer takes small steps. The returned model is refitted at the optimum from the prior mean, so it equals fit_latent_gp() with the optimized parameters.

$optimization has the fields of optimize_gp(), with optimize_likelihood and approximation in place of the noise fields, and at_bound, the parameters that ended at lower or upper (on the natural scale) or at the limits that optimize_gp() describes for changepoint locations and spectral frequencies.

Bounds and separable data

Every positive parameter is searched within [lower, upper]. If any parameter ends at a bound, a warning of class gaussianprocesses_bound_warning names it, because the likelihood may still increase beyond the bound and the bound, not the data, then determines the estimate. With binary responses this happens when the classes are separable: larger signal variances let the latent function separate them more sharply, and the approximate marginal likelihood keeps increasing up to very large variances. On twenty separable points on a line it peaks only near a variance of 1800, with a nearly flat maximum, and a smaller upper ends the search at the bound with this warning. Setting upper to the largest plausible variance makes the limit explicit.

Stability

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

Examples

x <- seq(0, 10, length.out = 40)
counts <- c(1, 0, 2, 1, 3, 2, 4, 5, 3, 6, 7, 5, 8, 6, 9, 7, 6, 8, 5, 6,
  4, 5, 3, 4, 2, 3, 2, 1, 2, 1, 1, 0, 1, 2, 1, 3, 2, 4, 3, 5)

model <- optimize_latent_gp(x, counts, rbf_kernel(), poisson_likelihood(),
  n_starts = 1)
kernel_parameters(model$kernel, flatten = TRUE)
#>     variance length_scale 
#>     1.939215     2.654603