Merge branch 'master-dev'

Author: theogf

Project.toml

@@ -10,10 +10,11 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 
 [compat]
-julia = "1.0"
+Distances = "0.8.2"
 PDMats = "0.9.9"
 SpecialFunctions = "0.7.2"
-Distances = "0.8.2"
+StatsFuns = "0.8"
+julia = "1.0"
 
 [extras]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

README.md

@@ -34,9 +34,13 @@ The aim is to make the API as model-agnostic as possible while still being user-
   <img src="docs/src/assets/heatmap_combination.png" width=400px>
 </p>
 
-## Objectives (by priority)
-- AD Compatibility (Zygote, ForwardDiff)
-- Toeplitz Matrices
+## Packages goals (by priority)
+- Ensure AD Compatibility (Zygote, ForwardDiff)
+- Toeplitz Matrices compatibility
 - BLAS backend
 
 Directly inspired by the [MLKernels](https://github.com/trthatcher/MLKernels.jl) package.
+
+## Issues/Contributing
+
+If you notice a problem or would like to contribute by adding more kernel functions or features please [submit an issue](https://github.com/theogf/KernelFunctions.jl/issues).

docs/make.jl

@@ -6,7 +6,8 @@ makedocs(
     format = Documenter.HTML(),
     modules = [KernelFunctions],
     pages = ["Home"=>"index.md",
-            "User Guide" => "userguide.md",
+             "User Guide" => "userguide.md",
+             "Examples"=>"example.md",
              "Kernel Functions"=>"kernels.md",
              "Transform"=>"transform.md",
              "Metrics"=>"metrics.md",

docs/src/api.md

@@ -18,7 +18,7 @@ KernelFunctions
 
 ```@docs
 SqExponentialKernel
-Exponential
+ExponentialKernel
 GammaExponentialKernel
 ExponentiatedKernel
 MaternKernel
@@ -28,6 +28,9 @@ LinearKernel
 PolynomialKernel
 RationalQuadraticKernel
 GammaRationalQuadraticKernel
+ZeroKernel
+ConstantKernel
+WhiteKernel
 ```
 
 ## Kernel Combinations
@@ -45,6 +48,7 @@ IdentityTransform
 ScaleTransform
 LowRankTransform
 FunctionTransform
+ChainTransform
 ```
 
 ## Functions
@@ -54,6 +58,7 @@ kernelmatrix
 kernelmatrix!
 kerneldiagmatrix
 kerneldiagmatrix!
+kernelpdmat
 transform
 ```
 

docs/src/example.md

@@ -0,0 +1,21 @@
+# Examples (WIP)
+
+Here are a few examples of known complex kernels and how to do them. Or how to use kernels in a certain context
+
+## Kernel Ridge Regression
+
+Make a simple example of kernel ridge regression
+
+## Gaussian Process Regression
+
+Make a simple example of gaussian process regression
+
+## Deep Kernel Learning
+
+Put a Flux neural net in front of the kernel
+cf. Wilson paper
+
+## Kernel Selection
+
+Create a large collection of kernels and optimize the weights
+cf AISTATS 2018 paper

docs/src/index.md

@@ -2,4 +2,16 @@
 
 Model agnostic kernel functions compatible with automatic differentiation
 
-*** In Construction ***
+**KernelFunctions.jl** is a general purpose kernel package.
+It aims at providing a flexible framework for creating kernels and manipulating them.
+The main goals of this package compared to its predecessors/concurrents in [MLKernels.jl](https://github.com/trthatcher/MLKernels.jl), [Stheno.jl](https://github.com/willtebbutt/Stheno.jl), [GaussianProcesses.jl](https://github.com/STOR-i/GaussianProcesses.jl) and [AugmentedGaussianProcesses.jl](https://github.com/theogf/AugmentedGaussianProcesses.jl) are:
+- **Automatic Differentation** compatibility: all kernel functions should be differentiable via packages like [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) or [Zygote.jl](https://github.com/FluxML/Zygote.jl)
+- **Flexibility**: operations between kernels should be fluid and easy without breaking.
+- **Plug-and-play**: including the kernels before/after other steps should be straightforward.
+
+The methodology of how kernels are computed is quite simple and is done in three phases :
+- A `Transform` object is applied sample-wise on every sample
+- The pairwise matrix is computed using [Distances.jl](https://github.com/JuliaStats/Distances.jl) by using a `Metric` proper to each kernel
+- The `Kernel` function is applied element-wise on the pairwise matrix
+
+For a quick introduction on how to use it go to [User guide](@ref)

docs/src/kernels.md

@@ -1,34 +1,93 @@
 ```@meta
-CurrentModule = KernelFunctions
+  CurrentModule = KernelFunctions
 ```
 
 ## Exponential Kernels
 
-```@docs
-ExponentialKernel
-SqExponentialKernel
-GammaExponentialKernel
+### Exponential Kernel
+
+The [Exponential Kernel](@ref ExponentialKernel) is defined as
+```math
+  k(x,x') = \exp\left(-|x-x'|\right)
+```
+
+### Square Exponential Kernel
+
+The [Square Exponential Kernel](@ref KernelFunctions.SqExponentialKernel) is defined as 
+```math
+  k(x,x') = \exp\left(-\|x-x'\|^2\right)
+```
+
+### Gamma Exponential Kernel
+
+```math
+  k(x,x';\gamma) = \exp\left(-\|x-x'\|^{2\gamma}\right)
 ```
 
 ## Matern Kernels
 
-```@docs
-MaternKernel
-Matern32Kernel
-Matern52Kernel
+### Matern Kernel
+
+```math
+  k(x,x';\nu) = \frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}|x-x'|\right)K_\nu\left(\sqrt{2\nu}|x-x'|\right)
+```
+
+### Matern 3/2 Kernel
+
+```math
+  k(x,x') = \left(1+\sqrt{3}|x-x'|\right)\exp\left(\sqrt{3}|x-x'|\right)
+```
+
+### Matern 5/2 Kernel
+
+```math
+  k(x,x') = \left(1+\sqrt{5}|x-x'|+\frac{5}{2}\|x-x'\|^2\right)\exp\left(\sqrt{5}|x-x'|\right)
+```
+
+## Rational Quadratic
+
+### Rational Quadratic Kernel
+
+```math
+  k(x,x';\alpha) = \left(1+\frac{\|x-x'\|^2}{\alpha}\right)^{-\alpha}
+```
+
+### Gamma Rational Quadratic Kernel
+
+```math
+  k(x,x';\alpha,\gamma) = \left(1+\frac{\|x-x'\|^{2\gamma}}{\alpha}\right)^{-\alpha}
 ```
 
 ## Polynomial Kernels
 
-```@docs
-LinearKernel
-PolynomialKernel
+### LinearKernel
+
+```math
+  k(x,x';c) = \langle x,x'\rangle + c
+```
+
+### PolynomialKernel
+
+```math
+  k(x,x';c,d) = \left(\langle x,x'\rangle + c\right)^d
 ```
 
 ## Constant Kernels
 
-```@docs
-ConstantKernel
-WhiteKernel
-ZeroKernel
+### ConstantKernel
+
+```math
+  k(x,x';c) = c
+```
+
+### WhiteKernel
+
+```math
+  k(x,x') = \delta(x-x')
+```
+
+### ZeroKernel
+
+```math
+  k(x,x') = 0
 ```

docs/src/theory.md

@@ -0,0 +1 @@
+See [Uncyclopedia article](https://en.wikipedia.org/wiki/Positive-definite_kernel)

docs/src/transform.md

@@ -7,6 +7,9 @@ You can also create a pipeline of `Transform` via `TransformChain`. For example
 One apply a transformation on a matrix or a vector via `transform(t::Transform,v::AbstractVecOrMat)`
 
 ## Transforms :
+```@meta
+CurrentModule = KernelFunctions
+```
 
 ```@docs
   IdentityTransform

docs/src/userguide.md

@@ -1 +1,28 @@
-# Building kernel and matrices easily!
+# User guide
+
+## Kernel creation
+
+To create a kernel chose one of the kernels proposed, see [Kernels](@ref), or create your own, see [Creating Kernels](@ref)
+For example to create a square exponential kernel
+```julia
+  k = SqExponentialKernel()
+```
+All kernels can take as argument a `Transform` object (see [Transform](@ref)) which is directly going to act on the inputs before it's processes.
+But it's also possible to simply give a scalar or a vector if all you are interested in is to modify the lengthscale, respectively for all dimensions or independently for each dimension.
+
+## Kernel matrix creation
+
+Matrix are created via the `kernelmatrix` function or `kerneldiagmatrix`.
+An important argument to give is the dimensionality of the input `obsdim`. It tells if the matrix is of the type `# samples X # features` (`obsdim`=1) or `# features X # samples`(`obsdim`=2) (similarly to [Distances.jl](https://github.com/JuliaStats/Distances.jl))
+For example:
+```julia
+  k = SqExponentialKernel()
+  A = rand(10,5)
+  kernelmatrix(k,A,obsdim=1) # Return a 10x10 matrix
+  kernelmatrix(k,A,obsdim=2) # Return a 5x5 matrix
+```
+
+## Kernel manipulation
+
+One can create combinations of kernels via `KernelSum` and `KernelProduct` or using simple operators `+` and `*`.
+For

src/KernelFunctions.jl

@@ -1,6 +1,9 @@
 module KernelFunctions
 
-export kernel, kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!, kappa
+export kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!, kappa, kernelpdmat
+export get_params, set_params!
+
+
 export Kernel
 export ConstantKernel, WhiteKernel, ZeroKernel
 export SqExponentialKernel, ExponentialKernel, GammaExponentialKernel

src/generic.jl

@@ -1,10 +1,9 @@
 @inline metric(κ::Kernel) = κ.metric
 
 ## Allows to iterate over kernels
-Base.length(::Kernel) = 1 #TODO Add test
-
-Base.iterate(k::Kernel) = (k,nothing) #TODO Add test
-Base.iterate(k::Kernel, ::Any) = nothing #TODO Add test
+Base.length(::Kernel) = 1
+Base.iterate(k::Kernel) = (k,nothing)
+Base.iterate(k::Kernel, ::Any) = nothing
 
 ### 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]
@@ -24,8 +23,19 @@ end
 ## Constructors for kernels without parameters
 for kernel in [:ExponentialKernel,:SqExponentialKernel,:Matern32Kernel,:Matern52Kernel,:ExponentiatedKernel]
     @eval begin
-        $kernel(ρ::T=1.0) where {T<:Real} =   $kernel{T,ScaleTransform{T}}(ScaleTransform(ρ))
+        $kernel(ρ::T=1.0) where {T<:Real} =   $kernel{T,ScaleTransform{Base.RefValue{T}}}(ScaleTransform(ρ))
         $kernel(ρ::A) where {A<:AbstractVector{<:Real}} = $kernel{eltype(A),ScaleTransform{A}}(ScaleTransform(ρ))
         $kernel(t::Tr) where {Tr<:Transform} = $kernel{eltype(t),Tr}(t)
     end
 end
+
+function set!(k::Kernel,x)
+    @error "Setting parameters to this kernel is either not possible or has not been implemented"
+end
+
+set_params!(k::Kernel{T,<:ScaleTransform{<:Base.RefValue{<:Tρ}}},ρ::AbstractVector{<:Tρ}) where {T,Tρ<:Real} = set!(k.transform,ρ[1])
+set_params!(k::Kernel{T,<:ScaleTransform{<:AbstractVector{<:Tρ}}},ρ::AbstractVector{<:Tρ}) where {T,Tρ<:Real} = set!(k.transform,ρ)
+set_params!(k::Kernel{T,<:LowRankTransform{<:AbstractMatrix{<:Tm}}},m::AbstractMatrix{<:Tm}) where {T,Tm<:Real} = set!(k.transform,m)
+
+get_params(k::Kernel{T,<:ScaleTransform{<:Base.RefValue{<:Tρ}}}) where {T,Tρ} = [k.transform.s[]]
+get_params(k::Kernel{T,<:ScaleTransform{<:AbstractVector{<:Tρ}}}) where {T,Tρ} = k.transform.s

src/kernels/constant.jl

@@ -1,8 +1,7 @@
 """
-    ZeroKernel([tr=IdentityTransform()])
-
-    Create a kernel that always return a zero kernel matrix
+ZeroKernel([tr=IdentityTransform()])
 
+Create a kernel always returning zero
 """
 struct ZeroKernel{T,Tr} <: Kernel{T,Tr}
     transform::Tr
@@ -19,12 +18,12 @@ end
 @inline kappa(κ::ZeroKernel,d::T) where {T<:Real} = zero(T)
 
 """
-    WhiteKernel([tr=IdentityTransform()])
+`WhiteKernel([tr=IdentityTransform()])`
 
 ```
     κ(x,y) = δ(x,y)
 ```
-    Kernel function working as an equivalent to add white noise.
+Kernel function working as an equivalent to add white noise.
 """
 struct WhiteKernel{T,Tr} <: Kernel{T,Tr}
     transform::Tr
@@ -41,12 +40,11 @@ end
 @inline kappa(κ::WhiteKernel,δₓₓ::Real) = δₓₓ
 
 """
-    ConstantKernel([tr=IdentityTransform(),[c=1.0]])
-
+`ConstantKernel([tr=IdentityTransform(),[c=1.0]])`
 ```
     κ(x,y) = c
 ```
-    Kernel function always returning a constant value `c`
+Kernel function always returning a constant value `c`
 """
 struct ConstantKernel{T,Tr,Tc<:Real} <: Kernel{T,Tr}
     transform::Tr