RBF

With \(f = \exp(-r^2/2)\), \(f' = -f/2\), so

\[ g = \sigma^2 f = k, \qquad \frac{\partial k}{\partial \log \ell_d} = k \, s_d. \]

Matérn 1/2

With \(f = e^{-r}\) and \(r = \sqrt{r^2}\), \(f'(r^2) = -e^{-r} / (2r)\), so

\[ g = \frac{\sigma^2 e^{-r}}{r} = \frac{k}{r}. \]

The factor is singular at \(r = 0\), but the product is not: \(s_d \le r^2\) gives \(g\, s_d \le k\, r \to 0\). Coincident inputs therefore have zero length-scale derivative.

Matérn 3/2

With \(f = (1 + \sqrt{3} r) e^{-\sqrt{3} r}\), \(\mathrm{d}f/\mathrm{d}r = -3 r e^{-\sqrt{3} r}\), so \(f'(r^2) = -\tfrac{3}{2} e^{-\sqrt{3} r}\) and

\[ g = 3 \sigma^2 e^{-\sqrt{3} r} = \frac{3k}{1 + \sqrt{3} r}. \]

Matérn 5/2

With \(f = (1 + \sqrt{5} r + \tfrac{5}{3} r^2) e^{-\sqrt{5} r}\), \(\mathrm{d}f/\mathrm{d}r = -\tfrac{5}{3} r (1 + \sqrt{5} r) e^{-\sqrt{5} r}\), so

\[ g = \frac{5}{3} \sigma^2 (1 + \sqrt{5} r) e^{-\sqrt{5} r} = \frac{5}{3}\, k\, \frac{1 + \sqrt{5} r}{1 + \sqrt{5} r + \tfrac{5}{3} r^2}. \]

Rational quadratic

Write \(u = 1 + r^2 / (2\alpha)\), so \(f = u^{-\alpha}\). Then \(f'(r^2) = -\tfrac{1}{2} u^{-\alpha - 1}\) and

\[ g = \sigma^2 u^{-\alpha - 1} = \frac{k}{u}. \]

For the shape parameter, \(\log k = \log \sigma^2 - \alpha \log u\) and \(\partial u / \partial \alpha = -r^2 / (2\alpha^2)\), which give

\[ \frac{\partial k}{\partial \log \alpha} = \alpha\, k\, \frac{\partial \log k}{\partial \alpha} = k \left( \frac{r^2}{2u} - \alpha \log u \right). \]

The implementation evaluates \(\log u\) as log1p(r^2 / (2 * alpha)), which stays accurate when \(r^2 / (2\alpha)\) is small.

Derivation

With centred observations \(\tilde y = y - m(X)\), the mean held fixed, and \(C = K(\theta) + \Sigma_\varepsilon\), the log marginal likelihood is

\[ \log p(y \mid X, \theta) = -\frac{1}{2} \tilde y^\top C^{-1} \tilde y -\frac{1}{2} \log |C| -\frac{n}{2} \log 2\pi. \]

For a symmetric positive-definite \(C\) that depends on a parameter \(\eta\),

\[ \frac{\partial C^{-1}}{\partial \eta} = -C^{-1} \frac{\partial C}{\partial \eta} C^{-1}, \qquad \frac{\partial \log |C|}{\partial \eta} = \operatorname{tr}\!\left( C^{-1} \frac{\partial C}{\partial \eta} \right). \]

With \(\alpha = C^{-1} \tilde y\), the first identity turns the derivative of the quadratic term into \(\tfrac{1}{2} \alpha^\top (\partial C / \partial \eta) \alpha = \tfrac{1}{2} \operatorname{tr}(\alpha\alpha^\top \partial C / \partial \eta)\), and the second gives the determinant term. Together, for \(\eta_j = \log \theta_j\),

\[ \frac{\partial \log p(y \mid X, \theta)}{\partial \eta_j} = \frac{1}{2} \operatorname{tr}\!\left( \left(\alpha\alpha^\top - C^{-1}\right) \frac{\partial C}{\partial \eta_j} \right). \]

For kernel parameters \(\partial C / \partial \eta_j = \partial K / \partial \eta_j\), as derived above. For a homoscedastic noise variance \(\sigma_n^2\), \(C = K + \sigma_n^2 I\) gives \(\partial C / \partial \log \sigma_n^2 = \sigma_n^2 I\), so

\[ \frac{\partial \log p}{\partial \log \sigma_n^2} = \frac{\sigma_n^2}{2} \left( \alpha^\top \alpha - \operatorname{tr}(C^{-1}) \right). \]

Known observation-specific noise variances are data, not parameters, so they have no derivative. If the Cholesky factorization needed numerical jitter \(\delta\), the factorized matrix is \(C + \delta I\). The jitter is a constant, so the gradient is that of the jitter-stabilized likelihood reported by log_marginal_likelihood().

Computation

Each evaluation reuses the Cholesky factor \(C = R^\top R\) computed for the likelihood; \(\alpha\) already comes from two triangular solves.

• The quadratic term \(\alpha^\top (\partial C / \partial \eta_j) \alpha\) costs \(O(n^2)\) per parameter.
• Because both matrices are symmetric, the trace term is \(\operatorname{tr}(C^{-1} \partial C / \partial \eta_j) = \sum_{a,b} (C^{-1})_{ab} (\partial C / \partial \eta_j)_{ab}\). Every entry of \(C^{-1}\) enters, so the entries are computed once per evaluation from the Cholesky factor with chol2inv(). This is the Cholesky solve of \(C X = I\), the same computation as cho_solve(L, I) in scikit-learn. \(C\) is never passed to solve() or to a general-purpose inverse, and the precision matrix is never used to solve linear systems. Each parameter then costs another \(O(n^2)\).

The alternative that never forms \(C^{-1}\) solves \(C X_j = \partial C / \partial \eta_j\) with the same factor for every parameter and takes \(\operatorname{tr}(X_j)\). It gives identical traces (the package tests check this) but costs \(O(n^3)\) per parameter. With \(n = 800\) and seven parameters it took about seven GP fits per gradient on a local machine, against less than half a fit through the precision matrix.

An analytical gradient therefore costs about one GP fit plus the precision matrix, whatever the number of parameters. Central finite differences need two fits per parameter. inst/benchmarks/optimizer-gradients.R measures both methods end to end.

x_train <- seq(-2, 2, length.out = 20)
model <- fit_gp(
  x_train,
  sin(2 * x_train),
  kernel = rbf_kernel(variance = 1, length_scale = 0.6),
  noise_variance = 0.02
)

log_marginal_likelihood_gradient(model)
#>       variance   length_scale noise_variance 
#>      -2.405276       9.770035      -5.635096

Derivative processes

Differentiation is a linear operation, so the partial derivatives of a Gaussian process \(f \sim \mathcal{GP}(m, k)\) are jointly Gaussian with the process itself. Their covariances are derivatives of the kernel:

\[ \operatorname{cov}\!\left(\frac{\partial f(x)}{\partial x_d}, f(x')\right) = \frac{\partial k(x, x')}{\partial x_d}, \qquad \operatorname{cov}\!\left( \frac{\partial f(x)}{\partial x_d}, \frac{\partial f(x')}{\partial x'_e} \right) = \frac{\partial^2 k(x, x')}{\partial x_d\, \partial x'_e}. \]

The derivative process exists in mean square exactly when the mixed second derivative exists at \(x = x'\). For a stationary kernel this means that the kernel is twice differentiable at zero distance. kernel_input_gradient() computes the first matrices, and the second ones give the variance of the gradient.

Smoothness

The number of mean-square derivatives of each process is a property of its kernel. Sums, products, and scaled kernels are as smooth as their least smooth component.

Asking for derivatives of a process that has none raises an error of class gaussianprocesses_smoothness_error. Observation noise belongs in noise_variance, not in a white-noise kernel component, so that the latent function stays differentiable.

The Matérn 3/2 process is differentiable once: its gradient process exists, but that process is not itself differentiable.

Radial kernels

For a radial kernel \(k = \kappa(r)\), with \(r^2 = \sum_d (x_d - x'_d)^2 / \ell_d^2\) as before and

\[ u_d = \frac{x_d - x'_d}{\ell_d^2}, \qquad \frac{\partial r}{\partial x_d} = \frac{u_d}{r}, \qquad \frac{\partial u_d}{\partial x'_e} = -\frac{\delta_{de}}{\ell_d^2}, \]

the chain rule gives

\[ \frac{\partial k}{\partial x_d} = a\, u_d, \qquad \frac{\partial^2 k}{\partial x_d\, \partial x'_e} = -b\, u_d u_e - a\, \frac{\delta_{de}}{\ell_d^2}, \]

with

\[ a = \frac{\kappa'(r)}{r}, \qquad b = \frac{1}{r^2}\left(\kappa''(r) - \frac{\kappa'(r)}{r}\right). \]

Both coefficients have closed forms in which \(r = 0\) needs no limit.

Here \(w = 1 + r^2 / (2\alpha)\). The Matérn 3/2 coefficient \(b\) is singular at \(r = 0\), but it multiplies \(u_d u_e = O(r^2)\), so the product is zero there. At coincident inputs every cross derivative reduces to \(-a(0) \delta_{de} / \ell_d^2\). This is the prior variance of \(\partial f / \partial x_d\): \(\sigma^2 / \ell_d^2\) for the RBF and rational-quadratic kernels, \(3\sigma^2 / \ell_d^2\) for Matérn 3/2, and \(\tfrac{5}{3} \sigma^2 / \ell_d^2\) for Matérn 5/2.

The Matérn 1/2 kernel has \(a = -\sigma^2 e^{-r} / r\), which is unbounded at \(r = 0\): its process is continuous but not differentiable.

Periodic and linear kernels

For the periodic kernel, with \(a_d = \pi (x_d - x'_d) / p\) as before and

\[ g_d = -\frac{2\pi}{p\, \ell_d^2} \sin 2 a_d, \]

\[ \frac{\partial k}{\partial x_d} = k\, g_d, \qquad \frac{\partial^2 k}{\partial x_d\, \partial x'_e} = -k\, g_d g_e + \delta_{de}\, k\, \frac{4 \pi^2}{p^2 \ell_d^2} \cos 2 a_d. \]

The linear kernel \(k = \sigma^2 x^\top x'\) has \(\partial k / \partial x_d = \sigma^2 x'_d\) and \(\partial^2 k / \partial x_d\, \partial x'_e = \sigma^2 \delta_{de}\).

Composite kernels

Sums and scaled kernels are differentiated term by term. For a product \(k = k_1 k_2\), the product rule gives

\[ \frac{\partial^2 k}{\partial x_d\, \partial x'_e} = \frac{\partial^2 k_1}{\partial x_d\, \partial x'_e} k_2 + \frac{\partial k_1}{\partial x_d} \frac{\partial k_2}{\partial x'_e} + \frac{\partial k_1}{\partial x'_e} \frac{\partial k_2}{\partial x_d} + k_1 \frac{\partial^2 k_2}{\partial x_d\, \partial x'_e}, \]

so a product needs the derivatives of each factor with respect to both arguments. Longer products apply the rule pairwise.

kernel <- product_kernel(
  matern52_kernel(variance = 1.2, length_scale = c(0.6, 1.4)),
  linear_kernel(variance = 0.5)
)
x <- cbind(c(-0.5, 0, 0.8), c(1, 0.3, -0.2))

gradient <- kernel_input_gradient(kernel, x)

step <- 1e-6
shift <- cbind(step, rep(0, nrow(x)))
numerical <- (evaluate_kernel(kernel, x + shift, x) -
  evaluate_kernel(kernel, x - shift, x)) / (2 * step)

max(abs(gradient[[1]] - numerical))
#> [1] 6.413786e-11

The posterior gradient

For an exact model with factorized covariance \(C = R^\top R\) and \(\alpha = C^{-1}(y - m(X))\), conditioning the joint Gaussian on the observations gives

\[ \mathrm{E}\!\left[\frac{\partial f(x_*)}{\partial x_{*,d}} \,\middle|\, y\right] = \frac{\partial m(x_*)}{\partial x_{*,d}} + \frac{\partial k(x_*, X)}{\partial x_{*,d}}\, \alpha, \]

\[ \operatorname{Var}\!\left[\frac{\partial f(x_*)}{\partial x_{*,d}} \,\middle|\, y\right] = \left.\frac{\partial^2 k(x_*, x')}{\partial x_{*,d}\, \partial x'_d}\right|_{x' = x_*} - v_d^\top v_d, \qquad v_d = R^{-\top} \frac{\partial k(X, x_*)}{\partial x_{*,d}}. \]

The zero and constant mean functions have zero gradient. predict_gradient_gp() computes these quantities in blocks of prediction inputs, as predict_gp() does, so memory use does not grow with the number of prediction inputs.

x_train <- seq(0, 2 * pi, length.out = 30)
model <- fit_gp(
  x_train,
  sin(x_train),
  kernel = rbf_kernel(length_scale = 1.2),
  noise_variance = 1e-4
)

x_new <- seq(0, 2 * pi, length.out = 7)
gradient <- predict_gradient_gp(model, x_new)

data.frame(
  x = round(x_new, 3),
  estimate = gradient$mean[, 1],
  sd = gradient$sd[, 1],
  derivative = cos(x_new)
)
#>       x   estimate         sd derivative
#> 1 0.000  0.9759805 0.06870568        1.0
#> 2 1.047  0.4979566 0.01447450        0.5
#> 3 2.094 -0.4993927 0.01318788       -0.5
#> 4 3.142 -0.9995817 0.01274587       -1.0
#> 5 4.189 -0.4993927 0.01318788       -0.5
#> 6 5.236  0.4979566 0.01447450        0.5
#> 7 6.283  0.9759805 0.06870568        1.0

The estimates follow \(\cos x\), the derivative of the data-generating function. The uncertainty is largest at the ends of the data, where differences between neighbouring observations constrain the slope least.

A model whose kernel has no derivatives says so:

rough <- fit_gp(x_train, sin(x_train), matern12_kernel(), noise_variance = 1e-4)

tryCatch(
  predict_gradient_gp(rough, 1),
  gaussianprocesses_smoothness_error = conditionMessage
)
#> [1] "The kernel is not mean-square differentiable, so its derivatives with respect to the inputs do not exist: Matern-1/2. Use smoother components, and observation noise instead of a white-noise kernel."

Verification

The test suite checks the input derivatives of every differentiable kernel, and of sums, products, scaled and nested kernels, in one and three dimensions with ARD length scales. Each check compares with central differences at input pairs a distance \(r \in \{0, 10^{-8}, 10^{-4}, 1\}\) apart, to a relative error of \(10^{-6}\). Kernels without a Matérn 3/2 component agree to better than \(10^{-8}\). The Matérn 3/2 cross derivative at coincident inputs converges only at first order in the step, because \(\kappa\) has a \(|r|^3\) term, and agrees to about \(2.5 \times 10^{-7}\).

The tests also check that

• the joint covariance of values and gradients is positive semidefinite, including at repeated inputs;
• the posterior gradient mean equals a central difference of the predict_gp() mean, to \(10^{-6}\);
• the posterior gradient variance equals the limit, by Richardson extrapolation, of the variance of \((f(x_* + h e_d) - f(x_* - h e_d)) / 2h\) computed from the full posterior covariance, to \(10^{-5}\);
• the results do not depend on the block size.

The joint model

The same covariances let a model observe derivatives. In engineering and physics, gradients are often cheap: an adjoint solver returns them with each function evaluation. Because differentiation is linear, values at \(X\) and partial derivatives at \(X'\) are jointly Gaussian, and

\[ \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'_e} 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 observation has its own direction \(d\) and the blocks are the input derivatives above. Every observation has its own noise variance in \(\Sigma\). Writing this matrix as \(C\), conditioning,

\[ \mathrm{E}[f(x_*) \mid y] = m(x_*) + k_*^\top C^{-1}(y - \mu), \qquad \operatorname{Var}[f(x_*) \mid y] = k(x_*, x_*) - k_*^\top C^{-1} k_*, \]

the log marginal likelihood, and its gradient

\[ \frac{1}{2} \operatorname{tr}\!\left( \left(\alpha\alpha^\top - C^{-1}\right) \frac{\partial C}{\partial \eta_j} \right) \]

keep their form. Only \(C\), the cross-covariances \(k_*\), and \(\partial C / \partial \eta_j\) change. The last needs the hyperparameter derivatives of the input derivatives. Every differentiable kernel supplies them in closed form; for a radial kernel with length scale \(\ell_c\) they follow from \(\partial a / \partial (r^2) = b / 2\) and \(\partial b / \partial (r^2) = c / 2\) with \(c = b'(r) / r\), for example

\[ \frac{\partial (a\, u_d)}{\partial \log \ell_c} = -s_c\, b\, u_d - 2 a\, u_d\, \delta_{dc}, \qquad s_c = \frac{(x_c - x'_c)^2}{\ell_c^2}. \]

Sums, scaled kernels, and selections of input columns are linear, and the product rule for \(k = k_1 k_2\) is bilinear in the two factors, so a hyperparameter derivative of a product is the product rule applied to one differentiated factor and one undifferentiated one.

fit_gp() and optimize_gp() take the type of each observation in derivative: 0 for a value and \(d\) for \(\partial f / \partial x_d\). Parametric means contribute \(\partial h(x)^\top \beta / \partial x_d\) to derivative rows; a constant mean contributes 0.

Example

Observe \(f(x) = \sin x\) at five points, then add its derivative \(\cos x\) at the same points:

x_observed <- c(-2, -1, 0, 1, 2)
kernel <- rbf_kernel(length_scale = 0.8)

values_only <- fit_gp(
  x_observed,
  sin(x_observed),
  kernel = kernel,
  noise_variance = 1e-4
)
with_gradients <- fit_gp(
  c(x_observed, x_observed),
  c(sin(x_observed), cos(x_observed)),
  kernel = kernel,
  noise_variance = 1e-4,
  derivative = rep(c(0, 1), each = 5)
)
with_gradients
#> Exact Gaussian-process regression model
#>   observations: 10 (5 values, 5 derivatives)
#>   input dimensions: 1
#>   mean: ZeroMean()
#>   noise structure: homoscedastic
#>   noise variance: 1e-04
#>   numerical jitter: 0
#>   kernel:
#>     RBF(variance=1, length_scale=0.8)

The joint covariance is a valid covariance matrix. The fitted model stores its Cholesky factor \(R\), with \(R^\top R = C\), and the smallest eigenvalue of \(C\) is positive:

joint <- crossprod(with_gradients$cholesky)
min(eigen(joint, symmetric = TRUE, only.values = TRUE)$values)
#> [1] 0.004783267

Between the observations, the gradients remove most of the posterior uncertainty:

grid <- seq(-2.5, 2.5, length.out = 201)
before <- predict_gp(values_only, grid)
after <- predict_gp(with_gradients, grid)

plot(
  grid,
  before$latent_sd,
  type = "l",
  lty = 2,
  ylim = c(0, max(before$latent_sd)),
  xlab = "x",
  ylab = "posterior standard deviation"
)
lines(grid, after$latent_sd)
points(x_observed, rep(0, 5), pch = 19)
legend(
  "top",
  legend = c("values", "values and gradients"),
  lty = c(2, 1),
  bty = "n"
)

between <- abs(grid - round(grid)) > 0.4 & abs(grid) < 2
c(
  values = max(before$latent_sd[between]),
  values_and_gradients = max(after$latent_sd[between])
)
#>               values values_and_gradients 
#>           0.22944994           0.01335559

The posterior gradient is constrained too. At the observed points it reproduces \(\cos x\) up to the noise:

predict_gradient_gp(with_gradients, x_observed)$mean[, 1] - cos(x_observed)
#> [1] 1.012559e-05 2.584515e-05 2.948141e-05 2.584515e-05 1.012559e-05

Derivative blocks scale with \(1 / \ell^2\), so they can be conditioned differently from the values. gp_numerical_diagnostics() reports each observation type separately:

gp_numerical_diagnostics(with_gradients)$observation_types
#>            type observations covariance_scale reciprocal_condition_number
#> 1         value            5           1.0001                   0.1095987
#> 2 derivative x1            5           1.5626                   0.1157423

Verification

The tests check that conditioning on a derivative equals the limit \(h \to 0\) of conditioning on difference quotients \((f(x + h e_d) - f(x - h e_d)) / 2h\), for the posterior mean, covariance, and log marginal likelihood, by Richardson extrapolation to \(10^{-6}\). They also check that the log-marginal-likelihood gradient, including the derivative blocks, matches central differences to \(10^{-6}\) for single, ARD, periodic, product, sum, scaled, and column-restricted kernels, and that the hyperparameter derivatives of every kernel’s input derivatives match central differences.