Merge branch 'master' into remove-transform

Author: theogf

.github/workflows/TagBot.yml

@@ -0,0 +1,11 @@
+name: TagBot
+on:
+  schedule:
+    - cron: 0 * * * *
+jobs:
+  TagBot:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: JuliaRegistries/TagBot@v1
+        with:
+          token: ${{ secrets.GITHUB_TOKEN }}

Project.toml

@@ -7,25 +7,29 @@ Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 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"
 ZygoteRules = "700de1a5-db45-46bc-99cf-38207098b444"
 
 [compat]
 Compat = "2.2, 3"
 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"
 ZygoteRules = "0.2"
 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

@@ -61,6 +61,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
 export transform
 export params, duplicate, set! # Helpers
 
@@ -16,12 +16,14 @@ export TransformedKernel, ScaledKernel
 
 export Transform, SelectTransform, ChainTransform, ScaleTransform, LowRankTransform, IdentityTransform, FunctionTransform
 
+export NystromFact, nystrom
+
 using Compat
+using Requires
 using Distances, LinearAlgebra
 using SpecialFunctions: logabsgamma, besselk
 using ZygoteRules: @adjoint
 using StatsFuns: logtwo
-using PDMats: PDMat
 using InteractiveUtils: subtypes
 
 const defaultobs = 2
@@ -44,12 +46,17 @@ end
 include("kernels/transformedkernel.jl")
 include("kernels/scaledkernel.jl")
 include("matrix/kernelmatrix.jl")
-include("matrix/kernelpdmat.jl")
 include("kernels/kernelsum.jl")
 include("kernels/kernelproduct.jl")
+include("approximations/nystrom.jl")
 
 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/approximations/nystrom.jl

@@ -0,0 +1,102 @@
+# Following the algorithm by William and Seeger, 2001
+# Cs is equivalent to X_mm and C to X_mn
+
+function sampleindex(X::AbstractMatrix, r::Real; obsdim::Integer=defaultobs)
+    0 < r <= 1 || throw(ArgumentError("Sample rate `r` must be in range (0,1]"))
+    n = size(X, obsdim)
+    m = ceil(Int, n*r)
+    S = StatsBase.sample(1:n, m; replace=false, ordered=true)
+    return S
+end
+
+function nystrom_sample(k::Kernel, X::AbstractMatrix, S::Vector{<:Integer}; obsdim::Integer=defaultobs)
+    obsdim ∈ [1, 2] || throw(ArgumentError("`obsdim` should be 1 or 2 (see docs of kernelmatrix))"))
+    Xₘ = obsdim == 1 ? X[S, :] : X[:, S]
+    C = k(Xₘ, X; obsdim=obsdim)
+    Cs = C[:, S]
+    return (C, Cs)
+end
+
+function nystrom_pinv!(Cs::Matrix{T}, tol::T=eps(T)*size(Cs,1)) where {T<:Real}
+    # Compute eigendecomposition of sampled component of K
+    QΛQᵀ = LinearAlgebra.eigen!(LinearAlgebra.Symmetric(Cs))
+
+    # Solve for D = Λ^(-1/2) (pseudo inverse - use tolerance from before factorization)
+    D = QΛQᵀ.values
+    λ_tol = maximum(D)*tol
+
+    for i in eachindex(D)
+        @inbounds D[i] = abs(D[i]) <= λ_tol ? zero(T) : one(T)/sqrt(D[i])
+    end
+
+    # Scale eigenvectors by D
+    Q = QΛQᵀ.vectors
+    QD = LinearAlgebra.rmul!(Q, LinearAlgebra.Diagonal(D))  # Scales column i of Q by D[i]
+
+    # W := (QD)(QD)ᵀ = (QΛQᵀ)^(-1)  (pseudo inverse)
+    W = QD*QD'
+
+    # Symmetrize W
+    return LinearAlgebra.copytri!(W, 'U')
+end
+
+@doc raw"""
+    NystromFact
+
+Type for storing a Nystrom factorization. The factorization contains two fields: `W` and
+`C`, two matrices satisfying:
+```math
+\mathbf{K} \approx \mathbf{C}^{\intercal}\mathbf{W}\mathbf{C}
+```
+"""
+struct NystromFact{T<:Real}
+    W::Matrix{T}
+    C::Matrix{T}
+end
+
+function NystromFact(W::Matrix{<:Real}, C::Matrix{<:Real})
+    T = Base.promote_eltypeof(W, C)
+    return NystromFact(convert(Matrix{T}, W), convert(Matrix{T}, C))
+end
+
+@doc raw"""
+    nystrom(k::Kernel, X::Matrix, S::Vector; obsdim::Int=defaultobs)
+
+Computes a factorization of Nystrom approximation of the square kernel matrix of data
+matrix `X` with respect to kernel `k`. Returns a `NystromFact` struct which stores a
+Nystrom factorization satisfying:
+```math
+\mathbf{K} \approx \mathbf{C}^{\intercal}\mathbf{W}\mathbf{C}
+```
+"""
+function nystrom(k::Kernel, X::AbstractMatrix, S::Vector{<:Integer}; obsdim::Int=defaultobs)
+    C, Cs = nystrom_sample(k, X, S; obsdim=obsdim)
+    W = nystrom_pinv!(Cs)
+    return NystromFact(W, C)
+end
+
+@doc raw"""
+    nystrom(k::Kernel, X::Matrix, r::Real; obsdim::Int=defaultobs)
+
+Computes a factorization of Nystrom approximation of the square kernel matrix of data
+matrix `X` with respect to kernel `k` using a sample ratio of `r`.
+Returns a `NystromFact` struct which stores a Nystrom factorization satisfying:
+```math
+\mathbf{K} \approx \mathbf{C}^{\intercal}\mathbf{W}\mathbf{C}
+```
+"""
+function nystrom(k::Kernel, X::AbstractMatrix, r::Real; obsdim::Int=defaultobs)
+    S = sampleindex(X, r; obsdim=obsdim)
+    return nystrom(k, X, S; obsdim=obsdim)
+end
+
+"""
+    nystrom(CᵀWC::NystromFact)
+
+Compute the approximate kernel matrix based on the Nystrom factorization.
+"""
+function kernelmatrix(CᵀWC::NystromFact{<:Real})
+    W = CᵀWC.W
+    C = CᵀWC.C
+    return C'*W*C
+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/runtests.jl

@@ -5,6 +5,7 @@ using Random
 
 @testset "KernelFunctions" begin
 include("test_kernelmatrix.jl")
+include("test_approximations.jl")
 include("test_constructors.jl")
 # include("test_AD.jl")
 include("test_transform.jl")

test/test_approximations.jl

@@ -0,0 +1,15 @@
+using Distances, LinearAlgebra
+using Test
+using KernelFunctions
+
+dims = [10,5]
+X = rand(dims...)
+k = SqExponentialKernel()
+@testset "Kernel Matrix Approximations" begin
+    @testset "Nystrom" begin
+        for obsdim in [1, 2]
+            @test kernelmatrix(k, X; obsdim=obsdim) ≈ kernelmatrix(nystrom(k, X, 1.0; obsdim=obsdim))
+            @test kernelmatrix(k, X; obsdim=obsdim) ≈ kernelmatrix(nystrom(k, X, collect(1:dims[obsdim]); obsdim=obsdim))
+        end
+    end
+end

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...)
@@ -84,4 +84,11 @@ kt = transform(SqExponentialKernel(),s)
             # @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