Update of docs/README

Author: theogf

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
@@ -45,6 +45,7 @@ IdentityTransform
 ScaleTransform
 LowRankTransform
 FunctionTransform
+ChainTransform
 ```
 
 ## Functions
@@ -54,6 +55,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/theory.md

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

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,6 @@
 module KernelFunctions
 
-export kernel, kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!, kappa
+export kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!, kappa, kernelpdmat
 export get_params, set_params!
 
 

src/kernels/exponential.jl

@@ -4,7 +4,7 @@
 The squared exponential kernel is an isotropic Mercer kernel given by the formula:
 
 ```
-    κ(x,y) = exp(-‖x-y‖²)
+    κ(x,y) = exp(-ρ²‖x-y‖²)
 ```
 
 See also [`ExponentialKernel`](@ref) for a
@@ -19,8 +19,9 @@ struct SqExponentialKernel{T,Tr} <: Kernel{T,Tr}
 end
 
 @inline kappa(κ::SqExponentialKernel, d²::Real) = exp(-d²)
+@inline iskroncompatible(::SqExponentialKernel) = true
 
-### Aliases
+## Aliases ##
 const RBFKernel = SqExponentialKernel
 const GaussianKernel = SqExponentialKernel
 
@@ -30,7 +31,7 @@ const GaussianKernel = SqExponentialKernel
 The exponential kernel is an isotropic Mercer kernel given by the formula:
 
 ```
-    κ(x,y) = exp(-‖x-y‖)
+    κ(x,y) = exp(-ρ‖x-y‖)
 ```
 
 """
@@ -43,8 +44,9 @@ struct ExponentialKernel{T,Tr} <: Kernel{T,Tr}
 end
 
 @inline kappa(κ::ExponentialKernel, d::Real) = exp(-d)
+@inline iskroncompatible(::ExponentialKernel) = true
 
-### Aliases
+## Alias ##
 const LaplacianKernel = ExponentialKernel
 
 """
@@ -53,7 +55,7 @@ const LaplacianKernel = ExponentialKernel
 The γ-exponential kernel is an isotropic Mercer kernel given by the formula:
 
 ```
-    κ(x,y) = exp(-‖x-y‖^2γ)
+    κ(x,y) = exp(-ρ^(2γ)‖x-y‖^(2γ))
 ```
 """
 struct GammaExponentialKernel{T,Tr,Tᵧ<:Real} <: Kernel{T,Tr}
@@ -81,3 +83,4 @@ function GammaExponentialKernel(t::Tr,gamma::T₁=2.0) where {Tr<:Transform,T₁
 end
 
 @inline kappa(κ::GammaExponentialKernel, d²::Real) = exp(-d²^κ.γ)
+@inline iskroncompatible(::GammaExponentialKernel) = true

src/matrix/kernelkroeneckermat.jl

@@ -3,12 +3,25 @@ function kernelkronmat(
     X::AbstractVector,
     dims::Int
     )
-    @assert iskroncompatible(κ) "The kernel chosed is not compatible for kroenecker matrices"
-    K = kernelmatrix(κ,reshape(X,:,1),obsdim=1)
-
+    @assert iskroncompatible(κ) "The kernel chosed is not compatible for kroenecker matrices (see `iskroncompatible()`)"
+    k = kernelmatrix(κ,reshape(X,:,1),obsdim=1)
+    K = kron()
 end
 
+function kernelkronmat(
+    κ::Kernel,
+    X::AbstractVector{<:AbstractVector};
+    obsdim::Int=defaultobs
+    )
+    @assert iskroncompatible(κ) "The kernel chosed is not compatible for kroenecker matrices"
+    Ks = kernelmatrix.(κ,X,obsdim=obsdim)
+    K = kron(Ks)
+end
 
-function iskroncompatible(κ::Kernel)
 
-end
+"""
+    To be compatible with kroenecker constructions the kernel must satisfy
+    the property : for x,x' ∈ ℜᴰ
+    k(x,x') = ∏ᵢᴰ k(xᵢ,x'ᵢ)
+"""
+@inline iskroncompatible(κ::Kernel) = false # Default return for kernels

src/matrix/kernelmatrix.jl

@@ -1,51 +1,47 @@
 """
 ```
-    kernelmatrix!(K::Matrix, κ::Kernel, X::Matrix; obsdim::Integer=2, symmetrize::Bool=true)
+    kernelmatrix!(K::Matrix, κ::Kernel, X::Matrix; obsdim::Integer=2)
+    kernelmatrix!(K::Matrix, κ::Kernel, X::Matrix, Y::Matrix; obsdim::Integer=2)
 ```
 In-place version of `kernelmatrix` where pre-allocated matrix `K` will be overwritten with the kernel matrix.
 """
+kernelmatrix!
+
+
 function kernelmatrix!(
         K::Matrix,
         κ::Kernel,
         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
         map!(x->kappa(κ,x),K,pairwise(metric(κ),transform(κ,X,obsdim),dims=obsdim))
 end
 
-"""
-```
-    kernelmatrix!(K::Matrix, κ::Kernel, X::Matrix, Y::Matrix; obsdim::Integer=2)
-```
-In-place version of `kernelmatrix` where pre-allocated matrix `K` will be overwritten with the kernel matrix.
-"""
 function kernelmatrix!(
         K::AbstractMatrix,
         κ::Kernel,
         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,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
         map!(x->kappa(κ,x),K,pairwise(metric(κ),transform(κ,X,obsdim),transform(κ,Y,obsdim),dims=obsdim))
 end
 
-"""
-```
-    kernel(κ::Kernel, x, y; obsdim=2)
-```
-Apply the kernel `κ` to `x` and `y`.
-"""
-function kernel(κ::Kernel, x::Real, y::Real)
+## Apply kernel on two reals ##
+function _kernel(κ::Kernel, x::Real, y::Real)
     kernel(κ, [x], [y])
 end
 
-function kernel(
+## Apply kernel on two vectors ##
+function _kernel(
         κ::Kernel,
         x::AbstractVector,
         y::AbstractVector;
@@ -57,35 +53,31 @@ end
 
 """
 ```
-    kernelmatrix(κ::Kernel, X::Matrix ; obsdim::Int=2, symmetrize::Bool=true)
+    kernelmatrix(κ::Kernel, X::Matrix ; obsdim::Int=2)
+    kernelmatrix(κ::Kernel, X::Matrix, Y::Matrix; obsdim::Int=2)
 ```
-Calculate the kernel matrix of `X` with respect to kernel `κ`.
-`obsdim=1` means the matrix `X` has size #samples x #dimension
-`obsdim=2` means the matrix `X` has size #dimension x #samples
+Calculate the kernel matrix of `X` (and `Y`) with respect to kernel `κ`.
+`obsdim=1` means the matrix `X` (and `Y`) has size #samples x #dimension
+`obsdim=2` means the matrix `X` (and `Y`) has size #dimension x #samples
 """
+kernelmatrix
+
+
 function kernelmatrix(
         κ::Kernel,
         X::AbstractMatrix;
-        obsdim::Int = defaultobs,
-        symmetrize::Bool = true
+        obsdim::Int = defaultobs
     )
     K = map(x->kappa(κ,x),pairwise(metric(κ),transform(κ,X,obsdim),dims=obsdim))
 end
 
-"""
-```
-    kernelmatrix(κ::Kernel, X::Matrix, Y::Matrix; obsdim::Int=2)
-```
-Calculate the base matrix of `X` and `Y` with respect to kernel `κ`.
-`obsdim=1` means the matrices `X` and `Y` have sizes #samples x #dimension
-`obsdim=2` means the matrices `X` and `Y` have size #dimension x #samples
-"""
 function kernelmatrix(
         κ::Kernel,
         X::AbstractMatrix,
         Y::AbstractMatrix;
         obsdim=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
@@ -107,6 +99,7 @@ function kerneldiagmatrix(
         X::AbstractMatrix;
         obsdim::Int = defaultobs
         )
+        @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))"
         if obsdim == 1
             [@views kernel(κ,X[i,:],X[i,:]) for i in 1:size(X,obsdim)]
         elseif obsdim == 2
@@ -126,6 +119,7 @@ function kerneldiagmatrix!(
         X::AbstractMatrix;
         obsdim::Int = defaultobs
         )
+        @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))"
         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/matrix/kernelpdmat.jl

@@ -1,17 +1,19 @@
 """
-    Guarantees to return a positive-definite matrix in the form of a `PDMat` matrix with the cholesky decomposition precomputed
+    Compute a positive-definite matrix in the form of a `PDMat` matrix see [PDMats.jl]() with the cholesky decomposition precomputed
+    The algorithm recursively tries to add recursively a diagonal nugget until positive definiteness is achieved or that the noise is too big
 """
 function kernelpdmat(
         κ::Kernel,
         X::AbstractMatrix;
         obsdim::Int = defaultobs
         )
     K = kernelmatrix(κ,X,obsdim=obsdim)
+    Kmax =maximum(K)
     α = eps(eltype(K))
-    while !isposdef(K+αI) && α < 0.01*maximum(K)
+    while !isposdef(K+αI) && α < 0.01*Kmax
         α *= 2.0
     end
-    if α >= 0.01*maximum(K)
+    if α >= 0.01*Kmax
         @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)