Merge pull request #83 from theogf/general_kernelmatrix

[WIP] Rework on kernelmatrix to work with Vectors and more complex kernels

Author: theogf

src/KernelFunctions.jl

@@ -10,6 +10,7 @@ export duplicate, set! # Helpers
 
 export Kernel
 export ConstantKernel, WhiteKernel, EyeKernel, ZeroKernel
+export CosineKernel
 export SqExponentialKernel, RBFKernel, GaussianKernel, SEKernel
 export LaplacianKernel, ExponentialKernel, GammaExponentialKernel
 export ExponentiatedKernel
@@ -43,6 +44,7 @@ Abstract type defining a slice-wise transformation on an input matrix
 abstract type Transform end
 abstract type Kernel end
 abstract type BaseKernel <: Kernel end
+abstract type SimpleKernel <: BaseKernel end
 
 include("utils.jl")
 include("distances/dotproduct.jl")

src/basekernels/constant.jl

@@ -7,7 +7,7 @@ Create a kernel that always returning zero
 ```
 The output type depends of `x` and `y`
 """
-struct ZeroKernel <: BaseKernel end
+struct ZeroKernel <: SimpleKernel end
 
 kappa(κ::ZeroKernel, d::T) where {T<:Real} = zero(T)
 
@@ -24,7 +24,7 @@ Base.show(io::IO, ::ZeroKernel) = print(io, "Zero Kernel")
 ```
 Kernel function working as an equivalent to add white noise. Can also be called via `EyeKernel()`
 """
-struct WhiteKernel <: BaseKernel end
+struct WhiteKernel <: SimpleKernel end
 
 """
     EyeKernel()
@@ -48,7 +48,7 @@ Kernel function always returning a constant value `c`
     κ(x,y) = c
 ```
 """
-struct ConstantKernel{Tc<:Real} <: BaseKernel
+struct ConstantKernel{Tc<:Real} <: SimpleKernel
     c::Vector{Tc}
     function ConstantKernel(;c::T=1.0) where {T<:Real}
         new{T}([c])

src/basekernels/cosine.jl

@@ -6,7 +6,7 @@ The cosine kernel is a stationary kernel for a sinusoidal given by
     κ(x,y) = cos( π * (x-y) )
 ```
 """
-struct CosineKernel <: BaseKernel end
+struct CosineKernel <: SimpleKernel end
 
 kappa(κ::CosineKernel, d::Real) = cospi(d)
 metric(::CosineKernel) = Euclidean()

src/basekernels/exponential.jl

@@ -9,7 +9,7 @@ Can also be called via `SEKernel`, `GaussianKernel` or `SEKernel`.
 See also [`ExponentialKernel`](@ref) for a
 related form of the kernel or [`GammaExponentialKernel`](@ref) for a generalization.
 """
-struct SqExponentialKernel <: BaseKernel end
+struct SqExponentialKernel <: SimpleKernel end
 
 kappa(κ::SqExponentialKernel, d²::Real) = exp(-d²)
 iskroncompatible(::SqExponentialKernel) = true
@@ -30,7 +30,7 @@ The exponential kernel is a Mercer kernel given by the formula:
     κ(x,y) = exp(-‖x-y‖)
 ```
 """
-struct ExponentialKernel <: BaseKernel end
+struct ExponentialKernel <: SimpleKernel end
 
 kappa(κ::ExponentialKernel, d::Real) = exp(-d)
 iskroncompatible(::ExponentialKernel) = true
@@ -51,7 +51,7 @@ The γ-exponential kernel is an isotropic Mercer kernel given by the formula:
 Where `γ > 0`, (the keyword `γ` can be replaced by `gamma`)
 For `γ = 1`, see `SqExponentialKernel` and `γ = 0.5`, see `ExponentialKernel`
 """
-struct GammaExponentialKernel{Tγ<:Real} <: BaseKernel
+struct GammaExponentialKernel{Tγ<:Real} <: SimpleKernel
     γ::Vector{Tγ}
     function GammaExponentialKernel(; gamma::T=2.0, γ::T=gamma) where {T<:Real}
         @check_args(GammaExponentialKernel, γ, γ >= zero(T), "γ > 0")

src/basekernels/exponentiated.jl

@@ -6,7 +6,7 @@ The exponentiated kernel is a Mercer kernel given by:
     κ(x,y) = exp(xᵀy)
 ```
 """
-struct ExponentiatedKernel <: BaseKernel end
+struct ExponentiatedKernel <: SimpleKernel end
 
 kappa(κ::ExponentiatedKernel, xᵀy::Real) = exp(xᵀy)
 metric(::ExponentiatedKernel) = DotProduct()

src/basekernels/fbm.jl

@@ -66,17 +66,6 @@ function kernelmatrix!(
     return K
 end
 
-## Apply kernel on two vectors ##
-function _kernel(
-        κ::FBMKernel,
-        x::AbstractVector,
-        y::AbstractVector;
-        obsdim::Int = defaultobs
-    )
-    @assert length(x) == length(y) "x and y don't have the same dimension!"
-    return kappa(κ, x, y)
-end
-
 function kappa(κ::FBMKernel, x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
     modX = sum(abs2, x)
     modY = sum(abs2, y)

src/basekernels/maha.jl

@@ -8,7 +8,7 @@ Mahalanobis distance-based kernel given by
 where the matrix P is the metric.
 
 """
-struct MahalanobisKernel{T<:Real, A<:AbstractMatrix{T}} <: BaseKernel
+struct MahalanobisKernel{T<:Real, A<:AbstractMatrix{T}} <: SimpleKernel
     P::A
     function MahalanobisKernel(P::AbstractMatrix{T}) where {T<:Real}
         LinearAlgebra.checksquare(P)

src/basekernels/matern.jl

@@ -7,7 +7,7 @@ The matern kernel is a Mercer kernel given by the formula:
 ```
 For `ν=n+1/2, n=0,1,2,...` it can be simplified and you should instead use [`ExponentialKernel`](@ref) for `n=0`, [`Matern32Kernel`](@ref), for `n=1`, [`Matern52Kernel`](@ref) for `n=2` and [`SqExponentialKernel`](@ref) for `n=∞`.
 """
-struct MaternKernel{Tν<:Real} <: BaseKernel
+struct MaternKernel{Tν<:Real} <: SimpleKernel
     ν::Vector{Tν}
     function MaternKernel(;nu::T=1.5, ν::T=nu) where {T<:Real}
         @check_args(MaternKernel, ν, ν > zero(T), "ν > 0")
@@ -37,7 +37,7 @@ The matern 3/2 kernel is a Mercer kernel given by the formula:
     κ(x,y) = (1+√(3)‖x-y‖)exp(-√(3)‖x-y‖)
 ```
 """
-struct Matern32Kernel <: BaseKernel end
+struct Matern32Kernel <: SimpleKernel end
 
 kappa(κ::Matern32Kernel, d::Real) = (1 + sqrt(3) * d) * exp(-sqrt(3) * d)
 metric(::Matern32Kernel) = Euclidean()
@@ -52,7 +52,7 @@ The matern 5/2 kernel is a Mercer kernel given by the formula:
     κ(x,y) = (1+√(5)‖x-y‖ + 5/3‖x-y‖^2)exp(-√(5)‖x-y‖)
 ```
 """
-struct Matern52Kernel <: BaseKernel end
+struct Matern52Kernel <: SimpleKernel end
 
 kappa(κ::Matern52Kernel, d::Real) = (1 + sqrt(5) * d + 5 * d^2 / 3) * exp(-sqrt(5) * d)
 metric(::Matern52Kernel) = Euclidean()

src/basekernels/periodic.jl

@@ -8,7 +8,7 @@ Periodic Kernel as described in http://www.inference.org.uk/mackay/gpB.pdf eq. 4
     κ(x,y) = exp( - 0.5 sum_i(sin (π(x_i - y_i))/r_i))
 ```
 """
-struct PeriodicKernel{T} <: BaseKernel
+struct PeriodicKernel{T} <: SimpleKernel
     r::Vector{T}
     function PeriodicKernel(; r::AbstractVector{T} = ones(Float64, 1)) where {T<:Real}
         @assert all(r .> 0)

src/basekernels/piecewisepolynomial.jl

@@ -10,12 +10,14 @@ processes are hence v times  mean-square differentiable. The kernel function is:
 where `r` is the Mahalanobis distance mahalanobis(x,y) with `maha` as the metric.
 
 """
-struct PiecewisePolynomialKernel{V, A<:AbstractMatrix{<:Real}} <: BaseKernel
+struct PiecewisePolynomialKernel{V, A<:AbstractMatrix{<:Real}} <: SimpleKernel
     maha::A
+    j::Int
     function PiecewisePolynomialKernel{V}(maha::AbstractMatrix{<:Real}) where V
         V in (0, 1, 2, 3) || error("Invalid paramter v=$(V). Should be 0, 1, 2 or 3.")
         LinearAlgebra.checksquare(maha)
-        return new{V,typeof(maha)}(maha)
+        j = div(size(maha, 1), 2) + V + 1
+        return new{V,typeof(maha)}(maha, j)
     end
 end
 
@@ -29,78 +31,7 @@ _f(κ::PiecewisePolynomialKernel{2}, r, j) = 1 + (j + 2) * r + (j^2 + 4 * j + 3)
 _f(κ::PiecewisePolynomialKernel{3}, r, j) = 1 + (j + 3) * r +
     (6 * j^2 + 36j + 45) / 15 * r.^2 + (j^3 + 9 * j^2 + 23j + 15) / 15 * r.^3
 
-function _piecewisepolynomial(κ::PiecewisePolynomialKernel{V}, r, j) where V
-    return max(1 - r, 0)^(j + V) * _f(κ, r, j)
-end
-
-function kappa(
-    κ::PiecewisePolynomialKernel{V},
-    x::AbstractVector{<:Real},
-    y::AbstractVector{<:Real},
-) where {V}
-    r = evaluate(metric(κ), x, y)
-    j = div(size(x, 2), 1) + V + 1
-    return _piecewisepolynomial(κ, r, j)
-end
-
-function _kernel(
-    κ::PiecewisePolynomialKernel,
-    x::AbstractVector,
-    y::AbstractVector;
-    obsdim::Int = defaultobs,
-)
-    @assert length(x) == length(y) "x and y don't have the same dimension!"
-    return kappa(κ,x,y)
-end
-
-function kernelmatrix(
-    κ::PiecewisePolynomialKernel{V},
-    X::AbstractMatrix;
-    obsdim::Int = defaultobs
-) where {V}
-    j = div(size(X, feature_dim(obsdim)), 2) + V + 1
-    return map(r->_piecewisepolynomial(κ, r, j), pairwise(metric(κ), X; dims=obsdim))
-end
-
-function _kernelmatrix(κ::PiecewisePolynomialKernel{V}, X, Y, obsdim) where {V}
-    j = div(size(X, feature_dim(obsdim)), 2) + V + 1
-    return map(r->_piecewisepolynomial(κ, r, j), pairwise(metric(κ), X, Y; dims=obsdim))
-end
-
-function kernelmatrix!(
-    K::AbstractMatrix,
-    κ::PiecewisePolynomialKernel{V},
-    X::AbstractMatrix;
-    obsdim::Int = defaultobs
-) where {V}
-    @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
-    j = div(size(X, feature_dim(obsdim)), 2) + V + 1
-    return map!(r->_piecewisepolynomial(κ,r,j), K, pairwise(metric(κ), X; dims=obsdim))
-end
-
-function kernelmatrix!(
-    K::AbstractMatrix,
-    κ::PiecewisePolynomialKernel{V},
-    X::AbstractMatrix,
-    Y::AbstractMatrix;
-    obsdim::Int = defaultobs,
-) where {V}
-    @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))"
-    if !check_dims(K, X, Y, feature_dim(obsdim), obsdim)
-        throw(DimensionMismatch(
-            "Dimensions $(size(K)) of the target array K are not consistent with X " *
-            "($(size(X))) and Y ($(size(Y)))",
-        ))
-    end
-    j = div(size(X, feature_dim(obsdim)), 2) + V + 1
-    return map!(r->_piecewisepolynomial(κ,r,j), K, pairwise(metric(κ), X, Y; dims=obsdim))
-end
+kappa(κ::PiecewisePolynomialKernel{V}, r) where V = max(1 - r, 0)^(κ.j + V) * _f(κ, r, κ.j)
 
 metric(κ::PiecewisePolynomialKernel) = Mahalanobis(κ.maha)