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,37 @@
+"""
+
+    GibbsKernel(x, y)
+
+# Definition
+
+The Gibbs kernel is non-stationary generalisation of the Squared-Exponential
+kernel. The lengthscale parameter ``\\ell`` becomes a function of
+position ``\\ell(x)``.
+
+``\\ell(x) = 1.`` then you recover the standard Squared-Exponential kernel
+with constant lengthscale.
+
+```math
+k(x, y) = \\sqrt{ \\left( \\frac{2 \\ell(x) \\ell(y)}{\\ell(x)^2 + \\ell(y)^2} \\right) } \\quad \\rm{exp} \\left( - \\frac{(x - y)^2}{\\ell(x)^2 + \\ell(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
+"""
+
+#function (k::GibbsKernel)(x, y, ell)
+function GibbsKernel(x::Any, y::Any, ell::Function)
+    fac = ell(x)^2 + ell(y)^2
+    term1 = 2.0 * ell(x) * ell(y) / fac
+    term2 = - (x-y)^2 / fac
+    return term1 * exp(term2)
+end

test/kernels/gibbskernel.jl

@@ -0,0 +1,13 @@
+using Test
+using KernelFunctions
+
+@testset "gibbskernel" begin
+    x = 1.0
+    y = 2.0
+
+    k = SqExponentialKernel()
+
+    ell(x) = 1.0
+
+    @test GibbsKernel(x, y, ell) == k(x, y)
+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