Add general TensorProduct kernel

src/KernelFunctions.jl

@@ -19,6 +19,7 @@ export MahalanobisKernel, GaborKernel, PiecewisePolynomialKernel
 export PeriodicKernel
 export KernelSum, KernelProduct
 export TransformedKernel, ScaledKernel
+export TensorProduct
 
 export Transform, SelectTransform, ChainTransform, ScaleTransform, LowRankTransform, IdentityTransform, FunctionTransform
 
@@ -56,6 +57,7 @@ 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("generic.jl")
 

src/kernels/tensorproduct.jl

@@ -0,0 +1,175 @@
+"""
+    TensorProduct(kernels...)
+
+Create a tensor product of kernels.
+"""
+struct TensorProduct{K} <: Kernel
+    kernels::K
+end
+
+function TensorProduct(kernel::Kernel, kernels::Kernel...)
+    return TensorProduct((kernel, kernels...))
+end
+
+Base.length(kernel::TensorProduct) = length(kernel.kernels)
+
+(kernel::TensorProduct)(x, y) = kappa(kernel, x, y)
+function kappa(kernel::TensorProduct, x, y)
+    return prod(kappa(k, xi, yi) for (k, xi, yi) in zip(kernel.kernels, x, y))
+end
+
+# TODO: General implementation of `kernelmatrix` and `kerneldiagmatrix`
+# Default implementation assumes 1D observations
+
+function kernelmatrix!(
+    K::AbstractMatrix,
+    kernel::TensorProduct,
+    X::AbstractMatrix;
+    obsdim::Int = defaultobs
+)
+    obsdim ∈ (1, 2) || "obsdim should be 1 or 2 (see docs of kernelmatrix))"
+
+    featuredim = feature_dim(obsdim)
+    if !check_dims(K, X, X, featuredim, obsdim)
+        throw(DimensionMismatch("Dimensions of the target array K $(size(K)) are not consistent with X $(size(X))"))
+    end
+
+    size(X, featuredim) == length(kernel) ||
+        error("number of kernels and groups of features are not consistent")
+
+    kernelmatrix!(K, kernel.kernels[1], selectdim(X, featuredim, 1))
+    for (k, Xi) in Iterators.drop(zip(kernel.kernels, eachslice(X; dims = featuredim)), 1)
+        K .*= kernelmatrix(k, Xi)
+    end
+
+    return K
+end
+
+function kernelmatrix!(
+    K::AbstractMatrix,
+    kernel::TensorProduct,
+    X::AbstractMatrix,
+    Y::AbstractMatrix;
+    obsdim::Int = defaultobs
+)
+    obsdim ∈ (1, 2) || error("obsdim should be 1 or 2 (see docs of kernelmatrix))")
+
+    featuredim = feature_dim(obsdim)
+    if !check_dims(K, X, Y, featuredim, obsdim)
+        throw(DimensionMismatch("Dimensions $(size(K)) of the target array K are not consistent with X ($(size(X))) and Y ($(size(Y)))"))
+    end
+
+    size(X, featuredim) == length(kernel) ||
+        error("number of kernels and groups of features are not consistent")
+
+    kernelmatrix!(K, kernel.kernels[1], selectdim(X, featuredim, 1),
+                  selectdim(Y, featuredim, 1))
+    for (k, Xi, Yi) in Iterators.drop(zip(kernel.kernels,
+                                          eachslice(X; dims = featuredim),
+                                          eachslice(Y; dims = featuredim)), 1)
+        K .*= kernelmatrix(k, Xi, Yi)
+    end
+
+    return K
+end
+
+# mapreduce with multiple iterators requires Julia 1.2 or later.
+
+function kernelmatrix(
+    kernel::TensorProduct,
+    X::AbstractMatrix;
+    obsdim::Int = defaultobs
+)
+    obsdim ∈ (1, 2) || error("obsdim should be 1 or 2 (see docs of kernelmatrix))")
+
+    featuredim = feature_dim(obsdim)
+    if !check_dims(X, X, featuredim, obsdim)
+        throw(DimensionMismatch("Dimensions of the target array K $(size(K)) are not consistent with X $(size(X))"))
+    end
+
+    size(X, featuredim) == length(kernel) ||
+        error("number of kernels and groups of features are not consistent")
+
+    return mapreduce((x, y) -> x .* y,
+                     zip(kernel.kernels, eachslice(X; dims = featuredim))) do (k, Xi)
+        kernelmatrix(k, Xi)
+    end
+end
+
+function kernelmatrix(
+    kernel::TensorProduct,
+    X::AbstractMatrix,
+    Y::AbstractMatrix;
+    obsdim::Int = defaultobs
+)
+    @assert obsdim ∈ (1, 2) || error("obsdim should be 1 or 2 (see docs of kernelmatrix))")
+
+    featuredim = feature_dim(obsdim)
+    if !check_dims(X, Y, featuredim, obsdim)
+        throw(DimensionMismatch("Dimensions $(size(K)) of the target array K are not consistent with X ($(size(X))) and Y ($(size(Y)))"))
+    end
+
+    size(X, featuredim) == length(kernel) ||
+        error("number of kernels and groups of features are not consistent")
+
+    return mapreduce((x, y) -> x .* y,
+                     zip(kernel.kernels,
+                         eachslice(X; dims = featuredim),
+                         eachslice(Y; dims = featuredim))) do (k, Xi, Yi)
+        kernelmatrix(k, Xi, Yi)
+    end
+end
+
+function kerneldiagmatrix!(
+    K::AbstractVector,
+    kernel::TensorProduct,
+    X::AbstractMatrix;
+    obsdim::Int = defaultobs
+)
+    obsdim ∈ (1, 2) || error("obsdim should be 1 or 2 (see docs of kernelmatrix))")
+    if length(K) != size(X, obsdim)
+        throw(DimensionMismatch("Dimensions of the target array K $(size(K)) are not consistent with X $(size(X))"))
+    end
+
+    featuredim = feature_dim(obsdim)
+    size(X, featuredim) == length(kernel) ||
+        error("number of kernels and groups of features are not consistent")
+
+    kerneldiagmatrix!(K, kernel.kernels[1], selectdim(X, featuredim, 1))
+    for (k, Xi) in Iterators.drop(zip(kernel.kernels, eachslice(X; dims = featuredim)), 1)
+        K .*= kerneldiagmatrix(k, Xi)
+    end
+
+    return K
+end
+
+function kerneldiagmatrix(
+    kernel::TensorProduct,
+    X::AbstractMatrix;
+    obsdim::Int = defaultobs
+)
+    obsdim ∈ (1,2) || error("obsdim should be 1 or 2 (see docs of kernelmatrix))")
+
+    featuredim = feature_dim(obsdim)
+    size(X, featuredim) == length(kernel) ||
+        error("number of kernels and groups of features are not consistent")
+
+    return mapreduce((x, y) -> x .* y,
+                     zip(kernel.kernels, eachslice(X; dims = featuredim))) do (k, Xi)
+         kerneldiagmatrix(k, Xi)
+    end
+end
+
+Base.show(io::IO, kernel::TensorProduct) = printshifted(io, kernel, 0)
+
+function printshifted(io::IO, kernel::TensorProduct, shift::Int)
+    print(io, "Tensor product of ", length(kernel), " kernels:")
+    for k in kernel.kernels
+        print(io, "\n")
+        for _ in 1:(shift + 1)
+            print(io, "\t")
+        end
+        print(io, "- ")
+        printshifted(io, k, shift + 2)
+    end
+end

src/matrix/kernelmatrix.jl

@@ -1,11 +1,19 @@
 """
-    kernelmatrix!(K::Matrix, κ::Kernel, X::Matrix; obsdim::Integer = 2)
-    kernelmatrix!(K::Matrix, κ::Kernel, X::Matrix, Y::Matrix; obsdim::Integer = 2)
+    kernelmatrix!(K::AbstractMatrix, κ::Kernel, X; obsdim::Integer = 2)
+    kernelmatrix!(K::AbstractMatrix, κ::Kernel, X, Y; obsdim::Integer = 2)
 
 In-place version of [`kernelmatrix`](@ref) where pre-allocated matrix `K` will be overwritten with the kernel matrix.
 """
 kernelmatrix!
 
+function kernelmatrix!(
+    K::AbstractMatrix,
+    kernel::Kernel,
+    X::AbstractVector;
+    obsdim::Int = defaultobs
+)
+    return kernelmatrix!(K, kernel, reshape(X, 1, :); obsdim = 2)
+end
 
 function kernelmatrix!(
         K::AbstractMatrix,
@@ -23,6 +31,16 @@ end
 kernelmatrix!(K::AbstractMatrix, κ::TransformedKernel, X::AbstractMatrix; obsdim::Int = defaultobs) =
         kernelmatrix!(K, kernel(κ), apply(κ.transform, X, obsdim = obsdim), obsdim = obsdim)
 
+function kernelmatrix!(
+    K::AbstractMatrix,
+    kernel::Kernel,
+    X::AbstractVector,
+    Y::AbstractVector;
+    obsdim::Int = defaultobs
+)
+    return kernelmatrix!(K, kernel, reshape(X, 1, :), reshape(Y, 1, :); obsdim = 2)
+end
+
 function kernelmatrix!(
         K::AbstractMatrix,
         κ::Kernel,
@@ -60,8 +78,8 @@ _kernel(κ::TransformedKernel, x::AbstractVector, y::AbstractVector; obsdim::Int
         _kernel(kernel(κ), apply(κ.transform, x), apply(κ.transform, y), obsdim = obsdim)
 
 """
-    kernelmatrix(κ::Kernel, X::Matrix; obsdim::Int = 2)
-    kernelmatrix(κ::Kernel, X::Matrix, Y::Matrix; obsdim::Int = 2)
+    kernelmatrix(κ::Kernel, X; obsdim::Int = 2)
+    kernelmatrix(κ::Kernel, X, Y; obsdim::Int = 2)
 
 Calculate the kernel matrix of `X` (and `Y`) with respect to kernel `κ`.
 `obsdim = 1` means the matrix `X` (and `Y`) has size #samples x #dimension
@@ -88,6 +106,15 @@ end
 kernelmatrix(κ::TransformedKernel, X::AbstractMatrix; obsdim::Int = defaultobs) =
         kernelmatrix(kernel(κ), apply(κ.transform, X, obsdim = obsdim), obsdim = obsdim)
 
+function kernelmatrix(
+    kernel::Kernel,
+    X::AbstractVector{<:Real},
+    Y::AbstractVector{<:Real};
+    obsdim::Int = defaultobs
+)
+    return kernelmatrix(kernel, reshape(X, 1, :), reshape(Y, 1, :); obsdim = 2)
+end
+
 function kernelmatrix(
         κ::Kernel,
         X::AbstractMatrix,
@@ -107,12 +134,20 @@ kernelmatrix(κ::TransformedKernel, X::AbstractMatrix, Y::AbstractMatrix; obsdim
         kernelmatrix(kernel(κ), apply(κ.transform, X, obsdim = obsdim), apply(κ.transform, Y, obsdim = obsdim), obsdim = obsdim)
 
 """
-    kerneldiagmatrix(κ::Kernel, X::Matrix; obsdim::Int = 2)
+    kerneldiagmatrix(κ::Kernel, X; obsdim::Int = 2)
 
 Calculate the diagonal matrix of `X` with respect to kernel `κ`
 `obsdim = 1` means the matrix `X` has size #samples x #dimension
 `obsdim = 2` means the matrix `X` has size #dimension x #samples
 """
+function kerneldiagmatrix(
+    kernel::Kernel,
+    X::AbstractVector;
+    obsdim::Int = defaultobs
+)
+    return kerneldiagmatrix(kernel, reshape(X, 1, :); obsdim = 2)
+end
+
 function kerneldiagmatrix(
         κ::Kernel,
         X::AbstractMatrix;
@@ -127,10 +162,19 @@ function kerneldiagmatrix(
 end
 
 """
-    kerneldiagmatrix!(K::AbstractVector,κ::Kernel, X::Matrix; obsdim::Int = 2)
+    kerneldiagmatrix!(K::AbstractVector, κ::Kernel, X; obsdim::Int = 2)
 
 In place version of [`kerneldiagmatrix`](@ref)
 """
+function kerneldiagmatrix!(
+    K::AbstractVector,
+    kernel::Kernel,
+    X::AbstractVector;
+    obsdim::Int = defaultobs
+)
+    return kerneldiagmatrix!(K, kernel, reshape(X, 1, :); obsdim = 2)
+end
+
 function kerneldiagmatrix!(
         K::AbstractVector,
         κ::Kernel,

test/kernels/tensorproduct.jl

@@ -0,0 +1,76 @@
+@testset "tensorproduct" begin
+    rng = MersenneTwister(123456)
+    u1 = rand(rng, 10)
+    u2 = rand(rng, 10)
+    v1 = rand(rng, 5)
+    v2 = rand(rng, 5)
+
+    # kernels
+    k1 = SqExponentialKernel()
+    k2 = ExponentialKernel()
+    kernel1 = TensorProduct(k1, k2)
+    kernel2 = TensorProduct([k1, k2])
+
+    @test kernel1.kernels === (k1, k2) === TensorProduct((k1, k2)).kernels
+
+    @testset "kappa" begin
+        for (x, y) in (((v1, u1), (v2, u2)), ([v1, u1], [v2, u2]))
+            val = k1(x[1], y[1]) * k2(x[2], y[2])
+
+            @test kernel1(x, y) == kernel2(x, y) == val
+            @test KernelFunctions.kappa(kernel1, x, y) ==
+                KernelFunctions.kappa(kernel2, x, y) == val
+        end
+    end
+
+    @testset "kernelmatrix" begin
+        X = rand(2, 10)
+        Y = rand(2, 10)
+        trueX = kernelmatrix(k1, X[1, :]) .* kernelmatrix(k2, X[2, :])
+        trueXY = kernelmatrix(k1, X[1, :], Y[1, :]) .* kernelmatrix(k2, X[2, :], Y[2, :])
+        tmp = Matrix{Float64}(undef, 10, 10)
+
+        for kernel in (kernel1, kernel2)
+            @test kernelmatrix(kernel, X) == trueX
+            @test kernelmatrix(kernel, X'; obsdim = 1) == trueX
+
+            @test kernelmatrix(kernel, X, Y) == trueXY
+            @test kernelmatrix(kernel, X', Y'; obsdim = 1) == trueXY
+
+            fill!(tmp, 0)
+            kernelmatrix!(tmp, kernel, X)
+            @test tmp == trueX
+
+            fill!(tmp, 0)
+            kernelmatrix!(tmp, kernel, X'; obsdim = 1)
+            @test tmp == trueX
+
+            fill!(tmp, 0)
+            kernelmatrix!(tmp, kernel, X, Y)
+            @test tmp == trueXY
+
+            fill!(tmp, 0)
+            kernelmatrix!(tmp, kernel, X', Y'; obsdim = 1)
+            @test tmp == trueXY
+        end
+    end
+
+    @testset "kerneldiagmatrix" begin
+        X = rand(2, 10)
+        trueval = ones(10)
+        tmp = Vector{Float64}(undef, 10)
+
+        for kernel in (kernel1, kernel2)
+            @test kerneldiagmatrix(kernel, X) == trueval
+            @test kerneldiagmatrix(kernel, X'; obsdim = 1) == trueval
+
+            fill!(tmp, 0)
+            kerneldiagmatrix!(tmp, kernel, X)
+            @test tmp == trueval
+
+            fill!(tmp, 0)
+            kerneldiagmatrix!(tmp, kernel, X'; obsdim = 1)
+            @test tmp == trueval
+        end
+    end
+end

test/runtests.jl

@@ -77,6 +77,7 @@ using KernelFunctions: metric
         include(joinpath("kernels", "piecewisepolynomial.jl"))
         include(joinpath("kernels", "rationalquad.jl"))
         include(joinpath("kernels", "scaledkernel.jl"))
+        include(joinpath("kernels", "tensorproduct.jl"))
         include(joinpath("kernels", "transformedkernel.jl"))
 
         # Legacy tests that don't correspond to anything meaningful in src. Unclear how