Gaussian Process inference for regression: can we derive new evaluation from just the new input points? I'm getting lost on the inference part of Gaussian process regression:
We start with a dataset $X \in \mathbb{R}^{n \times d}$ with corresponding observations $f = f(X)$ and some test data points $X_* \in \mathbb{R}^{n_* \times d}$ for which we want to infer function values $f_* = f(X_*)$.
The $GP$ defines a joint distribution for $p(f, f_*|X, X_*)$
$$\begin{bmatrix}\textbf{f}\\ \textbf{f}_*\end{bmatrix} \sim \mathcal{N}\left(\begin{bmatrix}\mu \\ \mu_*\end{bmatrix}, \, \begin{bmatrix}K & K_*\\ K_*^T & K_{**}\end{bmatrix}\right)$$
with ${\mu} = m(X)$, $\mu_* = m(X_*)$, $K = K(X,X) \, , \, K_* = K(X, X_*) \, , \, K_{**} = K(X_*, X_*)$.
then we apply rules of Multivariate normal conditioning to infer
$$p(f_* | f, X,X_*) \sim \mathcal{N}(\mu_*+K_*^TK^{-1}(f-\mu) \, , \,K_{**}-K_*^TK^{-1}K_*)$$
Until here all makes sense since any subset of uncountable family of gaussians has still a joint gaussian distribution, but:
Are information on the right hand side all known to us, i.e. deriving directly from the test set $X_*$ for which the evaluation $f(X_*)$ are unknown?
What I don't understand, for example, is that by definition:
$$\mu_* = m(X_*) = \mathbb{E}[f(X_*)]$$
and while for known tuples $(x, f(x))$ this actually makes sense, instead we don't know the values $f(X_*)$! What am I missing here?
many thanks,
James
 A: It's the same as with any other predictive model, say linear regression. At training time, you know both the inputs $X$ and the outputs $\mathbf{f}$, so you can learn from them. At prediction time, new inputs arrive $X_*$ and you make predictions for $\mathbf{f}_*$ by plugging it into the model. In the Gaussian process, you calculate the mean $\boldsymbol{\mu}_*$ and covariance $\boldsymbol{\Sigma}_*$ as functions of the data. You can think of the kernels $K$ as a kind of distance matrices, where the prediction is made by finding a value of $\mathbf{f}_*$ such that is close to what we observed in the training data, where the distances inform you as well about the uncertainty (large uncertainty for points unlike the ones you saw in training data). Given $\boldsymbol{\mu}_*$ you can treat it as a point prediction for $\mathbf{f}_*$ (as we do in linear regression), you can use it together with $\boldsymbol{\Sigma}_*$ to draw prediction interval over it, or you can use both $\boldsymbol{\mu}_*$ and $\boldsymbol{\Sigma}_*$, plug them into a Gaussian random number generator to draw samples of the possible values of $\mathbf{f}_*$. The left-hand side values are the random variables following the distribution defined on the right-hand side.
Here and here you can find nice tutorials including code examples in Python and here is a simple notebook in Julia I prepared some time ago. I recommend looking at the code or implementing it by yourself. The math in Gaussian process literature can get dense and it is often easier to understand when you code it.
