What is the shape of the decision surface of a Gaussian Process classifier?

I've got a binary classification problem, which I am trying to solve using a generative classifier. If I use Gaussian Discriminant Analysis, and fit two Gaussian distributions to my two classes, the decision surface between them will be a quadratic function (derivable analytically).

I guess a closed form solution will not be possible for Gaussian Processes, but I am still very interested in what can be said about the shape of the decision surface of GP classifiers in general? Will they be polynomial? Can something be said about their degree? Rasmussen & Williams briefly mention the expectation of the decision boundary function on page 45 of their book (http://www.gaussianprocess.org/gpml/chapters/RW3.pdf), but they never solve it...

I am using a standard Squared Exponential Kernel, but am interested in the decision surface shape arising from other kernels as well, if they are easier to derive.

As you can see in the example I crafted below, the probability surface in the case of the squared exponential (Gaussian) Kernel as covariance function for Gaussian processes looks like a smooth density.

A good read about covariance functions in the context of Gaussian processes is Chapter 4 - Covariance Functions [1]

So, I'll be using R, specifically the kernlab package for this example. There's a nice gausspr function that accepts different kernel types as parameters. I'll be using the iris dataset as a binary classification problem with only two dimensions.

library(kernlab)
data(iris)

#let's use only two variables plus the target variable so we can accurately plot it
data = iris[, 3:5]
data$Species[data$Species == "virginica"] = "versicolor"
data$Species = factor(data$Species)
levels(data$Species) = c("setosa", "virginica OR versicolor") #The fitting, you can change kernel and parameters, check the kernlab manual fit = gausspr(Species~., data = data, kernel = "rbfdot") pred = predict(fit, data) N = 250L #integer that gives the number of unique values in each dimension of the grid grid = expand.grid( Petal.Length = seq(min(data$Petal.Length), max(data$Petal.Length), length.out = N), Petal.Width = seq(min(data$Petal.Width), max(data$Petal.Width), length.out = N) ) pred.grid = predict(fit, grid, type = "probabilities")[, 1, drop = FALSE] #This maps the predictions a matrix representing the dimensions of data pred.grid = matrix(pred.grid, ncol = N) #The color part is thanks to http://www.r-bloggers.com/how-to-correctly-set-color-in-the-image-function/ collist<-c("#053061","#2166AC","#4393C3","#92C5DE","#D1E5F0","#F7F7F7","#FDDBC7","#F4A582","#D6604D","#B2182B","#67001F") ColorRamp<-colorRampPalette(collist)(100L) tiff(filename = "Rplot_rbfdot.tiff") image(unique(grid$Petal.Length), unique(grid$Petal.Width), pred.grid, useRaster = TRUE, col = ColorRamp, ylab = "Petal.Width", xlab = "Petal.Length", main = "kernel = \"rbfdot\"" ) points(data[,1:2], pch = c(16,17)[as.numeric(pred)], col = adjustcolor("black", alpha = 0.5)) contour(unique(grid$Petal.Length), unique(grid$Petal.Width), pred.grid, add = TRUE, levels = c(.4,.5,.6), labcex = 1, lwd = 1.75 ) legend("topleft", legend = levels(data$Species), pch = c(16,17), bg = "white")
dev.off()


These plots were produced with the above code.

[1] C. E. Rasmussen & C. K. I. Williams, Gaussian Processes for Machine Learning, the MIT Press, 2006, ISBN 026218253X. © 2006 Massachusetts Institute of Technology. www.GaussianProcess.org/gpml