Add i-times integrated Wiener Kernel

src/KernelFunctions.jl

@@ -9,7 +9,7 @@ export transform
 export duplicate, set! # Helpers
 
 export Kernel
-export ConstantKernel, WhiteKernel, EyeKernel, ZeroKernel
+export ConstantKernel, WhiteKernel, EyeKernel, ZeroKernel, WienerKernel
 export SqExponentialKernel, ExponentialKernel, GammaExponentialKernel
 export ExponentiatedKernel
 export MaternKernel, Matern32Kernel, Matern52Kernel
@@ -46,7 +46,8 @@ include("distances/dotproduct.jl")
 include("distances/delta.jl")
 include("transform/transform.jl")
 
-for k in ["exponential","matern","polynomial","constant","rationalquad","exponentiated","cosine","maha","fbm","gabor"]
+for k in ["exponential","matern","polynomial","constant","rationalquad","exponentiated",
+    "cosine","maha","fbm","gabor","wiener"]
     include(joinpath("kernels",k*".jl"))
 end
 include("kernels/transformedkernel.jl")

src/kernels/wiener.jl

@@ -0,0 +1,158 @@
+"""
+    WienerKernel{i}()
+
+i-times integrated Wiener process kernel function given by
+```julia
+    κ(x,y) =  kᵢ(x,y)
+```
+
+For i=-1, this is just the white noise covariance, see WhiteKernel.\\
+For i= 0, this is the Wiener process covariance,\\
+for i= 1, this is the integrated Wiener process covariance (velocity),\\
+for i= 2, this is the twice-integrated Wiener process covariance (accel.),\\
+for i= 3, this is the thrice-integrated Wiener process covariance. where `kᵢ` is given by\\
+
+```julia
+    k₋₁(x,y) =  δ(x,y)
+    i >= 0, kᵢ(x,y) = 1/ai * min(x,y)^(2i + 1) + bi * min(x,y)^(i+1) * |x-y| * ri(x,y),
+    with the coefficients ai, bi and the residual ri(x,y) defined as follows:
+        i = 0, ai =   1, bi = 0
+        i = 1, ai =   3, bi = 1/  2, ri(x,y) = 1
+        i = 2, ai =  20, bi = 1/ 12, ri(x,y) = x + y - 1/2 * min(x,y)
+        i = 3, ai = 252, bi = 1/720, ri(x,y) = 5 * max(x,y)² + 2xz + 3min(x,y)²
+```
+
+**References:**\\
+See the paper *Probabilistic ODE Solvers with Runge-Kutta Means* by Schober, Duvenaud and Hennig, NIPS, 2014, for more details.
+
+"""
+struct WienerKernel{I} <: BaseKernel
+    function WienerKernel{I}() where I
+        I in (-1, 0, 1, 2, 3) || error("Invalid paramter i=$(I). Should be -1, 0, 1, 2 or 3.")
+        if I==-1
+            return WhiteKernel()
+        end
+        new{I}()
+    end
+end
+
+function WienerKernel(;i=0)
+    return WienerKernel{i}()
+end
+
+function _wiener(κ::WienerKernel{I},x,y) where I
+    X = sqrt(sum(abs2.(x)))
+    Y = sqrt(sum(abs2.(y)))
+    minXY = min(X,Y)
+    if I==0
+        return         minXY^(2I + 1)
+    elseif I==1
+        return 1/3   * minXY^(2I + 1) + 1/2   * minXY^(I+1) * euclidean(x,y)
+    elseif I==2
+        return 1/20  * minXY^(2I + 1) + 1/12  * minXY^(I+1) * euclidean(x,y) * (X + Y - 1/2 * minXY)
+    elseif I==3
+        return 1/252 * minXY^(2I + 1) + 1/720 * minXY^(I+1) * euclidean(x,y) * (5*max(X,Y)^2 + 2*X*Y + 3 * minXY^2)
+    else
+        error("Invalid I")
+    end
+end
+
+function kappa(κ::WienerKernel, x,y)
+    return _wiener(κ, x, y)
+end
+
+(κ::WienerKernel)(x::Real, y::Real) = kappa(κ,x,y)
+
+function _kernel(
+    κ::WienerKernel,
+    x::AbstractVector,
+    y::AbstractVector;
+    obsdim::Int = defaultobs
+)
+    @assert length(x) == length(y) "x and y don't have the same dimension!"
+    kappa(κ,x,y)
+end
+
+function kernelmatrix!(
+    K::AbstractMatrix,
+    κ::WienerKernel,
+    X::AbstractMatrix,
+    Y::AbstractMatrix;
+    obsdim::Int = defaultobs
+)
+    @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))"
+    if !check_dims(K,X,X,feature_dim(obsdim),obsdim)
+        throw(DimensionMismatch("Dimensions of the target array K $(size(K)) are not consistent with X $(size(X))"))
+    end
+
+    if obsdim == 1
+        for j = 1:size(K,2)
+            for i = 1:size(K,1)
+                @inbounds @views K[i,j] = _kernel(κ, X[i,:],Y[j,:])
+            end
+        end
+    else
+        for j = 1:size(K,2)
+            for i = 1:size(K,1)
+                @inbounds @views K[i,j] = _kernel(κ, X[:,i],Y[:,j])
+            end
+        end
+    end
+    return K
+end
+
+function kernelmatrix(
+    κ::WienerKernel,
+    X::AbstractMatrix,
+    Y::AbstractMatrix;
+    obsdim::Int = defaultobs
+)
+    @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))"
+    if !check_dims(X,Y,feature_dim(obsdim),obsdim)
+        throw(DimensionMismatch("X $(size(X)) and Y $(size(Y)) do not have the same number of features on the dimension : $(feature_dim(obsdim))"))
+    end
+    if obsdim == 1
+        outdim = size(X,1)
+    else
+        outdim = size(X,2)
+    end
+    K = zeros(outdim,outdim)
+    if obsdim == 1
+        for j = 1:size(K,2)
+            for i = 1:size(K,1)
+                @inbounds @views K[i,j] = _kernel(κ, X[i,:],Y[j,:])
+            end
+        end
+    else
+        for j = 1:size(K,2)
+            for i = 1:size(K,1)
+                @inbounds @views K[i,j] = _kernel(κ, X[:,i],Y[:,j])
+            end
+        end
+    end
+    return K
+end
+
+function kernelmatrix!(
+    K::AbstractMatrix,
+    κ::WienerKernel,
+    X::AbstractMatrix;
+    obsdim::Int = defaultobs
+)
+    @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))"
+    if !check_dims(K,X,X,feature_dim(obsdim),obsdim)
+        throw(DimensionMismatch("Dimensions of the target array K $(size(K)) are not consistent with X $(size(X))"))
+    end
+    kernelmatrix!(K,κ,X,X;obsdim=obsdim)
+    return K
+end
+
+function kernelmatrix(
+    κ::WienerKernel,
+    X::AbstractMatrix;
+    obsdim::Int = defaultobs
+)
+    return kernelmatrix(κ,X,X;obsdim=obsdim)
+end
+
+Base.show(io::IO, κ::WienerKernel{I}) where I = print(io, "Wiener Kernel $(I)-times integrated")

test/kernels/wiener.jl

@@ -0,0 +1,46 @@
+@testset "wiener" begin
+    k_1 = WienerKernel(i=-1)
+    @test typeof(k_1) <: WhiteKernel
+
+    k0 = WienerKernel()
+    @test typeof(k0) <: WienerKernel{0}
+
+    k1 = WienerKernel(i=1)
+    @test typeof(k1) <: WienerKernel{1}
+
+    k2 = WienerKernel(i=2)
+    @test typeof(k2) <: WienerKernel{2}
+
+    k3 = WienerKernel(i=3)
+    @test typeof(k3) <: WienerKernel{3}
+
+    @test_throws ErrorException WienerKernel(i=4)
+    @test_throws ErrorException WienerKernel(i=-2)
+
+    v1 = rand(3); v2 = rand(3)
+    @test k0(v1,v2) ≈ kappa(k0,v1,v2)
+    @test k1(v1,v2) ≈ kappa(k1,v1,v2)
+    @test k2(v1,v2) ≈ kappa(k2,v1,v2)
+    @test k3(v1,v2) ≈ kappa(k3,v1,v2)
+
+    # kernelmatrix tests
+    m1 = rand(3,4)
+    m2 = rand(3,4)
+    @test kernelmatrix(k0, m1, m1) ≈ kernelmatrix(k0, m1) atol=1e-5
+    @test kernelmatrix(k0, m1, m2) ≈ k0(m1, m2) atol=1e-5
+
+    K = zeros(4,4)
+    kernelmatrix!(K,k0,m1,m2)
+    @test K ≈ kernelmatrix(k0, m1, m2) atol=1e-5
+
+    V = zeros(4)
+    kerneldiagmatrix!(V,k0,m1)
+    @test V ≈ kerneldiagmatrix(k0,m1) atol=1e-5
+
+    x1 = rand()
+    x2 = rand()
+    @test kernelmatrix(k0, x1*ones(1,1), x2*ones(1,1))[1] ≈ k0(x1, x2) atol=1e-5
+    @test kernelmatrix(k1, x1*ones(1,1), x2*ones(1,1))[1] ≈ k1(x1, x2) atol=1e-5
+    @test kernelmatrix(k2, x1*ones(1,1), x2*ones(1,1))[1] ≈ k2(x1, x2) atol=1e-5
+    @test kernelmatrix(k3, x1*ones(1,1), x2*ones(1,1))[1] ≈ k3(x1, x2) atol=1e-5
+end

test/runtests.jl

@@ -63,6 +63,7 @@ using KernelFunctions: metric
 
     @testset "kernels" begin
         include(joinpath("kernels", "constant.jl"))
+        include(joinpath("kernels", "wiener.jl"))
         include(joinpath("kernels", "cosine.jl"))
         include(joinpath("kernels", "exponential.jl"))
         include(joinpath("kernels", "exponentiated.jl"))