Some improvement and additions to MO kernels

docs/src/api.md

@@ -80,6 +80,11 @@ type enables specialised implementations of e.g. [`kernelmatrix`](@ref) for
 
 To find out more about the background, read this [review of kernels for vector-valued functions](https://arxiv.org/pdf/1106.6251.pdf).
 
+If you are interested in the matrix-kernel interpretation, Kernelfunction.jl provides a convenience function that computes the resulting kernel when applied to a pair of inputs
+```@docs
+matrixkernel
+```
+
 ## Generic Utilities
 
 KernelFunctions also provides miscellaneous utility functions.

src/KernelFunctions.jl

@@ -39,6 +39,7 @@ export ColVecs, RowVecs
 export MOInput, prepare_isotopic_multi_output_data, prepare_heterotopic_multi_output_data
 export IndependentMOKernel,
     LatentFactorMOKernel, IntrinsicCoregionMOKernel, LinearMixingModelKernel
+export matrixkernel
 
 # Reexports
 export tensor, ⊗, compose

src/matrix/kernelkroneckermat.jl

@@ -25,3 +25,37 @@ end
     k(x,x') = ∏ᵢᴰ k(xᵢ,x'ᵢ)
 """
 @inline iskroncompatible(κ::Kernel) = false # Default return for kernels
+
+@inline ismatrixkroncompatible(κ::MOKernel) = false # Default return for kernels
+@inline ismatrixkroncompatible(κ::IndependentMOKernel) = true
+@inline ismatrixkroncompatible(κ::IntrinsicCoregionMOKernel) = true
+
+function _mo_kernelmatrix_kron(::MOInputIsotopicByFeatures, K, B)
+    return Kronecker.kronecker(K, B)
+end
+
+function _mo_kernelmatrix_kron(::MOInputIsotopicByOutputs, K, B)
+    return Kronecker.kronecker(B, K)
+end
+
+function kernelkronmat(
+    k::IndependentMOKernel, x::MOI, y::MOI
+) where {MOI<:IsotopicMOInputsUnion}
+    @assert x.out_dim == y.out_dim
+    Ktmp = kernelmatrix(k.kernel, x.x, y.x)
+    mtype = eltype(Ktmp)
+    return _mo_kernelmatrix_kron(Ktmp, Eye{mtype}(x.out_dim), x)
+end
+
+function kernelkronmat(
+    k::IntrinsicCoregionMOKernel, x::MOI, y::MOI
+) where {MOI<:IsotopicMOInputsUnion}
+    @assert x.out_dim == y.out_dim
+    Ktmp = kernelmatrix(k.kernel, x.x, y.x)
+    return _mo_kernelmatrix_kron(Ktmp, k.B, x)
+end
+
+function kernelkronmat(k::MOK, x::MOI) where {MOI<:IsotopicMOInputsUnion,MOK<:MOKernel}
+    @assert ismatrixkroncompatible(κ) "The chosen kernel is not compatible for Kronecker matrices"
+    return kernelkronmat(k, x, x)
+end

src/mokernels/independent.jl

@@ -24,17 +24,43 @@ struct IndependentMOKernel{Tkernel<:Kernel} <: MOKernel
 end
 
 function (κ::IndependentMOKernel)((x, px)::Tuple{Any,Int}, (y, py)::Tuple{Any,Int})
-    if px == py
-        return κ.kernel(x, y)
-    else
-        return 0.0
-    end
+    return κ.kernel(x, y) * (px == py)
+end
+
+function _kronkernelmatrix(Ktmp, B, ::MOInputIsotopicByFeatures)
+    return kron(Ktmp, B)
 end
 
-function kernelmatrix(k::IndependentMOKernel, x::MOInput, y::MOInput)
+function _kronkernelmatrix(Ktmp, B, ::MOInputIsotopicByOutputs)
+    return kron(B, Ktmp)
+end
+
+function kernelmatrix(
+    k::IndependentMOKernel, x::MOI, y::MOI
+) where {MOI<:IsotopicMOInputsUnion}
     @assert x.out_dim == y.out_dim
-    temp = k.kernel.(x.x, permutedims(y.x))
-    return cat((temp for _ in 1:(y.out_dim))...; dims=(1, 2))
+    Ktmp = kernelmatrix(k.kernel, x.x, y.x)
+    mtype = eltype(Ktmp)
+    return _kronkernelmatrix(Ktmp, Eye{mtype}(x.out_dim), x)
+end
+
+if VERSION >= v"1.6"
+    function _kronkernelmatrix!(K, Ktmp, B, ::MOInputIsotopicByFeatures)
+        return kron!(K, Ktmp, B)
+    end
+
+    function _kronkernelmatrix!(K, Ktmp, B, ::MOInputIsotopicByOutputs)
+        return kron!(K, B, Ktmp)
+    end
+
+    function kernelmatrix!(
+        K::AbstractMatrix, k::IndependentMOKernel, x::MOI, y::MOI
+    ) where {MOI<:IsotopicMOInputsUnion}
+        @assert x.out_dim == y.out_dim
+        Ktmp = kernelmatrix(k.kernel, x.x, y.x)
+        mtype = eltype(Ktmp)
+        return _kronkernelmatrix!(K, Ktmp, Matrix{mtype}(I, x.out_dim, x.out_dim), x)
+    end
 end
 
 function Base.show(io::IO, k::IndependentMOKernel)

src/mokernels/intrinsiccoregion.jl

@@ -34,10 +34,36 @@ function IntrinsicCoregionMOKernel(; kernel::Kernel, B::AbstractMatrix)
     return IntrinsicCoregionMOKernel{typeof(kernel),typeof(B)}(kernel, B)
 end
 
+function IntrinsicCoregionMOKernel(kernel::Kernel, B::AbstractMatrix)
+    return IntrinsicCoregionMOKernel{typeof(kernel),typeof(B)}(kernel, B)
+end
+
 function (k::IntrinsicCoregionMOKernel)((x, px)::Tuple{Any,Int}, (y, py)::Tuple{Any,Int})
     return k.B[px, py] * k.kernel(x, y)
 end
 
+function matrixkernel(k::IntrinsicCoregionMOKernel, x, y)
+    return matrixkernel(k, x, y, size(k.B, 1))
+end
+
+function kernelmatrix(
+    k::IntrinsicCoregionMOKernel, x::MOI, y::MOI
+) where {MOI<:IsotopicMOInputsUnion}
+    @assert x.out_dim == y.out_dim
+    Ktmp = kernelmatrix(k.kernel, x.x, y.x)
+    return _kronkernelmatrix(Ktmp, k.B, x)
+end
+
+if VERSION >= v"1.6"
+    function kernelmatrix!(
+        K::AbstractMatrix, k::IntrinsicCoregionMOKernel, x::MOI, y::MOI
+    ) where {MOI<:IsotopicMOInputsUnion}
+        @assert x.out_dim == y.out_dim
+        Ktmp = kernelmatrix(k.kernel, x.x, y.x)
+        return _kronkernelmatrix!(K, Ktmp, k.B, x)
+    end
+end
+
 function Base.show(io::IO, k::IntrinsicCoregionMOKernel)
     return print(
         io, "Intrinsic Coregion Kernel: ", k.kernel, " with ", size(k.B, 1), " outputs"

src/mokernels/lmm.jl

@@ -1,7 +1,8 @@
 @doc raw"""
-    LinearMixingModelKernel(g, e::MOKernel, A::AbstractMatrix)
+    LinearMixingModelKernel(k::Kernel, H::AbstractMatrix)
+    LinearMixingModelKernel(Tk::AbstractVector{<:Kernel},Th::AbstractMatrix)
 
-Kernel associated with the linear mixing model.
+Kernel associated with the linear mixing model, taking a vector of `m` kernels and a `m × p` matrix H for a function with `p` outputs. Also accepts a single kernel `k` for use across all `m` basis vectors. 
 
 # Definition
 
@@ -20,6 +21,10 @@ mixing matrix of ``m`` basis vectors spanning the output space.
 struct LinearMixingModelKernel{Tk<:AbstractVector{<:Kernel},Th<:AbstractMatrix} <: MOKernel
     K::Tk
     H::Th
+    function LinearMixingModelKernel(Tk::AbstractVector{<:Kernel}, H::AbstractMatrix)
+        @assert length(Tk) == size(H, 1) "Number of kernels and number of rows in H must match"
+        return new{typeof(Tk),typeof(H)}(Tk, H)
+    end
 end
 
 function LinearMixingModelKernel(k::Kernel, H::AbstractMatrix)
@@ -32,6 +37,10 @@ function (κ::LinearMixingModelKernel)((x, px)::Tuple{Any,Int}, (y, py)::Tuple{A
     return sum(κ.H[i, px] * κ.K[i](x, y) * κ.H[i, py] for i in 1:length(κ.K))
 end
 
+function matrixkernel(k::LinearMixingModelKernel, x, y)
+    return matrixkernel(k, x, y, size(k.H, 2))
+end
+
 function Base.show(io::IO, k::LinearMixingModelKernel)
     return print(io, "Linear Mixing Model Multi-Output Kernel")
 end

src/mokernels/moinput.jl

@@ -199,6 +199,8 @@ function prepare_isotopic_multi_output_data(x::AbstractVector, y::RowVecs)
     return MOInputIsotopicByOutputs(x, size(y.X, 2)), vec(y.X)
 end
 
+IsotopicMOInputsUnion = Union{MOInputIsotopicByFeatures,MOInputIsotopicByOutputs}
+
 """
     prepare_heterotopic_multi_output_data(
         x::AbstractVector, y::AbstractVector{<:Real}, output_indices::AbstractVector{Int},

src/mokernels/mokernel.jl

@@ -4,3 +4,16 @@
 Abstract type for kernels with multiple outpus.
 """
 abstract type MOKernel <: Kernel end
+
+"""
+    matrixkernel(k::MOK, x, y)
+    matrixkernel(k::IndependentMOKernel, x, y, outputsize)
+
+Convenience function to compute the matrix kernel for two inputs `x` and `y`. The `outputsize` keyword is only required for the `IndependentMOKernel` to indicated the number of outputs. 
+"""
+function matrixkernel(k::MOKernel, x, y, outputsize)
+    @assert size(x) == size(y)
+    xMO = MOInputIsotopicByFeatures([x], outputsize)
+    yMO = MOInputIsotopicByFeatures([y], outputsize)
+    return kernelmatrix(k, xMO, yMO)
+end

src/mokernels/slfm.jl

@@ -1,5 +1,5 @@
 @doc raw"""
-    LatentFactorMOKernel(g, e::MOKernel, A::AbstractMatrix)
+    LatentFactorMOKernel(g::AbstractVector{<:Kernel}, e::MOKernel, A::AbstractMatrix)
 
 Kernel associated with the semiparametric latent factor model.
 
@@ -33,6 +33,10 @@ function (κ::LatentFactorMOKernel)((x, px)::Tuple{Any,Int}, (y, py)::Tuple{Any,
     return cov_f + κ.e((x, px), (y, py))
 end
 
+function matrixkernel(k::LatentFactorMOKernel, x, y)
+    return matrixkernel(k, x, y, size(k.A, 1))
+end
+
 function kernelmatrix(k::LatentFactorMOKernel, x::MOInput, y::MOInput)
     x.out_dim == y.out_dim || error("`x` and `y` should have the same output dimension")
     x.out_dim == size(k.A, 1) ||

temporary_script.jl

@@ -0,0 +1,78 @@
+# this script should be removed when the alternative MO kernelmatrix is accepted or rejected. 
+
+using KernelFunctions, BenchmarkTools
+using LinearAlgebra
+
+mrank = 1
+dims = (in=5, out=3)
+x = [rand(dims.in) for _ in 1:20]
+
+xMOF = KernelFunctions.MOInputIsotopicByFeatures(x, dims.out)
+xMOO = KernelFunctions.MOInputIsotopicByOutputs(x, dims.out)
+
+indk = IndependentMOKernel(GaussianKernel())
+
+Kind1 = kernelmatrix(indk, xMOF, xMOF)
+Kind2 = kernelmatrix2(indk, xMOF, xMOF)
+
+Kind1 ≈ Kind2
+# true
+
+@benchmark kernelmatrix($indk, $xMOF, $xMOF)
+# BenchmarkTools.Trial: 756 samples with 1 evaluation.
+#  Range (min … max):  6.186 ms …  11.470 ms  ┊ GC (min … max): 0.00% … 35.30%
+#  Time  (median):     6.423 ms               ┊ GC (median):    0.00%
+#  Time  (mean ± σ):   6.614 ms ± 806.792 μs  ┊ GC (mean ± σ):  2.76% ±  7.77%
+
+#    ▅▇█▇▃                                                       
+#   ▆██████▅▆▄▅▄▁▄▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▄▁▁▁▁▁▁▆██▇▇ ▇
+#   6.19 ms      Histogram: log(frequency) by time        10 ms <
+
+#  Memory estimate: 2.34 MiB, allocs estimate: 60060.
+
+@benchmark kernelmatrix2($indk, $xMOF, $xMOF)
+# BenchmarkTools.Trial: 10000 samples with 5 evaluations.
+#  Range (min … max):  6.162 μs … 341.226 μs  ┊ GC (min … max):  0.00% … 96.29%
+#  Time  (median):     6.947 μs               ┊ GC (median):     0.00%
+#  Time  (mean ± σ):   8.235 μs ±  16.802 μs  ┊ GC (mean ± σ):  12.37% ±  5.92%
+
+#    ▁▄▆▇██▇▆▅▄▃▂▁▁▁                                            ▂
+#   ▅██████████████████▇▇▆▇▅▆▆▆▅▅▆▅▅▃▄▅▅▅▁▁▃▁▃▅▅▆▆▆▇▇▇▇▇▇▇██▇▇▇ █
+#   6.16 μs      Histogram: log(frequency) by time      13.7 μs <
+
+#  Memory estimate: 34.89 KiB, allocs estimate: 7.
+
+A = randn(dims.out, mrank)
+B = A * transpose(A) + Diagonal(rand(dims.out))
+
+ickernel = IntrinsicCoregionMOKernel(GaussianKernel(), B)
+
+Kic1 = kernelmatrix(ickernel, xMOF, xMOF)
+Kic2 = kernelmatrix2(ickernel, xMOF, xMOF)
+
+Kic1 ≈ Kic2
+#true
+
+@benchmark kernelmatrix($ickernel, $xMOF, $xMOF)
+# BenchmarkTools.Trial: 1874 samples with 1 evaluation.
+#  Range (min … max):  2.522 ms …   5.424 ms  ┊ GC (min … max): 0.00% … 51.13%
+#  Time  (median):     2.601 ms               ┊ GC (median):    0.00%
+#  Time  (mean ± σ):   2.666 ms ± 369.030 μs  ┊ GC (mean ± σ):  2.01% ±  7.04%
+
+#   ▂█▆▁  ▁                                                      
+#   ███████▅▄▅▅▃▁▁▃▅▁▁▁▃▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇▇ █
+#   2.52 ms      Histogram: log(frequency) by time       570 ms <
+
+#  Memory estimate: 985.47 KiB, allocs estimate: 39639.
+
+@benchmark kernelmatrix2($ickernel, $xMOF, $xMOF)
+# BenchmarkTools.Trial: 10000 samples with 5 evaluations.
+#  Range (min … max):  6.152 μs … 322.002 μs  ┊ GC (min … max):  0.00% … 96.14%
+#  Time  (median):     6.676 μs               ┊ GC (median):     0.00%
+#  Time  (mean ± σ):   8.100 μs ±  16.959 μs  ┊ GC (mean ± σ):  12.53% ±  5.85%
+
+#    ▄▇██▇▅▄▃▂▁▁                                        ▁▁▁▁    ▂
+#   ▆██████████████▇▇▇▆▆▇▇▇▇▆▇▆▄▅▄▃▆▁▅▄▅▆▆▅▅▄▅▄▆▇▇█▇▇▇▇██████▇▇ █
+#   6.15 μs      Histogram: log(frequency) by time      13.5 μs <
+
+#  Memory estimate: 34.77 KiB, allocs estimate: 6.

test/mokernels/independent.jl

@@ -1,20 +1,41 @@
 @testset "independent" begin
-    x = MOInput([rand(5) for _ in 1:4], 3)
-    y = MOInput([rand(5) for _ in 1:4], 3)
+    outdim = 3
+    x = KernelFunctions.MOInputIsotopicByOutputs([rand(5) for _ in 1:4], outdim)
+    y = KernelFunctions.MOInputIsotopicByOutputs([rand(5) for _ in 1:4], outdim)
+    z = KernelFunctions.MOInputIsotopicByOutputs([rand(5) for _ in 1:2], outdim)
+
+    xIF = KernelFunctions.MOInputIsotopicByFeatures(x.x, outdim)
+    yIF = KernelFunctions.MOInputIsotopicByFeatures(y.x, outdim)
+    zIF = KernelFunctions.MOInputIsotopicByFeatures(z.x, outdim)
 
     k = IndependentMOKernel(GaussianKernel())
     @test k isa IndependentMOKernel
     @test k isa MOKernel
     @test k isa Kernel
     @test k.kernel isa Kernel
-    @test k(x[2], y[2]) isa Real
 
     @test kernelmatrix(k, x, y) == kernelmatrix(k, collect(x), collect(y))
-    @test kernelmatrix(k, x, x) == kernelmatrix(k, x)
 
-    x1 = MOInput(rand(5), 3) # Single dim input
-    @test k(x1[1], x1[1]) isa Real
-    @test kernelmatrix(k, x1) isa Matrix
+    ## accuracy
+    KernelFunctions.TestUtils.test_interface(k, x, y, z)
+    KernelFunctions.TestUtils.test_interface(k, xIF, yIF, zIF)
+
+    # in-place
+    if VERSION >= v"1.6"
+        K = zeros(12, 12)
+        kernelmatrix!(K, k, x, y)
+        @test K ≈ k.(x, permutedims(y))
+
+        K = zeros(12, 12)
+        kernelmatrix!(K, k, xIF, yIF)
+        @test K ≈ k.(xIF, permutedims(yIF))
+    end
+
+    # type stability (maybe move to test_interface?)
+    x2 = MOInput(rand(Float32, 4), 2)
+    @test k(x2[1], x2[2]) isa Float32
+    @test k(x2[1], x2[1]) isa Float32
+    @test eltype(typeof(kernelmatrix(k, x2))) <: Float32
 
     @test string(k) ==
           "Independent Multi-Output Kernel\n" *