Added KernelSum and KernelProduct and removed some useless parametrization

Author: theogf

README.md

@@ -1,8 +1,12 @@
 [![Build Status](https://travis-ci.org/theogf/KernelFunctions.jl.svg?branch=master)](https://travis-ci.org/theogf/AugmentedGaussianProcesses.jl)
 [![Coverage Status](https://coveralls.io/repos/github/theogf/KernelFunctions.jl/badge.svg?branch=master)](https://coveralls.io/github/theogf/KernelFunctions.jl?branch=master)
 [![Documentation](https://img.shields.io/badge/docs-dev-blue.svg)](https://theogf.github.io/KernelFunctions.jl/dev/)
-# KernelFunctions.jl (WIP)
-Julia Package for kernel functions for machine learning
+# KernelFunctions.jl
+## Kernel functions for machine learning
+
+KernelFunctions.jl provide a flexible and complete framework for kernel functions, pretransforming the input data.
+
+The aim is to make the API as model-agnostic as possible while still being user-friendly.
 
 ## Objectives (by priority)
 - ARD Kernels

src/KernelFunctions.jl

@@ -8,6 +8,7 @@ export ExponentiatedKernel
 export MaternKernel, Matern32Kernel, Matern52Kernel
 export LinearKernel, PolynomialKernel
 export RationalQuadraticKernel, GammaRationalQuadraticKernel
+export KernelSum, KernelProduct
 
 
 
@@ -32,6 +33,8 @@ for k in kernels
     include(joinpath("kernels",k*".jl"))
 end
 include("kernelmatrix.jl")
+include("kernels/kernelsum.jl")
+include("kernels/kernelproduct.jl")
 
 include("generic.jl")
 

src/generic.jl

@@ -1,5 +1,11 @@
 @inline metric(κ::Kernel) = κ.metric
 
+## Allows to iterate over kernels
+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]
     @eval begin

src/kernelmatrix.jl

@@ -5,11 +5,11 @@
 In-place version of `kernelmatrix` where pre-allocated matrix `K` will be overwritten with the kernel matrix.
 """
 function kernelmatrix!(
-        K::Matrix{T₁},
-        κ::Kernel{T},
-        X::AbstractMatrix{T₂};
+        K::Matrix,
+        κ::Kernel,
+        X::AbstractMatrix;
         obsdim::Int = defaultobs
-        ) where {T,T₁<:Real,T₂<:Real}
+        )
         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
@@ -23,12 +23,12 @@ end
 In-place version of `kernelmatrix` where pre-allocated matrix `K` will be overwritten with the kernel matrix.
 """
 function kernelmatrix!(
-        K::AbstractMatrix{T₁},
-        κ::Kernel{T},
-        X::AbstractMatrix{T₂},
-        Y::AbstractMatrix{T₃};
+        K::AbstractMatrix,
+        κ::Kernel,
+        X::AbstractMatrix,
+        Y::AbstractMatrix;
         obsdim::Int = defaultobs
-        ) where {T,T₁,T₂,T₃}
+        )
         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
@@ -41,16 +41,16 @@ end
 ```
 Apply the kernel `κ` to `x` and `y`.
 """
-function kernel(κ::Kernel{T}, x::Real, y::Real) where {T}
-    kernel(κ, [T(x)], [T(y)])
+function kernel(κ::Kernel, x::Real, y::Real)
+    kernel(κ, [x], [y])
 end
 
 function kernel(
-        κ::Kernel{T},
-        x::AbstractArray{T₁},
-        y::AbstractArray{T₂};
+        κ::Kernel,
+        x::AbstractVector,
+        y::AbstractVector;
         obsdim::Int = defaultobs
-    ) where {T,T₁<:Real,T₂<:Real}
+    )
     @assert length(x) == length(y) "x and y don't have the same dimension!"
     kappa(κ, evaluate(metric(κ),transform(κ,x),transform(κ,y)))
 end
@@ -64,11 +64,11 @@ Calculate the kernel matrix of `X` with respect to kernel `κ`.
 `obsdim=2` means the matrix `X` has size #dimension x #samples
 """
 function kernelmatrix(
-        κ::Kernel{T,<:Transform},
+        κ::Kernel,
         X::AbstractMatrix;
         obsdim::Int = defaultobs,
         symmetrize::Bool = true
-    ) where {T}
+    )
     K = map(x->kappa(κ,x),pairwise(metric(κ),transform(κ,X,obsdim),dims=obsdim))
 end
 
@@ -81,18 +81,19 @@ Calculate the base matrix of `X` and `Y` with respect to kernel `κ`.
 `obsdim=2` means the matrices `X` and `Y` have size #dimension x #samples
 """
 function kernelmatrix(
-        κ::Kernel{T},
-        X::AbstractMatrix{T₁},
-        Y::AbstractMatrix{T₂};
+        κ::Kernel,
+        X::AbstractMatrix,
+        Y::AbstractMatrix;
         obsdim=defaultobs
-    ) where {T,T₁<:Real,T₂<:Real}
+    )
     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))"))
     end
-    K = map(x->kappa(κ,x),pairwise(metric(κ),transform(κ,X,obsdim),transform(κ,Y,obsdim),dims=obsdim))
-    return K
+    _kernelmatrix(κ,X,Y,obsdim)
 end
 
+@inline _kernelmatrix(κ,X,Y,obsdim) = map(x->kappa(κ,x),pairwise(metric(κ),transform(κ,X,obsdim),transform(κ,Y,obsdim),dims=obsdim))
+
 """
 ```
     kerneldiagmatrix(κ::Kernel, X::Matrix; obsdim::Int=2)
@@ -102,14 +103,14 @@ Calculate the diagonal matrix of `X` with respect to kernel `κ`
 `obsdim=2` means the matrix `X` has size #dimension x #samples
 """
 function kerneldiagmatrix(
-        κ::Kernel{T},
-        X::AbstractMatrix{T₁};
+        κ::Kernel,
+        X::AbstractMatrix;
         obsdim::Int = defaultobs
-        ) where {T,T₁}
+        )
         if obsdim == 1
             [@views kernel(κ,X[i,:],X[i,:]) for i in 1:size(X,obsdim)]
         elseif obsdim == 2
-            [@views kernel(κ,X[i,:],X[i,:]) for i in 1:size(X,obsdim)]
+            [@views kernel(κ,X[:,i],X[:,i]) for i in 1:size(X,obsdim)]
         end
 end
 

src/kernels/kernelproduct.jl

@@ -1,3 +1,41 @@
 struct KernelProduct{T,Tr} <: Kernel{T,Tr}
+    kernels::Vector{Kernel}
+end
+
+function KernelProduct(kernels::AbstractVector{<:Kernel})
+    KernelProduct{eltype(kernels),Transform}(kernels)
+end
+
+Base.:*(k1::Kernel,k2::Kernel) = KernelProduct([k1,k2])
+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)
+
+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(κ))
+end
+
+function kernelmatrix(
+    κ::KernelProduct,
+    X::AbstractMatrix,
+    Y::AbstractMatrix;
+    obsdim::Int=defaultobs)
+    reduce(hadamard,_kernelmatrix(κ.kernels[i],X,Y,obsdim) for i in 1:length(κ))
+end
 
+function kerneldiagmatrix(
+    κ::KernelProduct,
+    X::AbstractMatrix;
+    obsdim::Int=defaultobs)
+    reduce(hadamard,kerneldiagmatrix(κ.kernels[i],X,obsdim=obsdim) for i in 1:length(κ))
 end

src/kernels/kernelsum.jl

@@ -1,8 +1,46 @@
 struct KernelSum{T,Tr} <: Kernel{T,Tr}
     kernels::Vector{Kernel}
     weights::Vector{Real}
+    function KernelSum{T,Tr}(kernels::AbstractVector{<:Kernel},weights::AbstractVector{<:Real}) where {T,Tr}
+        new{T,Tr}(kernels,weights)
+    end
 end
 
-function Base.:+(k1::Kernel,k2::Kernel)
-    KernelSum([k1,k2],[1.0,1.0])
+
+function KernelSum(kernels::AbstractVector{<:Kernel}; weights::AbstractVector{<:Real}=ones(Float64,length(kernels)))
+    @assert length(kernels)==length(weights) "Weights and kernel vector should be of the same length"
+    @assert all(weights.>=0) "All weights should be positive"
+    KernelSum{eltype(kernels),Transform}(kernels,weights)
+end
+
+Base.:+(k1::Kernel,k2::Kernel) = KernelSum([k1,k2],weights=[1.0,1.0])
+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.length(k::KernelSum) = length(k.kernels)
+metric(k::KernelSum) = metric.(k.kernels)
+transform(k::KernelSum) = transform.(k.kernels)
+transform(k::KernelSum,x::AbstractVecOrMat) = transform.(k.kernels,[x])
+transform(k::KernelSum,x::AbstractVecOrMat,obsdim::Int) = transform.(k.kernels,[x],obsdim)
+
+function kernelmatrix(
+    κ::KernelSum,
+    X::AbstractMatrix;
+    obsdim::Int=defaultobs)
+    sum(κ.weights[i]*kernelmatrix(κ.kernels[i],X,obsdim=obsdim) for i in 1:length(κ))
+end
+
+function kernelmatrix(
+    κ::KernelSum,
+    X::AbstractMatrix,
+    Y::AbstractMatrix;
+    obsdim::Int=defaultobs)
+    sum(κ.weights[i]*_kernelmatrix(κ.kernels[i],X,Y,obsdim) for i in 1:length(κ))
+end
+
+function kerneldiagmatrix(
+    κ::KernelSum,
+    X::AbstractMatrix;
+    obsdim::Int=defaultobs)
+    sum(kerneldiagmatrix(κ.kernels[i],X,obsdim=obsdim) for i in 1:length(κ))
 end

test/test_kernelmatrix.jl

@@ -23,12 +23,34 @@ k = SqExponentialKernel()
     end
     @testset "Kernel matrix" begin
         for obsdim in [1,2]
-            @test kernelmatrix(k,A,B,obsdim=obsdim) == kappa.([k],pairwise(KernelFunctions.metric(k),A,B,dims=obsdim))
-            @test kernelmatrix(k,A,obsdim=obsdim) == kappa.([k],pairwise(KernelFunctions.metric(k),A,dims=obsdim))
+            @test kernelmatrix(k,A,B,obsdim=obsdim) == kappa.(k,pairwise(KernelFunctions.metric(k),A,B,dims=obsdim))
+            @test kernelmatrix(k,A,obsdim=obsdim) == kappa.(k,pairwise(KernelFunctions.metric(k),A,dims=obsdim))
+            @test kerneldiagmatrix(k,A,obsdim=obsdim) == diag(kernelmatrix(k,A,obsdim=obsdim))
             @test k(A,B,obsdim=obsdim) == kernelmatrix(k,A,B,obsdim=obsdim)
             @test k(A,obsdim=obsdim) == kernelmatrix(k,A,obsdim=obsdim)
             @test kernel(k,1.0,2.0) == kernel(k,[1.0],[2.0])
             @test_throws DimensionMismatch kernelmatrix(k,A,C,obsdim=obsdim)
         end
     end
+    @testset "KernelSum" begin
+        k1 = SqExponentialKernel()
+        k2 = LinearKernel()
+        k3 =
+        ks = k1 + k2
+        w1 = 0.4; w2 = 1.2;
+        ks2 = KernelSum([k1,k2],weights=[w1,w2])
+        @test all(kernelmatrix(ks,A) .== kernelmatrix(k1,A) + kernelmatrix(k2,A))
+        @test all(kernelmatrix(ks,A,B) .== kernelmatrix(k1,A,B) + kernelmatrix(k2,A,B))
+        @test all(kerneldiagmatrix(ks,A) .== kerneldiagmatrix(k1,A) + kerneldiagmatrix(k2,A))
+        @test all(kernelmatrix(ks2,A) .== w1*kernelmatrix(k1,A) + w2*kernelmatrix(k2,A))
+    end
+    @testset "KernelProduct" begin
+        k1 = SqExponentialKernel()
+        k2 = LinearKernel()
+        kp = k1 * k2
+        @test all(kernelmatrix(kp,A) .== kernelmatrix(k1,A) .* kernelmatrix(k2,A))
+        @test all(kernelmatrix(kp,A,B) .== kernelmatrix(k1,A,B) .* kernelmatrix(k2,A,B))
+        @test all(kerneldiagmatrix(kp,A) .== kerneldiagmatrix(k1,A) .* kerneldiagmatrix(k2,A))
+        @test all(kernelmatrix(kp,A) .== kernelmatrix(k1,A) .* kernelmatrix(k2,A))
+    end
 end