Is my Gaussian Process implementation working correctly?

Question

I'm using the standard kernel: $K(x_1, x_2) = \exp(-\frac{1}{2} (\frac{(x_1-x_2)}{l})^2)$, with $l=1.0$.

I picked some arbitrary function, and then picked 15 points from which to learn from. Here's the graph I get when I plot those points, the real underlying function, and the learned conditional distributions confidence intervals.

I'm surprised that the confidence intervals are so "wobbly" and don't hug the real function more closely. Is this the normal behaviour? Am I getting these results because I'm not actually learning the value $l$ in the kernel and just plugging in $l=1.0$? Or is there possibly some bug in my code?

Answer 1

I cooked this up really quick, so beware. But the results appear to be consistent with yours. If you think about it, the non-convex shapes make sense because they're relative to the predicted mean, and because the mean is estimated very near zero and sigma is small, this model will tend towards zero in between the data. Your plot omits the predicted means, which is why the behavior appears to be counter-intuitive. Gaussian process models have the property that they will interpolate exactly in the absence of signal noise, so that's why the predictions have the U shape.

Playing with the sigma value will alter this effect, because it changes the behavior of the mean function between the observed data. Larger values of sigma produce dramatically different results. Below is a plot for sigma = 3.0.

The Gaussian process's interpolation property is why the confidence bands shrink to 0 where the data are observed. This is true for these models and all GP models which do not have observation noise, and omit a noise parameter. (To add a noise term, set eta to a positive value.)

gaussian_kernel <- function(x, y, sigma){
    d_sq <- sum( (x - y) ^ 2)
    out <- exp(- 0.5 * d_sq / (sigma ^ 2))
    return(out)
}

get_K <- function(X, Y = "None", sigma=1.0){
    check <- all(Y == "None")
    if (check){
        size <- nrow(X)
        K <- matrix(1.0, nrow=size, ncol=size)
        for(i in 1:(size - 1) ){
            for(j in ((i + 1):size)){
                new_k <- gaussian_kernel(X[i,], X[j,], sigma=sigma)
                K[i,j] <- new_k
                K[j,i] <- new_k
            }
        }
    }
    else{
        size1 <- nrow(X)
        size2 <- nrow(Y)
        K <- matrix(0.0, size1, size2)
        for(i in 1:size1){
            for(j in 1:size2){
                K[i,j] <- gaussian_kernel(X[i,], Y[j,], sigma=sigma)
            }
        }
    }
    return(K)
}

eta <- 0.0

x <- matrix(seq(-60,60,length.out=15), ncol=1)
y <- 3 * sin(0.7 * x)
mu <- mean(y)

x_tilde <- matrix(seq(-60,60,length.out=1024)), ncol=1)

K <- get_K(x, sigma=1)
K_star <- get_K(x_tilde, x, sigma=1)
K_star_star <- get_K(x_tilde, sigma=1)

eye <- diag(nrow(x))

K_qr <- qr(K + eta * eye)
K_R <- qr.R(K_qr)
K_Q <- qr.Q(K_qr)


y_hat <- mu + K_star %*% backsolve(K_R, crossprod(K_Q, y - mu))
cov <- K_star_star - K_star %*% backsolve(K_R, t(K_Q) %*% t(K_star))

y_ci <- 2 * sqrt(diag(cov))

y_ci_upper <- y_hat + y_ci
y_ci_lower <- y_hat - y_ci

png("cv_gp.png", height=6, width=6, units="in", res=300)
    plot(x, y, xlim=c(-60,60), ylim=c(-4,4))
    x_poly <- c(rev(x_tilde), x_tilde)
    y_poly <- c(rev(y_ci_lower), y_ci_upper)
    polygon(x_poly, y_poly, col="grey")
    lines(x_tilde, y_hat, col="red", lwd=2)
    points(x,y,lwd=2)
dev.off()

To make these diagrams even more interesting, try randomly sampling the observed data.