# Can you use a gaussian process to model the smoothness of residuals?

I see a lot of use of Gaussian Processes for regression - fitting a GP model to data points, with a prior specifying the smoothness of the function, and using it to predict new values. However, I'm wondering if it's possible to use a GP to aid in calculating the likelihood of an explicit (parametric) model when you can expect the residuals to be smooth. Say you have data points $D$, an underlying model $X$, a forward model $f(X)$, and covariance $\Sigma$. In that case the likelihood could be written:

$$L(\theta) = \mathcal{N}(D - f(X),\Sigma)$$

Where $\mathcal{N}$ is the multivariate normal. In my case (processed image data), $D$ is an image that was produced by projecting a model, adding random noise, and blurring with some PSF. $X$ is the "true" shape, and $f(X)$ is a the blurred image of that shape. Because of the blurring I expect errors to have smoothness - this is why I don't write the likelihood as just a product of independent likelihoods for each data point (pixel). Higher error covariance means more tolerance of disagreement between $f(X)$ and $D$, which makes sense for more distorted data.

I have a few questions:

• Is this the right way to write the likelihood, where the covariance is used to describe the smoothness of the error? Am I right to think of the residuals as a Gaussian Process?
• How do you estimate the covariance matrix?
• Do you have to invert the covariance in order to calculate the likelihood? For a large image this could be impossible. Assuming the pixels only co-vary locally, can you speed this up?

Edit: What I've done so far.

Let $\mathbf{y} = \mathbf{D}-\mathbf{f(X)}$ be the RV for the residuals. Say that y are observations of a GP with zero mean and the following covariance (let's say square exponential from Rasmussen p. 19):

$$K(\mathbf{x},\mathbf{x}') = \sigma_f^2\exp \left(-\frac{1}{2l^2}|\mathbf{x}-\mathbf{x}'|^2\right )+\sigma_n\delta_{xx'}$$

In that case the distribution for y is: $$\mathbf{y} \sim \mathcal{N}(\mathbf{0},K+\sigma_n I)$$ So the posterior probability is possible to calculate as long as you can invert $K$. So now in my case I'd alternate between optimizing the hyperparameters of $K$ to better estimate the covariance of the residuals, and altering my underlying model X to reduce the residuals (which I guess will also bring the points closer to the GP prior since it's zero mean?). So where I'm stuck is: is the posterior probability for $\mathbf{y}$ the same as the one in my first equation? Or am I totally lost?

(And: is there any way to speed up / avoid the calculation of $\mathbf{y}^T(K+\sigma_n I)^{-1}\mathbf{y}$? Remember I have millions of residuals...)

• Maybe I'll come back and write a full answer later, but: yes, it's reasonable to model residuals with a GP. This basically amounts to using the parametric model as your mean function. In terms of speeding up inference: since you're working with images, the points are probably aligned on a grid, so Kronecker-based inference methods are likely to work quite well. – Dougal Aug 22 '15 at 6:39
• @Dougal thanks for your comment. Thinking of the parametric model as the mean function does make sense - I guess I hadn't realized this was what the mean function was for. Is the likelihood of the GP what I want to maximize? I mean, I haven't seen elsewhere a justification for altering the mean as well as the hyperparameters to improve this likelihood. Also, I hadn't heard of Kronecker methods before but they sound really cool. Do they work well in practice? – cgreen Aug 23 '15 at 20:29
• The mean function is just another hyperparameter! Depending on how the parametric model works it might be easier to fit it first and then fit the residuals, but if you can optimize the combined likelihood, I'm not aware of any bad properties that come out of that. – Dougal Aug 24 '15 at 0:05