Add Spectral Mixture Kernel

src/KernelFunctions.jl

@@ -28,6 +28,9 @@ export Transform, SelectTransform, ChainTransform, ScaleTransform, LinearTransfo
 
 export NystromFact, nystrom
 
+export spectral_mixture_kernel, spectral_mixture_product_kernel
+
+
 using Compat
 using Requires
 using Distances, LinearAlgebra

src/basekernels/sm.jl

@@ -0,0 +1,108 @@
+"""
+    spectral_mixture_kernel(
+        h::Kernel=SqExponentialKernel(),
+        αs::AbstractVector{<:Real},
+        γs::AbstractMatrix{<:Real},
+        ωs::AbstractMatrix{<:Real},
+    )
+
+where αs are the weights of dimension (A, ), γs is the covariance matrix of
+dimension (D, A) and ωs are the mean vectors and is of dimension (D, A).
+Here, D is input dimension and A is the number of spectral components.
+
+`h` is the kernel, which defaults to [`SqExponentialKernel`](@ref) if not specified.
+
+Generalised Spectral Mixture kernel function. This family of functions is  dense
+in the family of stationary real-valued kernels with respect to the pointwise convergence.[1]
+
+```math
+   κ(x, y) = αs' (h(-(γs' * t)^2) .* cos(π * ωs' * t), t = x - y
+```
+
+# References:
+    [1] Generalized Spectral Kernels, by Yves-Laurent Kom Samo and Stephen J. Roberts
+    [2] SM: Gaussian Process Kernels for Pattern Discovery and Extrapolation,
+            ICML, 2013, by Andrew Gordon Wilson and Ryan Prescott Adams,
+    [3] Covariance kernels for fast automatic pattern discovery and extrapolation
+        with Gaussian processes, Andrew Gordon Wilson, PhD Thesis, January 2014.
+        http://www.cs.cmu.edu/~andrewgw/andrewgwthesis.pdf
+    [4] http://www.cs.cmu.edu/~andrewgw/pattern/.
+
+"""
+function spectral_mixture_kernel(
+    h::Kernel,
+    αs::AbstractVector{<:Real},
+    γs::AbstractMatrix{<:Real},
+    ωs::AbstractMatrix{<:Real},
+)
+    if !(size(αs, 1) == size(γs, 2) == size(ωs, 2))
+        throw(DimensionMismatch("The dimensions of αs, γs, ans ωs do not match"))
+    end
+    if size(γs) != size(ωs)
+        throw(DimensionMismatch("The dimensions of γs ans ωs do not match"))
+    end
+
+    return sum(zip(αs, eachcol(γs), eachcol(ωs))) do (α, γ, ω)
+        a = TransformedKernel(h, LinearTransform(γ'))
+        b = TransformedKernel(CosineKernel(), LinearTransform(ω'))
+        return α * a * b
+    end
+end
+
+function spectral_mixture_kernel(
+    αs::AbstractVector{<:Real},
+    γs::AbstractMatrix{<:Real},
+    ωs::AbstractMatrix{<:Real}
+)
+    spectral_mixture_kernel(SqExponentialKernel(), αs, γs, ωs)
+end
+
+"""
+    spectral_mixture_product_kernel(
+        h::Kernel=SqExponentialKernel(),
+        αs::AbstractMatrix{<:Real},
+        γs::AbstractMatrix{<:Real},
+        ωs::AbstractMatrix{<:Real},
+    )
+
+where αs are the weights of dimension (D, A), γs is the covariance matrix of
+dimension (D, A) and ωs are the mean vectors and is of dimension (D, A).
+Here, D is input dimension and A is the number of spectral components.
+
+Spectral Mixture Product Kernel. With enough components A, the SMP kernel
+can model any product kernel to arbitrary precision, and is flexible even
+with a small number of components [1]
+
+
+`h` is the kernel, which defaults to [`SqExponentialKernel`](@ref) if not specified.
+
+```math
+   κ(x, y) = Πᵢ₌₁ᴷ Σ(αsᵢᵀ .* (h(-(γsᵢᵀ * tᵢ)²) .* cos(ωsᵢᵀ * tᵢ))), tᵢ = xᵢ - yᵢ
+```
+
+# References:
+    [1] GPatt: Fast Multidimensional Pattern Extrapolation with GPs,
+        arXiv 1310.5288, 2013, by Andrew Gordon Wilson, Elad Gilboa,
+        Arye Nehorai and John P. Cunningham
+"""
+function spectral_mixture_product_kernel(
+    h::Kernel,
+    αs::AbstractMatrix{<:Real},
+    γs::AbstractMatrix{<:Real},
+    ωs::AbstractMatrix{<:Real},
+)
+    if !(size(αs) == size(γs) == size(ωs))
+        throw(DimensionMismatch("The dimensions of αs, γs, ans ωs do not match"))
+    end
+    return TensorProduct(spectral_mixture_kernel(h, α, reshape(γ, 1, :), reshape(ω, 1, :))
+               for (α, γ, ω) in zip(eachrow(αs), eachrow(γs), eachrow(ωs)))
+end
+
+function spectral_mixture_product_kernel(
+    αs::AbstractMatrix{<:Real},
+    γs::AbstractMatrix{<:Real},
+    ωs::AbstractMatrix{<:Real}
+)
+    spectral_mixture_product_kernel(SqExponentialKernel(), αs, γs, ωs)
+end
+

test/basekernels/sm.jl

@@ -0,0 +1,24 @@
+@testset "sm" begin
+    v1 = rand(5)
+    v2 = rand(5)
+
+    αs₁ = rand(3)
+    αs₂ = rand(5, 3)
+    γs = rand(5, 3)
+    ωs = rand(5, 3)
+
+    k1 = spectral_mixture_kernel(αs₁, γs, ωs)
+    k2 = spectral_mixture_product_kernel(αs₂, γs, ωs)
+
+    t = v1 - v2
+
+    @test k1(v1, v2) ≈ sum(αs₁ .* exp.(-(t' * γs)'.^2) .*
+                           cospi.((t' * ωs)')) atol=1e-5
+
+    @test k2(v1, v2) ≈ prod(sum(αs₂[i,:]' .* exp.(-(γs[i,:]' * t[i]).^2) .*
+                                cospi.(ωs[i,:]' * t[i])) for i in 1:length(t)) atol=1e-5
+
+    @test_throws DimensionMismatch spectral_mixture_kernel(rand(5) ,rand(4,3), rand(4,3))
+    @test_throws DimensionMismatch spectral_mixture_kernel(rand(3) ,rand(4,3), rand(5,3))
+    @test_throws DimensionMismatch spectral_mixture_product_kernel(rand(5,3) ,rand(4,3), rand(5,3))
+end

test/runtests.jl

@@ -76,6 +76,7 @@ using KernelFunctions: metric, kappa
         include(joinpath("basekernels", "polynomial.jl"))
         include(joinpath("basekernels", "piecewisepolynomial.jl"))
         include(joinpath("basekernels", "rationalquad.jl"))
+        include(joinpath("basekernels", "sm.jl"))
         include(joinpath("basekernels", "wiener.jl"))
     end