Why do sampled functions from a Gaussian Process perfectly interpolate the data? I've been studying GPs with Rasmussens book as well as a few of my old probability favorites.
I am confused by the results we're getting (doing also a lot of practice) when sampling from our Gaussian Process.
Conditional on the data $ D_n = [{X_n, Y_n}] $ for $n$ datapoints the predictive posterior distribution for a new datapoint $x$  is $GP(x) = Y(x) \sim N(\mu(x), \Sigma(x))$ Where $\mu(x)$ is linearly dependent on the $y_n$.
I do find most of it intuitive, except for the "knots" see the screenshot from Rasmussen's book
.
At first, I thought that it was just illustrative but I was able to get the same "shapes" (though not as pretty as my plotting skills are limited) with arbitrary data, though also not consistently, just most of the time, at times I don't get the knots everywhere, especially with many points.
If I understand correctly since the covariance at those points is 1, for any realization, at those points $Y(x_1) = Y(x_2)$ if $x_1 = x_2$ But only for that sampled function.
Why are all the sampled functions getting the same value?
I do understand that this intuitively is the point to the whole exercise.
We sample from possible functions in the $GP$ that map $x$, to $y$. I am just confused as on how to see it from the math.
This is probably quite the silly question and I am missing something here, I apologize in advance and thank you for your patience!

After playing around with a lot of algebra I've come to a conclusion though i don't feel confident enough to accept my answer I thought I may provide reasoning.
If we have two gaussian from the GP, $f(x_1) \sim N(\mu_1, \Sigma_{11})$ and $f(x_2) \sim N(\mu_2, \Sigma_{22})$
The conditional distribution
$$ f(x_1) | f(x_2) = y \sim N(\mu^*, \Sigma^*)$$
has mean function $$\mu^* = \mu_1 + \Sigma_{12}\Sigma_{22}^{-1}(y-\mu_2)$$
and Covariance $$\Sigma^* = \Sigma_{11} - \Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}$$
Now $\Sigma^*$ turns out to be $0$ if $x1 = x2$
because $\Sigma_{12}$ should equal $\Sigma_{22}$ and by the same reason we get
$$\mu^* = \mu_1 + (y-\mu_2)$$
but $\mu_1 = \mu_2$ therefore we have $\mu^* = y$
since the variance = 0, we get a perfectly interpolated point.
Feel free to correct me, this took a long think and I am not confident on it haha.
