Gaussian Process inference

\gdef\delim#1#2#3{\mathopen{}\mathclose{\left#1 #2 \right#3}} \gdef\p#1{\delim({#1})} \gdef\ps#1{\delim\{{#1}\}} \gdef\N{\mathcal N} \gdef\vec#1{\bm #1} \gdef\mat#1{\mathrm #1} \gdef\k{\mathrm k} \gdef\u{\mathrm u} \gdef\c{\mathrm c} \gdef\T{\mathrm T}

Given \vec μ and \mat Σ such that \vec y is multivariate normal distributed.

P\p{\vec y} ∼ \N\p{\vec μ, \mat Σ}

Given a draw of standard normal individually independent values \vec n, we can transform this to a draw from \N\p{\vec μ, \mat Σ}:

\vec y = \vec μ + \mat L \cdot \vec n

Where \mat L is the Cholesky factor such that \mat L ⋅ \mat L^{\T} = \mat Σ.

Computing the Cholesky decompostion

https://en.wikipedia.org/wiki/Cholesky_decomposition#The_Cholesky_algorithm

https://en.wikipedia.org/wiki/Cholesky_decomposition#The_Cholesky%E2%80%93Banachiewicz_and_Cholesky%E2%80%93Crout_algorithms

do i = 1, size(A,1)
    L(i,i) = sqrt(A(i,i) - dot_product(L(i,1:i-1), L(i,1:i-1)))
    L(i+1:,i) = (A(i+1:,i) - matmul(conjg(L(i,1:i-1)), L(i+1:,1:i-1))) / L(i,i)
end do

Inversion using Cholesky

Inference

Supposed we know a subset of values \vec y_{\mathrm k} and we want to use this to infer about the remaining unknown values \vec y_{\mathrm u}.

P\p{\vec y_{\u} \vert \vec y_{\k}} ∼ \N\p{\vec μ_{\c}, \mat Σ_{\c}}

This is again normally distributed with

\begin{aligned} \vec μ_{\c} &= \vec μ_{\u} + Σ_{\k\u}⋅Σ_{\k\k}^{-1}⋅ \p{\vec y_{\k} - \vec μ_{\k}} \\ \mat Σ_{\c} &= \mat Σ_{\u\u} - Σ_{\u\k}⋅Σ_{\k\k}^{-1}⋅ \mat Σ_{\k\u} \end{aligned}

Where {}_\k represents known indices and {}_\u the unknown indices. To avoid inverting we want to re-use \mat L.

Note that the second matrix is also kind of an update:

\mat Σ_{\c} = \mat Σ_{\u\u} - Σ_{\u\k} ⋅ \p{\mat L_{\k\k} ⋅ \mat L_{\k\k}^\T}^{-1} ⋅ \mat Σ_{\u\k}^\T

Efficient computation

Suppose we know \mat L. Can we use this to quickly compute \mat L_{\k\k}? This amounts to computing the Cholesky factor of a subset of the matrix. This can be done using repeated single-row deletion update:

https://en.wikipedia.org/wiki/Cholesky_decomposition#Adding_and_removing_rows_and_columns

Can we also compute \mat L_{\c} efficiently?