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)
 

src/basekernels/polynomial.jl

@@ -7,7 +7,7 @@ The linear kernel is a Mercer kernel given by
 ```
 Where `c` is a real number
 """
-struct LinearKernel{Tc<:Real} <: BaseKernel
+struct LinearKernel{Tc<:Real} <: SimpleKernel
     c::Vector{Tc}
     function LinearKernel(;c::T=0.0) where {T}
         new{T}([c])
@@ -28,7 +28,7 @@ The polynomial kernel is a Mercer kernel given by
 ```
 Where `c` is a real number, and `d` is a shape parameter bigger than 1. For `d = 1` see [`LinearKernel`](@ref)
 """
-struct PolynomialKernel{Td<:Real, Tc<:Real} <: BaseKernel
+struct PolynomialKernel{Td<:Real, Tc<:Real} <: SimpleKernel
     d::Vector{Td}
     c::Vector{Tc}
     function PolynomialKernel(; d::Td=2.0, c::Tc=0.0) where {Td<:Real, Tc<:Real}

src/basekernels/rationalquad.jl

@@ -7,7 +7,7 @@ The rational-quadratic kernel is a Mercer kernel given by the formula:
 ```
 where `α` is a shape parameter of the Euclidean distance. Check [`GammaRationalQuadraticKernel`](@ref) for a generalization.
 """
-struct RationalQuadraticKernel{Tα<:Real} <: BaseKernel
+struct RationalQuadraticKernel{Tα<:Real} <: SimpleKernel
     α::Vector{Tα}
     function RationalQuadraticKernel(;alpha::T=2.0, α::T=alpha) where {T}
         @check_args(RationalQuadraticKernel, α, α > zero(T), "α > 1")
@@ -28,7 +28,7 @@ The Gamma-rational-quadratic kernel is an isotropic Mercer kernel given by the f
 ```
 where `α` is a shape parameter of the Euclidean distance and `γ` is another shape parameter.
 """
-struct GammaRationalQuadraticKernel{Tα<:Real, Tγ<:Real} <: BaseKernel
+struct GammaRationalQuadraticKernel{Tα<:Real, Tγ<:Real} <: SimpleKernel
     α::Vector{Tα}
     γ::Vector{Tγ}
     function GammaRationalQuadraticKernel(;alpha::Tα=2.0, gamma::Tγ=2.0, α::Tα=alpha, γ::Tγ=gamma) where {Tα<:Real, Tγ<:Real}

src/generic.jl

@@ -22,7 +22,7 @@ end
 
 for k in concretetypes(Kernel, [])
     @eval begin
-        @inline (κ::$k)(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}) = kappa(κ, x, y)
+        @inline (κ::$k)(x, y) = kappa(κ, x, y)
         @inline (κ::$k)(X::AbstractMatrix{T}, Y::AbstractMatrix{T}; obsdim::Integer=defaultobs) where {T} = kernelmatrix(κ, X, Y, obsdim=obsdim)
         @inline (κ::$k)(X::AbstractMatrix{T}; obsdim::Integer=defaultobs) where {T} = kernelmatrix(κ, X, obsdim=obsdim)
     end

src/kernels/kernelproduct.jl

@@ -29,23 +29,26 @@ hadamard(x,y) = x.*y
 function kernelmatrix(
     κ::KernelProduct,
     X::AbstractMatrix;
-    obsdim::Int=defaultobs)
-    reduce(hadamard,kernelmatrix(κ.kernels[i],X,obsdim=obsdim) for i in 1:length(κ))
+    obsdim::Int=defaultobs,
+)
+    reduce(hadamard, kernelmatrix(κ.kernels[i], X, obsdim = obsdim) for i in 1:length(κ))
 end
 
 function kernelmatrix(
     κ::KernelProduct,
     X::AbstractMatrix,
     Y::AbstractMatrix;
-    obsdim::Int=defaultobs)
-    reduce(hadamard,_kernelmatrix(κ.kernels[i],X,Y,obsdim) for i in 1:length(κ))
+    obsdim::Int=defaultobs,
+)
+    reduce(hadamard, kernelmatrix(κ.kernels[i], X, Y, obsdim = obsdim) for i in 1:length(κ))
 end
 
 function kerneldiagmatrix(
     κ::KernelProduct,
     X::AbstractMatrix;
-    obsdim::Int=defaultobs) #TODO Add test
-    reduce(hadamard,kerneldiagmatrix(κ.kernels[i],X,obsdim=obsdim) for i in 1:length(κ))
+    obsdim::Int=defaultobs,
+) #TODO Add test
+    reduce(hadamard, kerneldiagmatrix(κ.kernels[i], X, obsdim = obsdim) for i in 1:length(κ))
 end
 
 function Base.show(io::IO, κ::KernelProduct)

src/kernels/kernelsum.jl

@@ -58,7 +58,7 @@ function kernelmatrix(
     Y::AbstractMatrix;
     obsdim::Int = defaultobs,
 )
-    sum(κ.weights[i] * _kernelmatrix(κ.kernels[i], X, Y, obsdim) for i in 1:length(κ))
+    sum(κ.weights[i] * kernelmatrix(κ.kernels[i], X, Y, obsdim = obsdim) for i in 1:length(κ))
 end
 
 function kerneldiagmatrix(

src/kernels/scaledkernel.jl

@@ -15,6 +15,8 @@ end
 
 kappa(k::ScaledKernel, x) = first(k.σ²) * kappa(k.kernel, x)
 
+kappa(k::ScaledKernel, x, y) = first(k.σ²) * kappa(k.kernel, x, y)
+
 metric(k::ScaledKernel) = metric(k.kernel)
 
 Base.:*(w::Real, k::Kernel) = ScaledKernel(k, w)

src/kernels/tensorproduct.jl

@@ -22,7 +22,6 @@ end
 
 Base.length(kernel::TensorProduct) = length(kernel.kernels)
 
-(kernel::TensorProduct)(x, y) = kappa(kernel, x, y)
 function kappa(kernel::TensorProduct, x, y)
     return prod(kappa(k, xi, yi) for (k, xi, yi) in zip(kernel.kernels, x, y))
 end
@@ -97,7 +96,7 @@ function kernelmatrix(
     obsdim ∈ (1, 2) || error("obsdim should be 1 or 2 (see docs of kernelmatrix))")
 
     featuredim = feature_dim(obsdim)
-    if !check_dims(X, X, featuredim, obsdim)
+    if !check_dims(X, X, featuredim)
         throw(DimensionMismatch("Dimensions of the target array K $(size(K)) are not " *
                                 "consistent with X $(size(X))"))
     end
@@ -120,7 +119,7 @@ function kernelmatrix(
     obsdim ∈ (1, 2) || error("obsdim should be 1 or 2 (see docs of kernelmatrix))")
 
     featuredim = feature_dim(obsdim)
-    if !check_dims(X, Y, featuredim, obsdim)
+    if !check_dims(X, Y, featuredim)
         throw(DimensionMismatch("Dimensions $(size(K)) of the target array K are not " *
                                 "consistent with X ($(size(X))) and Y ($(size(Y)))"))
     end

src/kernels/transformedkernel.jl

@@ -39,3 +39,17 @@ function printshifted(io::IO, κ::TransformedKernel, shift::Int)
     printshifted(io, κ.kernel, shift)
     print(io,"\n" * ("\t" ^ (shift + 1)) * "- $(κ.transform)")
 end
+
+# Kernel matrix operations
+
+kernelmatrix!(K::AbstractMatrix, κ::TransformedKernel, X::AbstractMatrix; obsdim::Int = defaultobs) =
+    kernelmatrix!(K, kernel(κ), apply(κ.transform, X, obsdim = obsdim), obsdim = obsdim)
+
+kernelmatrix!(K::AbstractMatrix, κ::TransformedKernel, X::AbstractMatrix, Y::AbstractMatrix; obsdim::Int = defaultobs) =
+    kernelmatrix!(K, kernel(κ), apply(κ.transform, X, obsdim = obsdim), apply(κ.transform, Y, obsdim = obsdim), obsdim = obsdim)
+
+kernelmatrix(κ::TransformedKernel, X::AbstractMatrix; obsdim::Int = defaultobs) =
+    kernelmatrix(kernel(κ), apply(κ.transform, X, obsdim = obsdim), obsdim = obsdim)
+    
+kernelmatrix(κ::TransformedKernel, X::AbstractMatrix, Y::AbstractMatrix; obsdim::Int = defaultobs) =
+    kernelmatrix(kernel(κ), apply(κ.transform, X, obsdim = obsdim), apply(κ.transform, Y, obsdim = obsdim), obsdim = obsdim)