Add NormalizedKernel (#274)

Co-authored-by: Théo Galy-Fajou <[email protected]>
Co-authored-by: David Widmann <[email protected]>

Author: 3 people

Project.toml

@@ -1,12 +1,13 @@
 name = "KernelFunctions"
 uuid = "ec8451be-7e33-11e9-00cf-bbf324bd1392"
-version = "0.9.2"
+version = "0.9.3"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 CompositionsBase = "a33af91c-f02d-484b-be07-31d278c5ca2b"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Functors = "d9f16b24-f501-4c13-a1f2-28368ffc5196"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -23,6 +24,7 @@ ChainRulesCore = "0.9"
 Compat = "3.7"
 CompositionsBase = "0.1"
 Distances = "0.10"
+FillArrays = "0.10, 0.11"
 Functors = "0.1"
 Requires = "1.0.1"
 SpecialFunctions = "0.8, 0.9, 0.10, 1"

docs/src/kernels.md

@@ -123,6 +123,7 @@ ScaledKernel
 KernelSum
 KernelProduct
 KernelTensorProduct
+NormalizedKernel
 ```
 
 ## Multi-output Kernels

src/KernelFunctions.jl

@@ -17,7 +17,7 @@ export RationalQuadraticKernel, GammaRationalQuadraticKernel
 export GaborKernel, PiecewisePolynomialKernel
 export PeriodicKernel, NeuralNetworkKernel
 export KernelSum, KernelProduct, KernelTensorProduct
-export TransformedKernel, ScaledKernel
+export TransformedKernel, ScaledKernel, NormalizedKernel
 
 export Transform,
     SelectTransform,
@@ -53,6 +53,7 @@ using ZygoteRules: ZygoteRules
 using StatsFuns: logtwo, twoπ
 using StatsBase
 using TensorCore
+using FillArrays
 
 abstract type Kernel end
 abstract type SimpleKernel <: Kernel end
@@ -89,6 +90,7 @@ include(joinpath("basekernels", "wiener.jl"))
 
 include(joinpath("kernels", "transformedkernel.jl"))
 include(joinpath("kernels", "scaledkernel.jl"))
+include(joinpath("kernels", "normalizedkernel.jl"))
 include(joinpath("matrix", "kernelmatrix.jl"))
 include(joinpath("kernels", "kernelsum.jl"))
 include(joinpath("kernels", "kernelproduct.jl"))

src/kernels/normalizedkernel.jl

@@ -0,0 +1,81 @@
+"""
+    NormalizedKernel(k::Kernel)
+
+A normalized kernel derived from `k`.
+
+# Definition
+
+For inputs ``x, x'``, the normalized kernel ``\\widetilde{k}`` derived from
+kernel ``k`` is defined as
+```math
+\\widetilde{k}(x, x'; k) = \\frac{k(x, x')}{\\sqrt{k(x, x) k(x', x')}}.
+```
+"""
+struct NormalizedKernel{Tk<:Kernel} <: Kernel
+    kernel::Tk
+end
+
+@functor NormalizedKernel
+
+(κ::NormalizedKernel)(x, y) = κ.kernel(x, y) / sqrt(κ.kernel(x, x) * κ.kernel(y, y))
+
+function kernelmatrix(κ::NormalizedKernel, x::AbstractVector, y::AbstractVector)
+    return kernelmatrix(κ.kernel, x, y) ./
+           sqrt.(
+        kernelmatrix_diag(κ.kernel, x) .* permutedims(kernelmatrix_diag(κ.kernel, y))
+    )
+end
+
+function kernelmatrix(κ::NormalizedKernel, x::AbstractVector)
+    x_diag = kernelmatrix_diag(κ.kernel, x)
+    return kernelmatrix(κ.kernel, x) ./ sqrt.(x_diag .* permutedims(x_diag))
+end
+
+function kernelmatrix_diag(κ::NormalizedKernel, x::AbstractVector)
+    first_x = first(x)
+    return Ones{typeof(κ(first_x, first_x))}(length(x))
+end
+
+function kernelmatrix_diag(κ::NormalizedKernel, x::AbstractVector, y::AbstractVector)
+    return kernelmatrix_diag(κ.kernel, x, y) ./
+           sqrt.(kernelmatrix_diag(κ.kernel, x) .* kernelmatrix_diag(κ.kernel, y))
+end
+
+function kernelmatrix!(
+    K::AbstractMatrix, κ::NormalizedKernel, x::AbstractVector, y::AbstractVector
+)
+    kernelmatrix!(K, κ.kernel, x, y)
+    K ./=
+        sqrt.(kernelmatrix_diag(κ.kernel, x) .* permutedims(kernelmatrix_diag(κ.kernel, y)))
+    return K
+end
+
+function kernelmatrix!(K::AbstractMatrix, κ::NormalizedKernel, x::AbstractVector)
+    kernelmatrix!(K, κ.kernel, x)
+    x_diag = kernelmatrix_diag(κ.kernel, x)
+    K ./= sqrt.(x_diag .* permutedims(x_diag))
+    return K
+end
+
+function kernelmatrix_diag!(
+    K::AbstractVector, κ::NormalizedKernel, x::AbstractVector, y::AbstractVector
+)
+    kernelmatrix_diag!(K, κ.kernel, x, y)
+    K ./= sqrt.(kernelmatrix_diag(κ.kernel, x) .* kernelmatrix_diag(κ.kernel, y))
+    return K
+end
+
+function kernelmatrix_diag!(K::AbstractVector, κ::NormalizedKernel, x::AbstractVector)
+    first_x = first(x)
+    return fill!(K, κ(first_x, first_x))
+end
+
+Base.show(io::IO, κ::NormalizedKernel) = printshifted(io, κ, 0)
+
+function printshifted(io::IO, κ::NormalizedKernel, shift::Int)
+    println(io, "Normalized Kernel:")
+    for _ in 1:(shift + 1)
+        print(io, "\t")
+    end
+    return printshifted(io, κ.kernel, shift + 1)
+end

test/kernels/normalizedkernel.jl

@@ -0,0 +1,16 @@
+@testset "normalizedkernel" begin
+    rng = MersenneTwister(123456)
+    x = randn(rng)
+    y = randn(rng)
+
+    k = 4 * SqExponentialKernel()
+    kn = NormalizedKernel(k)
+    @test kn(x, y) == k(x, y) / sqrt(k(x, x) * k(y, y))
+    @test kn(x, x) ≈ one(x) atol = 1e-5
+
+    # Standardised tests.
+    TestUtils.test_interface(kn, Float64)
+    test_ADs(x -> NormalizedKernel(exp(x[1]) * SqExponentialKernel()), rand(1))
+
+    test_params(kn, k)
+end

test/runtests.jl

@@ -123,6 +123,7 @@ include("test_utils.jl")
         include(joinpath("kernels", "overloads.jl"))
         include(joinpath("kernels", "scaledkernel.jl"))
         include(joinpath("kernels", "transformedkernel.jl"))
+        include(joinpath("kernels", "normalizedkernel.jl"))
     end
     @info "Ran tests on Kernel"