When I was taking classes on machine learning, I learned about SVMs and the kernel trick to get non-linear SVMs. The last slide usually said something along the lines of "you can also kernelize other methods" without giving anymore hints as to which methods this refers to that could be rewritten in terms of inner products and thus kernelized. So it came as a surprise to me when I began reading Gaussian Processes for Machine Learning and learned that not only is Bayesian linear regression (BLR) one of those methods but kernelizing it gets you Gaussian Processes (GPs).¹

Even after Rasmussen introduced the kernelization of BLR as the weight-view of GPs, I was still sceptical and expected to read that this is kernel-BLR and one can generalize it further to get GPs. But no, kernelized BLR and GPs are actually identical. To really convince myself of this fact, I filled in the details in the kernelization process that Rasmussen omits.

In Bayesian linear regression we have a point $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \v{x} \in \R^{D}$ and a noisy observation $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} y \in \R$ and we assume that the observation was generated by the model $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} y = \v{w}\T \v{x} + \epsilon$. Here $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \epsilon \sim \N(0, \sigma^2)$ is independent Gaussian noise and $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \v{w} \sim \N(\v{0}, \m{\Sigma}_p)$ is a latent variable that describes the linear relationship between $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \v{x}$ and $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} y$. Making a prediction $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \pred{y}$ at a new point $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \pred{\v{x}}$ given $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} n$ known data points $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \data = (\obs{\m{X}}, \obs{\v{y}}) \in (\R^{D \times n}, \R^n)$ in the Bayesian setting means to marginalize over the latent variables, i.e.

Some quite laborious transformations (or a high-level argument) will get you the result that $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \pred{y}$ follows another Gaussian distribution

where $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \m{A} = \sigma^{-2}\obs{\m{X}}\obs{\m{X}}\T + \m{\Sigma}_{p}\inv$.

To kernelize the predictive distribution and thus BLR, we have to rewrite both $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \pred{\mu}$ and $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \pred{\sigma}$ using only inner products between data points. Let's begin with the mean.

Now we pull out $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \m{\Sigma_p}$ using $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} (\m{A} + \m{B}\inv)\inv = \m{B}(\m{A}\m{B} + \m{I})\inv$.

Next, we move $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \m{X}$ to form inner products using the push-through identity $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} (\m{I} + \m{A}\m{B})\inv\m{A} = \m{A}(\m{I} + \m{B}\m{A})\inv$ and push $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \sigma$ inside.

At this point we can introduce the kernel $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} k(\v{x}, \v{y}) = \langle \phi(\v{x}), \phi(\v{y}) \rangle$ with the feature mapping $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \phi(\v{x}) = \m{L}\T\v{x}$ where $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \m{L}\m{L}\T = \m{\Sigma}_{p}$ is the Cholesky decomposition of $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \m{\Sigma}_{p}$. Expanding the kernel, we see that $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} k(\v{x}, \v{y}) = (\m{L}\T\v{x})\T(\m{L}\T\v{y}) = \v{x}\m{L}\m{L}\T\v{y} = \v{x}\m{\Sigma}_{p}\v{y}$ and get the final kernelized mean

Kernelizing the variance begins in the same way with pulling out $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \m{\Sigma}_{p}$.

This looks quite similar to what we had before but we are missing an $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \m{X}$ on the right side that prevents us from forming inner products. Instead we need to invoke the Woodbury matrix identity

with $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \m{A} = \m{I}$, $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \m{U} = \obs{\m{X}}$, $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \m{C} = \sigma^{-2}\m{I}$ and $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \m{V} = \obs{\m{X}}\T\m{\Sigma}_{p}$.

At this point, the expression admits the same kernel formulation as the mean did before and we get

Compare the equations for $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \pred{\mu}$ and $\newcommand{\nth}[2]{#1^{\left(#2\right)}} \newcommand{\el}[2]{#1_{#2}} \newcommand{\col}[2]{#1_{:,#2}} \newcommand{\row}[2]{#1_{#2,:}} \renewcommand{\v}[1]{{\bm{#1}}} \newcommand{\m}[1]{{\bm{#1}}} \newcommand{\s}[1]{{\mathcal{#1}}} \renewcommand{\set}[1]{{\left\{#1\right\}}} \def\d{{\mathrm{d}}} \def\R{{\mathbb{R}}} \def\Rplus{{\pos{\R}}} \def\N{{\mathcal{N}}} \def\E{\operatorname*{\mathbb{E}}} \def\T{^{\intercal}} \def\inv{^{-1}} \def\sign{\operatorname{sign}} \def\Tr{\operatorname{Tr}} \def\ceil#1{\lceil #1 \rceil} \def\floor#1{\lfloor #1 \rfloor} \def\sp{\operatorname{sp}} \def\p{\mathrm{p}} \def\data{\mathcal{D}} \newcommand{\pred}[1]{\hat{#1}} \newcommand{\obs}[1]{#1} \pred{\sigma}$ with equations (2.22)-(2.24) in GPML for GP prediction and you will find them equivalent except for additional prediction noise on our part.

It is known that GPs are equivalent to deep neural networks in the infinite-width limit². Speaking informally, this shows that "Bayesianization" and kernelizing can generalize linear regression into something as powerful as deep learning. Goes to show that I underestimated the increase in "power" that kernelizing a method brings.

¹ The distill journal has a nice visual introduction to GPs. ↩

² Apparently, GPs are also equivalent to spline smoothing but that looks to be more of a theoretical result. ↩