Variance of kernel ridge regression vs Gaussian process

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.

Mixing your notation ($\mu_p, V_p$) with that of Rasmussen, the predictive equation for Gaussian Process Regression (GPR) (whose GP has a mean function of $0$) is
$\mathbb{f}_* | X, \mathbb{y}, X_* \sim \mathcal{N}(\mu_p, V_p)$
and for GPR you've identified the conditional predictive mean/expectation $\mu_p = \text{E}[\mathbb{f}_* | X, \mathbb{y}, X_*]$ and conditional covariance $V_p = \text{Cov}(\mathbb{f}_* | X, \mathbb{y}, X_*)$
And you've noticed that (in this case) the predictive output of the Kernel Ridge Regression (KRR) procedure is the same as the conditional predictive mean/expectation $\text{E}[\mathbb{f}_* | X, \mathbb{y}, X_*]$ from the GPR's predictive distribution (i.e. output).
But the variance of this conditional predictive mean/expectation $\text{Var} \Big[ \text{E}[\mathbb{f}_* | X, \mathbb{y}, X_*] \Big]$ (which is what you are looking for) is not necessarilty the same thing as $\text{Cov}(\mathbb{f}_* | X, \mathbb{y}, X_*)$.