Implemented Kronecker matrices with Kronecker.jl

Author: theogf

Project.toml

@@ -5,8 +5,10 @@ version = "0.2.4"
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+Kronecker = "2c470bb0-bcc8-11e8-3dad-c9649493f05e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
+Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
@@ -15,6 +17,7 @@ ZygoteRules = "700de1a5-db45-46bc-99cf-38207098b444"
 [compat]
 Compat = "2.2, 3.2"
 Distances = "0.8"
+Kronecker = "0.3.1"
 PDMats = "0.9"
 SpecialFunctions = "0.8, 0.9, 0.10"
 StatsFuns = "0.8, 0.9"

src/KernelFunctions.jl

@@ -1,6 +1,6 @@
 module KernelFunctions
 
-export kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!, kappa, kernelpdmat # Main matrix functions
+export kernelmatrix, kernelmatrix!, kerneldiagmatrix, kerneldiagmatrix!, kappa # Main matrix functions
 export params, duplicate, set! # Helpers
 
 export Kernel
@@ -17,6 +17,7 @@ export Transform, SelectTransform, ChainTransform, ScaleTransform, LowRankTransf
 export NystromFact, nystrom
 
 using Compat
+using Requires
 using Distances, LinearAlgebra
 using SpecialFunctions: logabsgamma, besselk
 using ZygoteRules: @adjoint
@@ -41,7 +42,6 @@ for k in ["exponential","matern","polynomial","constant","rationalquad","exponen
     include(joinpath("kernels",k*".jl"))
 end
 include("matrix/kernelmatrix.jl")
-include("matrix/kernelpdmat.jl")
 include("kernels/kernelsum.jl")
 include("kernels/kernelproduct.jl")
 include("approximations/nystrom.jl")
@@ -50,4 +50,9 @@ include("generic.jl")
 
 include("zygote_adjoints.jl")
 
+function __init__()
+    @require Kronecker="2c470bb0-bcc8-11e8-3dad-c9649493f05e" include("matrix/kernelkroneckermat.jl")
+    @require PDMats="90014a1f-27ba-587c-ab20-58faa44d9150" include("matrix/kernelpdmat.jl")
+end
+
 end

src/matrix/kernelkroeneckermat.jl renamed to src/matrix/kernelkroneckermat.jl

@@ -1,21 +1,25 @@
+using .Kronecker
+
+export kernelkronmat
+
 function kernelkronmat(
     κ::Kernel,
     X::AbstractVector,
     dims::Int
     )
-    @assert iskroncompatible(κ) "The kernel chosed is not compatible for kroenecker matrices (see `iskroncompatible()`)"
+    @assert iskroncompatible(κ) "The chosen kernel is not compatible for kroenecker matrices (see `iskroncompatible()`)"
     k = kernelmatrix(κ,reshape(X,:,1),obsdim=1)
-    K = kron()
+    kronecker(k,dims)
 end
 
 function kernelkronmat(
     κ::Kernel,
     X::AbstractVector{<:AbstractVector};
     obsdim::Int=defaultobs
     )
-    @assert iskroncompatible(κ) "The kernel chosed is not compatible for kroenecker matrices"
+    @assert iskroncompatible(κ) "The chosen kernel is not compatible for kroenecker matrices"
     Ks = kernelmatrix.(κ,X,obsdim=obsdim)
-    K = kron(Ks)
+    K = reduce(⊗,Ks)
 end
 
 

src/matrix/kernelpdmat.jl

@@ -1,3 +1,5 @@
+export kernelpdmat
+
 """
     Compute a positive-definite matrix in the form of a `PDMat` matrix see [PDMats.jl]() with the cholesky decomposition precomputed
     The algorithm recursively tries to add recursively a diagonal nugget until positive definiteness is achieved or that the noise is too big

test/test_kernelmatrix.jl

@@ -2,7 +2,7 @@ using Distances, LinearAlgebra
 using Test
 using KernelFunctions
 using PDMats
-
+using Kronecker
 dims = [10,5]
 
 A = rand(dims...)
@@ -67,4 +67,11 @@ k = SqExponentialKernel()
             # @test_throws ErrorException kernelpdmat(k,ones(100,100),obsdim=obsdim)
         end
     end
+    @testset "Kronecker" begin
+        x = range(0,1,length=10)
+        X = vcat(collect.(Iterators.product(x,x))'...)
+        @test all(collect(kernelkronmat(k,collect(x),2)).≈kernelmatrix(k,X,obsdim=1))
+        @test all(collect(kernelkronmat(k,[x,x])).≈kernelmatrix(k,X,obsdim=1))
+        @test_throws AssertionError kernelkronmat(LinearKernel(),collect(x),2)
+    end
 end