Add NeuralNetOneKernel

src/KernelFunctions.jl

@@ -18,7 +18,7 @@ export MaternKernel, Matern32Kernel, Matern52Kernel
 export LinearKernel, PolynomialKernel
 export RationalQuadraticKernel, GammaRationalQuadraticKernel
 export MahalanobisKernel, GaborKernel, PiecewisePolynomialKernel
-export PeriodicKernel
+export PeriodicKernel, NeuralNetworkKernel
 export KernelSum, KernelProduct
 export TransformedKernel, ScaledKernel
 export TensorProduct

src/basekernels/nn.jl

@@ -0,0 +1,26 @@
+"""
+    NeuralNetworkKernel()
+
+Neural network kernel function.
+
+```math
+    κ(x, y) =  asin(x' * y / sqrt[(1 + x' * x) * (1 + y' * y)])
+```
+# Significance
+Neal (1996) pursued the limits of large models, and showed that a Bayesian neural network
+becomes a Gaussian process with a **neural network kernel** as the number of units
+approaches infinity. Here, we give the neural network kernel for single hidden layer
+Bayesian neural network with erf (Error Function) as activation function.
+
+# References:
+- [GPML Pg 105](http://www.gaussianprocess.org/gpml/chapters/RW4.pdf)
+- [Neal(1996)](https://www.cs.toronto.edu/~radford/bnn.book.html)
+- [Andrew Gordon's Thesis Pg 45](http://www.cs.cmu.edu/~andrewgw/andrewgwthesis.pdf)
+"""
+struct NeuralNetworkKernel <: BaseKernel end
+
+function (κ::NeuralNetworkKernel)(x, y)
+    return asin(dot(x, y) / sqrt((1 + sum(abs2, x)) * (1 + sum(abs2, y))))
+end
+
+Base.show(io::IO, κ::NeuralNetworkKernel) = print(io, "Neural Network Kernel")

test/basekernels/nn.jl

@@ -0,0 +1,47 @@
+@testset "nn" begin
+    using LinearAlgebra
+    k = NeuralNetworkKernel()
+    v1 = rand(3); v2 = rand(3)
+    @test k(v1,v2) ≈ asin(v1' * v2 / sqrt((1 + v1' * v1) * (1 + v2' * v2))) atol=1e-5
+
+    # kernelmatrix tests
+    m1 = rand(3,4)
+    m2 = rand(3,4)
+    @test kernelmatrix(k, m1, m1) ≈ kernelmatrix(k, m1) atol=1e-5
+    @test_broken kernelmatrix(k, m1, m2) ≈ k(m1, m2) atol=1e-5
+
+
+    x1 = rand()
+    x2 = rand()
+    @test kernelmatrix(k, x1*ones(1,1), x2*ones(1,1))[1] ≈ k(x1, x2) atol=1e-5
+
+    @test k(v1, v2) ≈ k(v1, v2) atol=1e-5
+    @test typeof(k(v1, v2)) <: Real
+
+    @test_broken size(k(m1, m2)) == (4, 4)
+    @test_broken size(k(m1)) == (4, 4)
+
+    A1 = ones(4, 4)
+    kernelmatrix!(A1, k, m1, m2)
+    @test A1 ≈ kernelmatrix(k, m1, m2) atol=1e-5
+
+    A2 = ones(4, 4)
+    kernelmatrix!(A2, k, m1)
+    @test A2 ≈ kernelmatrix(k, m1) atol=1e-5
+
+    @test size(kerneldiagmatrix(k, m1)) == (4,)
+    A3 = kernelmatrix(k, m1)
+    @test kerneldiagmatrix(k, m1) ≈ [A3[i, i] for i in 1:LinearAlgebra.checksquare(A3)] atol=1e-5
+
+    A4 = ones(4)
+    kerneldiagmatrix!(A4, k, m1)
+    @test kerneldiagmatrix(k, m1) ≈ A4 atol=1e-5
+
+    A5 = ones(4,4)
+    @test_throws AssertionError kernelmatrix!(A5, k, m1, m2, obsdim=3)
+    @test_throws AssertionError kernelmatrix!(A5, k, m1, obsdim=3)
+    @test_throws DimensionMismatch kernelmatrix!(A5, k, ones(4,3), ones(3,4))
+
+    @test k([x1], [x2]) ≈ k(x1, x2) atol=1e-5
+
+end

test/runtests.jl

@@ -71,6 +71,7 @@ using KernelFunctions: metric, kappa
         include(joinpath("basekernels", "gabor.jl"))
         include(joinpath("basekernels", "maha.jl"))
         include(joinpath("basekernels", "matern.jl"))
+        include(joinpath("basekernels", "nn.jl"))
         include(joinpath("basekernels", "periodic.jl"))
         include(joinpath("basekernels", "polynomial.jl"))
         include(joinpath("basekernels", "piecewisepolynomial.jl"))