I am looking to compute the covariance matrix of an inferred Gaussian process in R. Below I outline how I would do this manually, but I realize that the kernlab package provides the function gausspr, which would be preferable for many reasons. Unfortunately, while I see how to get the expected value (predict.gausspr) using the package, I do not see how to get the associated covariance matrix.

To provide a minimal working example for comparison, I loosely follow Rasmussen & Williams (2006) Chapter 2:

Consider we have the observed x,y points and x values where we desire predicted y values:

obs <- data.frame(x = c(-4, -3, -1,  0,  2),
                  y = c(-2,  0,  1,  2, -1))
x_predict <- seq(-5,5,len=50)

Use a radial basis kernel:

SE <- function(Xi,Xj, l=1) exp(-0.5 * (Xi - Xj) ^ 2 / l ^ 2)
cov <- function(X, Y) outer(X, Y, SE)

Calculate mean and covariance matrix:

sigma.n <- 0.3
cov_xx_inv <- solve(cov(obs$x, obs$x) + sigma.n^2 * diag(1, length(obs$x)))
Ef <- cov(x_predict, obs$x) %*% cov_xx_inv %*% obs$y
Cf <- cov(x_predict, x_predict) - cov(x_predict, obs$x)  %*% cov_xx_inv %*% cov(obs$x, x_predict)

kernlab approach

I think I see how I get the equivalent expected values in kernlab:

gp <- gausspr(obs$x, obs$y, kernel="rbfdot", kpar=list(sigma=1/(2*l^2)), fit=FALSE, scaled=FALSE, var=0.8)
Ef_kernlab <- predict(gp, x_predict)

but I don't see how I get the associated covariances? Can the covariance matrix be reconstructed by knowing alpha?

Ideally there is a solution that doesn't involve explicitly inverting a matrix again, since that inverse was already calculated in determining the alpha values during the fit done by gausspr...

There are many cases where it would be nice to have access to the resulting Gaussian process, such as in generating plots as in Rasmussen and Williams:

dat <- data.frame(x=x_predict, y=(Ef), ymin=(Ef-2*sqrt(diag(Cf))), ymax=(Ef+2*sqrt(diag(Cf))))    
ggplot(dat) +
  geom_ribbon(aes(x=x,y=y, ymin=ymin, ymax=ymax), fill="grey80") + # Var
  geom_line(aes(x=x,y=y), size=1) + #MEAN
  geom_point(data=obs,aes(x=x,y=y)) +  #OBSERVED DATA
  scale_y_continuous(lim=c(-3,3), name="output, f(x)") + xlab("input, x")


1 Answer 1


I think latest package has this feature. Using the option "variance.model = TRUE" in gausspr function can generate the variance at each point. Document for kernlab elaborates this.

  • 3
    $\begingroup$ Welcome to CV.stackexchange. Maybe you could demonstrate with the given example on question how this new version of gausspr function would provide the posterior covariance matrix returned by kernlab. $\endgroup$ May 8, 2013 at 20:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.