Fits a Gaussian-process regression model by conditioning a Gaussian prior on observed responses. The implementation uses a Cholesky factorization and triangular solves; it never forms an explicit matrix inverse.
Arguments
- x
Numeric vector or matrix of training inputs. Matrix rows are observations and columns are input dimensions.
- y
Numeric response vector with one value per training observation.
- kernel
A Gaussian-process kernel specification.
- noise_variance
Non-negative scalar observation-noise variance or a non-negative vector with one known variance per training observation.
- mean
Gaussian-process mean specification:
zero_mean(),constant_mean(), or a parametric mean such aslinear_mean()whose coefficients are fixed, estimated, or marginalized (mean_coefficients).- initial_jitter
Non-negative jitter tried first, relative to the covariance scale (see Numerical jitter).
- fallback_jitter
Positive relative jitter tried after the first failed factorization.
- jitter_multiplier
Multiplicative increase applied to jitter after each failed factorization.
- max_attempts
Maximum number of Cholesky attempts.
- symmetry_tolerance
Relative tolerance for checking that the covariance matrix is symmetric.
- derivative
Optional observation types, one per observation: 0 if
y[i]observes the function value \(f(x_i)\), and \(d\) if it observes the derivative \(\partial f(x_i) / \partial x_d\).NULL(the default) means every observation is a function value.
Details
A scalar noise_variance gives the usual homoscedastic model. A vector
adds a known diagonal observation-noise variance for each training point.
Observation noise and numerical jitter are always stored separately.
Derivative observations
Differentiation is linear, so derivatives of a Gaussian process are
jointly Gaussian with its values, and observing them is exact
conditioning. With values at \(X\) and derivatives at \(X'\),
$$\begin{pmatrix} y \\ y' \end{pmatrix} \sim N\left(
\begin{pmatrix} m(X) \\ \partial_d m(X') \end{pmatrix},
\begin{pmatrix} K(X, X) & \partial_{x'_d} K(X, X') \\
\partial_{x_d} K(X', X) & \partial_{x_d} \partial_{x'_e} K(X', X')
\end{pmatrix} + \Sigma \right),$$
where each derivative row uses its own direction \(d\). The blocks come
from the kernel's input derivatives (kernel_input_gradient()), so the
kernel must be mean-square differentiable; otherwise an error of class
gaussianprocesses_smoothness_error is raised. Parametric means
contribute their input derivatives. noise_variance may differ between
value and derivative observations.
Conditioning, prediction, log_marginal_likelihood(), its gradient, and
optimize_gp() work unchanged with this covariance;
predict_gp() and predict_gradient_gp() predict values and
derivatives at new inputs. Derivative blocks scale with
\(1 / \ell^2\), so gp_numerical_diagnostics() reports the
conditioning of each observation type separately.
Numerical jitter
If \(C = K + \Sigma_\varepsilon\) does not factorize, the factorization
is retried with \(C + \delta I\), where \(\delta = \epsilon s\) and
\(s\) is the covariance scale, the mean of the diagonal of \(C\). The
relative jitter \(\epsilon\) starts at initial_jitter (if positive),
then takes the values fallback_jitter, fallback_jitter * jitter_multiplier, and so on, for at most max_attempts attempts.
Because \(\delta\) is proportional to \(C\), results do not depend
on the units of y: multiplying y by \(c\) and every variance
parameter by \(c^2\) multiplies the jitter by \(c^2\). The fitted
model reports the absolute jitter \(\delta\) in jitter, and
gp_numerical_diagnostics() also reports \(\epsilon\) and \(s\).
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 derivative argument (derivative observations) is experimental.
See also
gaussianprocesses_model for the structure, versioning, and persistence of fitted models.
Examples
x <- seq(-2, 2, length.out = 15)
y <- sin(2 * x)
model <- fit_gp(x, y, kernel = rbf_kernel(length_scale = 0.6), noise_variance = 0.01)
model
#> Exact Gaussian-process regression model
#> observations: 15
#> input dimensions: 1
#> mean: ZeroMean()
#> noise structure: homoscedastic
#> noise variance: 0.01
#> numerical jitter: 0
#> kernel:
#> RBF(variance=1, length_scale=0.6)
#>
# Known, observation-specific noise variances.
fit_gp(
x,
y,
kernel = rbf_kernel(length_scale = 0.6),
noise_variance = seq(0.01, 0.1, length.out = 15)
)
#> Exact Gaussian-process regression model
#> observations: 15
#> input dimensions: 1
#> mean: ZeroMean()
#> noise structure: known_heteroscedastic
#> training noise variance range: [0.01, 0.1]
#> numerical jitter: 0
#> kernel:
#> RBF(variance=1, length_scale=0.6)
#>