Interpretation of inverse matrix in mean calculation in Gaussian Process
Mathematical Interpretation
Note that equation (2.23) is simply calculating the conditional distribution of equation (2.21) and then finding the mean. Your question reduces to:
&quot;Given normal variables $X$ and $Y$, why is $\mathbb{E}[Y|X] = \mu_y + Cov(X, Y)Cov(Y, Y)^{-1}(x - \mu_x)$? (note: in the book, $\mu_x = \mu_y = 0$)
Deriving the conditional probability mean is complicated (see The Bivariate Normal Distribution, page 3). A more intuitive look can be seen in the first graph in this page.
Here, the mean of $Y|X$ is linear in what value $X$ takes. The line starts at the intercept $\mu_y$, and increases with slope $\rho (\frac{\sigma_y}{\sigma_x}) = \frac{Cov(X,Y)}{\sigma_x \sigma_y}\left(\frac{\sigma_y}{\sigma_x}\right)=\frac{Cov(X,Y)}{\sigma_x^2} = \frac{Cov(X,Y)}{Var(x)} = \frac{Cov(X,Y)}{Cov(X,X)}$. So the mathematical interpretation of $Cov(X,Y)Cov(X,X)^{-1}$ is that it is the slope of the relationship between the mean of $Y|X$ and the value of $x$ that you are given. As you are given a higher value $x$, say $x + \delta$, then the mean of $Y|X$ raises by $Cov(X,Y)Cov(X,X)^{-1}\delta$.
Why is there even a $Cov(X,Y)Cov(X,X)^{-1}$ term there? For some reason, multiplying $Y|X$ by a $Cov(X,Y)Cov(X,X)^{-1}$ term makes $Y|X$ completely independent of $X$ (which makes sense as the definition of &quot;conditional probability&quot;, because you are already given a value of $X$). This is just a mathematical property, I don't know if there's an intuitive explanation as to why.
Human Interpretation
In case your post just want an intuition as to why there is a $Cov(X,X)^{-1}$ in the prediction of a Gaussian process (and ignoring the conditional probability fluff), I don't think there's a real basis for this, it would only be coincidental as the authors simply used the conditional probability mean formula, but I would guess $Cov(X,X)^{-1}$ somehow normalizes the values of covariance matrix $Cov(X_*, X)$.
For example, if the training set $X$ has a lot of outliers and therefore extremely high variance (e.g. all non-diagonal entries in millions), then it is very likely that $Cov(X_*, X)$ would also be extremely high as $X_*$ follows the same distribution as $X$ (unless each data in $X_*$ matches the exact same variance in $X$). It doesn't make sense to multiply $y$ by millions though, as $y$ is already a somewhat decent estimator/prior.
It makes more sense to normalize $Cov(X_*, X)$ by dividing it with the training data variance $Cov(X,X) = Var(X)$ so that the ratio $Cov(X_*, X)Cov(X, X)^{-1}$ tends to be more closer to 1 when $X_*$ follows the same distribution as $X$ (which should be the case). If the ratio is exactly 1, then $X_*$ has the exact same distribution as $X$, so you just return the prior estimate $y$. If the ratio is far away from $1$, then the test set distribution is wildly different than the training set distribution, so you return a number far away from $y$.