So recently I had to learn Gaussian Process (GP) fitting for something I was working on, and man, I was struggling so bad. Every resource just throws \(f_* = k(x_*,X)\,K(X,X)^{-1}y\) at you like it fell from the sky, and expects you to just accept it and move on. Nobody tells you where this formula actually comes from. I hate that. I wanted to build the whole thing from scratch — starting from literally nothing — and land on the exact formula Google gives you.
Turns out, the cleanest way to understand GP is to first re-derive plain old Linear Regression (OLS) in a slightly different way than how they taught us in class — not by minimizing squared error, but by just matching covariances directly. Once you see that, GP basically writes itself. Let's go.
Part 1: OLS, but the "matching" way
Everyone knows OLS as "minimize the squared error." Cool, but let's not do that today. Instead, assume the model:
where \(\varepsilon\) is just leftover noise that has nothing to do with \(X\) (we're just assuming this — it's a modeling choice, called independence of noise and regressor).
Now here's the trick: instead of doing calculus, let's just take covariance with \(X\) on both sides of the equation.
That second term, \(\text{Cov}(\varepsilon,X)\), is zero — because we assumed the noise is independent of \(X\). That's literally the only reason it vanishes, nothing magical.
Rename \(\text{Cov}(Y,X) = \sigma_{xy}\) and \(\text{Var}(X)=\sigma_x^2\), and boom:
No calculus, no \(\arg\min\), just "match the covariance and solve." Same answer as normal OLS, just a different (and honestly cleaner) route.
What about the leftover noise?
Same trick, but take Variance instead of Covariance:
So:
Now, if we're told \(X\) takes a specific known value \(x\), then \(a\cdot x\) is just a fixed number, and \(\varepsilon \sim N(0,\text{Var}(\varepsilon))\), so:
This is literally predicting a distribution for \(Y\) given \(x\) — sounds suspiciously like what GP does, right? Hold that thought.
Part 2: Rewriting this using kernel notation
Now let's just relabel things. Instead of writing \(\text{Var}(X)\), let's call it \(k(x,x)\), and instead of \(\text{Cov}(y,x)\), call it \(k(y,x)\). It's the exact same numbers, just fancier names:
And the prediction:
Nothing changed mathematically — I literally just renamed variance/covariance as "kernel" values. But hold on to this shape, because now let's look at GP.
Part 3: Enter GP
GP regression's mean prediction formula is:
Now compare this to what we just built:
Look closely — it's the same structure, just flipped. In our OLS version, the raw new input \(x_*\) gets multiplied directly, and \(y\) is hidden inside the coefficient. In GP, it's the exact opposite — \(y\) is multiplied directly, and \(x_*\) gets stuffed inside the kernel function.
That's genuinely one real difference already. But there's a bigger, cooler one: GP doesn't pass \(x_*\) through raw — it passes it through a non-linear function \(k(\cdot,\cdot)\), like the RBF kernel:
This looks like "bro that's just one exponential, what's the big deal." But that exponential secretly contains an infinite expansion:
So mathematically, using this kernel is equivalent to having transformed your input into an infinite-dimensional feature vector \([x, x^2, x^3, \dots]\), without ever writing that vector down. That's the whole "kernel trick" — massive expressive power, hidden inside one clean exponential formula.
Part 4: Can we do this trick to OLS too?
Obviously yes. If non-linearity is just "transform the input before doing linear stuff," nothing's stopping us from doing that in plain OLS too. Define a feature map:
Stack all training points' features into a matrix \(\Phi \in \mathbb{R}^{n \times d}\), and just run ordinary OLS on \(\Phi\) instead of \(X\):
And to predict, first transform the new point too:
This is genuinely cool because now you get non-linearity and you're still just doing plain linear algebra — no mysterious kernel magic needed.
So... is GP just OLS wearing a mask?
Not exactly — but the spirit is definitely borrowed from the same place, and there's actually a beautiful matrix identity connecting both formulas exactly (I derived it properly in a longer write-up, drop a comment if you want that). But conceptually here's the real difference, and honestly both have their own trade-offs:
| OLS + feature map | GP + kernel | |
|---|---|---|
| Matrix you invert | \(d \times d\) (feature dim) | \(n \times n\) (data size) |
| Expressiveness | controlled, whatever \(d\) you pick | can be infinite (like RBF) |
| Cost grows with | number of features | number of data points |
| Risk | underfitting if \(d\) too small | overfitting if kernel too flexible |
That last row is honestly the whole trade-off in one line. With OLS + feature maps, you decide how expressive the model gets, by picking \(d\) — you're in control, and if you don't need infinite flexibility (which, let's be real, most of the time you don't), you save yourself from a massive \(n \times n\) matrix inversion and from overfitting garbage.
With GP + kernels, you get insane expressive power basically for free, packed into one clean function — but you pay for it with an \(n \times n\) matrix that gets expensive fast, and a model so flexible it can happily overfit if you're not careful with your kernel's hyperparameters (like that \(\ell\) in the RBF kernel).
So nah, GP isn't "OLS behind a mask" — they're solving genuinely different problems (GP is fundamentally about a distribution over functions, with uncertainty built in; OLS just gives you one line). But the underlying trick — "sneak non-linearity in by transforming the input before doing linear stuff" — is the exact same clever idea in both. Once you see that, both formulas stop looking like they fell from the sky, and start looking like two flavors of the same underlying hack. And honestly, that's what made this whole thing click for me.