EulerPad

Practice / Probability for ML

Gaussian process conditioning

Ten problems on Gaussian process regression: the predictive mean and covariance from the joint Gaussian, two training points worked by hand, the posterior mean as a sum of kernels and why the variance only falls, the RBF kernel's derivatives, the gradient of the log marginal likelihood with respect to the lengthscale, signal variance and noise, the noise gradient through the eigenvalues, conditioning one observation at a time, the linear kernel as ridge regression, the lengthscale limits, and the derivative of the posterior mean with respect to the lengthscale, with worked solutions and the mistakes that leave out the noise, the log-determinant, or the second term of a derivative through an inverse.

Before you start

A Gaussian process is a prior over functions under which any finite set of function values is jointly Gaussian with a covariance given by a kernel, and Gaussian process regression is nothing more than the conditioning formula of the multivariate-Gaussian page applied to that joint: condition the unseen values on the observed ones. Every formula in the subject, the predictive mean and variance, the marginal likelihood, its hyperparameter gradients, comes from that one step plus the inverse-and-determinant calculus of the previous page. These ten problems do the step, run it on two data points by hand, read what the result says (a weighted sum of kernels, a variance that can only fall), differentiate the kernel and then the marginal likelihood with respect to the lengthscale, the signal variance and the noise, condition one observation at a time, show that a linear kernel recovers ridge regression, and take the lengthscale to its two limits. The five mistakes are each a missing term: the noise left out of the matrix that is inverted, noise added where none belongs, a marginal likelihood gradient without its log-determinant, a lengthscale derivative short a factor of 22, and a derivative of the posterior mean that holds the inverse fixed.

  • ff is a function with a Gaussian process prior with zero mean and kernel kk: for any inputs x1,…,xnx_1, \dots, x_n, the vector (f(x1),…,f(xn))(f(x_1), \dots, f(x_n)) is Gaussian with mean 00 and covariance matrix Kij=k(xi,xj)K_{ij} = k(x_i, x_j). Kernels here are positive definite, so every such matrix is symmetric positive semidefinite.
  • Training inputs X=(x1,…,xn)X = (x_1, \dots, x_n) with function values f=(f(x1),…,f(xn))⊤f = (f(x_1), \dots, f(x_n))^\top and noisy observations y=f+εy = f + \varepsilon, ε∼N(0,σ2I)\varepsilon \sim \mathcal N(0, \sigma^2I) independent of ff. Test inputs X∗=(x1∗,…,xm∗)X_* = (x_1^*, \dots, x_m^*) with values f∗∈Rmf_* \in \mathbb{R}^m.
  • Kernel blocks: K=k(X,X)K = k(X, X) is n×nn\times n, K∗=k(X,X∗)K_* = k(X, X_*) is n×mn\times m with (K∗)ij=k(xi,xj∗)(K_*)_{ij} = k(x_i, x_j^*), and K∗∗=k(X∗,X∗)K_{**} = k(X_*, X_*) is m×mm\times m. For a single test point, k∗=k(X,x∗)∈Rnk_* = k(X, x_*) \in \mathbb{R}^n and k∗∗=k(x∗,x∗)k_{**} = k(x_*, x_*). Ky=K+σ2IK_y = K + \sigma^2I and α=Ky−1y\alpha = K_y^{-1}y.
  • Gaussian conditioning (the multivariate-Gaussian page, Problem 9): if (xa,xb)(x_a, x_b) is jointly Gaussian with mean (μa,μb)(\mu_a, \mu_b) and covariance blocks Σaa\Sigma_{aa}, Σab\Sigma_{ab}, Σbb\Sigma_{bb}, then xa∣xb∼N(μa+ΣabΣbb−1(xb−μb), Σaa−ΣabΣbb−1Σba)x_a \mid x_b \sim \mathcal N\big(\mu_a + \Sigma_{ab}\Sigma_{bb}^{-1}(x_b - \mu_b),\ \Sigma_{aa} - \Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba}\big).
  • The RBF (squared-exponential) kernel on scalar or vector inputs is k(x,x′)=s2exp⁡(−∥x−x′∥2/(2ℓ2))k(x, x') = s^2\exp\big(-\|x - x'\|^2/(2\ell^2)\big) with signal variance s2s^2 and lengthscale ℓ>0\ell > 0. Its hyperparameters, with σ2\sigma^2, are written θ\theta. Numbers on this page use scalar inputs.
  • From the previous page: d(A−1)=−A−1(dA)A−1d(A^{-1}) = -A^{-1}(dA)A^{-1}, dlog⁡det⁡A=tr⁡(A−1dA)d\log\det A = \operatorname{tr}(A^{-1}dA), and tr⁡\operatorname{tr} is cyclic. tr⁡(A)\operatorname{tr}(A) of a symmetric matrix with eigenvalues λi\lambda_i is ∑iλi\sum_i\lambda_i.

Builds on: The multivariate Gaussian: gradients and identities, Matrix inverse and log-determinant derivatives

Problems

  1. ··

    Write the joint distribution of (f∗,y)(f_*, y) and condition on yy to obtain the predictive distribution f∗∣y∼N(μ∗,Σ∗)f_* \mid y \sim \mathcal N(\mu_*, \Sigma_*) with μ∗=K∗⊤Ky−1y\mu_* = K_*^\top K_y^{-1}y and Σ∗=K∗∗−K∗⊤Ky−1K∗\Sigma_* = K_{**} - K_*^\top K_y^{-1}K_*.

  2. ··

    Two training points x1=0x_1 = 0, x2=1x_2 = 1 with y=(1,2)⊤y = (1, 2)^\top, RBF kernel with s2=1s^2 = 1, ℓ=1\ell = 1, noise σ2=0.1\sigma^2 = 0.1. Compute the predictive mean and variance at x∗=0.5x_* = 0.5 to four decimal places.

  3. ··

    Show that the posterior mean at any x∗x_* is μ(x∗)=∑i=1nαik(xi,x∗)\mu(x_*) = \sum_{i=1}^n\alpha_ik(x_i, x_*) with α=Ky−1y\alpha = K_y^{-1}y not depending on x∗x_*; that the posterior variance never exceeds the prior variance; and that with σ2=0\sigma^2 = 0 the posterior mean passes through every observation with zero variance there.

  4. ·

    For the RBF kernel with r=∥x−x′∥r = \|x - x'\|, compute ∂k/∂ℓ\partial k/\partial\ell, ∂k/∂s2\partial k/\partial s^2 and, for a scalar test input, ∂k(xi,x∗)/∂x∗\partial k(x_i, x_*)/\partial x_*. Use the last to write the slope μ′(x∗)\mu'(x_*) of the posterior mean.

  5. ···

    The log marginal likelihood is log⁡p(y∣X,θ)=−12y⊤Ky−1y−12log⁡det⁡Ky−n2log⁡2π\log p(y \mid X, \theta) = -\tfrac12y^\top K_y^{-1}y - \tfrac12\log\det K_y - \tfrac n2\log2\pi. For a hyperparameter θ\theta with ∂Ky/∂θ\partial K_y/\partial\theta known, show that ∂∂θlog⁡p(y∣X,θ)=12tr⁡((αα⊤−Ky−1)∂Ky∂θ)\dfrac{\partial}{\partial\theta}\log p(y \mid X, \theta) = \tfrac12\operatorname{tr}\Big((\alpha\alpha^\top - K_y^{-1})\dfrac{\partial K_y}{\partial\theta}\Big), and write the two terms for θ=ℓ\theta = \ell.

  6. ··

    Let KK have eigenvalues λ1,…,λn≥0\lambda_1, \dots, \lambda_n \ge 0. Show that log⁡det⁡Ky=∑ilog⁡(λi+σ2)\log\det K_y = \sum_i\log(\lambda_i + \sigma^2) and tr⁡(Ky−1)=∑i1/(λi+σ2)\operatorname{tr}(K_y^{-1}) = \sum_i1/(\lambda_i + \sigma^2), hence that ∂log⁡p∂σ2=12(∥α∥2−tr⁡Ky−1)\dfrac{\partial\log p}{\partial\sigma^2} = \tfrac12\big(\|\alpha\|^2 - \operatorname{tr}K_y^{-1}\big), and state the condition at a maximum over σ2\sigma^2.

  7. ···

    After conditioning on nn observations, the posterior over (f∗,f(xn+1))(f_*, f(x_{n+1})) is Gaussian with mean (μ∗,μn+1)(\mu_*, \mu_{n+1}), covariance Σ∗\Sigma_* for f∗f_*, cross-covariance c=Cov⁡(f∗,f(xn+1))c = \operatorname{Cov}(f_*, f(x_{n+1})) and variance v=Var⁡(f(xn+1))v = \operatorname{Var}(f(x_{n+1})). A new observation yn+1=f(xn+1)+εn+1y_{n+1} = f(x_{n+1}) + \varepsilon_{n+1} arrives. Show that the updated posterior is μ∗′=μ∗+c (yn+1−μn+1)v+σ2\mu_*' = \mu_* + \dfrac{c\,(y_{n+1} - \mu_{n+1})}{v + \sigma^2} and Σ∗′=Σ∗−cc⊤v+σ2\Sigma_*' = \Sigma_* - \dfrac{cc^\top}{v + \sigma^2}, and say why this equals the batch posterior on all n+1n + 1 points.

  8. ··

    Take the linear kernel k(x,x′)=x⊤x′k(x, x') = x^\top x' with inputs as rows of X∈Rn×dX \in \mathbb{R}^{n\times d} and X∗∈Rm×dX_* \in \mathbb{R}^{m\times d}. Show that the posterior mean is X∗wX_*w with w=(X⊤X+σ2I)−1X⊤yw = (X^\top X + \sigma^2I)^{-1}X^\top y, the ridge solution with penalty σ2\sigma^2, and that the posterior covariance is σ2X∗(X⊤X+σ2I)−1X∗⊤\sigma^2X_*(X^\top X + \sigma^2I)^{-1}X_*^\top.

  9. ··

    With the RBF kernel, find the posterior mean and variance at a test point as ℓ→∞\ell \to \infty and as ℓ→0\ell \to 0 (for x∗x_* not equal to any training input).

  10. ···

    For a single test point, μ∗=k∗⊤Ky−1y\mu_* = k_*^\top K_y^{-1}y depends on ℓ\ell through both k∗k_* and KK. Compute dμ∗dℓ\dfrac{d\mu_*}{d\ell} in terms of α\alpha, ∂k∗/∂ℓ\partial k_*/\partial\ell and ∂K/∂ℓ\partial K/\partial\ell.

Answers

  1. f∗∣y∼N(K∗⊤Ky−1y, K∗∗−K∗⊤Ky−1K∗)f_* \mid y \sim \mathcal N\big(K_*^\top K_y^{-1}y,\ K_{**} - K_*^\top K_y^{-1}K_*\big)
  2. μ∗≈1.5514\mu_* \approx 1.5514 and Σ∗=1−0.9128≈0.0873\Sigma_* = 1 - 0.9128 \approx 0.0873
  3. μ(x∗)=∑iαik(xi,x∗)\mu(x_*) = \sum_i\alpha_ik(x_i, x_*); Var⁡(f(x∗)∣y)≤k∗∗\operatorname{Var}(f(x_*) \mid y) \le k_{**} everywhere; with σ2=0\sigma^2 = 0, μ(xj)=yj\mu(x_j) = y_j and the variance at xjx_j is 00
  4. ∂k/∂ℓ=k r2/ℓ3\partial k/\partial\ell = k\,r^2/\ell^3, ∂k/∂s2=k/s2\partial k/\partial s^2 = k/s^2, ∂k(xi,x∗)/∂x∗=−k(xi,x∗)(x∗−xi)/ℓ2\partial k(x_i, x_*)/\partial x_* = -k(x_i, x_*)(x_* - x_i)/\ell^2, and μ′(x∗)=−1ℓ2∑iαi(x∗−xi)k(xi,x∗)\mu'(x_*) = -\tfrac1{\ell^2}\sum_i\alpha_i(x_* - x_i)k(x_i, x_*)
  5. ∂log⁡p∂θ=12tr⁡((αα⊤−Ky−1)∂Ky∂θ)=12α⊤∂Ky∂θα−12tr⁡(Ky−1∂Ky∂θ)\dfrac{\partial\log p}{\partial\theta} = \tfrac12\operatorname{tr}\Big((\alpha\alpha^\top - K_y^{-1})\dfrac{\partial K_y}{\partial\theta}\Big) = \tfrac12\alpha^\top\dfrac{\partial K_y}{\partial\theta}\alpha - \tfrac12\operatorname{tr}\Big(K_y^{-1}\dfrac{\partial K_y}{\partial\theta}\Big)
  6. log⁡det⁡Ky=∑ilog⁡(λi+σ2)\log\det K_y = \sum_i\log(\lambda_i + \sigma^2), tr⁡Ky−1=∑i(λi+σ2)−1\operatorname{tr}K_y^{-1} = \sum_i(\lambda_i + \sigma^2)^{-1}, ∂log⁡p/∂σ2=12(∥α∥2−tr⁡Ky−1)\partial\log p/\partial\sigma^2 = \tfrac12(\|\alpha\|^2 - \operatorname{tr}K_y^{-1}), and at the optimum ∥α∥2=tr⁡Ky−1\|\alpha\|^2 = \operatorname{tr}K_y^{-1}
  7. μ∗′=μ∗+c (yn+1−μn+1)v+σ2\mu_*' = \mu_* + \dfrac{c\,(y_{n+1} - \mu_{n+1})}{v + \sigma^2}, Σ∗′=Σ∗−cc⊤v+σ2\Sigma_*' = \Sigma_* - \dfrac{cc^\top}{v + \sigma^2}, identical to the batch posterior on n+1n + 1 points
  8. μ∗=X∗(X⊤X+σ2I)−1X⊤y\mu_* = X_*(X^\top X + \sigma^2I)^{-1}X^\top y and Σ∗=σ2X∗(X⊤X+σ2I)−1X∗⊤\Sigma_* = \sigma^2X_*(X^\top X + \sigma^2I)^{-1}X_*^\top
  9. ℓ→∞\ell \to \infty: μ∗→ns2ns2+σ2yˉ\mu_* \to \dfrac{ns^2}{ns^2 + \sigma^2}\bar y and Σ∗→s2σ2ns2+σ2\Sigma_* \to \dfrac{s^2\sigma^2}{ns^2 + \sigma^2} everywhere; ℓ→0\ell \to 0: μ∗→0\mu_* \to 0 and Σ∗→s2\Sigma_* \to s^2 away from the data
  10. dμ∗dℓ=(∂k∗∂ℓ)⊤α−k∗⊤Ky−1∂K∂ℓα\dfrac{d\mu_*}{d\ell} = \Big(\dfrac{\partial k_*}{\partial\ell}\Big)^\top\alpha - k_*^\top K_y^{-1}\dfrac{\partial K}{\partial\ell}\alpha

Worked solutions

Problem 1

Write the joint distribution of (f∗,y)(f_*, y) and condition on yy to obtain the predictive distribution f∗∣y∼N(μ∗,Σ∗)f_* \mid y \sim \mathcal N(\mu_*, \Sigma_*) with μ∗=K∗⊤Ky−1y\mu_* = K_*^\top K_y^{-1}y and Σ∗=K∗∗−K∗⊤Ky−1K∗\Sigma_* = K_{**} - K_*^\top K_y^{-1}K_*.

  1. (f,f∗)(f, f_*) is jointly Gaussian with mean 00 and covariance (KK∗K∗⊤K∗∗)\begin{pmatrix}K & K_*\\ K_*^\top & K_{**}\end{pmatrix}.The prior applied to the n+mn + m inputs XX and X∗X_* together; the blocks are the kernel evaluated between the two sets.
  2. y=f+εy = f + \varepsilon is Gaussian with Cov⁡(y)=K+σ2I=Ky\operatorname{Cov}(y) = K + \sigma^2I = K_y.A sum of independent Gaussians is Gaussian and their covariances add; Cov⁡(ε)=σ2I\operatorname{Cov}(\varepsilon) = \sigma^2I.
  3. Cov⁡(f∗,y)=Cov⁡(f∗,f)+Cov⁡(f∗,ε)=K∗⊤+0\operatorname{Cov}(f_*, y) = \operatorname{Cov}(f_*, f) + \operatorname{Cov}(f_*, \varepsilon) = K_*^\top + 0.Covariance is bilinear, and ε\varepsilon is independent of everything in f∗f_*. The noise touches only the yy block.
  4. (f∗,y)∼N(0,(K∗∗K∗⊤K∗Ky))(f_*, y) \sim \mathcal N\Big(0, \begin{pmatrix}K_{**} & K_*^\top\\ K_* & K_y\end{pmatrix}\Big).Steps 1 to 3 assembled with xa=f∗x_a = f_* and xb=yx_b = y; (f∗,y)(f_*, y) is jointly Gaussian because it is a linear function of the Gaussian vector (f,f∗,ε)(f, f_*, \varepsilon) (the multivariate-Gaussian page, Problem 7).
  5. μ∗=0+K∗⊤Ky−1(y−0)\mu_* = 0 + K_*^\top K_y^{-1}(y - 0) and Σ∗=K∗∗−K∗⊤Ky−1K∗\Sigma_* = K_{**} - K_*^\top K_y^{-1}K_*.The conditioning formula with Σab=K∗⊤\Sigma_{ab} = K_*^\top, Σbb=Ky\Sigma_{bb} = K_y, Σaa=K∗∗\Sigma_{aa} = K_{**} and both means 00.
  6. f∗∣y∼N(K∗⊤Ky−1y, K∗∗−K∗⊤Ky−1K∗)f_* \mid y \sim \mathcal N\big(K_*^\top K_y^{-1}y,\ K_{**} - K_*^\top K_y^{-1}K_*\big)μ∗\mu_* is m×1m\times1 and Σ∗\Sigma_* is m×mm\times m. The noise appears once, in the matrix that is inverted, because the observations are noisy and the quantity predicted, f∗f_*, is not; to predict a noisy y∗y_* instead, add σ2I\sigma^2I to Σ∗\Sigma_* afterwards. The inverse of KyK_y is the whole cost, O(n3)O(n^3) once, after which each test point needs O(n)O(n) for the mean and O(n2)O(n^2) for the variance.

Problem 2

Two training points x1=0x_1 = 0, x2=1x_2 = 1 with y=(1,2)⊤y = (1, 2)^\top, RBF kernel with s2=1s^2 = 1, ℓ=1\ell = 1, noise σ2=0.1\sigma^2 = 0.1. Compute the predictive mean and variance at x∗=0.5x_* = 0.5 to four decimal places.

  1. K=(1e−1/2e−1/21)K = \begin{pmatrix}1 & e^{-1/2}\\ e^{-1/2} & 1\end{pmatrix}, and e−1/2≈0.6065e^{-1/2} \approx 0.6065, so Ky=(1.10.60650.60651.1)K_y = \begin{pmatrix}1.1 & 0.6065\\ 0.6065 & 1.1\end{pmatrix}.k(0,1)=exp⁡(−1/2)k(0, 1) = \exp(-1/2); k(x,x)=1k(x, x) = 1; add σ2=0.1\sigma^2 = 0.1 to the diagonal.
  2. k∗=(e−1/8,e−1/8)⊤≈(0.8825,0.8825)⊤k_* = (e^{-1/8}, e^{-1/8})^\top \approx (0.8825, 0.8825)^\top and k∗∗=1k_{**} = 1.Both training points are at distance 0.50.5 from x∗x_*: exp⁡(−0.25/2)=exp⁡(−1/8)\exp(-0.25/2) = \exp(-1/8).
  3. Ky−1=1a2−b2(a−b−ba)K_y^{-1} = \dfrac{1}{a^2 - b^2}\begin{pmatrix}a & -b\\ -b & a\end{pmatrix} with a=1.1a = 1.1, b=0.6065b = 0.6065, a2−b2≈1.21−0.3679=0.8421a^2 - b^2 \approx 1.21 - 0.3679 = 0.8421.The 2×22\times2 inverse of a symmetric matrix with equal diagonal: swap nothing, negate the off-diagonal, divide by the determinant.
  4. α=Ky−1y=10.8421(1.1−2(0.6065), −0.6065+2(1.1))⊤≈(−0.1343, 1.8922)⊤\alpha = K_y^{-1}y = \dfrac{1}{0.8421}\big(1.1 - 2(0.6065),\ -0.6065 + 2(1.1)\big)^\top \approx (-0.1343,\ 1.8922)^\top.Multiply out: 1.1−1.2131=−0.11311.1 - 1.2131 = -0.1131 and −0.6065+2.2=1.5935-0.6065 + 2.2 = 1.5935, each divided by 0.84210.8421.
  5. μ∗=k∗⊤α=0.8825 (−0.1343+1.8922)=0.8825×1.7580≈1.5514\mu_* = k_*^\top\alpha = 0.8825\,(-0.1343 + 1.8922) = 0.8825\times1.7580 \approx 1.5514.Both entries of k∗k_* are equal, so the dot product is 0.88250.8825 times the sum of α\alpha.
  6. k∗⊤Ky−1k∗=0.88252×(sum of all entries of Ky−1)=0.7788×2(a−b)a2−b2=0.7788×0.98700.8421≈0.7788×1.1720≈0.9128k_*^\top K_y^{-1}k_* = 0.8825^2\times(\text{sum of all entries of } K_y^{-1}) = 0.7788\times\dfrac{2(a - b)}{a^2 - b^2} = 0.7788\times\dfrac{0.9870}{0.8421} \approx 0.7788\times1.1720 \approx 0.9128.With k∗=c1k_* = c\mathbf{1}, k∗⊤Mk∗=c21⊤M1k_*^\top Mk_* = c^2\mathbf{1}^\top M\mathbf{1}, the sum of the entries; for step 3's inverse that is 2(a−b)/(a2−b2)2(a - b)/(a^2 - b^2).
  7. μ∗≈1.5514\mu_* \approx 1.5514 and Σ∗=1−0.9128≈0.0873\Sigma_* = 1 - 0.9128 \approx 0.0873A little above the midpoint 1.51.5 of the two observations, because the second observation has the larger weight α2\alpha_2 and the first a small negative one: α\alpha is not a set of averaging weights but the solution of Kyα=yK_y\alpha = y. The variance has fallen from the prior 11 to 0.08730.0873, near the noise level 0.10.1, as a point half a lengthscale from two observations should.

Problem 3

Show that the posterior mean at any x∗x_* is μ(x∗)=∑i=1nαik(xi,x∗)\mu(x_*) = \sum_{i=1}^n\alpha_ik(x_i, x_*) with α=Ky−1y\alpha = K_y^{-1}y not depending on x∗x_*; that the posterior variance never exceeds the prior variance; and that with σ2=0\sigma^2 = 0 the posterior mean passes through every observation with zero variance there.

  1. μ(x∗)=k∗⊤Ky−1y=k∗⊤α=∑i(k∗)iαi=∑iαik(xi,x∗)\mu(x_*) = k_*^\top K_y^{-1}y = k_*^\top\alpha = \sum_i(k_*)_i\alpha_i = \sum_i\alpha_ik(x_i, x_*).Problem 1 for one test point; α\alpha depends on the training data only, so it is computed once.
  2. K∗⊤Ky−1K∗=(Ky−1/2K∗)⊤(Ky−1/2K∗)K_*^\top K_y^{-1}K_* = (K_y^{-1/2}K_*)^\top(K_y^{-1/2}K_*) is positive semidefinite.KyK_y is positive definite (a positive semidefinite KK plus σ2I\sigma^2I), so it has a symmetric square root and Ky−1=Ky−1/2Ky−1/2K_y^{-1} = K_y^{-1/2}K_y^{-1/2}; a matrix of the form B⊤BB^\top B is positive semidefinite.
  3. So Σ∗=K∗∗−K∗⊤Ky−1K∗⪯K∗∗\Sigma_* = K_{**} - K_*^\top K_y^{-1}K_* \preceq K_{**}, and in particular Var⁡(f(x∗)∣y)=k∗∗−k∗⊤Ky−1k∗≤k∗∗\operatorname{Var}(f(x_*) \mid y) = k_{**} - k_*^\top K_y^{-1}k_* \le k_{**}.Subtracting a positive semidefinite matrix can only lower each diagonal entry. The data cannot increase uncertainty about ff anywhere.
  4. With σ2=0\sigma^2 = 0 and x∗=xjx_* = x_j: k∗=Kejk_* = Ke_j (column jj of KK) and k∗∗=Kjjk_{**} = K_{jj}.k(xi,xj)k(x_i, x_j) is entry (i,j)(i, j) of KK.
  5. μ(xj)=ej⊤K⊤K−1y=ej⊤y=yj\mu(x_j) = e_j^\top K^\top K^{-1}y = e_j^\top y = y_j and Var⁡=Kjj−ej⊤KK−1Kej=Kjj−Kjj=0\operatorname{Var} = K_{jj} - e_j^\top KK^{-1}Ke_j = K_{jj} - K_{jj} = 0.K⊤=KK^\top = K and KK−1=IKK^{-1} = I, assuming KK is invertible.
  6. μ(x∗)=∑iαik(xi,x∗)\mu(x_*) = \sum_i\alpha_ik(x_i, x_*); Var⁡(f(x∗)∣y)≤k∗∗\operatorname{Var}(f(x_*) \mid y) \le k_{**} everywhere; with σ2=0\sigma^2 = 0, μ(xj)=yj\mu(x_j) = y_j and the variance at xjx_j is 00The posterior mean is a fixed linear combination of nn copies of the kernel, one centred on each training input, so it inherits the kernel's shape: smooth for the RBF, returning to the prior mean 00 far from the data (Problem 9). The variance depends on the inputs XX and x∗x_* but not on yy: where the data are is what reduces uncertainty, not what they say. Noise-free conditioning interpolates; with σ2>0\sigma^2 > 0 the mean passes near, not through, the observations, and the variance at a training input is positive.

Problem 4

For the RBF kernel with r=∥x−x′∥r = \|x - x'\|, compute ∂k/∂ℓ\partial k/\partial\ell, ∂k/∂s2\partial k/\partial s^2 and, for a scalar test input, ∂k(xi,x∗)/∂x∗\partial k(x_i, x_*)/\partial x_*. Use the last to write the slope μ′(x∗)\mu'(x_*) of the posterior mean.

  1. ∂∂ℓ(−r22ℓ2)=−r22⋅(−2)ℓ−3=r2ℓ3\dfrac{\partial}{\partial\ell}\Big(-\dfrac{r^2}{2\ell^2}\Big) = -\dfrac{r^2}{2}\cdot(-2)\ell^{-3} = \dfrac{r^2}{\ell^3}.ddℓℓ−2=−2ℓ−3\tfrac{d}{d\ell}\ell^{-2} = -2\ell^{-3}.
  2. ∂k∂ℓ=k⋅r2ℓ3\dfrac{\partial k}{\partial\ell} = k\cdot\dfrac{r^2}{\ell^3}.Chain rule through the exponential: ddℓeu(ℓ)=euu′\tfrac{d}{d\ell}e^{u(\ell)} = e^{u}u', and s2eu=ks^2e^{u} = k.
  3. ∂k∂s2=exp⁡(−r22ℓ2)=ks2\dfrac{\partial k}{\partial s^2} = \exp\Big(-\dfrac{r^2}{2\ell^2}\Big) = \dfrac{k}{s^2}.kk is linear in s2s^2.
  4. ∂k(xi,x∗)∂x∗=k(xi,x∗)⋅∂∂x∗(−(x∗−xi)22ℓ2)=−k(xi,x∗) x∗−xiℓ2\dfrac{\partial k(x_i, x_*)}{\partial x_*} = k(x_i, x_*)\cdot\dfrac{\partial}{\partial x_*}\Big(-\dfrac{(x_* - x_i)^2}{2\ell^2}\Big) = -k(x_i, x_*)\,\dfrac{x_* - x_i}{\ell^2}.Chain rule; ddx∗(x∗−xi)2=2(x∗−xi)\tfrac{d}{dx_*}(x_* - x_i)^2 = 2(x_* - x_i) and the 22 cancels.
  5. μ′(x∗)=∑iαi∂k(xi,x∗)∂x∗=−1ℓ2∑iαi(x∗−xi) k(xi,x∗)\mu'(x_*) = \sum_i\alpha_i\dfrac{\partial k(x_i, x_*)}{\partial x_*} = -\dfrac{1}{\ell^2}\sum_i\alpha_i(x_* - x_i)\,k(x_i, x_*).Differentiate Problem 3's representer form term by term; α\alpha does not depend on x∗x_*.
  6. ∂k/∂ℓ=k r2/ℓ3\partial k/\partial\ell = k\,r^2/\ell^3, ∂k/∂s2=k/s2\partial k/\partial s^2 = k/s^2, ∂k(xi,x∗)/∂x∗=−k(xi,x∗)(x∗−xi)/ℓ2\partial k(x_i, x_*)/\partial x_* = -k(x_i, x_*)(x_* - x_i)/\ell^2, and μ′(x∗)=−1ℓ2∑iαi(x∗−xi)k(xi,x∗)\mu'(x_*) = -\tfrac1{\ell^2}\sum_i\alpha_i(x_* - x_i)k(x_i, x_*)The lengthscale derivative is largest at r≈2 ℓr \approx \sqrt2\,\ell and vanishes at r=0r = 0 (the diagonal of KK never depends on ℓ\ell) and at large rr; the signal-variance derivative is the kernel with s2s^2 stripped out. The slope formula is how a Gaussian process gives derivatives of its mean function for free, and the same ∂k/∂x∗\partial k/\partial x_* is the cross-covariance between ff and f′f'.

Problem 5

The log marginal likelihood is log⁡p(y∣X,θ)=−12y⊤Ky−1y−12log⁡det⁡Ky−n2log⁡2π\log p(y \mid X, \theta) = -\tfrac12y^\top K_y^{-1}y - \tfrac12\log\det K_y - \tfrac n2\log2\pi. For a hyperparameter θ\theta with ∂Ky/∂θ\partial K_y/\partial\theta known, show that ∂∂θlog⁡p(y∣X,θ)=12tr⁡((αα⊤−Ky−1)∂Ky∂θ)\dfrac{\partial}{\partial\theta}\log p(y \mid X, \theta) = \tfrac12\operatorname{tr}\Big((\alpha\alpha^\top - K_y^{-1})\dfrac{\partial K_y}{\partial\theta}\Big), and write the two terms for θ=ℓ\theta = \ell.

  1. ∂∂θ(y⊤Ky−1y)=y⊤∂Ky−1∂θy=−y⊤Ky−1∂Ky∂θKy−1y\dfrac{\partial}{\partial\theta}\big(y^\top K_y^{-1}y\big) = y^\top\dfrac{\partial K_y^{-1}}{\partial\theta}y = -y^\top K_y^{-1}\dfrac{\partial K_y}{\partial\theta}K_y^{-1}y.yy is data, so only the inverse moves; the previous page's Problem 4 gives ∂(A−1)/∂θ=−A−1(∂A/∂θ)A−1\partial(A^{-1})/\partial\theta = -A^{-1}(\partial A/\partial\theta)A^{-1}.
  2. =−α⊤∂Ky∂θα=−tr⁡(αα⊤∂Ky∂θ)= -\alpha^\top\dfrac{\partial K_y}{\partial\theta}\alpha = -\operatorname{tr}\Big(\alpha\alpha^\top\dfrac{\partial K_y}{\partial\theta}\Big).Ky−1y=αK_y^{-1}y = \alpha and y⊤Ky−1=α⊤y^\top K_y^{-1} = \alpha^\top since KyK_y is symmetric; a scalar is its own trace and the cyclic property moves α\alpha to the back.
  3. ∂∂θlog⁡det⁡Ky=tr⁡(Ky−1∂Ky∂θ)\dfrac{\partial}{\partial\theta}\log\det K_y = \operatorname{tr}\Big(K_y^{-1}\dfrac{\partial K_y}{\partial\theta}\Big).The previous page's Problem 4 again, for log⁡det⁡\log\det.
  4. ∂∂θlog⁡p=12tr⁡(αα⊤∂Ky∂θ)−12tr⁡(Ky−1∂Ky∂θ)=12tr⁡((αα⊤−Ky−1)∂Ky∂θ)\dfrac{\partial}{\partial\theta}\log p = \tfrac12\operatorname{tr}\Big(\alpha\alpha^\top\dfrac{\partial K_y}{\partial\theta}\Big) - \tfrac12\operatorname{tr}\Big(K_y^{-1}\dfrac{\partial K_y}{\partial\theta}\Big) = \tfrac12\operatorname{tr}\Big((\alpha\alpha^\top - K_y^{-1})\dfrac{\partial K_y}{\partial\theta}\Big).−12-\tfrac12 times step 2 and −12-\tfrac12 times step 3; the constant term has no θ\theta; the trace is linear.
  5. For θ=ℓ\theta = \ell: ∂Ky∂ℓ=∂K∂ℓ\dfrac{\partial K_y}{\partial\ell} = \dfrac{\partial K}{\partial\ell} with entries Kij (xi−xj)2/ℓ3K_{ij}\,(x_i - x_j)^2/\ell^3, and the two terms are 12α⊤∂K∂ℓα\tfrac12\alpha^\top\dfrac{\partial K}{\partial\ell}\alpha and −12tr⁡(Ky−1∂K∂ℓ)-\tfrac12\operatorname{tr}\Big(K_y^{-1}\dfrac{\partial K}{\partial\ell}\Big).σ2I\sigma^2I does not depend on ℓ\ell; Problem 4 entry by entry.
  6. ∂log⁡p∂θ=12tr⁡((αα⊤−Ky−1)∂Ky∂θ)=12α⊤∂Ky∂θα−12tr⁡(Ky−1∂Ky∂θ)\dfrac{\partial\log p}{\partial\theta} = \tfrac12\operatorname{tr}\Big((\alpha\alpha^\top - K_y^{-1})\dfrac{\partial K_y}{\partial\theta}\Big) = \tfrac12\alpha^\top\dfrac{\partial K_y}{\partial\theta}\alpha - \tfrac12\operatorname{tr}\Big(K_y^{-1}\dfrac{\partial K_y}{\partial\theta}\Big)The first term is the derivative of the data fit −12y⊤Ky−1y-\tfrac12y^\top K_y^{-1}y and the second of the complexity penalty −12log⁡det⁡Ky-\tfrac12\log\det K_y; the marginal likelihood trades them off automatically, which is why its maximum is a usable choice of hyperparameters without a validation set. The same ∂Ky/∂θ\partial K_y/\partial\theta with K/s2K/s^2 (Problem 4) gives the signal-variance gradient, and with II the noise gradient (Problem 6). Each gradient costs one n×nn\times n product once Ky−1K_y^{-1} and α\alpha are known.

Problem 6

Let KK have eigenvalues λ1,…,λn≥0\lambda_1, \dots, \lambda_n \ge 0. Show that log⁡det⁡Ky=∑ilog⁡(λi+σ2)\log\det K_y = \sum_i\log(\lambda_i + \sigma^2) and tr⁡(Ky−1)=∑i1/(λi+σ2)\operatorname{tr}(K_y^{-1}) = \sum_i1/(\lambda_i + \sigma^2), hence that ∂log⁡p∂σ2=12(∥α∥2−tr⁡Ky−1)\dfrac{\partial\log p}{\partial\sigma^2} = \tfrac12\big(\|\alpha\|^2 - \operatorname{tr}K_y^{-1}\big), and state the condition at a maximum over σ2\sigma^2.

  1. K=QΛQ⊤K = Q\Lambda Q^\top with QQ orthogonal, so Ky=Q(Λ+σ2I)Q⊤K_y = Q(\Lambda + \sigma^2I)Q^\top.KK is symmetric; I=QQ⊤I = QQ^\top puts the noise into the same basis, as on the previous page's Problem 4.
  2. det⁡Ky=det⁡(Λ+σ2I)=∏i(λi+σ2)\det K_y = \det(\Lambda + \sigma^2I) = \prod_i(\lambda_i + \sigma^2), so log⁡det⁡Ky=∑ilog⁡(λi+σ2)\log\det K_y = \sum_i\log(\lambda_i + \sigma^2).det⁡(QDQ⊤)=det⁡Qdet⁡Ddet⁡Q⊤=det⁡D\det(QDQ^\top) = \det Q\det D\det Q^\top = \det D since det⁡Q=±1\det Q = \pm1; a diagonal determinant is the product of the diagonal.
  3. Ky−1=Q(Λ+σ2I)−1Q⊤K_y^{-1} = Q(\Lambda + \sigma^2I)^{-1}Q^\top and tr⁡Ky−1=∑i1λi+σ2\operatorname{tr}K_y^{-1} = \sum_i\dfrac{1}{\lambda_i + \sigma^2}.Inverse of QDQ⊤QDQ^\top is QD−1Q⊤QD^{-1}Q^\top; the trace is cyclic, tr⁡(QD−1Q⊤)=tr⁡(D−1)\operatorname{tr}(QD^{-1}Q^\top) = \operatorname{tr}(D^{-1}).
  4. ∂Ky∂σ2=I\dfrac{\partial K_y}{\partial\sigma^2} = I, so Problem 5 gives ∂log⁡p∂σ2=12α⊤α−12tr⁡Ky−1\dfrac{\partial\log p}{\partial\sigma^2} = \tfrac12\alpha^\top\alpha - \tfrac12\operatorname{tr}K_y^{-1}.Substitute ∂Ky/∂θ=I\partial K_y/\partial\theta = I into both terms.
  5. At a maximum over σ2\sigma^2 (an interior one), ∥α∥2=tr⁡Ky−1=∑i1λi+σ2\|\alpha\|^2 = \operatorname{tr}K_y^{-1} = \sum_i\dfrac{1}{\lambda_i + \sigma^2}.Set the derivative to zero.
  6. log⁡det⁡Ky=∑ilog⁡(λi+σ2)\log\det K_y = \sum_i\log(\lambda_i + \sigma^2), tr⁡Ky−1=∑i(λi+σ2)−1\operatorname{tr}K_y^{-1} = \sum_i(\lambda_i + \sigma^2)^{-1}, ∂log⁡p/∂σ2=12(∥α∥2−tr⁡Ky−1)\partial\log p/\partial\sigma^2 = \tfrac12(\|\alpha\|^2 - \operatorname{tr}K_y^{-1}), and at the optimum ∥α∥2=tr⁡Ky−1\|\alpha\|^2 = \operatorname{tr}K_y^{-1}∥α∥2=y⊤Ky−2y\|\alpha\|^2 = y^\top K_y^{-2}y measures how hard the model is working to explain yy and tr⁡Ky−1\operatorname{tr}K_y^{-1} is what that quantity would be on average for data drawn from the model itself (E[y⊤Ky−2y]=tr⁡(Ky−2Ky)\mathbb E[y^\top K_y^{-2}y] = \operatorname{tr}(K_y^{-2}K_y) for y∼N(0,Ky)y \sim \mathcal N(0, K_y)); the noise is tuned until the two agree. The sum form shows the log-determinant is dominated by the small eigenvalues of KK, which σ2\sigma^2 lifts: the noise is also the regulariser that keeps KyK_y invertible, and the "jitter" added to a noise-free kernel matrix is exactly a small σ2\sigma^2.

Problem 7

After conditioning on nn observations, the posterior over (f∗,f(xn+1))(f_*, f(x_{n+1})) is Gaussian with mean (μ∗,μn+1)(\mu_*, \mu_{n+1}), covariance Σ∗\Sigma_* for f∗f_*, cross-covariance c=Cov⁡(f∗,f(xn+1))c = \operatorname{Cov}(f_*, f(x_{n+1})) and variance v=Var⁡(f(xn+1))v = \operatorname{Var}(f(x_{n+1})). A new observation yn+1=f(xn+1)+εn+1y_{n+1} = f(x_{n+1}) + \varepsilon_{n+1} arrives. Show that the updated posterior is μ∗′=μ∗+c (yn+1−μn+1)v+σ2\mu_*' = \mu_* + \dfrac{c\,(y_{n+1} - \mu_{n+1})}{v + \sigma^2} and Σ∗′=Σ∗−cc⊤v+σ2\Sigma_*' = \Sigma_* - \dfrac{cc^\top}{v + \sigma^2}, and say why this equals the batch posterior on all n+1n + 1 points.

  1. Given the first nn observations, (f∗,yn+1)(f_*, y_{n+1}) is jointly Gaussian with mean (μ∗,μn+1)(\mu_*, \mu_{n+1}), Cov⁡(f∗)=Σ∗\operatorname{Cov}(f_*) = \Sigma_*, Cov⁡(f∗,yn+1)=c\operatorname{Cov}(f_*, y_{n+1}) = c and Var⁡(yn+1)=v+σ2\operatorname{Var}(y_{n+1}) = v + \sigma^2.The posterior after nn points is a Gaussian (Problem 1), and yn+1y_{n+1} adds independent noise of variance σ2\sigma^2 to f(xn+1)f(x_{n+1}), which changes its variance but not its covariance with f∗f_* (as in Problem 1, step 3).
  2. Condition on yn+1y_{n+1}: μ∗′=μ∗+c (v+σ2)−1(yn+1−μn+1)\mu_*' = \mu_* + c\,(v + \sigma^2)^{-1}(y_{n+1} - \mu_{n+1}).The conditioning formula with Σab=c\Sigma_{ab} = c (m×1m\times1) and Σbb=v+σ2\Sigma_{bb} = v + \sigma^2, a scalar.
  3. Σ∗′=Σ∗−c (v+σ2)−1c⊤\Sigma_*' = \Sigma_* - c\,(v + \sigma^2)^{-1}c^\top.The same formula's covariance; cc⊤cc^\top is m×mm\times m of rank one.
  4. Conditioning on y1,…,yny_1, \dots, y_n and then on yn+1y_{n+1} is conditioning on all n+1n + 1 at once.p(f∗∣y1:n+1)∝p(f∗,yn+1∣y1:n)p(f_* \mid y_{1:n+1}) \propto p(f_*, y_{n+1} \mid y_{1:n}) by the definition of conditional probability, applied to the posterior after nn points as the new "prior"; the Gaussian family is closed under conditioning, so the sequential answer is the batch answer.
  5. μ∗′=μ∗+c (yn+1−μn+1)v+σ2\mu_*' = \mu_* + \dfrac{c\,(y_{n+1} - \mu_{n+1})}{v + \sigma^2}, Σ∗′=Σ∗−cc⊤v+σ2\Sigma_*' = \Sigma_* - \dfrac{cc^\top}{v + \sigma^2}, identical to the batch posterior on n+1n + 1 pointsThe update is the residual yn+1−μn+1y_{n+1} - \mu_{n+1} (what the new observation says beyond what was predicted) spread over the test points in proportion to their posterior covariance with the new input, divided by the predicted variance of the observation; the variance drops by a rank-one amount that again does not involve yn+1y_{n+1}. In the KyK_y picture this is a rank-one update of the inverse, the previous page's Sherman–Morrison formula, and it is the step a Kalman filter takes for one scalar measurement.

Problem 8

Take the linear kernel k(x,x′)=x⊤x′k(x, x') = x^\top x' with inputs as rows of X∈Rn×dX \in \mathbb{R}^{n\times d} and X∗∈Rm×dX_* \in \mathbb{R}^{m\times d}. Show that the posterior mean is X∗wX_*w with w=(X⊤X+σ2I)−1X⊤yw = (X^\top X + \sigma^2I)^{-1}X^\top y, the ridge solution with penalty σ2\sigma^2, and that the posterior covariance is σ2X∗(X⊤X+σ2I)−1X∗⊤\sigma^2X_*(X^\top X + \sigma^2I)^{-1}X_*^\top.

  1. K=XX⊤K = XX^\top, K∗=XX∗⊤K_* = XX_*^\top, K∗∗=X∗X∗⊤K_{**} = X_*X_*^\top.Entry (i,j)(i, j) of each is the dot product of the corresponding rows.
  2. X⊤(XX⊤+σ2I)=(X⊤X+σ2I)X⊤X^\top(XX^\top + \sigma^2I) = (X^\top X + \sigma^2I)X^\top.Multiply out: both sides are X⊤XX⊤+σ2X⊤X^\top XX^\top + \sigma^2X^\top.
  3. (X⊤X+σ2I)−1X⊤=X⊤(XX⊤+σ2I)−1(X^\top X + \sigma^2I)^{-1}X^\top = X^\top(XX^\top + \sigma^2I)^{-1}.Multiply step 2 by (X⊤X+σ2I)−1(X^\top X + \sigma^2I)^{-1} on the left and (XX⊤+σ2I)−1(XX^\top + \sigma^2I)^{-1} on the right; both matrices are positive definite, hence invertible. This is the push-through identity: a d×dd\times d inverse on one side, an n×nn\times n inverse on the other.
  4. μ∗=K∗⊤Ky−1y=X∗X⊤(XX⊤+σ2I)−1y=X∗(X⊤X+σ2I)−1X⊤y=X∗w\mu_* = K_*^\top K_y^{-1}y = X_*X^\top(XX^\top + \sigma^2I)^{-1}y = X_*(X^\top X + \sigma^2I)^{-1}X^\top y = X_*w.Problem 1, then step 3.
  5. Σ∗=X∗X∗⊤−X∗X⊤(XX⊤+σ2I)−1XX∗⊤=X∗(I−(X⊤X+σ2I)−1X⊤X)X∗⊤\Sigma_* = X_*X_*^\top - X_*X^\top(XX^\top + \sigma^2I)^{-1}XX_*^\top = X_*\big(I - (X^\top X + \sigma^2I)^{-1}X^\top X\big)X_*^\top.Step 3 applied inside the second term, then factor X∗X_* and X∗⊤X_*^\top out.
  6. I−(X⊤X+σ2I)−1X⊤X=(X⊤X+σ2I)−1(X⊤X+σ2I−X⊤X)=σ2(X⊤X+σ2I)−1I - (X^\top X + \sigma^2I)^{-1}X^\top X = (X^\top X + \sigma^2I)^{-1}\big(X^\top X + \sigma^2I - X^\top X\big) = \sigma^2(X^\top X + \sigma^2I)^{-1}.Write II as (X⊤X+σ2I)−1(X⊤X+σ2I)(X^\top X + \sigma^2I)^{-1}(X^\top X + \sigma^2I) and subtract.
  7. μ∗=X∗(X⊤X+σ2I)−1X⊤y\mu_* = X_*(X^\top X + \sigma^2I)^{-1}X^\top y and Σ∗=σ2X∗(X⊤X+σ2I)−1X∗⊤\Sigma_* = \sigma^2X_*(X^\top X + \sigma^2I)^{-1}X_*^\topGaussian process regression with a linear kernel is Bayesian linear regression with prior w∼N(0,I)w \sim \mathcal N(0, I): the mean is the ridge fit with λ=σ2\lambda = \sigma^2 (the regression page, in the units of 12∥Xw−y∥2+λ2∥w∥2\tfrac12\|Xw - y\|^2 + \tfrac\lambda2\|w\|^2), and the covariance is X∗X_* times the posterior covariance of ww, which is σ2(X⊤X+σ2I)−1\sigma^2(X^\top X + \sigma^2I)^{-1}. The two sides of step 3 are the two ways to compute it: in weight space (d×dd\times d) when d<nd < n, in function space (n×nn\times n) when the kernel has no finite dd, which is the point of kernels.

Problem 9

With the RBF kernel, find the posterior mean and variance at a test point as ℓ→∞\ell \to \infty and as ℓ→0\ell \to 0 (for x∗x_* not equal to any training input).

  1. As ℓ→∞\ell \to \infty, k(x,x′)→s2k(x, x') \to s^2 for every pair, so K→s211⊤K \to s^2\mathbf{1}\mathbf{1}^\top, k∗→s21k_* \to s^2\mathbf{1} and k∗∗=s2k_{**} = s^2.exp⁡(−r2/2ℓ2)→1\exp(-r^2/2\ell^2) \to 1.
  2. (s211⊤+σ2I)1=(ns2+σ2)1(s^2\mathbf{1}\mathbf{1}^\top + \sigma^2I)\mathbf{1} = (ns^2 + \sigma^2)\mathbf{1}, so Ky−11=1ns2+σ2K_y^{-1}\mathbf{1} = \dfrac{\mathbf{1}}{ns^2 + \sigma^2}.1⊤1=n\mathbf{1}^\top\mathbf{1} = n: 1\mathbf{1} is an eigenvector of KyK_y with eigenvalue ns2+σ2ns^2 + \sigma^2, and the inverse has the reciprocal eigenvalue on the same vector.
  3. μ∗=s21⊤Ky−1y=s21⊤yns2+σ2=ns2ns2+σ2 yˉ\mu_* = s^2\mathbf{1}^\top K_y^{-1}y = \dfrac{s^2\mathbf{1}^\top y}{ns^2 + \sigma^2} = \dfrac{ns^2}{ns^2 + \sigma^2}\,\bar y.1⊤Ky−1=(Ky−11)⊤\mathbf{1}^\top K_y^{-1} = (K_y^{-1}\mathbf{1})^\top by symmetry, and 1⊤y=nyˉ\mathbf{1}^\top y = n\bar y.
  4. Σ∗=s2−s41⊤Ky−11=s2−ns4ns2+σ2=s2σ2ns2+σ2\Sigma_* = s^2 - s^4\mathbf{1}^\top K_y^{-1}\mathbf{1} = s^2 - \dfrac{ns^4}{ns^2 + \sigma^2} = \dfrac{s^2\sigma^2}{ns^2 + \sigma^2}.Step 2 again; put over the common denominator.
  5. As ℓ→0\ell \to 0, k(x,x′)→0k(x, x') \to 0 for x≠x′x \neq x', so K→s2IK \to s^2I, k∗→0k_* \to 0 and k∗∗=s2k_{**} = s^2.exp⁡(−r2/2ℓ2)→0\exp(-r^2/2\ell^2) \to 0 for r>0r > 0, and is 11 at r=0r = 0.
  6. μ∗→0\mu_* \to 0 and Σ∗→s2\Sigma_* \to s^2.k∗=0k_* = 0 kills both the mean and the correction to the variance.
  7. ℓ→∞\ell \to \infty: μ∗→ns2ns2+σ2yˉ\mu_* \to \dfrac{ns^2}{ns^2 + \sigma^2}\bar y and Σ∗→s2σ2ns2+σ2\Sigma_* \to \dfrac{s^2\sigma^2}{ns^2 + \sigma^2} everywhere; ℓ→0\ell \to 0: μ∗→0\mu_* \to 0 and Σ∗→s2\Sigma_* \to s^2 away from the dataA very long lengthscale says the function is constant, so the model fits one number, the sample mean shrunk towards the prior mean 00 by the factor ns2/(ns2+σ2)ns^2/(ns^2 + \sigma^2), with the variance of a mean of nn noisy observations; a very short one says the observations tell nothing about any other point, so the posterior is the prior except on the data. The marginal likelihood of Problem 5 picks the ℓ\ell between these at which the data look most like a draw from the prior.

Problem 10

For a single test point, μ∗=k∗⊤Ky−1y\mu_* = k_*^\top K_y^{-1}y depends on ℓ\ell through both k∗k_* and KK. Compute dμ∗dℓ\dfrac{d\mu_*}{d\ell} in terms of α\alpha, ∂k∗/∂ℓ\partial k_*/\partial\ell and ∂K/∂ℓ\partial K/\partial\ell.

  1. dμ∗dℓ=(∂k∗∂ℓ)⊤Ky−1y+k∗⊤∂Ky−1∂ℓy\dfrac{d\mu_*}{d\ell} = \Big(\dfrac{\partial k_*}{\partial\ell}\Big)^\top K_y^{-1}y + k_*^\top\dfrac{\partial K_y^{-1}}{\partial\ell}y.Product rule on the three factors; yy is constant.
  2. ∂Ky−1∂ℓ=−Ky−1∂K∂ℓKy−1\dfrac{\partial K_y^{-1}}{\partial\ell} = -K_y^{-1}\dfrac{\partial K}{\partial\ell}K_y^{-1}.The previous page's Problem 4; ∂Ky/∂ℓ=∂K/∂ℓ\partial K_y/\partial\ell = \partial K/\partial\ell because the noise term has no ℓ\ell.
  3. dμ∗dℓ=(∂k∗∂ℓ)⊤α−k∗⊤Ky−1∂K∂ℓα\dfrac{d\mu_*}{d\ell} = \Big(\dfrac{\partial k_*}{\partial\ell}\Big)^\top\alpha - k_*^\top K_y^{-1}\dfrac{\partial K}{\partial\ell}\alpha.Ky−1y=αK_y^{-1}y = \alpha in both terms.
  4. The entries are (∂k∗/∂ℓ)i=k(xi,x∗)(xi−x∗)2/ℓ3(\partial k_*/\partial\ell)_i = k(x_i, x_*)(x_i - x_*)^2/\ell^3 and (∂K/∂ℓ)ij=Kij(xi−xj)2/ℓ3(\partial K/\partial\ell)_{ij} = K_{ij}(x_i - x_j)^2/\ell^3.Problem 4.
  5. dμ∗dℓ=(∂k∗∂ℓ)⊤α−k∗⊤Ky−1∂K∂ℓα\dfrac{d\mu_*}{d\ell} = \Big(\dfrac{\partial k_*}{\partial\ell}\Big)^\top\alpha - k_*^\top K_y^{-1}\dfrac{\partial K}{\partial\ell}\alphaThe first term is how the prediction changes because the test point's correlations with the data change; the second, through the inverse, is how it changes because the data's correlations with each other change, which reweights α\alpha. The second term needs one more solve with KyK_y (of ∂K∂ℓα\tfrac{\partial K}{\partial\ell}\alpha, or of k∗k_* by symmetry), and dropping it is the last mistake below. The same two-term structure gives the derivative of any prediction with respect to any hyperparameter, which is what makes the predictions themselves differentiable for downstream use.

Where this goes wrong

1. Inverting K instead of K + σ²I

The conditioning formula is written with the prior covariance, and KK is the prior covariance.

  1. (f,f∗)(f, f_*) has prior covariance with blocks KK, K∗K_*, K∗∗K_{**}Right so far: Problem 1, step 1.
  2. “Condition f∗f_* on the observed values using the prior covariance of what was observed.”The slip that causes the mistake: what was observed is y=f+εy = f + \varepsilon, whose covariance is K+σ2IK + \sigma^2I, not KK.
  3. μ∗=K∗⊤K−1y\mu_* = K_*^\top K^{-1}yThis is the σ2=0\sigma^2 = 0 posterior (Problem 3), which interpolates every noisy observation exactly, so the mean wiggles through the noise and the variance is 00 at the data, when the data are known to be noisy. It also needs KK invertible, which a kernel matrix with two nearby inputs is not, numerically; σ2I\sigma^2I is what makes KyK_y well conditioned (Problem 6). In Problem 2 the wrong mean is 1.64801.6480 in place of 1.55141.5514 and the variance 0.03050.0305 in place of 0.08730.0873.

2. Noise added to the cross-covariance

A convenient implementation folds the noise into the kernel as ky(x,x′)=k(x,x′)+σ2[x=x′]k_y(x, x') = k(x, x') + \sigma^2[x = x'] and uses kyk_y for every block.

  1. Ky=K+σ2IK_y = K + \sigma^2I is correct for the training blockRight so far: the diagonal of KK is where the observation noise lives.
  2. “Use the same noisy kernel kyk_y to build K∗K_*, so a test point equal to a training input gets the same treatment.”The shortcut that causes the mistake: the noise is a property of the observations, not of the function, and f∗f_* is a function value; Cov⁡(f(x∗),yi)=k(x∗,xi)\operatorname{Cov}(f(x_*), y_i) = k(x_*, x_i) with no σ2\sigma^2 even when x∗=xix_* = x_i (Problem 1, step 3).
  3. (K∗)ij=k(xi,xj∗)+σ2(K_*)_{ij} = k(x_i, x_j^*) + \sigma^2 when xj∗=xix_j^* = x_iAt a test point that coincides with a training input the cross-covariance is too large by σ2\sigma^2, so the posterior mean leans towards that observation's noisy value and the posterior variance comes out too small (it can go negative). The bug is invisible at test points away from the data and appears exactly where predictions are compared with the training targets. The noise belongs in KyK_y only, and in Σ∗+σ2I\Sigma_* + \sigma^2I if a noisy y∗y_* is being predicted.

3. Marginal likelihood gradient without the log-determinant

The data-fit term is the one with yy in it, and its derivative is the one that looks like a gradient of a loss.

  1. ∂∂θ(−12y⊤Ky−1y)=12α⊤∂Ky∂θα\dfrac{\partial}{\partial\theta}\Big(-\tfrac12y^\top K_y^{-1}y\Big) = \tfrac12\alpha^\top\dfrac{\partial K_y}{\partial\theta}\alphaRight so far: Problem 5, steps 1 and 2.
  2. “The log⁡det⁡\log\det term is a normalising constant.”The assumption that causes the mistake: it is constant in yy, not in θ\theta; it is how the marginal likelihood charges for a kernel that could explain anything.
  3. ∂log⁡p∂θ=12α⊤∂Ky∂θα\dfrac{\partial\log p}{\partial\theta} = \tfrac12\alpha^\top\dfrac{\partial K_y}{\partial\theta}\alphaThe term −12tr⁡(Ky−1∂Ky/∂θ)-\tfrac12\operatorname{tr}(K_y^{-1}\partial K_y/\partial\theta) is missing (Problem 5). For θ=σ2\theta = \sigma^2 the surviving term is 12∥α∥2>0\tfrac12\|\alpha\|^2 > 0, so the "gradient" always says to increase the noise; for s2s^2 it is 12α⊤Kα/s2>0\tfrac12\alpha^\top K\alpha/s^2 > 0, so it always says to increase the signal variance; the optimiser drives the hyperparameters off to infinity, where the data fit is best because everything is explained as a draw from a huge prior. Problem 6's balance ∥α∥2=tr⁡Ky−1\|\alpha\|^2 = \operatorname{tr}K_y^{-1} exists only because the second term is there.

4. Lengthscale derivative short a factor of 2

The exponent −r2/(2ℓ2)-r^2/(2\ell^2) has a 22 in it, and it is tempting to let it cancel something.

  1. k=s2exp⁡(−r2/(2ℓ2))k = s^2\exp\big(-r^2/(2\ell^2)\big)Right so far.
  2. “Differentiate the exponent: −r22⋅ddℓℓ−2=−r22⋅(−ℓ−3)-\tfrac{r^2}{2}\cdot\tfrac{d}{d\ell}\ell^{-2} = -\tfrac{r^2}{2}\cdot(-\ell^{-3}).”The slip that causes the mistake: ddℓℓ−2=−2ℓ−3\tfrac{d}{d\ell}\ell^{-2} = -2\ell^{-3}, and the 22 is exactly what the 12\tfrac12 in the exponent cancels.
  3. ∂k∂ℓ=k r22ℓ3\dfrac{\partial k}{\partial\ell} = k\,\dfrac{r^2}{2\ell^3}The derivative is k r2/ℓ3k\,r^2/\ell^3 (Problem 4): every entry of ∂K/∂ℓ\partial K/\partial\ell is half its true value, so the lengthscale gradient of the marginal likelihood (Problem 5) is halved. The stationary point is unchanged, because a zero gradient is still zero, so a gradient-based optimiser still finds the right ℓ\ell, only more slowly; a finite-difference check of ∂log⁡p/∂ℓ\partial\log p/\partial\ell catches it at once, and a derivative-based test of μ′(x∗)\mu'(x_*) does not, since that uses ∂k/∂x∗\partial k/\partial x_*.

5. Derivative of the posterior mean with the inverse held fixed

α=Ky−1y\alpha = K_y^{-1}y is computed once and stored, and it is easy to forget that it moves with the kernel.

  1. μ∗=k∗⊤α\mu_* = k_*^\top\alphaRight so far: Problem 3.
  2. “α\alpha is the fitted weight vector; differentiate the kernel features k∗k_* and leave the weights alone.”The analogy that causes the mistake: in a linear model the weights are parameters, independent of the features; here α=Ky−1y\alpha = K_y^{-1}y is a function of ℓ\ell through KK.
  3. dμ∗dℓ=(∂k∗∂ℓ)⊤α\dfrac{d\mu_*}{d\ell} = \Big(\dfrac{\partial k_*}{\partial\ell}\Big)^\top\alphaThe second term −k∗⊤Ky−1(∂K/∂ℓ)α-k_*^\top K_y^{-1}(\partial K/\partial\ell)\alpha (Problem 10) is missing: it is the derivative through the inverse, d(A−1)=−A−1(dA)A−1d(A^{-1}) = -A^{-1}(dA)A^{-1}, and it is of the same order as the first. A gradient check against finite differences in ℓ\ell fails. In the marginal likelihood (Problem 5) every dependence on ℓ\ell runs through KyK_y, so there the derivative through the inverse is not a correction but the whole data-fit term.

Print this set: gaussian-process-conditioning.pdf (problems, answers, and worked solutions on separate pages).

More practice sets

All practice sets

eulerpad.app/practice/gaussian-process-conditioning