Explicit Mathematical Description of Gaussian Process Regression with 2D inputs

Question

Rasmussen and Williams section 2.2 (page 16) gives a formula for the posterior distribution of test points, $f_{\star}$, of a Gaussian Process when conditioned on some training points, $f$ in Equation 2.19. At the end of this section they claim that extending the analysis to multidimensional inputs is "trivial" and I do not see this fact at all. My question is how do I do any of this stuff if the inputs are in 2D?

It is clear from the definition of the covariance matrix that it is always a 2D object, in light of this fact, the resulting operations with $h$-dimensional inputs make no sense to me.

They give the process for sampling from a multivariate Gaussian distribution in Appendix A.2: If $x \sim \mathcal{N}(\mathbf{m}, K)$, then $\mathbf{x} = \mathbf{m} + L\mathbf{u}$, where $\mathbf{u}$ is a vector the length of $\mathbf{x}$ with each term drawn independently from a standard normal distribution and $L$ is the Cholesky decomposition of the matrix $K$.

When $h=1$, the process is exactly as described in the text. However, if $h=2$, what do I do? The output is scalar, the input is 2D. I know the covariance matrix relates indexes and not locations, and this works as simple matrix math when $h=1$, because vectors are a handy way to manipulate all the points being evaluated at once. In short, the index in the vector directly corresponds to the index of the covariance function.

Getting back to $h=2$, evaluating $K$ is straightforward and the Cholesky decomposition is the same as in the $h=1$ case. I guess that the output is then given by $X = \mu + (L\mathbf{u}_1) \otimes (L\mathbf{u}_2)$ where $X$ is a matrix value of the outputs at all of the finite number of test locations, $\mathbf{\mu}_x$ is the prior mean evaluated at the input coordinates (remember the input coordinates are 2D, meaning $\mathbf{\mu}_x$ is also a matrix), and $\mathbf{u}_1$ and $\mathbf{u}_2$ are two different vectors of lengths. I have no idea if this is correct, but it seems reasonable.

Moving on in my question(s). Given a joint normal distribution as in Appendix A.2 (Equation A.5)

$$ \begin{bmatrix} \mathbf{x} \\ \mathbf{y} \end{bmatrix} \sim \mathcal{N} \left( \begin{bmatrix} \mathbf{\mu}_x \\ \mathbf{\mu}_y \end{bmatrix} , \begin{bmatrix} A & C \\ C^T & B \end{bmatrix} \right), $$

they define (Equation A.6)

$$ \mathbf{x}|\mathbf{y} \sim \mathcal{N} \left( \mathbf{\mu}_x + CB^{-1} (\mathbf{y} - \mathbf{\mu}_y) , A - C B^{-1}C^T \right) $$

Once again, when $h=1$, I have no problem evaluating these functions. When $h=2$, I have no idea what to do. For example, assume I have 10 test points and 3 training points. Then $A$ is a $10\times10$ matrix, $B$ is a $3\times3$ matrix and $C$ is a $10\times3$ matrix. Further, $(\mathbf{y} - \mathbf{\mu}_y)$ is a $3\times1$ vector.

Looking at the posterior mean, the dimensions make no sense: $$ \mathbf{\mu}_x + CB^{-1} (\mathbf{y} - \mathbf{\mu}_y) \\ 10\times10 + (10\times3)(3\times3)(3\times1) $$

where the last line shows the sizes of the various matrices. This is very clearly the addition of a $10\times10$ matrix and a $10\times1$ matrix, which makes no sense.

The sum total of this is to say that I am very clearly missing something about how to extend Gaussian Processes to multidimensional inputs. What am I missing?

--

Note that this question references a defunct website and uses functions in R to perform the sampling. I am looking for a mathematical description that I can use to understand the problem better and expand to even higher dimension.

I found this website showing a loop-based process in python, but because of the loops and rearrange calls used, it isn't clear to me what is going on.

Answer 1

I have figured out that my problem was that I did not appreciate how much the 1D-input case allows one to use the index of the position of a point in the input and the set index. As such, my notation was "colliding" and causing my confusion. In an effort to clarify my thinking, I'm going to answer my own questions.

For a univariate distribution $f \sim \mathcal{N}(m,\sigma^2)$, I can get a single random sample point with the following formula: $f_r = m + \sigma \mathcal{N}(0,1)$ where $\mathcal{N}(0,1)$ is the standard normal distribution, available in most (all?) programming languages. It can be verified that the expected mean value of a number of such samples is $m$ and the expected variance is $\sigma^2$.

Let's generalize that last computational procedure by assuming now that the parameters of the normal distribution are functions of some input variable $x$, which is itself a scalar. Now we have $f \sim \mathcal{N}(m(x),\sigma^2(x))$ and the sampling distribution is still univariate. In principle, there is no difficulty in computing any number of samples at any number of input values $x_j$, by repeating the sampling procedure $f_{rj} = m(x_j) + \sigma(x_j)$ however many times I need.

At this point, the only information I need is the definitions of $m(x)$ and $\sigma^2(x)$, which I will come to a little while. But first, I want generalize to a vector input: $\mathbf{x} = (x_1,x_2,...x_n)$. Where I define the dimension of the input to be $n$. The output is still a scalar given by the sampling distribution $f \sim \mathcal{N}(m(\mathbf{x}),\sigma^2(\mathbf{x}))$. I still know how to compute a random sample from this distribution: $f_r = m(\mathbf{x}) + \sigma(\mathbf{x})\mathcal{N}(0,1)$, once again assuming I know $m(\mathbf{x})$ and $\sigma(\mathbf{x})$.

In this particular case I have assumed that the variance is only a function of the input, it can also be a function of other outputs, which is the central assumption of Gaussian Processes: the outputs are correlated with each other through a covariance matrix $\Sigma$, whose elements are defined by (this is a very simple case I use for clarity, in principle the matrix only needs to be symmetric and positive definite)

$$ \Sigma_{pq} = \exp\left( -\frac{1}{2} |\mathbf{x}_p - \mathbf{x}_q|^2 \right). $$

Now, a single random output of the sampling distribution is given by $f_r \sim \mathcal{N}(m(\mathbf{x}),\Sigma)$. Informally, this choice for $\Sigma$ means that we think outputs whose inputs are near to each other should be more similar than outputs with inputs that are farther from each other.

The key concept to notice here is that the covariance of the scalar outputs is a 2D matrix whose elements are computed through pairs of vector inputs. If I only sample a single point, the $\Sigma$ matrix is a scalar ($\Sigma \rightarrow \sigma^2$), and I can compute the output value of that sample as follows

$$ f_r = m(\mathbf{x}) + \sigma \mathcal{N}(0,1) = m(\mathbf{x}) + \mathcal{N}(0,1), $$

where the last line follows from the definition of the covariance. If I want to generate more than one sample output, I assumed earlier that they will be correlated through $\Sigma$. I can represent computing $k$ of them at once using a compact notation $\mathbf{f} = \mathcal{N}(\mathbf{m}, \Sigma)$, where $\mathbf{f} = (f_1,f_2,...,f_k)$ and $\mathbf{m} = (m(\mathbf{x}_1),m(\mathbf{x}_2),...,m(\mathbf{x}_k))$. I compute a specific set of these sampled outputs with the formula $\mathbf{f} = \mathbf{m} + L \mathbf{u}_k$, where $I_k$ is the identity matrix of size $k \times k$, $L$ is the Cholesky decomposition of the matrix $\Sigma$ and $\mathbf{u}_k = \mathcal{N}(0,I_k)$ means "make a vector of $k$ independent samples from the standard normal distribution."

The important point here is that once again the vector operation is over the scalar outputs and not over the (possibly) vector inputs. The index $k$ is an index over outputs which each have an associated input vector, $\mathbf{x}_k$.

This is the mistake I was making, and I'll give two examples to show what I did wrong.

Scalar Input

Assume the inputs are scalars and take on the values $x \in (-1,0,1) \equiv (x_1,x_2,x_3)$. For which I can compute (rounded to save space):

$$ \Sigma = \begin{bmatrix} 1.00 & 0.61 & 0.14 \\ 0.61 & 1.00 & 0.61 \\ 0.14 & 0.61 & 1.00 \end{bmatrix}. $$

$\Sigma$ is a $3\times3$ matrix and the index of the input in the vector $x$ corresponds exactly to the index of the output because of this. But this is just a coincidence because of the low dimensionality of the input.

Vector Input

Now assume the inputs are vectors that come from all the permutations of $x \in (-1,0,1) \equiv (x_1,x_2,x_3)$ and $y \in (-1,0,1) \equiv (y_1,y_2,y_3) $. If I incorrectly use the indexes of the inputs (i.e. 1,2 and 3) there is still a perfectly reasonable way to interpret the computation of the covariance matrix, and I end up with another $3\times3$ matrix, just like I did for the scalar input case. But this is not correct!

There are now nine unique inputs whose indexes are unrelated to the input indexes, they are

$$ \mathbf{v} = \begin{bmatrix} -1 \\ -1 \end{bmatrix}, \begin{bmatrix} -1 \\ 0 \end{bmatrix}, \begin{bmatrix} -1 \\ 1 \end{bmatrix}, \begin{bmatrix} 0 \\ -1 \end{bmatrix}, \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} 0 \\ 1 \end{bmatrix}, \begin{bmatrix} 1 \\ -1 \end{bmatrix}, \begin{bmatrix} 1 \\ 0 \end{bmatrix}, \begin{bmatrix} 1 \\ 1 \end{bmatrix}. $$

Associated with these nine vector inputs are nine scalar outputs that have covariance matrix $\Sigma$, which is a $9\times9$ matrix. Let me assume the mean function is given by $m(\mathbf{v}) = v_1 * v_2$. I now have everything I need to produce random prior and posterior functions of two dimensions. The prior is computed as

$$ \mathbf{f}_{prior} = \mathbf{m} + L \mathbf{u}_9. $$

Where $\Sigma = L L^T$ is the Cholesky decomposition of the matrix $\Sigma$.

The posterior is computed by first conditioning the output on the some input values, let me assume I have 3 of them $(\mathbf{x}_1,y_1)$, $(\mathbf{x}_2,y_2)$ and $(\mathbf{x}_3,y_3)$. Using the marginalization properties of normal distributions I construct an expanded covariance matrix:

$$ \Sigma_p = \begin{bmatrix} A & C \\ C^T & B \end{bmatrix} $$

where $A$ is the original $9\times9$ matrix (I called it $\Sigma$ for the prior), $B$ is a $3\times3$ matrix I compute in a way similar to $A$, and $C$ is a $9\times3$ matrix that is computed using the 9 $\mathbf{v}$ and the 3 $\mathbf{x}$ vectors I already have. The new covariance matrix is $12\times12$ and I compute the value of the posterior function at the 9 $v$ using the formulas given in my question $\mathbf{f}_\mathbf{v}|\mathbf{y} = \mathbf{m}(\mathbf{v}) + CB^{-1}(\mathbf{y} - \mathbf{m}(\mathbf{x}))+ L_p \mathbf{u}_9$ and $L_p$ is the Cholesky decomposition of the matrix $K = A - CB^{-1}C^T = L_p L_p^T$.

This probably seems like a lot of words to elucidate this concept, but I reckoned with it for a long time and hopefully this will help someone get through that hurdle faster than me!