Fix SqExponential and GammaExponential + style (#158)

* Fix style

* Fix convention

* Add Gamma Exponential kernel reference

* Update src/matrix/kernelpdmat.jl

Co-authored-by: Théo Galy-Fajou <theo.galyfajou@gmail.com>

* Warn about breaking change

* Only warn when first loaded

Co-authored-by: Will Tebbutt <will.tebbutt@invenialabs.co.uk>
Co-authored-by: Théo Galy-Fajou <theo.galyfajou@gmail.com>

Author: 3 people

.gitignore

@@ -2,3 +2,4 @@
 *.cov
 Manifest.toml
 coverage/
+src/update_v0.8.0

Project.toml

@@ -1,6 +1,6 @@
 name = "KernelFunctions"
 uuid = "ec8451be-7e33-11e9-00cf-bbf324bd1392"
-version = "0.7.2"
+version = "0.8.0"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"

src/KernelFunctions.jl

@@ -4,6 +4,19 @@ KernelFunctions. [Github](https://github.com/JuliaGaussianProcesses/KernelFuncti
 """
 module KernelFunctions
 
+if !isfile(joinpath(@__DIR__, "update_v0.8.0"))
+    printstyled(
+        stdout,
+        """
+        WARNING: SqExponentialKernel changed convention in version 0.8.0.
+        This kernel now divides the squared distance by 2 to align with standard practice.
+        This warning will be removed in 0.9.0.
+        """;
+        color = Base.info_color(),
+    )
+    touch(joinpath(@__DIR__, "update_v0.8.0"))
+end
+
 export kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!
 export transform
 export duplicate, set! # Helpers
@@ -53,35 +66,51 @@ abstract type Kernel end
 abstract type SimpleKernel <: Kernel end
 
 include("utils.jl")
-include("distances/pairwise.jl")
-include("distances/dotproduct.jl")
-include("distances/delta.jl")
-include("distances/sinus.jl")
-include("transform/transform.jl")
-
-for f in readdir(joinpath(@__DIR__, "basekernels"))
-    endswith(f, ".jl") && include(joinpath("basekernels", f))
-end
-
-include("kernels/transformedkernel.jl")
-include("kernels/scaledkernel.jl")
-include("matrix/kernelmatrix.jl")
-include("kernels/kernelsum.jl")
-include("kernels/kernelproduct.jl")
-include("kernels/tensorproduct.jl")
-include("approximations/nystrom.jl")
+include(joinpath("distances", "pairwise.jl"))
+include(joinpath("distances", "dotproduct.jl"))
+include(joinpath("distances", "delta.jl"))
+include(joinpath("distances", "sinus.jl"))
+include(joinpath("transform", "transform.jl"))
+
+include(joinpath("basekernels", "constant.jl"))
+include(joinpath("basekernels", "cosine.jl"))
+include(joinpath("basekernels", "exponential.jl"))
+include(joinpath("basekernels", "exponentiated.jl"))
+include(joinpath("basekernels", "fbm.jl"))
+include(joinpath("basekernels", "gabor.jl"))
+include(joinpath("basekernels", "maha.jl"))
+include(joinpath("basekernels", "matern.jl"))
+include(joinpath("basekernels", "nn.jl"))
+include(joinpath("basekernels", "periodic.jl"))
+include(joinpath("basekernels", "piecewisepolynomial.jl"))
+include(joinpath("basekernels", "polynomial.jl"))
+include(joinpath("basekernels", "rationalquad.jl"))
+include(joinpath("basekernels", "sm.jl"))
+include(joinpath("basekernels", "wiener.jl"))
+
+include(joinpath("kernels", "transformedkernel.jl"))
+include(joinpath("kernels", "scaledkernel.jl"))
+include(joinpath("matrix", "kernelmatrix.jl"))
+include(joinpath("kernels", "kernelsum.jl"))
+include(joinpath("kernels", "kernelproduct.jl"))
+include(joinpath("kernels", "tensorproduct.jl"))
+include(joinpath("approximations", "nystrom.jl"))
 include("generic.jl")
 
-include("mokernels/mokernel.jl")
-include("mokernels/moinput.jl")
-include("mokernels/independent.jl")
-include("mokernels/slfm.jl")
+include(joinpath("mokernels", "mokernel.jl"))
+include(joinpath("mokernels", "moinput.jl"))
+include(joinpath("mokernels", "independent.jl"))
+include(joinpath("mokernels", "slfm.jl"))
 
 include("zygote_adjoints.jl")
 
 function __init__()
-    @require Kronecker="2c470bb0-bcc8-11e8-3dad-c9649493f05e" include("matrix/kernelkroneckermat.jl")
-    @require PDMats="90014a1f-27ba-587c-ab20-58faa44d9150" include("matrix/kernelpdmat.jl")
+    @require Kronecker="2c470bb0-bcc8-11e8-3dad-c9649493f05e" begin
+        include(joinpath("matrix", "kernelkroneckermat.jl"))
+    end
+    @require PDMats="90014a1f-27ba-587c-ab20-58faa44d9150" begin
+        include(joinpath("matrix", "kernelpdmat.jl"))
+    end
 end
 
 end

src/basekernels/exponential.jl

@@ -11,7 +11,7 @@ related form of the kernel or [`GammaExponentialKernel`](@ref) for a generalizat
 """
 struct SqExponentialKernel <: SimpleKernel end
 
-kappa(κ::SqExponentialKernel, d²::Real) = exp(-d²)
+kappa(κ::SqExponentialKernel, d²::Real) = exp(-d² / 2)
 
 metric(::SqExponentialKernel) = SqEuclidean()
 
@@ -48,12 +48,15 @@ const LaplacianKernel = ExponentialKernel
 """
     GammaExponentialKernel(; γ = 2.0)
 
-The γ-exponential kernel is an isotropic Mercer kernel given by the formula:
+The γ-exponential kernel [1] is an isotropic Mercer kernel given by the formula:
 ```
-    κ(x,y) = exp(-‖x-y‖^(2γ))
+    κ(x,y) = exp(-‖x-y‖^γ)
 ```
 Where `γ > 0`, (the keyword `γ` can be replaced by `gamma`)
-For `γ = 1`, see `SqExponentialKernel` and `γ = 0.5`, see `ExponentialKernel`
+For `γ = 2`, see `SqExponentialKernel` and `γ = 1`, see `ExponentialKernel`.
+
+[1] - Gaussian Processes for Machine Learning, Carl Edward Rasmussen and Christopher K. I.
+    Williams, MIT Press, 2006.
 """
 struct GammaExponentialKernel{Tγ<:Real} <: SimpleKernel
     γ::Vector{Tγ}
@@ -65,10 +68,12 @@ end
 
 @functor GammaExponentialKernel
 
-kappa(κ::GammaExponentialKernel, d²::Real) = exp(-d²^first(κ.γ))
+kappa(κ::GammaExponentialKernel, d::Real) = exp(-d^first(κ.γ))
 
-metric(::GammaExponentialKernel) = SqEuclidean()
+metric(::GammaExponentialKernel) = Euclidean()
 
 iskroncompatible(::GammaExponentialKernel) = true
 
-Base.show(io::IO, κ::GammaExponentialKernel) = print(io, "Gamma Exponential Kernel (γ = ", first(κ.γ), ")")
+function Base.show(io::IO, κ::GammaExponentialKernel)
+    print(io, "Gamma Exponential Kernel (γ = ", first(κ.γ), ")")
+end

src/basekernels/gabor.jl

@@ -11,7 +11,7 @@ struct GaborKernel{K<:Kernel} <: Kernel
     kernel::K
     function GaborKernel(;ell=nothing, p=nothing)
         k = _gabor(ell=ell, p=p)
-        new{typeof(k)}(k)
+        return new{typeof(k)}(k)
     end
 end
 
@@ -57,24 +57,17 @@ end
 
 Base.show(io::IO, κ::GaborKernel) = print(io, "Gabor Kernel (ell = ", κ.ell, ", p = ", κ.p, ")")
 
-function kernelmatrix(
-    κ::GaborKernel,
-    X::AbstractMatrix;
-    obsdim::Int=defaultobs)
-    kernelmatrix(κ.kernel, X; obsdim=obsdim)
+function kernelmatrix(κ::GaborKernel, X::AbstractMatrix; obsdim::Int=defaultobs)
+    return kernelmatrix(κ.kernel, X; obsdim=obsdim)
 end
 
 function kernelmatrix(
-    κ::GaborKernel,
-    X::AbstractMatrix,
-    Y::AbstractMatrix;
-    obsdim::Int=defaultobs)
-    kernelmatrix(κ.kernel, X, Y; obsdim=obsdim)
+    κ::GaborKernel, X::AbstractMatrix, Y::AbstractMatrix;
+    obsdim::Int=defaultobs,
+)
+    return kernelmatrix(κ.kernel, X, Y; obsdim=obsdim)
 end
 
-function kerneldiagmatrix(
-    κ::GaborKernel,
-    X::AbstractMatrix;
-    obsdim::Int=defaultobs) #TODO Add test
-    kerneldiagmatrix(κ.kernel, X; obsdim=obsdim)
+function kerneldiagmatrix(κ::GaborKernel, X::AbstractMatrix; obsdim::Int=defaultobs) #TODO Add test
+    return kerneldiagmatrix(κ.kernel, X; obsdim=obsdim)
 end

src/matrix/kernelkroneckermat.jl

@@ -2,22 +2,16 @@ using .Kronecker
 
 export kernelkronmat
 
-function kernelkronmat(
-    κ::Kernel,
-    X::AbstractVector,
-    dims::Int
-    )
+function kernelkronmat(κ::Kernel, X::AbstractVector, dims::Int)
     @assert iskroncompatible(κ) "The chosen kernel is not compatible for kroenecker matrices (see [`iskroncompatible`](@ref))"
     k = kernelmatrix(κ, X)
     kronecker(k, dims)
 end
 
 function kernelkronmat(
-    κ::Kernel,
-    X::AbstractVector{<:AbstractVector};
-    obsdim::Int=defaultobs
-    )
-    @assert iskroncompatible(κ) "The chosen kernel is not compatible for kroenecker matrices"
+    κ::Kernel, X::AbstractVector{<:AbstractVector}; obsdim::Int=defaultobs,
+)
+    @assert iskroncompatible(κ) "The chosen kernel is not compatible for Kronecker matrices"
     Ks = kernelmatrix.(κ, X)
     K = reduce(⊗, Ks)
 end

src/matrix/kernelmatrix.jl

@@ -74,10 +74,7 @@ function kernelmatrix!(K::AbstractMatrix, κ::SimpleKernel, x::AbstractVector)
 end
 
 function kernelmatrix!(
-    K::AbstractMatrix,
-    κ::SimpleKernel,
-    x::AbstractVector,
-    y::AbstractVector,
+    K::AbstractMatrix, κ::SimpleKernel, x::AbstractVector, y::AbstractVector,
 )
     validate_inplace_dims(K, x, y)
     pairwise!(K, metric(κ), x, y)
@@ -102,19 +99,13 @@ end
 const defaultobs = 2
 
 function kernelmatrix!(
-    K::AbstractMatrix,
-    κ::Kernel,
-    X::AbstractMatrix;
-    obsdim::Int = defaultobs
+    K::AbstractMatrix, κ::Kernel, X::AbstractMatrix; obsdim::Int = defaultobs
 )
     return kernelmatrix!(K, κ, vec_of_vecs(X; obsdim=obsdim))
 end
 
 function kernelmatrix!(
-    K::AbstractMatrix,
-    κ::Kernel,
-    X::AbstractMatrix,
-    Y::AbstractMatrix;
+    K::AbstractMatrix, κ::Kernel, X::AbstractMatrix, Y::AbstractMatrix;
     obsdim::Int = defaultobs
 )
     x = vec_of_vecs(X; obsdim=obsdim)
@@ -133,10 +124,7 @@ function kernelmatrix(κ::Kernel, X::AbstractMatrix, Y::AbstractMatrix; obsdim=d
 end
 
 function kerneldiagmatrix!(
-    K::AbstractVector,
-    κ::Kernel,
-    X::AbstractMatrix;
-    obsdim::Int = defaultobs
+    K::AbstractVector, κ::Kernel, X::AbstractMatrix; obsdim::Int = defaultobs
 )
     return kerneldiagmatrix!(K, κ, vec_of_vecs(X; obsdim=obsdim))
 end

src/matrix/kernelpdmat.jl

@@ -3,24 +3,28 @@ using .PDMats: PDMat
 export kernelpdmat
 
 """
-    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
+    Compute a positive-definite matrix in the form of a `PDMat` matrix see [PDMats.jl](https://github.com/JuliaStats/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)
+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*Kmax
+    while !isposdef(K + α * I) && α < 0.01 * Kmax
         α *= 2.0
     end
-    if α >= 0.01*Kmax
-        throw(ErrorException("Adding noise on the diagonal was not sufficient to build a positive-definite matrix:\n\t- Check that your kernel parameters are not extreme\n\t- Check that your data is sufficiently sparse\n\t- Maybe use a different kernel"))
+    if α >= 0.01 * Kmax
+        error(
+            "Adding noise on the diagonal was not sufficient to build a positive-definite" *
+            " matrix:\n\t- Check that your kernel parameters are not extreme\n\t- Check" *
+            " that your data is sufficiently sparse\n\t- Maybe use a different kernel",
+        )
     end
-    return PDMat(K+α*I)
+    return PDMat(K + α * I)
 end
 
-kernelpdmat(κ::Kernel,X::AbstractVector{<:Real};obsdim=defaultobs) = kernelpdmat(κ,reshape(X,1,:),obsdim=2)
+function kernelpdmat(κ::Kernel, X::AbstractVector{<:Real}; obsdim=defaultobs)
+    return kernelpdmat(κ, reshape(X, 1, :); obsdim=2)
+end

test/basekernels/exponential.jl

@@ -5,8 +5,8 @@
     v2 = rand(rng, 3)
     @testset "SqExponentialKernel" begin
         k = SqExponentialKernel()
-        @test kappa(k,x) ≈ exp(-x)
-        @test k(v1,v2) ≈ exp(-norm(v1-v2)^2)
+        @test kappa(k,x) ≈ exp(-x / 2)
+        @test k(v1,v2) ≈ exp(-norm(v1-v2)^2 / 2)
         @test kappa(SqExponentialKernel(),x) == kappa(k,x)
         @test metric(SqExponentialKernel()) == SqEuclidean()
         @test RBFKernel == SqExponentialKernel
@@ -30,18 +30,20 @@
     @testset "GammaExponentialKernel" begin
         γ = 2.0
         k = GammaExponentialKernel(γ=γ)
-        @test kappa(k,x) ≈ exp(-(x)^(γ))
-        @test k(v1,v2) ≈ exp(-norm(v1-v2)^(2γ))
-        @test kappa(GammaExponentialKernel(),x) == kappa(k,x)
+        @test k(v1, v2) ≈ exp(-norm(v1 - v2)^γ)
+        @test kappa(GammaExponentialKernel(), x) == kappa(k, x)
         @test GammaExponentialKernel(gamma=γ).γ == [γ]
-        @test metric(GammaExponentialKernel()) == SqEuclidean()
-        @test metric(GammaExponentialKernel(γ=2.0)) == SqEuclidean()
+        @test metric(GammaExponentialKernel()) == Euclidean()
+        @test metric(GammaExponentialKernel(γ=2.0)) == Euclidean()
         @test repr(k) == "Gamma Exponential Kernel (γ = $(γ))"
         @test KernelFunctions.iskroncompatible(k) == true
         test_ADs(γ -> GammaExponentialKernel(gamma=first(γ)), [γ])
         test_params(k, ([γ],))
         #Coherence :
-        @test GammaExponentialKernel(γ=1.0)(v1,v2) ≈ SqExponentialKernel()(v1,v2)
-        @test GammaExponentialKernel(γ=0.5)(v1,v2) ≈ ExponentialKernel()(v1,v2)
+        @test isapprox(
+            GammaExponentialKernel(γ=2.0)(sqrt(0.5) * v1, sqrt(0.5) * v2),
+            SqExponentialKernel()(v1,v2),
+        )
+        @test GammaExponentialKernel(γ=1.0)(v1, v2) ≈ ExponentialKernel()(v1, v2)
     end
 end

test/basekernels/gabor.jl

@@ -6,7 +6,7 @@
     @test k.ell ≈ ell atol=1e-5
     @test k.p ≈ p atol=1e-5
 
-    k_manual = exp(-sqeuclidean(v1, v2) / (k.ell^2)) * cospi(euclidean(v1, v2) / k.p)
+    k_manual = exp(-sqeuclidean(v1, v2) / (2 * k.ell^2)) * cospi(euclidean(v1, v2) / k.p)
     @test k(v1,v2) ≈ k_manual atol=1e-5
 
     lhs_manual = transform(SqExponentialKernel(), 1/k.ell)(v1,v2)