Merge branch 'master-dev'

Author: theogf

Project.toml

@@ -6,6 +6,7 @@ version = "0.2.0"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 IRTools = "7869d1d1-7146-5819-86e3-90919afe41df"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

README.md

@@ -8,12 +8,35 @@ KernelFunctions.jl provide a flexible and complete framework for kernel function
 
 The aim is to make the API as model-agnostic as possible while still being user-friendly.
 
+## Examples
+
+```julia
+  X = reshape(collect(range(-3.0,3.0,length=100)),:,1)
+  # Set simple scaling of the data
+  k₁ = SqExponentialKernel(1.0)
+  K₁ = kernelmatrix(k,X,obsdim=1)
+
+  # Set a function transformation on the data
+  k₂ = MaternKernel(FunctionTransform(x->sin.(x)))
+  K₂ = kernelmatrix(k,X,obsdim=1)
+
+  # Set a matrix premultiplication on the data
+  k₃ = PolynomialKernel(LowRankTransform(randn(4,1)),0.0,2.0)
+  K₃ = kernelmatrix(k,X,obsdim=1)
+
+  # Add and sum kernels
+  k₄ = 0.5*SqExponentialKernel()*LinearKernel(0.5) + 0.4*k₂
+  K₄ = kernelmatrix(k,X,obsdim=1)
+
+  heatmap([K₁,K₂,K₃,K₄],yflip=false,colorbar=false)
+```
+<p align=center>
+  <img src="docs/src/assets/heatmap_combination.png" width=400px>
+</p>
+
 ## Objectives (by priority)
-- ARD Kernels
-- AD Compatible (Zygote, ForwardDiff, ReverseDiff)
-- Kernel sum and product
+- AD Compatibility (Zygote, ForwardDiff)
 - Toeplitz Matrices
 - BLAS backend
 
-
-Directly inspired by the [MLKernels](https://github.com/trthatcher/MLKernels.jl) package
+Directly inspired by the [MLKernels](https://github.com/trthatcher/MLKernels.jl) package.

docs/.gitignore

@@ -1,2 +1,4 @@
 build/
 site/
+
+#Temp to avoid to many changes

docs/create_kernel_plots.jl

@@ -10,6 +10,33 @@ x₀ = 0.0; l=0.1
 n_grid = 101
 fill(x₀,n_grid,1)
 xrange = reshape(collect(range(-3,3,length=n_grid)),:,1)
+
+k = SqExponentialKernel(1.0)
+K1 = kernelmatrix(k,xrange,obsdim=1)
+p = heatmap(K1,yflip=true,colorbar=false,framestyle=:none,background_color=RGBA(0.0,0.0,0.0,0.0))
+savefig(joinpath(@__DIR__,"src","assets","heatmap_sqexp.png"))
+
+
+k = Matern32Kernel(FunctionTransform(x->(sin.(x)).^2))
+K2 = kernelmatrix(k,xrange,obsdim=1)
+p = heatmap(K2,yflip=true,colorbar=false,framestyle=:none,background_color=RGBA(0.0,0.0,0.0,0.0))
+savefig(joinpath(@__DIR__,"src","assets","heatmap_matern.png"))
+
+
+k = PolynomialKernel(LowRankTransform(randn(3,1)),2.0,0.0)
+K3 = kernelmatrix(k,xrange,obsdim=1)
+p = heatmap(K3,yflip=true,colorbar=false,framestyle=:none,background_color=RGBA(0.0,0.0,0.0,0.0))
+savefig(joinpath(@__DIR__,"src","assets","heatmap_poly.png"))
+
+k = 0.5*SqExponentialKernel()*LinearKernel(0.5) + 0.4*Matern32Kernel(FunctionTransform(x->sin.(x)))
+K4 = kernelmatrix(k,xrange,obsdim=1)
+p = heatmap(K4,yflip=true,colorbar=false,framestyle=:none,background_color=RGBA(0.0,0.0,0.0,0.0))
+savefig(joinpath(@__DIR__,"src","assets","heatmap_prodsum.png"))
+
+plot(heatmap.([K1,K2,K3,K4],yflip=true,colorbar=false)...,layout=(2,2))
+savefig(joinpath(@__DIR__,"src","assets","heatmap_combination.png"))
+
+
 for k in [SqExponentialKernel,ExponentialKernel]
     K = kernelmatrix(k(),xrange,obsdim=1)
     v = rand(MvNormal(K+1e-7I))

docs/src/assets/heatmap_combination.png

@@ -16,6 +16,7 @@ using Distances, LinearAlgebra
 using Zygote: @adjoint
 using SpecialFunctions: lgamma, besselk
 using StatsFuns: logtwo
+using PDMats
 
 const defaultobs = 2
 
@@ -32,7 +33,7 @@ kernels = ["exponential","matern","polynomial","constant","rationalquad","expone
 for k in kernels
     include(joinpath("kernels",k*".jl"))
 end
-include("kernelmatrix.jl")
+include("matrix/kernelmatrix.jl")
 include("kernels/kernelsum.jl")
 include("kernels/kernelproduct.jl")
 

docs/src/assets/heatmap_matern.png

@@ -1,15 +1,15 @@
 @inline metric(κ::Kernel) = κ.metric
 
 ## Allows to iterate over kernels
-Base.length(::Kernel) = 1
+Base.length(::Kernel) = 1 #TODO Add test
 
-Base.iterate(k::Kernel) = (k,nothing)
-Base.iterate(k::Kernel, ::Any) = nothing
+Base.iterate(k::Kernel) = (k,nothing) #TODO Add test
+Base.iterate(k::Kernel, ::Any) = nothing #TODO Add test
 
 ### Syntactic sugar for creating matrices and using kernel functions
 for k in [:ExponentialKernel,:SqExponentialKernel,:GammaExponentialKernel,:MaternKernel,:Matern32Kernel,:Matern52Kernel,:LinearKernel,:PolynomialKernel,:ExponentiatedKernel,:ZeroKernel,:WhiteKernel,:ConstantKernel,:RationalQuadraticKernel,:GammaRationalQuadraticKernel]
     @eval begin
-        @inline (κ::$k)(d::Real) = kappa(κ,d)
+        @inline (κ::$k)(d::Real) = kappa(κ,d) #TODO Add test
         @inline (κ::$k)(x::AbstractVector{<:Real},y::AbstractVector{<:Real}) = kappa(κ,evaluate(κ.metric,transform(κ,x),transform(κ,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)

docs/src/assets/heatmap_poly.png

@@ -1,5 +1,5 @@
 """
-    ZeroKernel()
+    ZeroKernel([tr=IdentityTransform()])
 
     Create a kernel that always return a zero kernel matrix
 
@@ -19,7 +19,7 @@ end
 @inline kappa(κ::ZeroKernel,d::T) where {T<:Real} = zero(T)
 
 """
-    WhiteKernel()
+    WhiteKernel([tr=IdentityTransform()])
 
 ```
     κ(x,y) = δ(x,y)
@@ -41,7 +41,7 @@ end
 @inline kappa(κ::WhiteKernel,δₓₓ::Real) = δₓₓ
 
 """
-    ConstantKernel([c=1.0])
+    ConstantKernel([tr=IdentityTransform(),[c=1.0]])
 
 ```
     κ(x,y) = c

docs/src/assets/heatmap_prodsum.png

@@ -55,19 +55,6 @@ The γ-exponential kernel is an isotropic Mercer kernel given by the formula:
 ```
     κ(x,y) = exp(-‖x-y‖^2γ)
 ```
-
-# Examples
-
-```jldoctest; setup = :(using KernelFunctions)
-julia> GammaExponentialKernel()
-GammaExponentialKernel{Float64,Float64,Float64}(1.0,2.0)
-
-julia> GammaExponentialKernel(2.0f0,3.0)
-GammaExponentialKernel{Float32,Float32,Float64}(2.0,3.0)
-
-julia> GammaExponentialKernel([2.0,3.0],2f0)
-GammaExponentialKernel{Float64,Array{Float64},Float32}([2.0,3.0],2.0)
-```
 """
 struct GammaExponentialKernel{T,Tr,Tᵧ<:Real} <: Kernel{T,Tr}
     transform::Tr

docs/src/assets/heatmap_sqexp.png

@@ -1,5 +1,5 @@
 """
-    ExponentiatedKernel([α=1])
+    ExponentiatedKernel([ρ=1])
 
     The exponentiated kernel is a Mercer kernel given by:
 

src/KernelFunctions.jl

@@ -1,3 +1,11 @@
+"""
+    KernelProduct(kernels::Array{Kernel})
+Create a multiplication of kernels.
+One can also use the operator `*`
+```
+    kernelmatrix(SqExponentialKernel()*LinearKernel(),X) == kernelmatrix(SqExponentialKernel(),X).*kernelmatrix(LinearKernel(),X)
+```
+"""
 struct KernelProduct{T,Tr} <: Kernel{T,Tr}
     kernels::Vector{Kernel}
 end
@@ -7,14 +15,15 @@ function KernelProduct(kernels::AbstractVector{<:Kernel})
 end
 
 Base.:*(k1::Kernel,k2::Kernel) = KernelProduct([k1,k2])
+Base.:*(k1::KernelProduct,k2::KernelProduct) = KernelProduct(vcat(k1.kernels,k2.kernels)) #TODO Add test
 Base.:*(k::Kernel,kp::KernelProduct) = KernelProduct(vcat(k,kp.kernels))
 Base.:*(kp::KernelProduct,k::Kernel) = KernelProduct(vcat(kp.kernels,k))
 
 Base.length(k::KernelProduct) = length(k.kernels)
-metric(k::KernelProduct) = getmetric.(k.kernels)
-transform(k::KernelProduct) = transform.(k.kernels)
-transform(k::KernelProduct,x::AbstractVecOrMat) = transform.(k.kernels,[x])
-transform(k::KernelProduct,x::AbstractVecOrMat,obsdim::Int) = transform.(k.kernels,[x],obsdim)
+metric(k::KernelProduct) = getmetric.(k.kernels) #TODO Add test
+transform(k::KernelProduct) = transform.(k.kernels) #TODO Add test
+transform(k::KernelProduct,x::AbstractVecOrMat) = transform.(k.kernels,[x]) #TODO Add test
+transform(k::KernelProduct,x::AbstractVecOrMat,obsdim::Int) = transform.(k.kernels,[x],obsdim) #TODO Add test
 
 hadamard(x,y) = x.*y
 
@@ -36,6 +45,6 @@ end
 function kerneldiagmatrix(
     κ::KernelProduct,
     X::AbstractMatrix;
-    obsdim::Int=defaultobs)
+    obsdim::Int=defaultobs) #TODO Add test
     reduce(hadamard,kerneldiagmatrix(κ.kernels[i],X,obsdim=obsdim) for i in 1:length(κ))
 end

src/generic.jl

@@ -1,3 +1,11 @@
+"""
+    KernelSum(kernels::Array{Kernel};weights::Array{Real}=ones(length(kernels)))
+Create a positive weighted sum of kernels.
+One can also use the operator `+`
+```
+    kernelmatrix(SqExponentialKernel()+LinearKernel(),X) == kernelmatrix(SqExponentialKernel(),X).+kernelmatrix(LinearKernel(),X)
+```
+"""
 struct KernelSum{T,Tr} <: Kernel{T,Tr}
     kernels::Vector{Kernel}
     weights::Vector{Real}
@@ -14,8 +22,12 @@ function KernelSum(kernels::AbstractVector{<:Kernel}; weights::AbstractVector{<:
 end
 
 Base.:+(k1::Kernel,k2::Kernel) = KernelSum([k1,k2],weights=[1.0,1.0])
+Base.:+(k1::KernelSum,k2::KernelSum) = KernelSum(vcat(k1.kernels,k2.kernels),weights=vcat(k1.weights,k2.weights))
 Base.:+(k::Kernel,ks::KernelSum) = KernelSum(vcat(k,ks.kernels),weights=vcat(1.0,ks.weights))
 Base.:+(ks::KernelSum,k::Kernel) = KernelSum(vcat(ks.kernels,k),weights=vcat(ks.weights,1.0))
+Base.:*(w::Real,k::Kernel) = KernelSum([k],weights=[w]) #TODO add tests
+Base.:*(w::Real,k::KernelSum) = KernelSum(k.kernels,weights=w*k.weights) #TODO add tests
+
 
 Base.length(k::KernelSum) = length(k.kernels)
 metric(k::KernelSum) = metric.(k.kernels)

src/kernels/constant.jl

@@ -38,7 +38,7 @@ end
 @inline kappa(κ::MaternKernel, d::Real) = iszero(d) ? one(d) : exp((1.0-κ.ν)*logtwo-lgamma(κ.ν) + κ.ν*log(sqrt(2κ.ν)*d)+log(besselk(κ.ν,sqrt(2κ.ν)*d)))
 
 """
-    Matern32Kernel(ρ=1.0)
+    Matern32Kernel([ρ=1.0])
 
 The matern 3/2 kernel is an isotropic Mercer kernel given by the formula:
 
@@ -59,7 +59,7 @@ end
 @inline kappa(κ::Matern32Kernel, d::T) where {T<:Real} = (1+sqrt(3)*d)*exp(-sqrt(3)*d)
 
 """
-    Matern52Kernel(ρ=1.0)
+    Matern52Kernel([ρ=1.0])
 
 The matern 5/2 kernel is an isotropic Mercer kernel given by the formula:
 

src/kernels/exponential.jl

@@ -87,7 +87,7 @@ function kernelmatrix(
         obsdim=defaultobs
     )
     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 obsdim : $(feature_dim(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
     _kernelmatrix(κ,X,Y,obsdim)
 end
@@ -114,12 +114,18 @@ function kerneldiagmatrix(
         end
 end
 
+"""
+```
+    kerneldiagmatrix!(K::AbstractVector,κ::Kernel, X::Matrix; obsdim::Int=2)
+```
+In place version of `kerneldiagmatrix`
+"""
 function kerneldiagmatrix!(
-        K::AbstractVector{T₁},
-        κ::Kernel{T},
-        X::AbstractMatrix{T₂};
+        K::AbstractVector,
+        κ::Kernel,
+        X::AbstractMatrix;
         obsdim::Int = defaultobs
-        ) where {T,T₁,T₂}
+        )
         if length(K) != size(X,obsdim)
             throw(DimensionMismatch("Dimensions of the target array K $(size(K)) are not consistent with X $(size(X))"))
         end

src/kernels/exponentiated.jl

@@ -0,0 +1,18 @@
+"""
+    Guarantees to return a positive-definite matrix in the form of a `PDMat` matrix with the cholesky decomposition precomputed
+"""
+function kernelpdmat(
+        κ::Kernel,
+        X::AbstractMatrix;
+        obsdim::Int = defaultobs
+        )
+    K = kernelmatrix(κ,X,obsdim=obsdim)
+    α = eps(eltype(K))
+    while !isposdef(K+αI) && α < 0.01*maximum(K)
+        α *= 2.0
+    end
+    if α >= 0.01*maximum(K)
+        @error "Adding noise on the diagonal was not sufficient to build a positive-definite matrix:\n - Check that your kernel parameters are not extreme\n - Check that your data is sufficiently sparse\n - Maybe use a different kernel"
+    end
+    return PDMat(K+αI)
+end

src/kernels/kernelproduct.jl

@@ -1,6 +1,9 @@
 """
     FunctionTransform
-
+```
+    f(x) = abs.(x)
+    tr = FunctionTransform(f)
+```
     Take a function `f` as an argument which is going to act on each vector individually.
     Make sure that `f` is supposed to act on a vector by eventually using broadcasting
     For example `f(x)=sin(x)` -> `f(x)=sin.(x)`

src/kernels/kernelsum.jl

@@ -1,17 +1,27 @@
+"""
+    LowRankTransform
+    ```
+        P = rand(10,5)
+        tr = LowRankTransform(P)
+    ```
+    Apply the low-rank projection realised by the matrix `P`
+    The second dimension of `P` must match the number of features of the target.
+"""
 struct LowRankTransform{T<:AbstractMatrix{<:Real}} <: Transform
     proj::T
 end
 
 Base.size(tr::LowRankTransform,i::Int) = size(tr.proj,i)
-Base.size(tr::LowRankTransform) = size(tr.proj)
+Base.size(tr::LowRankTransform) = size(tr.proj) #  TODO Add test
 
 function transform(t::LowRankTransform,X::AbstractMatrix{<:Real},obsdim::Int=defaultobs)
     @boundscheck size(t,2) != size(X,feature_dim(obsdim)) ?
         throw(DimensionMismatch("The projection matrix has size $(size(t)) and cannot be used on X with dimensions $(size(X))")) : nothing
     @inbounds _transform(t,X,obsdim)
 end
-function transform(t::LowRankTransform,x::AbstractVector{<:Real})
-    @assert size(t,2) == length(x) "Vector has wrong dimensions"
+
+function transform(t::LowRankTransform,x::AbstractVector{<:Real},obsdim::Int=defaultobs) #TODO Add test
+    @assert size(t,2) == length(x) "Vector has wrong dimensions $(length(x)) compared to projection matrix"
     t.proj*X
 end
 

src/kernels/matern.jl

@@ -1,5 +1,13 @@
 """
     Scale Transform
+```
+    l = 2.0
+    tr = ScaleTransform(l)
+    v = rand(3)
+    tr = ScaleTransform(v)
+```
+    Multiply every element of the matrix by `l` for a scalar
+    Multiply every vector of observation by `v` element-wise for a vector
 """
 struct ScaleTransform{T<:Union{Real,AbstractVector{<:Real}}} <: Transform
     s::T
@@ -10,7 +18,7 @@ function ScaleTransform(s::T=1.0) where {T<:Real}
     ScaleTransform{T}(s)
 end
 
-function ScaleTransform(s::T,dims::Integer) where {T<:Real}
+function ScaleTransform(s::T,dims::Integer) where {T<:Real} # TODO Add test
     @check_args(ScaleTransform, s, s > zero(T), "s > 0")
     ScaleTransform{Vector{T}}(fill(s,dims))
 end
@@ -20,12 +28,12 @@ function ScaleTransform(s::A) where {A<:AbstractVector{<:Real}}
     ScaleTransform{A}(s)
 end
 
-dim(str::ScaleTransform{<:Real}) = 1
+dim(str::ScaleTransform{<:Real}) = 1 #TODO Add test
 dim(str::ScaleTransform{<:AbstractVector{<:Real}}) = length(str.s)
 
 function transform(t::ScaleTransform{<:AbstractVector{<:Real}},X::AbstractMatrix{<:Real},obsdim::Int)
     @boundscheck if dim(t) != size(X,!Bool(obsdim-1)+1)
-        throw(DimensionMismatch("Array has size $(size(X,!Bool(obsdim-1)+1)) on dimension $(!Bool(obsdim-1)+1)) which does not match the length of the scale transform length , $(dim(t))."))
+        throw(DimensionMismatch("Array has size $(size(X,!Bool(obsdim-1)+1)) on dimension $(!Bool(obsdim-1)+1)) which does not match the length of the scale transform length , $(dim(t)).")) #TODO Add test
     end
     _transform(t,X,obsdim)
 end