Merge pull request #36 from theogf/kronecker

Implemented Kronecker matrices with Kronecker.jl

Author: theogf

Project.toml

@@ -6,7 +6,7 @@ version = "0.2.4"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 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,7 +15,7 @@ ZygoteRules = "700de1a5-db45-46bc-99cf-38207098b444"
 [compat]
 Compat = "2.2, 3.2"
 Distances = "0.8"
-PDMats = "0.9"
+Requires = "1.0.1"
 SpecialFunctions = "0.8, 0.9, 0.10"
 StatsBase = "0.32"
 StatsFuns = "0.8, 0.9"
@@ -24,9 +24,11 @@ julia = "1.0"
 
 [extras]
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
+Kronecker = "2c470bb0-bcc8-11e8-3dad-c9649493f05e"
+PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Random", "Test", "FiniteDifferences", "Zygote"]
+test = ["Random", "Test", "FiniteDifferences", "Zygote", "PDMats", "Kronecker"]

docs/src/api.md

@@ -59,6 +59,7 @@ kernelmatrix!
 kerneldiagmatrix
 kerneldiagmatrix!
 kernelpdmat
+kernelkronmat
 transform
 ```
 

docs/src/userguide.md

@@ -22,6 +22,10 @@ For example:
   kernelmatrix(k,A,obsdim=2) # Return a 5x5 matrix
 ```
 
+We also support specific kernel matrices outputs:
+- For a positive-definite matrix object`PDMat` from [`PDMats.jl`](https://github.com/JuliaStats/PDMats.jl). Call `kernelpdmat(k,A,obsdim=1)`, it will create a matrix and in case of bad conditionning will add some diagonal noise until the matrix is considered PSD, it will then return a `PDMat` object. For this method to work in your code you need to include `using PDMats` first
+- For a Kronecker matrix, we rely on [`Kronecker.jl`](https://github.com/MichielStock/Kronecker.jl). We give two methods : `kernelkronmat(k,[x,y,z])` where `x` `y` and `z` are vectors which will return a `KroneckerProduct`, and `kernelkronmat(k,x,dims)` where `x` is a vector and dims and the number of features. Make sure that `k` is a vector compatible with such constructions (with `iskroncompatible`). Both method will return a . For those methods to work in your code you need to include `using Kronecker` first
+
 ## Kernel manipulation
 
 One can create combinations of kernels via `KernelSum` and `KernelProduct` or using simple operators `+` and `*`.

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,12 +17,12 @@ export Transform, SelectTransform, ChainTransform, ScaleTransform, LowRankTransf
 export NystromFact, nystrom
 
 using Compat
+using Requires
 using Distances, LinearAlgebra
 using SpecialFunctions: logabsgamma, besselk
 using ZygoteRules: @adjoint
 using StatsFuns: logtwo
 using StatsBase
-using PDMats: PDMat
 
 const defaultobs = 2
 
@@ -41,7 +41,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 +49,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,7 @@
+using PDMats: PDMat
+
+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