Consider the model
$ \mathbf{y} = f(\mathrm{X}) + \epsilon $.
Here $\mathrm{X}$ is a fixed $n \times d$ data matrix, and $\epsilon \sim \mathcal{N}(0, \sigma^2 I)$ is iid Gaussian noise. Assume that $\sigma^2$ is known.
First, consider modeling this using a Gaussian process i.e. $f \sim \mathcal{GP}(0, k)$. Then it can be shown that for a new point $x_\ast$, the predictive distribution is Gaussian with mean and variance given by
$ \mu_p = k(\mathrm{X}, x_\ast)^T(k(\mathrm{X}, \mathrm{X}) + \sigma^2 I)^{-1}\mathbf{y} $,
$ V_p = k(x_\ast, x_\ast) - k(\mathrm{X}, x_\ast)^T(k(\mathrm{X}, \mathrm{X}) + \sigma^2 I)^{-1}k(\mathrm{X}, x_\ast) $,
respectively. Now, consider modeling the same data using kernel ridge regression (with regularization parameter $\lambda$). In this case, we estimate $f$ (assumed to be in the RKHS corresponding to kernel $k$), and get predictions given by
$ \hat{f}(x_\ast) = k(\mathrm{X}, x_\ast)^T(k(\mathrm{X}, \mathrm{X}) + \lambda I)^{-1}\mathbf{y} $,
which is of course the same as the posterior Gaussian process mean (with $\lambda = \sigma^2$), because the two models are just different ways of looking at the same thing.
Now, here is where my confusion arises. Based on this equivalence, it seems to me that the variance of the ridge prediction should match the posterior Gaussian process variance. But this does not seem to be the case. We have,
$\mathbb{V}[\hat{f}(x_\ast)] = k(\mathrm{X}, x_\ast)^T(k(\mathrm{X}, \mathrm{X}) + \lambda I)^{-1} \mathbb{V}[\mathbf{y}] (k(\mathrm{X}, \mathrm{X}) + \lambda I)^{-1}k(\mathrm{X}, x_\ast) = \sigma^2 k(\mathrm{X}, x_\ast)^T(k(\mathrm{X}, \mathrm{X}) + \lambda I)^{-2}k(\mathrm{X}, x_\ast) $.
Using the Woodbury identity, this can be re-written as
$ \frac{\sigma^2}{\lambda}(\phi(x_\ast)^T(\phi(\mathrm{X})^T\phi(\mathrm{X}) + \lambda I)^{-1} \phi(\mathrm{X})^T\phi(\mathrm{X}) \phi(x_\ast) - \phi(x_\ast)^T(\phi(\mathrm{X})^T\phi(\mathrm{X}) + \lambda I)^{-1}\phi(\mathrm{X})^T\phi(\mathrm{X})\phi(\mathrm{X})^T(k(\mathrm{X}, \mathrm{X} + \lambda I)^{-1}\phi(\mathrm{X})\phi(x_\ast)) $,
where $\phi$ is the feature map corresponding to $k$. This is similar to the posterior Gaussian process variance, but not equal (with $\lambda = \sigma^2$). We can get approximate equality by taking $(\phi(\mathrm{X})^T\phi(\mathrm{X}) + \lambda I)^{-1}\phi(\mathrm{X})^T\phi(\mathrm{X}) \approx I$, but it is not clear to me why there isn't a strict equality like in the case of the mean.