How do I access or compute the posterior covariance matrix returned by kernlab::gausspr R function?

Question

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:

library(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:

require(ggplot2)
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")

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.