Why functions sampled from a linear kernel Gaussian Process are guaranteed to be a linear function?

Question

It's well known that a linear kernel Gaussian Process regression is equivalent to Bayesian Linear Regression, because the functions sampled from a linear kernel GP is bound to be a linear function. However, I have trouble to understand why this is the case. Just in case you are not sure, you can see it for your self with this nice visualization tool. Just choose the kernel type to be linear.

My confusion can be further boiled down to the following question.

Let's suppose a simple covariance matrix as below, which is derived from a linear kernel.

[[  0.   0.   0.]
 [  0.  25.  50.]
 [  0.  50. 100.]]

If you draw a sample Y from a 3D multivariate Gaussian with this covariance matrix, and plot the vector Y against X=[1,2,3], the 3 points will always be on a straight line. I just can't figure out why this is the case. I hope someone can help give an intuitive explanation.

Answer 1

Observe that your covariance is degenerate, as it has rank 1. Sampling from this distribution does not give you as much "randomness" as expected as it can be written as the push-forward of a random variable on a lower dimensional space. Let's make this explicit! Notice that we can write your covariance matrix as $$\boldsymbol{\Sigma} = \begin{pmatrix}0 & 0 & 0 \\ 0 & 25 & 50 \\ 0 & 50 & 100 \end{pmatrix} = \begin{pmatrix}0 & 5 & 10 \end{pmatrix}\begin{pmatrix} 0 \\ 5 \\ 10 \end{pmatrix} = \boldsymbol{v}\boldsymbol{v}^{T}$$ which reveals the rank 1 structure of $\boldsymbol{\Sigma}=\boldsymbol{v}\boldsymbol{v}^{T}$. Notice that this structure directly comes from the linearity of the kernel! Denote the samples as $\boldsymbol{y} \sim \mathcal{N}\left(\boldsymbol{0}, \boldsymbol{\Sigma}\right) \in \mathbb{R}^{3}$. The trick is now rewriting this as $$\boldsymbol{y} \stackrel{(d)}{=} z\boldsymbol{v}$$ where $z \sim \mathcal{N}(0,1)$ (one dimensional!). This is due to the fact that scaling $z$ with $v_i$ in each component is still Gaussian and we can check the statistics $$\mathbb{E}[z\boldsymbol{v}] = 0 \hspace{3mm} \text{ and } \hspace{3mm} \text{cov}(z\boldsymbol{v}) = \boldsymbol{v} \text{ cov}(z)\boldsymbol{v}^{T} = \boldsymbol{v}\boldsymbol{v}^{T} = \boldsymbol{\Sigma}$$ which shows $\boldsymbol{y} \stackrel{(d)}{=} z\boldsymbol{v}$. But this means we can also sample $\boldsymbol{y}$ as $z\boldsymbol{v}$ which always is a random multiple of $\boldsymbol{v}$ (again, $z$ is one dimensional) and hence reveals that the three points will always lie on the line given by $\boldsymbol{v}$!

Essentially, having a low rank covariance implies that we only will access low-dimensional randomness although the random variable itself seems to live in a higher dimensional space. This is why all your points end up on a line!

Another, maybe helpful approach looks at the diagonalization of your covariance: $$\boldsymbol{\Sigma} = \boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{T}$$ Due to rank 1, $\Lambda_{11} \not = 0$ while $\Lambda_{22}=\Lambda_{33}=0$. Since Gaussians are invariant under rotations, we can essentially use this new rotated coordinate system where your covariance $\boldsymbol{\Sigma}$ becomes the diagonal covariance $\boldsymbol{\Lambda}$. But here you see that the second and third entry of your random vector are fixed since they have no variance ($\Lambda_{22}=\Lambda_{33}=0$)!

You can check out this answer for a more general explanation of this effect.

Answer 2

Refer to Drawing values from the distribution section in Wikipedia. Let me use that to prove the linearity of functions with basic linear algebra.

We will use the following fact,

If we can write the covariance matrix $\Sigma$ of a multivariate normal distribution as $\Sigma=AA^T$ for any real $A$, then $\mathbf{y} = \boldsymbol{\mu} + A \mathbf{z}$ is a valid function drawn from $\mathcal{N}(\boldsymbol{\mu},\Sigma)$, where $\mathbf{z}$ is a function drawn from $\mathcal{N}(\mathbf{o},I)$ (standard multivariate normal distribution).

We want to show that any $\mathbf{y}$ follows $m\mathbf{x}+b$ (linear) form, where $m$ and $b$ are slope and offset respectively.

According to the distill article you refer to in the question (also in general), the linear kernel is given as follows, $$ K(x,x') = \sigma^2(x-c)(x'-c) + \sigma_b^2 $$ Writing it in covariance matrix form, $$ \Sigma = K(\mathbf{x},\mathbf{x}) = \begin{bmatrix} \sigma^2(x_1-c)^2+\sigma_b^2 & \sigma^2(x_1-c)(x_2-c)+\sigma_b^2 & \cdots\\ \sigma^2(x_2-c)(x_1-c)+\sigma_b^2 & \sigma^2(x_2-c)^2+\sigma_b^2 & \cdots\\ \cdots& \cdots & \cdots \end{bmatrix} $$

If we want to write it in $\Sigma = AA^T$ form, $A$ would be, $$ A= \begin{bmatrix} \sigma(x_1-c) & \sigma_b & 0 & \cdots\\ \sigma(x_2-c) & \sigma_b & 0 & \cdots\\ \cdots& \cdots & \cdots & \cdots\\ \sigma(x_n-c)& \sigma_b & 0 & \cdots \end{bmatrix} $$ Now, for GP, we take $\boldsymbol{\mu}=\mathbf{o}$, so, $\mathbf{y}=A\mathbf{z}$ is a valid function from $\mathcal{N}(\mathbf{o},\Sigma)$. $$ \begin{bmatrix} y_1\\y_2\\\cdots\\y_n \end{bmatrix} =A\mathbf{z}= \begin{bmatrix} \sigma(x_1-c) & \sigma_b & 0 & \cdots\\ \sigma(x_2-c) & \sigma_b & 0 & \cdots\\ \cdots& \cdots & \cdots & \cdots\\ \sigma(x_n-c)& \sigma_b & 0 & \cdots \end{bmatrix} \begin{bmatrix} z_1\\z_2\\\cdots\\z_n \end{bmatrix}=z_1\sigma \begin{bmatrix} x_1\\x_2\\\cdots\\x_n \end{bmatrix}-(z_1\sigma c - z_2\sigma_b) $$

As you can see, $\mathbf{y}$ follows $m\mathbf{x}+b$ form, where slope $m=z_1\sigma$ and offset $b=z_2\sigma_b - z_1\sigma c$.

Hence proved :)

(let me know in comments if any step requires more clarification).