Fit a latent Gaussian-process model

Fits a Gaussian process $f \sim GP(m, k)$ observed through a likelihood $p(y_i \mid f(x_i))$, such as Bernoulli responses for classification or Poisson counts, with the Laplace approximation to the posterior of $f$. With gaussian_likelihood() the approximation is exact and the model equals fit_gp().

Usage

fit_latent_gp(
  x,
  y,
  kernel,
  likelihood,
  mean = zero_mean(),
  exposure = NULL,
  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)
)

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: A Gaussian-process kernel specification.
likelihood: A likelihood specification (gp_likelihoods). It must be log-concave.
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()).
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

An object of class gaussianprocesses_latent_model with fields x, y (responses in the internal coding), response_levels (the levels of a factor response), exposure, kernel, mean, likelihood, method, mode ($\hat f$), alpha ($a$), curvature (the diagonal of $W$), b_cholesky ($R$ with $R^\top R = B + \delta I$), jitter ($\delta$), log_marginal_likelihood, laplace, n_observations, n_features, schema_version, and package_version. Other fields are internal.

Details

The posterior of $f = f(X)$ is proportional to $p(y \mid f) N(f \mid m, K)$. The Laplace approximation replaces it by $N(\hat f, (K^{-1} + W)^{-1})$, where $\hat f$ maximizes $$\Psi(f) = \log p(y \mid f) - \tfrac12 (f - m)^\top K^{-1} (f - m)$$ and $W = -\nabla\nabla \log p(y \mid \hat f)$ is diagonal. For a log-concave likelihood $\Psi$ is strictly concave, so the mode is unique; other likelihoods raise an error of class gaussianprocesses_log_concavity_error.

Newton's method finds the mode (Rasmussen and Williams, 2006, Algorithm 3.1) in the parameterization $f = m + K a$. Each iteration factorizes $$B = I + W^{1/2} K W^{1/2},$$ whose eigenvalues are at least 1, so it is well conditioned even when $K$ is nearly singular; the factorization follows the jitter policy of fit_gp(). A step-halving line search on $\Psi$ prevents overshooting. The iterations stop when the Newton decrement $$\lambda^2 = g^\top (K^{-1} + W)^{-1} g, \qquad g = \nabla \Psi(f),$$ which is twice the increase of $\Psi$ that the next Newton step predicts, is at most convergence_tolerance. $\lambda$ is about the distance to the mode measured in posterior standard deviations, so the default of $10^{-14}$ places the mode within about $10^{-7}$ standard deviations. $laplace records the number of iterations and step halvings, the final decrement, the largest component of $g$, the last change of $\Psi$, and whether the iterations converged; if they did not, a warning of class gaussianprocesses_convergence_warning is raised.

The approximate log marginal likelihood is $$\log q(y \mid X, \theta) = \log p(y \mid \hat f) - \tfrac12 (\hat f - m)^\top K^{-1} (\hat f - m) - \tfrac12 \log |B|,$$ reported by log_marginal_likelihood() and logLik(), with its gradient in log_marginal_likelihood_gradient().

A mean with fixed coefficients contributes $m(x)$. Coefficients with a Gaussian prior $N(b, B_\beta)$ are integrated out: the prior becomes $m(x) = h(x)^\top b$ with covariance $k(x, x') + h(x)^\top B_\beta h(x')$. Coefficients estimated by maximum likelihood (estimate_coefficients("ml"), the default of linear_mean()) maximize $\log q(y \mid X, \theta, \beta)$ at the given hyperparameters, by BFGS with the analytical gradient $H^\top a + [(I + K W)^{-1} H]^\top s$, where $s$ is the derivative of $-\frac12 \log |B|$ in $\hat f$, starting from the generalized linear model without the latent process. The estimate is plugged in, like the hyperparameters; coef() returns it and $mean_fit records it. A constant mean estimated this way is the usual mean level of a count or classification model. Restricted likelihoods and vague priors are not supported.

Accuracy of the Laplace approximation

For a few observations the exact marginal likelihood $\log p(y \mid X)$ can be computed by tensor-product Gauss-Hermite quadrature. With an RBF kernel (variance 2, length scale 1) and the first $n$ of the inputs $x = 0, 0.5, 2$ and of the responses below, the error $\log q - \log p$ of the Laplace approximation is:

Likelihood	Responses	$n = 1$	$n = 2$	$n = 3$
logit	1, 1, 0	-0.0171	-0.0192	-0.0387
probit	1, 1, 0	-0.0196	-0.0298	-0.0450
Poisson	0, 2, 7	0.0125	-0.0101	-0.0197
Poisson, mean 5	120, 150, 200	-0.000681	-0.00118	-0.00158

For binary responses and small counts the error is a few hundredths and grows with $n$; the approximation can over- or underestimate. For counts in the hundreds, where the likelihood is nearly Gaussian in $f$, it is ten times smaller. The test suite reproduces these values.

Stability

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

References

Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press. Chapters 3 and 5.

Examples

x <- seq(-3, 3, length.out = 30)
y <- as.numeric(sin(x) + 0.3 * cos(5 * x) > 0)

model <- fit_latent_gp(x, y, rbf_kernel(variance = 4), bernoulli_likelihood())
model
#> Latent Gaussian-process model (Laplace approximation)
#>   observations: 30
#>   input dimensions: 1
#>   likelihood: BernoulliLikelihood(link=logit)
#>   mean: ZeroMean()
#>   Newton iterations: 5 (converged)
#>   approximate log marginal likelihood: -10.63073
#>   numerical jitter: 0
#>   kernel:
#>     RBF(variance=4, length_scale=1)
#> 
predict_latent_gp(model, c(-2, 0, 2))$probability
#> [1] 0.0832661 0.6497020 0.8880144

Likelihood	Responses	\(n = 1\)	\(n = 2\)	\(n = 3\)
logit	1, 1, 0	-0.0171	-0.0192	-0.0387
probit	1, 1, 0	-0.0196	-0.0298	-0.0450
Poisson	0, 2, 7	0.0125	-0.0101	-0.0197
Poisson, mean 5	120, 150, 200	-0.000681	-0.00118	-0.00158