CS 6241: Numerics for Data Science

Computing with kernels

The basic tasks

In the past two lectures, we have discussed various interpretations of kernels and how they are used for interpolation and regression. In today’s lecture, we focus on the numerical methods used to fit kernel-based models and to use them for prediction and uncertainty estimation. More specifically, we consider four problems:

• For a fixed kernel, how do we fit the parameters for a kernel-based model? That is, how do we compute $c$ such that $$\overline{f}(x)=\sum _{j=1}^n{c}_{j}k(x,x_j)$$ approximately satisfies $\hat{f}(x_i)\approx y_i$, where $(x_i,y_i)$ are our observed data points?

• Many kernels involve hyper-parameters, such as length scale and smoothness parameters. How should we determine these hyper-parameters?

• How do we quickly evaluate the predicted values $\hat{f}(x)$ at new points $x$?

• How do we quickly evaluate measures that quantify uncertainty, such as the predictive variance?

Because fitting model parameters and hyper-parameters is more costly than evaluating the predicted mean field or predictive variance, we will focus most of our attention on the fitting tasks. We will also focus primarily on the probabilistic interpretation of kernels in terms of Gaussian processes. Because we have in mind downstream tasks in which a user can adaptively obtain new data (as opposed to fitting to a fixed data set), we will also discuss incremental fitting.

Learning at small $n$

We initially consider the case where there are few enough data points that it makes sense to use standard direct factorization methods to solve the fitting problem. These standard factorizations require $O(n^3)$ time, but if we use well-tuned libraries, the constants need not be large. Using a standard laptop, we can work with data sets of size greater than $n=1000$ as a matter of routine. The direct factorization approach is also useful as a starting point for scalable algorithms more appropriate to large $n$.

Cholesky and kernel-only fitting

The Cholesky factorization of a positive definite kernel matrix $K$ is $$K=R^{T}R$$ where $R$ is upper triangular. The fitting system in the case of a positive definite kernel with no tail is $$Kc=y,$$ and the standard solution algorithm is a Cholesky factorization of $K$ followed by two triangular solves $$c=R^{-1}(R^{-T}y).$$ The Cholesky factorization costs $O(n^3)$; the triangular solves each cost $O(n^2)$. Hence, the cost is dominated by the cost of the Cholesky factorization.

For computation, it is useful to think of the decomposition in block terms $$\left[\begin{array}{cc}{K}_{11}& K_{12}\\ K_{12}^T& K_{22}\end{array}\right]=\left[\begin{array}{cc}{R}_{11}^T& 0\\ R_{12}^T& R_{22}^T\end{array}\right]\left[\begin{array}{cc}{R}_{11}& R_{12}\\ 0& R_{22}\end{array}\right].$$ Rearranging block-by-block, we have $$\begin{array}{rlrl}{K}_{11}& =R_{11}^T{R}_{11}& R_{11}^T{R}_{11}& =K_{11}\\ K_{12}& =R_{11}^T{R}_{12}& R_{12}& =R_{11}^{-T}{K}_{12}\\ K_{22}& =R_{12}^T{R}_{12}+R_{22}^T{R}_{22}& R_{22}^T{R}_{22}& =K_{22}-R_{12}^T{R}_{12}\end{array}$$ That is, we can think of the decomposition in terms of a decomposition of the leading submatrix of $K$, a triangular solve to get an off-diagonal block, and then decomposition of the trailing Schur complement $$S_{22}=R_{22}^T{R}_{22}=K_{22}-R_{12}^T{R}_{12}=K_{22}-K_{21}{K}_{11}^{-1}{K}_{12}.$$ The Schur decomposition has a meaning independent of the Cholesky factorization, as we can see by solving the system $$\begin{gathered} \left[\begin{array}{cc}{R}_{11}^T& 0\\ R_{12}^T& R_{22}\end{array}\right]\left[\begin{array}{cc}{R}_{11}& R_{12}\\ 0& R_{22}\end{array}\right] \\ \left[\begin{array}{c}{X}_1\\ X_2\end{array}\right]=\left[\begin{array}{c}0\\ I\end{array}\right]. \end{gathered}$$ Block forward substitution gives $$\left[\begin{array}{cc}{R}_{11}& R_{12}\\ 0& R_{22}\end{array}\right]\left[\begin{array}{c}{X}_1\\ X_2\end{array}\right]=\left[\begin{array}{c}0\\ R_{22}^{-T}\end{array}\right]$$ and block back substitution then yields $$\left[\begin{array}{c}{X}_1\\ X_2\end{array}\right]=\left[\begin{array}{c}-R_{11}^{-1}{R}_{12}{R}_{22}^{-1}{R}_{22}^{-T}\\ R_{22}^{-1}{R}_{22}^{-T}\end{array}\right]=\left[\begin{array}{c}-S_{11}^{-1}{K}_{12}{S}_{22}^{-1}\\ S_{22}^{-1}\end{array}\right].$$ That is, the inverse of the Schur complement is a submatrix of the inverse of the original matrix $K$: $S_{22}^{-1}={\left[\begin{array}{c}{K}^{-1}\end{array}\right]}_{22}$. In the Gaussian process setting, if $K$ is the covariance matrix, then $K^{-1}$ is the precision matrix, and $S_{22}^{-1}$ is the precision matrix of the posterior for the second block of variables conditioned on observations of data points from the first block.

In introductory numerical methods courses, one typically first encounters the (scalar¹) right-looking version of the Cholesky algorithm:

  % Overwrite K with the Cholesky factor R
  for j = 1:n
    K(j,j) = sqrt(K(j,j));
    K(j,j+1:end) = K(j,j+1:end) / K(j,j);
    K(j+1:end,j) = 0;
    K(j+1:end,j+1:end) = K(j+1:end,j+1:end) - K(j,j+1:end)'*K(j,j+1:end);
  end

At step $j$ of the loop, we have computed the first $j-1$ rows of $R$, and the trailing submatrix contains the Schur complement. At the end of the step, we perform a rank 1 update to get a new (smaller) Schur complement matrix.

An alternate organization, the left-looking version, defers the Schur complement update until it is needed:

  % Overwrite K with the Cholesky factor R
  K(1,1) = sqrt(K(1,1));
  K(2:end,1) = 0;
  for j = 2:n
    K(1:j-1,j) = K(1:j-1,1:j-1) \ K(1:j-1,j)
    K(j,1:j-1) = 0;
    K(j,j) = sqrt(K(j,j) - K(1:j-1,j)'*K(1:j-1,j));
  end

Both organizations of Cholesky do the same operations, but in a different order. Both require $O(n^3)$ time; the dominant cost per step for the right-looking variant is the rank-1 update of the Schur complement, where the dominant cost per step for the left-looking variant is the triangular solve.

The advantage of the left-looking factorization is that it can be applied incrementally. That is, suppose that $K_{11}$ corresponds to the kernel matrix associated with an initial sample of data points, and we have computed $K_{11}=R_{11}^T{R}_{11}$. Then to add a second set of data points, we can do a left-looking block update $$\begin{array}{rl}{R}_{12}& =R_{11}^{-T}{K}_{12}\\ R_{22}^T{R}_{22}& =K_{22}-R_{12}^T{R}_{12}\end{array}$$ The cost of this update is $O((n_1^2+n_2^2)n_2)$, which is significantly less than the $O((n_1+n_2{)}^3)$ cost of recomputing the factorization from scratch in the case that $n_2\ll n_1$.

Fitting with a tail

We now consider fitting in the presence of a tail term. Polynomial tails are typically used with conditionally positive definite kernel functions, but they can also be used with positive definite kernels — and indeed they usually are used in some geostatistical applications. One can also incorporate non-polynomial terms into the tail if they are useful to the model.

The fitting problem with a tail looks like $$\hat{f}(x)=\sum _i{c}_{i}k(x,x_i)+\sum _j{d}_j{p}_j(x)$$ where the coefficients $d$ satisfy the discrete orthogonality condition $$\sum _i{p}_j(x_i)d_i=0$$ for each basis function $p_j(x)$. This gives us the linear system $$\left[\begin{array}{cc}K& P\\ P^T& 0\end{array}\right]\left[\begin{array}{c}c\\ d\end{array}\right]=\left[\begin{array}{c}y\\ 0\end{array}\right].$$ We can also see this linear system from the perspective of constrained optimization: we are minimizing a quadratic objective (associated with $K$) subject to linear constraints (associated with $P$). Note that this formulation is only well posed if $P$ is full rank, a condition sometimes known as unisolvency of the interpolation points.

One way to deal with the tail term is the “null space” approach; that is, rather than adding equations that enforce the discrete orthogonality constraint, we find a coordinate system in which the discrete orthogonality constraint is automatic. Specifically, suppose we write the full QR decomposition of $P$ as $$P=\left[\begin{array}{cc}{Q}_1& Q_2\end{array}\right]\left[\begin{array}{c}{R}_1\\ 0\end{array}\right].$$ Then the constraint $P^{T}c=0$ can be rewritten as $c=Q_{2}w$, giving us $$Q_2^{T}K{Q}_{2}w=Q_2^{T}y$$ where $Q_2^{T}K{Q}_2$ is generally symmetric and positive definite even for a conditionally positive definite kernel. Once we have computed $w$ (and from there $c$), we can compute $d$ by the relation $$R_{1}d=Q_1^T(y-Kc).$$

An alternate approach is to partition the data points into two groups, the first of which is unisolvent and has the same size as the dimension of the tail. Then we can write any $c$ satisfying $P^{T}c=0$ in the form $$\left[\begin{array}{c}{c}_1\\ c_2\end{array}\right]=\left[\begin{array}{c}-W^T\\ I\end{array}\right]c_2$$ where $W=P_2{P}_1^{-1}$. Substituting this into the constrained optimization of $c^{T}Kc$ gives the reduced problem $${\tilde{K}}_{22}{c}_2=(K_{22}-W{K}_{12}-K_{21}{W}^T+W{K}_{11}W)c_2=y_2,$$ from which we can recover the remaining coefficients by solving $$\begin{array}{rl}{c}_1& =-W^T{c}_2\\ P_{1}d& =y_1-K_{11}{c}_1-K_{12}{c}_2.\end{array}$$ This reduction is formally equivalent to Gaussian elimination and substitution on the system $$\left[\begin{array}{ccc}0& P_1^T& P_2^T\\ P_1& K_{11}& K_{12}\\ P_2& K_{21}& K_{22}\end{array}\right]\left[\begin{array}{c}d\\ c_1\\ c_2\end{array}\right]=\left[\begin{array}{c}0\\ y_1\\ y_2\end{array}\right].$$ As with the left-looking factorization described in the previous section, we can form $W$, $K$, ${\tilde{K}}_{22}$, and the Cholesky factorization of ${\tilde{K}}_{22}$ incrementally as new data points are added.

Likelihoods and gradients

Recall from the last lecture that the log likelihood function for a (mean zero) Gaussian process is $$\mathcal{L}=-\frac{1}{2}{y}^T{K}^{-1}y-\frac{1}{2}\log detK-\frac{n}{2}\log (2\pi ).$$ Given a Cholesky factorization $K=R^{T}R$, we can rewrite this as $$\mathcal{L}=-\frac{1}{2}\parallel R^{-T}y{\parallel }^2-\log detR-\frac{n}{2}\log (2\pi ),$$ and note that $$\log detR=\sum _i\log r_{ii}.$$ The cost of evaluating the log likelihood is dominated by the cost of the initial Cholesky factorization.

In order to optimize kernel hyper-parameters via maximum likelihood, we would also like to compute the gradient (and maybe the Hessian) with respect to the hyper-parameters. Recall from last time that we computed the derivative $$\delta \mathcal{L}=\frac{1}{2}{c}^T[\delta K]c-\frac{1}{2}\mathrm{tr}(K^{-1}\delta K)$$ where $Kc=y$. Unfortunately, I know of no tricks to exactly compute $\mathrm{tr}(K^{-1}\delta K)$ for arbitrary $\delta K$ without in time less than $O(n^3)$ without exploiting additional structure beyond what we have discussed so far.

Simply computing gradients of the log likelihood is sufficient for gradient descent or quasi-Newton methods such as BFGS. However, if the number of hyper-parameters is not too great, we may also decide to compute second derivatives and compute a true Newton iteration. Let $\theta$ be the vector of hyper-parameters and use $[f{]}_{,j}$ to denote the partial derivative of an expression $f$ with respect to ${\theta }_j$; then $$\begin{array}{rl}{\left[y^T{K}^{-1}y\right]}_{,ij}& ={[-y^T{K}^{-1}{K}_{,i}{K}^{-1}y]}_{,j}\\ & =2y^T{K}^{-1}{K}_{,i}{K}^{-1}{K}_{,j}{K}^{-1}y-y^T{K}^{-1}{K}_{,ij}{K}^{-1}y\\ & =2c^T{K}_{,i}{K}^{-1}{K}_{,j}c-c^T{K}_{,ij}c\\ {[\log detK]}_{,ij}& ={\left[\mathrm{tr}\left(K^{-1}{K}_{,i}\right)\right]}_{,j}\\ & =\mathrm{tr}(K^{-1}{K}_{,ij}-K^{-1}{K}_{,i}{K}^{-1}{K}_{,j})\\ {\mathcal{L}}_{,ij}& =\frac{1}{2}{c}^T{K}_{,ij}c-c^T{K}_{,i}{K}^{-1}{K}_{,j}c-\frac{1}{2}\mathrm{tr}(K^{-1}{K}_{,ij}-K^{-1}{K}_{,i}{K}^{-1}{K}_{,j}).\end{array}$$ Barring tricks that take advantage of further structure in the kernel matrix $K$ or its derivatives, computing these second partials has the same $O(n^3)$ complexity as computing the first derivatives.

Optimizing the nugget

An important special case is optimization with only a single hyper-parameter, the noise variance term or “nugget” term² That is, we seek to maximize $$\mathcal{L}(\eta )=-\frac{1}{2}{y}^T(K+\eta I{)}^{-1}y-\frac{1}{2}\log det(K+\eta I)-\frac{n}{2}\log (2\pi ).$$ In this case, rather than doing a Cholesky factorization, we may choose to compute an eigenvalue decomposition $K=Q\mathrm{\Lambda }{Q}^T$ in order to obtain $$\mathcal{L}(\eta )=-\frac{1}{2}{\hat{y}}^T(\mathrm{\Lambda }+\eta I{)}^{-1}\hat{y}-\frac{1}{2}\log det(\mathrm{\Lambda }+\eta I)-\frac{n}{2}\log (2\pi )$$ where $\hat{y}=Q^{T}y$. In this case, the stationary points satisfy a rational equation in $\eta$: $$\sum _{j=1}^n(\frac{{\hat{y}}_j^2}{({\lambda }_j+\eta {)}^2}-\frac{1}{2({\lambda }_j+\eta )})=0.$$ We can run Newton on this equation in $O(n)$ time per step.

Smooth kernels and low-rank structure

The cost of parameter fitting for a general kernel matrix is $O(n^3)$, and the cost of hyper-parameter fitting by a gradient-based method applied to maximization of the log-likelihood involves an additional $O(n^3)$ cost per hyper-parameter per step. However, for smooth kernels, the kernel matrix (without a nugget term) may be effectively very low rank. This case is of enough interest that we give it special treatment here. More specifically, we assume a kernel matrix of the form $$\tilde{K}=K+\eta I$$ where most of the eigenvalues of $K$ are much less than $\eta$. In this case, we can solve the fitting problem and the computation of the log-likelihood and gradients in much less than $O(n^3)$ time.

Pivoted Cholesky and kernel approximation

We begin with an outer product factorization of the kernel matrix $K$ without a regularizing nugget. Specifically, we consider the pivoted Cholesky factorization $${\mathrm{\Pi }}^{T}K\mathrm{\Pi }=R^{T}R$$ where the diagonal elements of $R$ appear in descending order and the off-diagonal elements in each row of $R$ are dominated by the diagonal element. The pivoted Cholesky factorization algorithm looks very much like the standard Cholesky algorithm, except before adding row $j$ to the $R$ factor, we first swap rows and columns $j$ and $j_{max}$ of $K$, where $j_{max}$ is the index of the largest entry in the Schur complement matrix.

The approximation error associated with truncating the pivoted Cholesky factorization is just the Schur complement matrix; that is (suppressing the permutation momentarily), $$\left[\begin{array}{cc}{K}_{11}& K_{12}\\ K_{21}& K_{22}\end{array}\right]-\left[\begin{array}{c}{R}_{11}^T\\ R_{12}^T\end{array}\right]\left[\begin{array}{cc}{R}_{11}& R_{12}\end{array}\right]=\left[\begin{array}{cc}0& 0\\ 0& K_{22}-R_{12}^T{R}_{12}\end{array}\right]$$ Typically, we would stop the factorization when all the diagonal elements in the Schur complement are less than some tolerance. Because the Schur complement is positive definite, the sum of the diagonal elements in the Schur complement is the same as the nuclear norm of the Schur complement, so controlling the size of the largest diagonal element also controls the approximation error in the nuclear norm (which dominates both the Frobenius norm and the spectral norm).

Part of what makes pivoted Cholesky attractive is that we can run the algorithm without ever forming $K$ — we only need the rows of $K$ associated with the pivots selected up to the current step, plus a running computation of the diagonal of the Schur complement. If we have formed $j-1$ steps of the pivoted Cholesky factorization, we can extend the factorization by:

• Searching the Schur complement diagonal $d$ to find the pivot row $j_{max}$ and swapping points $j$ and $j_{max}$ and columns $j$ and $j_{max}$ of $R$ (and $d$).

• Forming the top row of the current Schur complement using the formula $$s_{jl}=k(x_j,x_l)-\sum _{i=1}^{j-1}{r}_{ij}{r}_{il}$$

• Extending the factorization with $s_{jj}=\sqrt{s_{jj}}$ and $r_{jl}=s_{jl}/r_{jj}$ for $l>j$.

• Updating the diagonal of the Schur complement via $d_{ll}^{\mathrm{new}}=d_{ll}-r_{jl}^2$.

Running pivoted Cholesky in this way for $r$ steps only requires $O(nr)$ kernel evaluations, $O(nr)$ storage, and $O(n{r}^2)$ time.

The effects of truncation error

We now consider the solution of the linear system $$(K+\eta I)c=y$$ where $\eta I$ is the nugget term. Suppose we have a low-rank factorization (computed by truncated pivoted Cholesky, for example) of the form $K\approx W{W}^T$. Substituting gives us the approximate system $$(W{W}^T+\eta I)\hat{c}=y.$$ How close are $\hat{c}$ and $c$? Let $K=W{W}^T+E$, i.e. $E$ is the part of $K$ discarded by truncating a decomposition of $K$. Then subtracting the approximate linear systems gives $$(K+\eta I)(c-\hat{c})=E\hat{c},$$ and norm bounds yield $$\parallel c-\hat{c}\parallel \le \parallel (K+\eta I{)}^{-1}E\parallel \parallel \hat{c}\parallel \le \frac{\parallel E\parallel }{\eta }\parallel \hat{c}\parallel .$$ where we have throughout used the matrix two-norm. The matrix two-norm is bounded by the nuclear norm; and if $W$ comes from a pivoted Cholesky factorization, we can compute the nuclear norm of $E$ as the sum of the diagonal elements in the truncated Schur complement term. Hence, if we use pivoted Cholesky with the termination criterion $$\parallel S{\parallel }_{\ast }\le \delta \eta ,$$ then we obtain the relative error bound $$\frac{\parallel c-\hat{c}\parallel }{\parallel \hat{c}\parallel }\le \delta$$ if all computations are done in exact arithmetic.

Solving with the kernel

To solve a system $(W{W}^T+\eta I)c=y$ efficiently, let us introduce a new variable $z=W^{T}c$; putting together the definition of $z$ with the equation for $c$ gives us $$\left[\begin{array}{cc}\eta I& W\\ W^T& -I\end{array}\right]\left[\begin{array}{c}c\\ z\end{array}\right]=\left[\begin{array}{c}y\\ 0\end{array}\right].$$ Of course, we can eliminate $z$ to get back to the original equation; but what happens if we instead eliminate $c$? Then we get the reduced system $$(-I-{\eta }^{-1}{W}^{T}W)z=-{\eta }^{-1}{W}^{T}y.$$ Multiplying through by $-\eta$ gives $$(W^{T}W+\eta I)z=W^{T}y,$$ Back substitution gives the formula $$c=\frac{1}{\eta }(y-W(W^{T}W+\eta I{)}^{-1}{W}^{T}y).$$ This is a special case of the famous Sherman-Morrison-Woodbury formula for the inverse of a matrix (in this case ${\sigma }^{2}I$) plus a low-rank modification.

While we can certainly work with the Sherman-Morrison-Woodbury formula directly, we can do a little better if we recognize that the equations in $z$ alone are the regularized normal equations for the optimization $$\text{minimize }\frac{1}{2}\parallel Wz-y{\parallel }^2+\frac{\eta }{2}\parallel z{\parallel }^2.$$ Furthermore, we have that $$c={\eta }^{-1}(y-Wz);$$ that is, $c$ is just the scaled residual of the least squares problem. In terms of an economy QR decomposition $$\left[\begin{array}{c}W\\ \sqrt{\eta }I\end{array}\right]=\left[\begin{array}{c}{Q}_1\\ Q_2\end{array}\right]R$$ we have $$c={\eta }^{-1}(y-Q_1(Q_1^{T}y)).$$ This computational approach (using an economy QR decomposition) enjoys a modest stability advantage compared to working directly with the Sherman-Morrison-Woodbury formula.

Evaluating the kernel approximant

We have so far treated the pivoted Cholesky approach as a way to approximate the kernel matrix. Another interpretation is that we have computed a semidefinite kernel function $$\hat{k}(x,y)=k_{x\mathcal{I}}{K}_{\mathcal{I}\mathcal{I}}^{-1}{k}_{\mathcal{I}x}$$ where $\mathcal{I}$ refers to the $r$ indices selected in the Cholesky factorization. These two interpretations lead to the same linear system to approximate the coefficient vector $c$, but they suggest two different approximations to the function: $$s(x)=\sum _{i=1}^{n}k(x,x_i)c_i\text{ or }\hat{s}(x)=\sum _{i=1}^n\hat{k}(x,x_i)c_i.$$ The two approximations agree for points $x$ that are sufficiently close to the training data, but they may differ when $x$ is farther away. The former approximation is the one that may have motivated us initially, but it is worth spending a moment on the second approximation.

Writing $\hat{s}$ in matrix terms, we have $$\begin{array}{rl}\hat{s}(x)& ={\hat{k}}_{xX}({\hat{K}}_{XX}+\eta I{)}^{-1}y\\ & =k_{xI}{K}_{\mathcal{I}\mathcal{I}}^{-1}{K}_{\mathcal{I}X}(K_{X\mathcal{I}}{K}_{\mathcal{I}\mathcal{I}}^{-1}{K}_{\mathcal{I}X}+\eta I{)}^{-1}y.\end{array}$$ Equivalently, we can write $$\hat{s}(x)=k_{xI}d$$ where (omitting some algebra) $d$ is the solution to the regularized least squares problem $$\text{minimize }\frac{1}{2}\parallel K_{XI}\hat{d}-y{\parallel }^2+\frac{\eta }{2}{d}^T{K}_{II}d.$$ If $K_{II}=R_{II}^T{R}_{II}$ is the leading part of the pivoted Cholesky factorization, we can rewrite this minimization as $$\text{minimize }\frac{1}{2}{\Vert \left[\begin{array}{c}{K}_{XI}\\ \sqrt{\eta }{R}_{II}\end{array}\right]d-\left[\begin{array}{c}y\\ 0\end{array}\right]\Vert }^2.$$ Alternately, we can think of the pivoted Cholesky factorization as inducing a feature map associated with the approximate kernel $\hat{k}$: $$\begin{array}{rlrl}\hat{k}(x,y)& =\psi (x{)}^T\psi (y),& \psi (x)& =R_{II}^{-T}{k}_{Ix}.\end{array}$$ Let $W$ denote the matrix of feature vectors (i.e. row $i$ is $\psi (x_i{)}^T$); in the pivoted Cholesky factor, this is just $R_{I,:}^T$. Then the approximation system involves solving $$\text{minimize }\frac{1}{2}\parallel Wz-y\parallel +\frac{\eta }{2}\parallel z{\parallel }^2,$$ which is the same minimization we saw in the previous subsection; and we can write the $\hat{s}(x)$ function as $$\hat{s}(x)=\psi (x{)}^{T}z.$$

Predictive variance

Treating $k$ as the covariance kernel of a Gaussian process, the predictive variance at a test point $x$ (conditioned on noisy data at points $X$) is $$v(x)=k(x,x)-k_{xX}(K_{XX}+\eta I{)}^{-1}{k}_{Xx}.$$ As in the previous section, we can again think about using low rank approximation in one of two ways: either we can use the low rank structure with the original kernel for the linear solve $(W{W}^T+\eta I{)}^{-1}{k}_{Xx}$; or we can think of the low rank factorization (via pivoted Cholesky) as defining the new kernel $\hat{k}$, for which the predictive variance is $$\hat{v}(x)=k_{x\mathcal{I}}(K_{\mathcal{I}\mathcal{I}}+{\eta }^{-1}{K}_{\mathcal{I}X}{K}_{X\mathcal{I}}{)}^{-1}{k}_{\mathcal{I}x}.$$ Using the Cholesky factorization $R_{\mathcal{I}\mathcal{I}}^T{R}_{\mathcal{I}\mathcal{I}}=K_{\mathcal{I}\mathcal{I}}$, we can rewrite $$(K_{\mathcal{I}\mathcal{I}}+{\eta }^{-1}{K}_{\mathcal{I}X}{K}_{X\mathcal{I}}{)}^{-1}=\eta {\tilde{Q}}_1^T{\tilde{Q}}_1$$ where ${\tilde{Q}}_1$ comes from the economy QR decomposition $$\left[\begin{array}{c}{K}_{X\mathcal{I}}\\ \sqrt{\eta }{R}_{\mathcal{I}\mathcal{I}}\end{array}\right]=\left[\begin{array}{c}{\tilde{Q}}_1\\ {\tilde{Q}}_2\end{array}\right]\tilde{R},$$ which we may have already formed in the process of computing the coefficients in $\hat{s}(x)=k_{x\mathcal{I}}d$.

Likelihoods and gradients

Recall the log likelihood function is $$\mathcal{L}=\frac{1}{2}{y}^{T}c-\frac{1}{2}\log det\tilde{K}-\frac{n}{2}\log (2\pi )$$ where $\tilde{K}=K+\eta I$ is the regularized kernel matrix. The first and last terms need no special treatment in the case when the kernel matrix has special structure, so long as we are able to rapidly compute the coefficient vector $c$. The only non-obvious computation is the log determinant.

Consider the regularized kernel matrix $\tilde{K}=W{W}^T+\eta I$ where $W\in {\mathbb{R}}^{n\times r}$, and suppose we write a full SVD for $W$: $$W=U\mathrm{\Sigma }{V}^T=\left[\begin{array}{cc}{U}_1& U_2\end{array}\right]\left[\begin{array}{c}{\mathrm{\Sigma }}_1\\ 0\end{array}\right]V^T$$ Then $$U^T\tilde{K}U=\left[\begin{array}{cc}{\mathrm{\Sigma }}^2+\eta I& 0\\ 0& \eta I_2\end{array}\right]$$ Hence, $$\log det(\tilde{K})=\log det(U^T\tilde{K}U)=\log det({\mathrm{\Sigma }}^2+\eta I)+(n-r)\log (\eta ).$$ If we use that ${\mathrm{\Sigma }}^2+\eta I$ is the eigenvalue matrix for $W^{T}W+\eta I=R^{T}R$, with $R$ computed by the economy QR we used for kernel solves, we have $$\log det(K)=2\sum _{j=1}^r\log (r_{jj})-(n-r)\log \eta .$$

What of the derivative of the log likelihood? That is, we would like a fast method to compute $$\delta L=\frac{1}{2}{c}^T(\delta \tilde{K})c-\frac{1}{2}\mathrm{tr}({\tilde{K}}^{-1}\delta \tilde{K}).$$ We can simplify the latter term by observing that $${\tilde{K}}^{-1}=\frac{1}{\eta }(I-Q_1{Q}_1^T)$$ where $$\left[\begin{array}{c}W\\ \sqrt{\eta I}\end{array}\right]=\left[\begin{array}{c}{Q}_1\\ Q_2\end{array}\right]R;$$ substitution and using properties of traces gives us $$\mathrm{tr}({\tilde{K}}^{-1}\delta \tilde{K})=\frac{1}{\eta }(\mathrm{tr}(\delta \tilde{K})-\mathrm{tr}(Q_1^T\delta \tilde{K}{Q}_1)).$$ If we work with the original kernel, it is difficult to do better than this.

Now suppose we differentiate the log-likelihood for the approximate kernel induced by the points $\mathcal{I}$, i.e. $${\tilde{K}}_{XX}=K_{X\mathcal{I}}{K}_{\mathcal{I}\mathcal{I}}^{-1}{K}_{\mathcal{I}X}+\eta I.$$ Differentiating this expression gives us $$\delta \tilde{K}=(\delta K_{:,\mathcal{I}})F^T+F(\delta K_{:,\mathcal{I}}{)}^T-F(\delta K_{\mathcal{I}\mathcal{I}})F^T+(\delta \eta )I$$ where $F=K_{:,\mathcal{I}}{K}_{\mathcal{I}\mathcal{I}}^{-1}$. Therefore $$\begin{array}{rl}\frac{1}{2}{c}^T(\delta K)c& =(\delta K_{:,\mathcal{I}}c{)}^T(Fc)-\frac{1}{2}(Fc{)}^T(\delta K_{II})(Fc)+\frac{1}{2}\parallel c{\parallel }^2\delta \eta \\ \frac{1}{2}\mathrm{tr}({\tilde{K}}^{-1}\delta \tilde{K})& =⟨\tilde{F},\delta K_{:,\mathcal{I}}{⟩}_F-\frac{1}{2}⟨\tilde{F},F\delta K_{\mathcal{I}\mathcal{I}}{⟩}_F+\delta \eta (n-\parallel Q_1{\parallel }_F^2)\\ \tilde{F}& \equiv \frac{1}{\eta }(I-Q_1{Q}_1^T)F\end{array}$$ and so we can compute derivatives of the log likelihood in $O(nr)$ space and $O(n{r}^2)$ time in this approximation.

Beyond low rank

For a number of kernels, the kernel matrix is not all that low rank — the singular values decay, but not so quickly that we would be happy with $O(n{r}^2)$ time algorithms. However, in low-dimensional spaces, one may still have low-rank submatrices, and there are fast factorization techniques based on these rank-structured matrices as well. There are also a variety of iterative solvers that take advantage of fast matrix-vector multiplies; for example, we can use the method of conjugate gradients (CG) to solve linear systems using only matrix-vector products with a kernel matrix. For these types of solvers, the low rank approximations using pivoted Cholesky as described above still give us a good preconditioner that can accelerate convergence.

Footnotes

1. The scalar version of the algorithm works one column at a time. Various blocked versions of the factorization algorithm update in blocks of several columns at a time, with a small (scalar) Cholesky factorization for the diagonal blocks. Block algorithms have the same complexity as the corresponding scalar algorithms, but achieve better performance on modern hardware because they make better use of the cache.

2. The term “nugget” comes from the use of Gaussian process regression in geostatistical applications. When trying to use Gaussian processes to predict the location of precious mineral deposits based on the mineral content in bore measurements, it is necessary to account for the possibility that a measurement may accidentally happen on a mineral nugget.