Matrixkernel convenience functions and related performance improvements #363

Crown421 · 2021-08-28T11:04:46Z

Summary
Following up from #354 , this pull request contains the matrixkernel convenience function, some additional tests, a small fix, and performance improvements that came up along the way.

Proposed changes

Adding matrixkernel
Include changes discussed in Some improvement and additions to MO kernels #354
Fix for helper method in fbm kernel
Performance improvement for lmm kernelmatrix
Fix specialized implementation for slfm to work with both MOInput types
Related tests

What alternatives have you considered?

st-- · 2021-09-02T09:19:06Z

src/basekernels/fbm.jl

@@ -43,6 +43,8 @@ end
 _fbm(modX, modY, modXY, h) = (modX^h + modY^h - modXY^h) / 2

 _mod(x::AbstractVector{<:Real}) = abs2.(x)
+_mod(x::AbstractVector{<:AbstractVector{<:Real}}) = sum.(abs2, x)
+# two lines above could be combined into the second (dispatching on general AbstractVectors), but this (somewhat) more performant


What's this new line needed for ?

Suggested change

# two lines above could be combined into the second (dispatching on general AbstractVectors), but this (somewhat) more performant

# two lines above could be combined into the second (dispatching on general AbstractVectors), but this is (somewhat) more performant

Many other kernels work on arrays of arrays, but the fbm kernels errors as it could not find a method for _mod. I can try to recreate the example that I let me to add the line.

That'd be great! Sounds like something that should be included in the test suite (maybe for all kernels?).
It also seems quite orthogonal to the matrixkernel addition, in which case it'd be helpful if you could move that out to a separate branch/PR as well :)

Apologies for the delay, the error appears when doing

k = FBMKernel() xs = [rand(5) for _ in 1:4] kernelmatrix(k, xs)

which produces:
ERROR: LoadError: MethodError: no method matching _mod(::Array{Array{Float64,1},1})
This is needed during one of the matrixkernel calls. I can segment the fix into a separate pull request, but it is needed for this one. It would be easier for me if it could just be added via this PR.

I know it's very tempting to just keep things together in a single PR from the author's point of view, but I do highly recommend splitting it up - then e.g. it's easier for different people to review the different (smaller) PRs, it's easier to say "ok I'll spend ten minutes reviewing a ten-line PR" instead of thinking "oh when will i ever find enough time to review a several hundred lines PR", and so on!

src/mokernels/mokernel.jl

test/mokernels/intrinsiccoregion.jl

Co-authored-by: st-- <st--@users.noreply.github.com>

st-- · 2021-09-24T08:38:33Z

docs/src/api.md

@@ -80,6 +80,13 @@ 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 provides a convenience function that computes the resulting kernel for a pair of inputs directly. 


Suggested change

If you are interested in the matrix-kernel interpretation, Kernelfunction provides a convenience function that computes the resulting kernel for a pair of inputs directly.

If you are interested in the matrix-kernel interpretation, KernelFunctions provides a convenience function that computes the resulting kernel for a pair of inputs directly.

st-- · 2021-09-24T08:39:24Z

docs/src/api.md

+matrixkernel
+```
+<!-- Add when exporting IsotopicOutputs -->
+<!-- One way to look at this is that applying `matrixkernel` pairwise to a list of inputs results in a block matrix, which when flattened is the same is the `kernelmatrix` computed when providing the same list of inputs as `MOInputIsotopicByOutputs`.  -->


"flattened" = "reshape from [N,N] array of [Q,Q] arrays to [NQ, NQ] array"? would it be worth being more explicit here?

st-- · 2021-09-24T08:41:24Z

test/basekernels/fbm.jl

@@ -15,6 +15,8 @@
    @test repr(k) == "Fractional Brownian Motion Kernel (h = $(h))"
    test_ADs(FBMKernel; ADs=[:ReverseDiff])

+    # ToDo: needs tests for _mod


not specifically for _mod itself (though of course an equivalence of _mod([1, 2, 3]) == _mod([[1], [2], [3]]) would be good, but also for FBMKernel on higher-dimensional inputs

st-- · 2021-09-24T08:43:04Z

src/basekernels/fbm.jl

@@ -43,6 +43,8 @@ end
 _fbm(modX, modY, modXY, h) = (modX^h + modY^h - modXY^h) / 2

 _mod(x::AbstractVector{<:Real}) = abs2.(x)
+_mod(x::AbstractVector{<:AbstractVector{<:Real}}) = sum.(abs2, x)
+# two lines above could be combined into the second (dispatching on general AbstractVectors), but this (somewhat) more performant


I know it's very tempting to just keep things together in a single PR from the author's point of view, but I do highly recommend splitting it up - then e.g. it's easier for different people to review the different (smaller) PRs, it's easier to say "ok I'll spend ten minutes reviewing a ten-line PR" instead of thinking "oh when will i ever find enough time to review a several hundred lines PR", and so on!

st-- · 2021-09-24T08:44:07Z

src/mokernels/lmm.jl

+function kernelmatrix2(k::LinearMixingModelKernel, X, Y)
+    K = [kernelmatrix(ki, X.x, Y.x) for ki in k.K]
+    L = size(k.H, 2)
+    return reduce(
+        hcat, [reduce(vcat, [sum(k.H[:, i] .* (K .* k.H[:, j])) for i in 1:L]) for j in 1:L]
+    )
+end


For ease of reviewing, it would be great if you could move the performance improvements into a separate PR!

st-- · 2021-09-24T08:47:25Z

src/mokernels/mokernel.jl

+function matrixkernel(k::MOKernel, x, y)
+    return throw(
+        ArgumentError(
+            "This kernel does not have a specific matrixkernel implementation, you can call `matrixkernel(k, x, y, out_dim)`",


Suggested change

"This kernel does not have a specific matrixkernel implementation, you can call `matrixkernel(k, x, y, out_dim)`",

"For a $(nameof(typeof(k)), you must explicitly specify the requested output dimension: call `matrixkernel(k, x, y, out_dim)`",

st-- · 2021-09-24T08:47:40Z

src/mokernels/mokernel.jl

+function _kernelmatrix_kron_helper(::MOInputIsotopicByFeatures, Kfeatures, Koutputs)
+    return kron(Kfeatures, Koutputs)
+end
+
+function _kernelmatrix_kron_helper(::MOInputIsotopicByOutputs, Kfeatures, Koutputs)
+    return kron(Koutputs, Kfeatures)
+end


should not be part of this PR?

st-- · 2021-09-24T08:47:58Z

src/mokernels/mokernel.jl

@@ -4,3 +4,31 @@
 Abstract type for kernels with multiple outpus.
 """
 abstract type MOKernel <: Kernel end
+
+"""
+    matrixkernel(k::MOKernel, x, y)


Suggested change

matrixkernel(k::MOKernel, x, y)

matrixkernel(k::MOKernel, x, y[, out_dim])

& explanation below?

st-- · 2021-09-24T08:49:02Z

src/mokernels/slfm.jl

@@ -45,7 +58,7 @@ function kernelmatrix(k::LatentFactorMOKernel, x::MOInput, y::MOInput)
    H = [gi.(x.x, permutedims(y.x)) for gi in k.g]

    # Weighted latent kernel matrix ((N*out_dim) x (N*out_dim))
-    W_H = sum(kron(Wi, Hi) for (Wi, Hi) in zip(W, H))
+    W_H = sum(_kernelmatrix_kron_helper(x, Hi, Wi) for (Wi, Hi) in zip(W, H))


Will this work if x and y have different types?

st-- · 2021-09-24T08:49:40Z

src/mokernels/slfm.jl

+# function matrixkernel(k::LatentFactorMOKernel, x, y)
+#     return matrixkernel(k, x, y, size(k.A, 1))
+# end
+


Suggested change

# function matrixkernel(k::LatentFactorMOKernel, x, y)

# return matrixkernel(k, x, y, size(k.A, 1))

# end

but would be great to check the implementation for correctness by comparing against this!

Crown421 · 2021-09-24T16:34:16Z

As suggested, I am splitting this again into a few smaller PRs, not in the least because in doing so smaller bits spawned more larger issues.

Crown421 added 10 commits August 26, 2021 13:48

Restore additions

45096ed

Further improvements

738348c

Added missing method for _mod

7ef2db5

Add comment

4af63ac

Performance improvement lmm kernel

68ef765

Specialized for slfm

7a9ba3d

Move helper, fix slfm for both MOInput types

de51078

Switch order

61aca2b

Formatter and fallback

47b64e4

Add tests

2de7f5f

Crown421 mentioned this pull request Aug 30, 2021

Add lazy kronecker product for matrix kernels, if Kronecker.jl is loaded #364

Merged

st-- reviewed Sep 2, 2021

View reviewed changes

Crown421 and others added 2 commits September 16, 2021 11:48

Apply suggestions from code review

03338a3

Co-authored-by: st-- <st--@users.noreply.github.com>

Merge branch 'master' into mo-mk

6a65de7

st-- reviewed Sep 24, 2021

View reviewed changes

Crown421 changed the title ~~Matrixkernel convenience functions and related performancec improvements~~ Matrixkernel convenience functions and related performance improvements Sep 24, 2021

Crown421 mentioned this pull request Sep 24, 2021

Testing (and fixing) handling of AbstractVector{AbstractVector{T}} inputs #370

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Matrixkernel convenience functions and related performance improvements #363

Matrixkernel convenience functions and related performance improvements #363

Crown421 commented Aug 28, 2021

st-- Sep 2, 2021

Crown421 Sep 2, 2021

st-- Sep 2, 2021

Crown421 Sep 16, 2021

st-- Sep 24, 2021

st-- Sep 24, 2021

st-- Sep 24, 2021

st-- Sep 24, 2021

st-- Sep 24, 2021

st-- Sep 24, 2021

st-- Sep 24, 2021

st-- Sep 24, 2021

st-- Sep 24, 2021

st-- Sep 24, 2021

st-- Sep 24, 2021

Crown421 commented Sep 24, 2021

	# two lines above could be combined into the second (dispatching on general AbstractVectors), but this (somewhat) more performant
	# two lines above could be combined into the second (dispatching on general AbstractVectors), but this is (somewhat) more performant

		@@ -80,6 +80,13 @@ 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 provides a convenience function that computes the resulting kernel for a pair of inputs directly.

	"This kernel does not have a specific matrixkernel implementation, you can call `matrixkernel(k, x, y, out_dim)`",
	"For a $(nameof(typeof(k)), you must explicitly specify the requested output dimension: call `matrixkernel(k, x, y, out_dim)`",

	matrixkernel(k::MOKernel, x, y)
	matrixkernel(k::MOKernel, x, y[, out_dim])

	# function matrixkernel(k::LatentFactorMOKernel, x, y)
	# return matrixkernel(k, x, y, size(k.A, 1))
	# end

Matrixkernel convenience functions and related performance improvements #363

Are you sure you want to change the base?

Matrixkernel convenience functions and related performance improvements #363

Conversation

Crown421 commented Aug 28, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Crown421 commented Sep 24, 2021