Explicit Mathematical Description of Gaussian Process Regression with 2D inputs - Cross Validated most recent 30 from stats.stackexchange.com 2022-01-20T17:58:59Z https://stats.stackexchange.com/feeds/question/374429 https://creativecommons.org/licenses/by-sa/4.0/rdf https://stats.stackexchange.com/q/374429 1 Explicit Mathematical Description of Gaussian Process Regression with 2D inputs Finncent Price https://stats.stackexchange.com/users/224527 2018-10-30T13:09:18Z 2018-10-31T15:01:22Z <p>Rasmussen and Williams section 2.2 (page 16) gives a formula for the posterior distribution of test points, <span class="math-container">$f_{\star}$</span>, of a Gaussian Process when conditioned on some training points, <span class="math-container">$f$</span> 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?</p> <p>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 <span class="math-container">$h$</span>-dimensional inputs make no sense to me.</p> <p>They give the process for sampling from a multivariate Gaussian distribution in Appendix A.2: If <span class="math-container">$x \sim \mathcal{N}(\mathbf{m}, K)$</span>, then <span class="math-container">$\mathbf{x} = \mathbf{m} + L\mathbf{u}$</span>, where <span class="math-container">$\mathbf{u}$</span> is a vector the length of <span class="math-container">$\mathbf{x}$</span> with each term drawn independently from a standard normal distribution and <span class="math-container">$L$</span> is the Cholesky decomposition of the matrix <span class="math-container">$K$</span>.</p> <p>When <span class="math-container">$h=1$</span>, the process is exactly as described in the text. However, if <span class="math-container">$h=2$</span>, what do I do? The output is scalar, the input is 2D. I know the covariance matrix relates <em>indexes</em> and not locations, and this works as simple matrix math when <span class="math-container">$h=1$</span>, 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.</p> <p>Getting back to <span class="math-container">$h=2$</span>, evaluating <span class="math-container">$K$</span> is straightforward and the Cholesky decomposition is the same as in the <span class="math-container">$h=1$</span> case. I guess that the output is then given by <span class="math-container">$X = \mu + (L\mathbf{u}_1) \otimes (L\mathbf{u}_2)$</span> where <span class="math-container">$X$</span> is a matrix value of the outputs at all of the finite number of test locations, <span class="math-container">$\mathbf{\mu}_x$</span> is the prior mean evaluated at the input coordinates (remember the input coordinates are 2D, meaning <span class="math-container">$\mathbf{\mu}_x$</span> is also a matrix), and <span class="math-container">$\mathbf{u}_1$</span> and <span class="math-container">$\mathbf{u}_2$</span> are two different vectors of lengths. <strong>I have no idea if this is correct, but it seems reasonable.</strong></p> <p>Moving on in my question(s). Given a joint normal distribution as in Appendix A.2 (Equation A.5)</p> <p><span class="math-container">$$\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 &amp; C \\ C^T &amp; B \end{bmatrix} \right),$$</span></p> <p>they define (Equation A.6)</p> <p><span class="math-container">$$\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)$$</span></p> <p>Once again, when <span class="math-container">$h=1$</span>, I have no problem evaluating these functions. When <span class="math-container">$h=2$</span>, I have no idea what to do. For example, assume I have 10 test points and 3 training points. Then <span class="math-container">$A$</span> is a <span class="math-container">$10\times10$</span> matrix, <span class="math-container">$B$</span> is a <span class="math-container">$3\times3$</span> matrix and <span class="math-container">$C$</span> is a <span class="math-container">$10\times3$</span> matrix. Further, <span class="math-container">$(\mathbf{y} - \mathbf{\mu}_y)$</span> is a <span class="math-container">$3\times1$</span> vector.</p> <p>Looking at the posterior mean, the dimensions make no sense: <span class="math-container">$$\mathbf{\mu}_x + CB^{-1} (\mathbf{y} - \mathbf{\mu}_y) \\ 10\times10 + (10\times3)(3\times3)(3\times1)$$</span></p> <p>where the last line shows the sizes of the various matrices. This is very clearly the addition of a <span class="math-container">$10\times10$</span> matrix and a <span class="math-container">$10\times1$</span> matrix, which makes no sense.</p> <p>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?</p> <p>--</p> <p>Note that <a href="https://stats.stackexchange.com/questions/320068/sample-from-gaussian-process-across-2d">this question</a> 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.</p> <p>I found <a href="https://www.science-emergence.com/Codes/Gaussian-Processes-for-regression-and-classification-2d-example-with-python/" rel="nofollow noreferrer">this website</a> 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.</p> https://stats.stackexchange.com/questions/374429/explicit-mathematical-description-of-gaussian-process-regression-with-2d-inputs/374640#374640 0 Answer by Finncent Price for Explicit Mathematical Description of Gaussian Process Regression with 2D inputs Finncent Price https://stats.stackexchange.com/users/224527 2018-10-31T13:47:29Z 2018-10-31T15:01:22Z <p>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 <strong><em>and</em></strong> 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.</p> <p>For a univariate distribution <span class="math-container">$f \sim \mathcal{N}(m,\sigma^2)$</span>, I can get a single random sample point with the following formula: <span class="math-container">$f_r = m + \sigma \mathcal{N}(0,1)$</span> where <span class="math-container">$\mathcal{N}(0,1)$</span> 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 <span class="math-container">$m$</span> and the expected variance is <span class="math-container">$\sigma^2$</span>.</p> <p>Let's generalize that last computational procedure by assuming now that the parameters of the normal distribution are functions of some input variable <span class="math-container">$x$</span>, which is itself a scalar. Now we have <span class="math-container">$f \sim \mathcal{N}(m(x),\sigma^2(x))$</span> 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 <span class="math-container">$x_j$</span>, by repeating the sampling procedure <span class="math-container">$f_{rj} = m(x_j) + \sigma(x_j)$</span> however many times I need.</p> <p>At this point, the only information I need is the definitions of <span class="math-container">$m(x)$</span> and <span class="math-container">$\sigma^2(x)$</span>, which I will come to a little while. But first, I want generalize to a vector input: <span class="math-container">$\mathbf{x} = (x_1,x_2,...x_n)$</span>. Where I define the dimension of the input to be <span class="math-container">$n$</span>. The output is still a scalar given by the sampling distribution <span class="math-container">$f \sim \mathcal{N}(m(\mathbf{x}),\sigma^2(\mathbf{x}))$</span>. I still know how to compute a random sample from this distribution: <span class="math-container">$f_r = m(\mathbf{x}) + \sigma(\mathbf{x})\mathcal{N}(0,1)$</span>, once again assuming I know <span class="math-container">$m(\mathbf{x})$</span> and <span class="math-container">$\sigma(\mathbf{x})$</span>.</p> <p>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 <span class="math-container">$\Sigma$</span>, 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)</p> <p><span class="math-container">$$\Sigma_{pq} = \exp\left( -\frac{1}{2} |\mathbf{x}_p - \mathbf{x}_q|^2 \right).$$</span></p> <p>Now, a single random output of the sampling distribution is given by <span class="math-container">$f_r \sim \mathcal{N}(m(\mathbf{x}),\Sigma)$</span>. Informally, this choice for <span class="math-container">$\Sigma$</span> 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. </p> <p>The key concept to notice here is that the covariance of the <em>scalar outputs</em> is a 2D matrix whose elements are computed through pairs of <em>vector inputs</em>. If I only sample a single point, the <span class="math-container">$\Sigma$</span> matrix is a scalar (<span class="math-container">$\Sigma \rightarrow \sigma^2$</span>), and I can compute the output value of that sample as follows</p> <p><span class="math-container">$$f_r = m(\mathbf{x}) + \sigma \mathcal{N}(0,1) = m(\mathbf{x}) + \mathcal{N}(0,1),$$</span></p> <p>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 <span class="math-container">$\Sigma$</span>. I can represent computing <span class="math-container">$k$</span> of them at once using a compact notation <span class="math-container">$\mathbf{f} = \mathcal{N}(\mathbf{m}, \Sigma)$</span>, where <span class="math-container">$\mathbf{f} = (f_1,f_2,...,f_k)$</span> and <span class="math-container">$\mathbf{m} = (m(\mathbf{x}_1),m(\mathbf{x}_2),...,m(\mathbf{x}_k))$</span>. I compute a specific set of these sampled outputs with the formula <span class="math-container">$\mathbf{f} = \mathbf{m} + L \mathbf{u}_k$</span>, where <span class="math-container">$I_k$</span> is the identity matrix of size <span class="math-container">$k \times k$</span>, <span class="math-container">$L$</span> is the <a href="https://en.wikipedia.org/wiki/Cholesky_decomposition" rel="nofollow noreferrer">Cholesky decomposition</a> of the matrix <span class="math-container">$\Sigma$</span> and <span class="math-container">$\mathbf{u}_k = \mathcal{N}(0,I_k)$</span> means "make a vector of <span class="math-container">$k$</span> independent samples from the standard normal distribution."</p> <p>The important point here is that once again the vector operation is over the <em>scalar outputs</em> and not over the (possibly) vector <em>inputs</em>. The index <span class="math-container">$k$</span> is an index over outputs which each have an associated input vector, <span class="math-container">$\mathbf{x}_k$</span>.</p> <p>This is the mistake I was making, and I'll give two examples to show what I did wrong. </p> <p><strong>Scalar Input</strong></p> <p>Assume the inputs are scalars and take on the values <span class="math-container">$x \in (-1,0,1) \equiv (x_1,x_2,x_3)$</span>. For which I can compute (rounded to save space):</p> <p><span class="math-container">$$\Sigma = \begin{bmatrix} 1.00 &amp; 0.61 &amp; 0.14 \\ 0.61 &amp; 1.00 &amp; 0.61 \\ 0.14 &amp; 0.61 &amp; 1.00 \end{bmatrix}.$$</span></p> <p><span class="math-container">$\Sigma$</span> is a <span class="math-container">$3\times3$</span> matrix and the index of the input in the vector <span class="math-container">$x$</span> 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.</p> <p><strong>Vector Input</strong></p> <p>Now assume the inputs are vectors that come from all the permutations of <span class="math-container">$x \in (-1,0,1) \equiv (x_1,x_2,x_3)$</span> and <span class="math-container">$y \in (-1,0,1) \equiv (y_1,y_2,y_3)$</span>. 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 <span class="math-container">$3\times3$</span> matrix, just like I did for the scalar input case. But this is not correct!</p> <p>There are now <em>nine</em> unique inputs whose indexes are unrelated to the input indexes, they are</p> <p><span class="math-container">$$\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}.$$</span></p> <p>Associated with these nine vector inputs are nine scalar outputs that have covariance matrix <span class="math-container">$\Sigma$</span>, which is a <span class="math-container">$9\times9$</span> matrix. Let me assume the mean function is given by <span class="math-container">$m(\mathbf{v}) = v_1 * v_2$</span>. I now have everything I need to produce random prior and posterior functions of two dimensions. The prior is computed as </p> <p><span class="math-container">$$\mathbf{f}_{prior} = \mathbf{m} + L \mathbf{u}_9.$$</span></p> <p>Where <span class="math-container">$\Sigma = L L^T$</span> is the Cholesky decomposition of the matrix <span class="math-container">$\Sigma$</span>.</p> <p>The posterior is computed by first conditioning the output on the some input values, let me assume I have 3 of them <span class="math-container">$(\mathbf{x}_1,y_1)$</span>, <span class="math-container">$(\mathbf{x}_2,y_2)$</span> and <span class="math-container">$(\mathbf{x}_3,y_3)$</span>. Using the marginalization properties of normal distributions I construct an expanded covariance matrix:</p> <p><span class="math-container">$$\Sigma_p = \begin{bmatrix} A &amp; C \\ C^T &amp; B \end{bmatrix}$$</span></p> <p>where <span class="math-container">$A$</span> is the original <span class="math-container">$9\times9$</span> matrix (I called it <span class="math-container">$\Sigma$</span> for the prior), <span class="math-container">$B$</span> is a <span class="math-container">$3\times3$</span> matrix I compute in a way similar to <span class="math-container">$A$</span>, and <span class="math-container">$C$</span> is a <span class="math-container">$9\times3$</span> matrix that is computed using the 9 <span class="math-container">$\mathbf{v}$</span> and the 3 <span class="math-container">$\mathbf{x}$</span> vectors I already have. The new covariance matrix is <span class="math-container">$12\times12$</span> and I compute the value of the posterior function at the 9 <span class="math-container">$v$</span> using the formulas given in my question <span class="math-container">$\mathbf{f}_\mathbf{v}|\mathbf{y} = \mathbf{m}(\mathbf{v}) + CB^{-1}(\mathbf{y} - \mathbf{m}(\mathbf{x}))+ L_p \mathbf{u}_9$</span> and <span class="math-container">$L_p$</span> is the Cholesky decomposition of the matrix <span class="math-container">$K = A - CB^{-1}C^T = L_p L_p^T$</span>.</p> <p>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!</p>