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

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.

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.

• Thanks for your answer. While I can see your line of reasoning proves sampling from the family of z*v will result in a covariance matrix as specified, I still have difficulty to understand why it is the only possible type of vector which can be drawn from a multivariate with such covariance matrix. To put it another way, why it is not possible for a sample from a multivariate with such covariance matrix is not on a linear line? Commented Nov 27, 2020 at 5:44
• You're welcome! I edited my answer and added another short paragraph. I hope this explains it more intuitively. A rank 1 covariance matrix essentially means that you can find a rotation in which the rotated system has a diagonal covariance but only the first eigenvalue is non-zero. This means that all other entries are deterministic (they have zero variance) and hence you land in a 1d subspace, given by a line. Commented Nov 27, 2020 at 8:29

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).