Adding Gibbskernel

src/KernelFunctions.jl

@@ -17,6 +17,7 @@ export PiecewisePolynomialKernel
 export PeriodicKernel, NeuralNetworkKernel
 export KernelSum, KernelProduct, KernelTensorProduct
 export TransformedKernel, ScaledKernel, NormalizedKernel
+export GibbsKernel
 
 export Transform,
     SelectTransform,
@@ -94,6 +95,7 @@ include("basekernels/rational.jl")
 include("basekernels/sm.jl")
 include("basekernels/wiener.jl")
 
+include("kernels/gibbskernel.jl")
 include("kernels/scaledkernel.jl")
 include("kernels/normalizedkernel.jl")
 include("matrix/kernelmatrix.jl")

src/kernels/gibbskernel.jl

@@ -0,0 +1,44 @@
+"""
+    GibbsKernel(x, y)
+
+# Definition
+
+The Gibbs kernel is non-stationary generalisation of the Squared-Exponential
+kernel. The lengthscale parameter ``l`` becomes a function of
+position ``l(x)``.
+
+``l(x) = 1.`` then you recover the standard Squared-Exponential kernel
+with constant lengthscale.
+
+```math
+k(x, y) = \\sqrt{ \\left(\\frac{2 l(x) l(y)}{l(x)^2 + l(y)^2} \\right) }
+\\quad \\rm{exp} \\left( - \\frac{(x - y)^2}{l(x)^2 + l(y)^2} \\right)
+```
+
+[1] - Mark N. Gibbs. "Bayesian Gaussian Processes for Regression and Classication."
+    PhD thesis, 1997
+
+[2] - Christopher J. Paciorek and Mark J. Schervish. "Nonstationary Covariance
+    Functions for Gaussian Process Regression". NEURIPS, 2003
+
+[3] - Sami Remes, Markus Heinonen, Samuel Kaski.
+    "Non-Stationary Spectral Kernels". arXiV:1705.08736, 2017
+
+[4] - Sami Remes, Markus Heinonen, Samuel Kaski.
+    "Neural Non-Stationary Spectral Kernel". arXiv:1811.10978, 2018
+"""
+
+struct GibbsKernel{T} <: Kernel
+    lengthscale::T
+end
+
+GibbsKernel(; lengthscale) = GibbsKernel(lengthscale)
+
+function (k::GibbsKernel)(x, y)
+    lengthscale = k.lengthscale
+    lx = lengthscale(x)
+    ly = lengthscale(y)
+    l = hypot(lx, ly)
+    kernel = (sqrt(2 * lx * ly) / l) * with_lengthscale(SqExponentialKernel(), l)
+    return kernel(x, y)
+end

test/kernels/gibbskernel.jl

@@ -0,0 +1,32 @@
+@testset "gibbskernel" begin
+    x = randn(2)
+    y = randn(2)
+
+    # generate random number of standard SqExponentialKernel lengthscale
+    # taking exp() to ensure ell is positive
+    ell(x) = exp(sum(sin, x))
+
+    # this is the gibbs lengthscale function.
+    # to check that we can recover the stationary SqExponentialKernel
+    # with a constant lengthscale we just set this function
+    # equal to a constant.
+    l_func(x) = ell(1)
+
+    # in order to compare to the Gibbs kernel we compute the
+    # equivalent lengthscale l(x)^2 + l(y)^2.
+    # See the denominator of the exponential term in the gibbs kernel.
+    lengthscale = hypot(l_func(x), l_func(y))
+
+    # create a SqExponentialKernel with a lengthscale
+    k = with_lengthscale(SqExponentialKernel(), lengthscale)
+
+    # create a gibbs kernel with our constant lengthscale function
+    k_gibbs = GibbsKernel(l_func)
+
+    # check they are equal
+    @test k_gibbs(x, y) == k(x, y)
+
+    k_gibbs = GibbsKernel(ell)
+    @test k_gibbs(x, y) ≈ sqrt((2 * ell(x) * ell(y)) / (ell(x)^2 + ell(y)^2)) * exp(- norm(x - y)^2 / (ell(x)^2 + ell(y)^2))
+
+end

test/runtests.jl

@@ -126,6 +126,7 @@ include("test_utils.jl")
             include("kernels/transformedkernel.jl")
             include("kernels/normalizedkernel.jl")
             include("kernels/neuralkernelnetwork.jl")
+            include("kernels/gibbskernel.jl")
         end
         @info "Ran tests on Kernel"
     end