Gradient of the Gaussian-process log marginal likelihood — log_marginal_likelihood

Computes the exact gradient of log_marginal_likelihood() with respect to the logarithm of every kernel parameter and, for models with a single homoscedastic noise variance, the logarithm of that noise variance.

Usage

log_marginal_likelihood_gradient(model)

Arguments

model: A fitted gaussianprocesses_model, gaussianprocesses_latent_model, or gaussianprocesses_state_space_model.

Value

A named numeric vector. Names follow kernel_parameters(model$kernel, flatten = TRUE), followed by noise_variance for homoscedastic models, or by likelihood.<name> for the parameters of the likelihood of a latent model.

Details

Let $C$ be the observation covariance and $\alpha = C^{-1}(y - m(X))$. For $\eta_j = \log \theta_j$, $$\frac{\partial \log p(y \mid X, \theta)}{\partial \eta_j} = \frac{1}{2} \mathrm{tr}\left((\alpha \alpha^\top - C^{-1}) \frac{\partial C}{\partial \eta_j}\right).$$ For kernel parameters, $\partial C / \partial \eta_j$ is given by kernel_gradient(); for the noise variance $\sigma_n^2$ it is $\sigma_n^2 I$.

The trace requires every entry of the precision matrix $C^{-1}$. They are computed once from the stored Cholesky factor with chol2inv(), the Cholesky solve of $C X = I$; $C$ is never passed to solve() or a general-purpose inverse, and the precision matrix is not used to solve linear systems. The alternative, two triangular solves per parameter, never forms $C^{-1}$; against an exact reference it was equally accurate for condition numbers up to $10^{12}$ and 4 to 18 times slower (inst/benchmarks/gradient-numerics.R). As in log_marginal_likelihood(), any numerical jitter used during fitting is part of the factorized covariance and is treated as a constant.

Models fitted with observation-specific noise variances have no noise parameter, so only kernel derivatives are returned for them.

For a latent model the gradient is that of the Laplace approximation, including the implicit dependence of the mode on the parameters (optimize_latent_gp()), in the optimizer coordinates of the kernel parameters and of the likelihood parameters.

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.

Examples

x <- seq(-2, 2, length.out = 15)
model <- fit_gp(
  x,
  sin(x),
  kernel = rbf_kernel(length_scale = 0.8),
  noise_variance = 0.01
)
log_marginal_likelihood_gradient(model)
#>       variance   length_scale noise_variance 
#>      -2.574669      10.342921      -3.859761