Some improvement and additions to MO kernels

src/mokernels/independent.jl

@@ -23,20 +23,39 @@ struct IndependentMOKernel{Tkernel<:Kernel} <: MOKernel
     kernel::Tkernel
 end
 
+# kernel function should be symmetric
+# would really like (κ::IndependentMOKernel)((x, px)::Tuple{T,Int}, (y, py)::Tuple{T,Int}) where T, but seems to cause autodiff problems
 function (κ::IndependentMOKernel)((x, px)::Tuple{Any,Int}, (y, py)::Tuple{Any,Int})
     if px == py
         return κ.kernel(x, y)
     else
-        return 0.0
+        retType = Base.return_types(κ.kernel, (typeof(x), typeof(y)))[1]
+        return zero(retType)
     end
 end
 
+# this function never gets called it seems
 function kernelmatrix(k::IndependentMOKernel, x::MOInput, y::MOInput)
     @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))
 end
 
+export kernelmatrix2
+function kernelmatrix2(k::IndependentMOKernel, x::MOInputIsotopicByFeatures, y::MOInputIsotopicByFeatures)
+    @assert x.out_dim == y.out_dim
+    Ktmp = kernelmatrix(k.kernel, x.x, y.x)
+    mtype = eltype(Ktmp)
+    kron(Ktmp, Matrix{mtype}(I, x.out_dim, x.out_dim))
+end
+
+function kernelmatrix2(k::IndependentMOKernel, x::MOInputIsotopicByOutputs, y::MOInputIsotopicByOutputs)
+    @assert x.out_dim == y.out_dim
+    Ktmp = kernelmatrix(k.kernel, x.x, y.x)
+    mtype = eltype(Ktmp)
+    kron(Matrix{mtype}(I, x.out_dim, x.out_dim), Ktmp)
+end
+
 function Base.show(io::IO, k::IndependentMOKernel)
     return print(io, string("Independent Multi-Output Kernel\n\t", string(k.kernel)))
 end

src/mokernels/intrinsiccoregion.jl

@@ -34,10 +34,37 @@ 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})
+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
 
+# convenience function
+function (k::IntrinsicCoregionMOKernel)(x::Vector{T}, y::Vector{T}) where T <: Real
+    @assert size(x) == size(y)
+    outputsize = size(k.B, 1)
+    xMO = MOInputIsotopicByFeatures([x], outputsize)
+    yMO = MOInputIsotopicByFeatures([y], outputsize)
+    kernelmatrix(k, xMO, yMO)
+end
+
+# kernelmatrix (test for now)
+export kernelmatrix2
+function kernelmatrix2(k::IntrinsicCoregionMOKernel, x::MOInputIsotopicByFeatures, y::MOInputIsotopicByFeatures)
+    @assert x.out_dim == y.out_dim == size(k.B, 1)
+    Ktmp = kernelmatrix(k.kernel, x.x, y.x)
+    kron(Ktmp, k.B)
+end
+
+function kernelmatrix2(k::IntrinsicCoregionMOKernel, x::MOInputIsotopicByOutputs, y::MOInputIsotopicByOutputs)
+    @assert x.out_dim == y.out_dim == size(k.B, 1)
+    Ktmp = kernelmatrix(k.kernel, x.x, y.x)
+    kron(k.B, Ktmp)
+end
+
 function Base.show(io::IO, k::IntrinsicCoregionMOKernel)
     return print(
         io, "Intrinsic Coregion Kernel: ", k.kernel, " with ", size(k.B, 1), " outputs"

temporary_script.jl

@@ -0,0 +1,79 @@
+# 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

@@ -16,6 +16,12 @@
     @test k(x1[1], x1[1]) isa Real
     @test kernelmatrix(k, x1) isa Matrix
 
+    # type stability
+    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" *
           "\tSquared Exponential Kernel (metric = Euclidean(0.0))"

test/mokernels/intrinsiccoregion.jl

@@ -7,16 +7,27 @@
     A = randn(dims.out, rank)
     B = A * transpose(A) + Diagonal(rand(dims.out))
 
-    X = [(rand(dims.in), rand(1:(dims.out))) for i in 1:(dims.obs)]
+    # X = [(rand(dims.in), rand(1:(dims.out))) for i in 1:(dims.obs)]
+    x = [rand(dims.in) for _ in 1:2]
+    X = KernelFunctions.MOInputIsotopicByFeatures(x, dims.out)
 
     kernel = SqExponentialKernel()
     icoregionkernel = IntrinsicCoregionMOKernel(; kernel=kernel, B=B)
 
+    icoregionkernel2 = IntrinsicCoregionMOKernel(kernel, B)
+    @test icoregionkernel == icoregionkernel2
+
     @test icoregionkernel.B == B
     @test icoregionkernel.kernel == kernel
     @test icoregionkernel(X[1], X[1]) ≈ B[X[1][2], X[1][2]] * kernel(X[1][1], X[1][1])
     @test icoregionkernel(X[1], X[end]) ≈ B[X[1][2], X[end][2]] * kernel(X[1][1], X[end][1])
 
+    # test convenience function using kronecker product
+    @test icoregionkernel(X.x[1], X.x[2]) ≈ icoregionkernel.kernel(X.x[1], X.x[2])*B
+
+    # kernelmatrix
+    @test kernelmatrix(icoregionkernel, X) ≈ kron(kernelmatrix(kernel, X.x), B)
+
     KernelFunctions.TestUtils.test_interface(
         icoregionkernel, Vector{Tuple{Float64,Int}}; dim_out=dims.out
     )