Add NormalizedKernel

docs/src/kernels.md

@@ -125,6 +125,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,
@@ -88,6 +88,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 fill(κ(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"