Is Gaussian process functional regression a truly Bayesian method (again)?

Question

This question has been asked before but I'd like to come back to it because I point a precise issue out.

Suppose we want to estimate a function $f\left( x \right)$ from data $D = \left( {\left( {{x_1},{y_1}} \right),...,\left( {{x_n},{y_n}} \right)} \right)$ with ${y_i} = f\left( {{x_i}} \right) + {\xi _i}$, ${\xi _i}\mathop \sim \limits^{{\text{i}}{\text{.i}}{\text{.d}}{\text{.}}} \mathcal{N}\left( {0,{\sigma ^2}} \right)$ by Gaussian process functional regression.

Let $X = \left( {{x_1},...,{x_n}} \right)$ and $Y = \left( {{y_1},...,{y_n}} \right)$.

The likelihood is $p\left( {\left. Y \right|X , f , \sigma } \right) \propto {\sigma ^{ - n}}\prod\limits_{i = 1}^n {{e^{ - \frac{{{{\left( {{y_i} - f({x_i})} \right)}^2}}}{{2{\sigma ^2}}}}}} $

We have a ${\text{GP}}\left( {m\left( x \right),k\left( {x,x'} \right)} \right)$ prior on $f$ with hyperparameters $m$ and $k$.

Generally speaking, we can have hyperhyperparameters ${\rm M}$ and ${\rm K}$ for $m\left( x \right)$ and $k\left( {x,x'} \right)$.

Therefore Bayes rule writes

$ p\left( {\left. {f , \sigma , m , k ,{\rm M} , {\rm K}} \right|X , Y} \right) \propto \\ {\sigma ^{ - n}}\prod\limits_{i = 1}^n {{e^{ - \frac{{{{\left( {{y_i} - {x_i}} \right)}^2}}}{{2{\sigma ^2}}}}}} {\text{GP}}\left( {m\left( x \right),k\left( {x,x'} \right)} \right)p\left( {\left. m \right|{\rm M}} \right)p\left( {\left. k \right|{\rm K}} \right)p\left( {\rm M} \right)p\left( {\rm K} \right)p\left( \sigma \right) \\ $

and that's all we should need.

The problem is that we have something more that does not appear in Bayes rule at this point: the ${\text{GP}}\left( {m\left( x \right),k\left( {x,x'} \right)} \right)$ prior is used to assign the probability distribution $\mathcal{N}\left( {m\left( X \right),k\left( {X,X} \right)} \right)$ of r.v. $\left. {f\left( X \right)} \right|X,m,k$ that is used in the update equations to get the posterior Gaussian process.

It seems there is no place for $f(X)|X, m , k$ in Bayes rule because we can only add hyperparameters and $f\left( X \right)$ is not such an hyperparameter but a function of the piece of data $X$. $f(X)|X, m ,k$ doesn't look like a prior because it is conditional on $X$ and depends on $X$, nor a likelihood because the likelihood is ${\left. X , Y \right|f , \sigma }$ nor a posterior because all of them are conditional on $X , Y$.

How to plug $f(X)|X, m , k$ in Bayes rule above and where, given we can only add hyperparameters? If we can't, where does ${\mathcal{N}}\left( {m\left( X \right),k\left( {X,X} \right)} \right)$ come from?

So my question is: can we find a set/logical conjunction of hyperparameters that makes the ${\mathcal{N}}\left( {m\left( X \right),k\left( {X,X} \right)} \right)$ multivariate Gaussian appears somewhere in Bayes rule?

Answer 1

The answer to the question:

where does $\mathcal{N}(m(X),k(X,X))$ come from ?

relies on two elementary pieces of information:

(1) Understanding how marginal posterior probabilities can be derived from Bayes' rule. Consider for example a model with two sets of unknown parameters $\theta$ and $\phi$, but with data $y$ that depends only on $\theta$, that is $P(y|\theta,\phi) = P(y|\theta)$. Applying Bayes' rule we have

$$ P(\theta,\phi | y ) \propto P(y|\theta) \pi(\theta,\phi)$$

Where $\pi(\theta,\phi)$ denotes the joint prior of $\theta$ and $\phi$. Now we can integrate both sides with respect to $\phi$ to obtain the marginal posterior probability of $\theta$: $$ P(\theta| y ) \equiv \int d\phi P(\theta,\phi | y ) \propto P(y|\theta) \int d\phi\pi(\theta,\phi) \equiv P(y|\theta) \pi(\theta)$$

Where $\pi(\theta) \equiv \int d\phi\pi(\theta,\phi)$ is the marginal prior of $\theta$.

In other words, the marginal posterior can simply be obtained by applying Bayes' rule to the marginal prior, when the likelihood depends only on a subset of the parameters:

$$ P(\theta| y ) \propto P(y|\theta) \pi(\theta).$$

(2) Understanding the marginal distributions of a Gaussian process: this in fact follows directly from the very Definition of a Gaussian process: we say that a stochastic process $\{f(x)|x \in \mathcal{X}\}$ is Gaussian if and only if every finite subset $(f(x_1),f(x_2),...,f(x_n))$ has a marginal multivariate normal distribution. The notation $f \sim \mathcal{GP}(m,k)$ is therefore equivalent to $(f(x_1),f(x_2),...,f(x_n)) \sim \mathcal N(\mathbf{m}, \mathbf{K})$ where $\mathbf{m} = (m(x_1),m(x_2),...,m(x_n))$ and $\mathbf{K}_{ij} = k(x_i,x_j)$.

Combining the above two pieces of information leads us to the answer in a straightforward manner. In GP Regression, we are estimating an unknown function $f$, that is, a quantity that has an infinite number of degrees of freedom. (If you don't like thinking about infinite quantities, you can instead just think about a very large number $N$, e.g. by restricting $f$ to some grid points. This doesn't change anything about the logic).

The observed data $y_i$ only depend on the finite subset $\{f(x_i)|i=1,...,n\}$. If we would like to find the marginal posterior of some other subset, e.g. $\{f(x_i)|i=n+1,...,m+n\}$, then we only need to know the marginal prior of $\mathbf{f}=(f(x_1),f(x_2),...,f(x_{n+m}))$, which is by definition multivariate normal. This is how we arrive at

$$ P(\mathbf{f} | \{y_i\} ) \propto \Pi_{i=1}^n P(y_i|f(x_i)) \times \mathcal N(\mathbf f | \mathbf{m} , \mathbf{K})$$

where $\mathbf{m}$ and $\mathbf K$ are derived from the mean and covariance functions $m(x)$ and $k(x,x')$ as explained above.

Notice that, if $P(y_i|f(x_i))$ is normal, then the marginal posterior probability of $\mathbf{f}$ is also multivariate normal, And this is true for any subset $\mathbf{f}$. But this is exactly the requirement for the posterior distribution of the entire function $f$ to be a Gaussian process itself!. We established, therefore, that in GP regression, the posterior probability is also a GP.

Answer 2

where does ${\mathcal{N}}\left( {m\left( X \right),k\left( {X,X} \right)} \right)$ come from?

This occurs as prior.

$$f(X)|X,k \sim {\mathcal{N}}\left( {m\left( X \right),k\left( {X,X^\prime} \right)} \right)$$

Describes a distribution of different functions $f(X)$.

It does occur in your Bayesian rule

$$ p\left( {\left. {f , \sigma , m , k ,{\rm M} , {\rm K}} \right|X , Y} \right) \propto \\ {\sigma ^{ - n}}\prod\limits_{i = 1}^n {{e^{ - \frac{{{{\left( {{y_i} - {x_i}} \right)}^2}}}{{2{\sigma ^2}}}}}} \underbrace{{\text{GP}}\left( {m\left( x \right),k\left( {x,x'} \right)} \right)}_{{\mathcal{N}}\left( {m\left( X \right),k\left( {X,X} \right)} \right)}p\left( {\left. m \right|{\rm M}} \right)p\left( {\left. k \right|{\rm K}} \right)p\left( {\rm M} \right)p\left( {\rm K} \right)p\left( \sigma \right) \\ $$