Add Fractional Brownian motion kernel

[sharanry]

Issue #44

I am trying to implement Fractional Brownian motion covariance given by

    κ(x,y) =  ( |x|²ʰ + |z|²ʰ - |x-z|²ʰ ) / 2

with Hurst index h from (0,1).

I don't think the current structure of kappa and metric can be used directly implement this. I am considering implementing a modified kappa function to serve such cases. Any thoughts on this?

[sharanry]

@willtebbutt @theogf @devmotion

[willtebbutt]

Yeah -- doesn't look doable with kappa and metric

The best way to go about implementing this will be to implement kernelmatrix and kerneldiagmatrix directly. The kappa and metric abstraction is just a convenient way of doing that anyway, so if its not applicable to a particular problem, you can always just implement kernelmatrix and kerneldiagmatrix directly. See eg. the KernelSum implementation for an example of how this would work.

[willtebbutt]

Just a few minor things. Not properly reviewed yet.

[willtebbutt]

LGTM. I've left a few style-related comments that could do with addressing before merging though.

[sharanry]

@willtebbutt Can we merge this now?

[devmotion]

I guess it would be good to implement the inplace variants kernelmatrix! and kerneldiagmatrix! as well.

[sharanry]

@theogf @devmotion
What is the purpose of kerneldiagmatrix and kerneldiagmatrix!?
We currently have implementations in general for

kerneldiagmatrix(κ::Kernel, X::Matrix; obsdim::Int=2)
kerneldiagmatrix!(K::AbstractVector,κ::Kernel, X::Matrix; obsdim::Int=2)

however we do not have implementations for

kerneldiagmatrix(κ::Kernel, X::Matrix,Y::Matrix; obsdim::Int=2)
kerneldiagmatrix!(K::AbstractVector,κ::Kernel, X::Matrix,Y::Matrix; obsdim::Int=2)

Any particular reason begin this?

[willtebbutt]

Any particular reason begin this?

This is definitely something that I want, because it's crucial for Stheno.jl (see the docs for the rationale). I don't think this PR is the right place for it though -- I would stick with implementing the current interface, and we can re-visit later.

edit:

What is the purpose of kerneldiagmatrix and kerneldiagmatrix!?

The reason for these is computational efficiency. If you want to compute posterior marginal statistics or construct certain pseudo-point approximations, you only wind up needing the diagonal elements of certain covariance matrices.

[sharanry]

@willtebbutt So are we going to implement kerneldiagmatrix and kerneldiagmatrix! for FBMKernel before merging or not?

[theogf]

kerneldiagmatrix(k, X) aims at computing the diagonal only of kernelmatrix(k, X). As @willtebbutt pointed out, this extremely useful for Gaussian Processes (and probably other stuff)
This would not make sense to do it for kerneldiagmatrix(k, X, Y) because the diagonal of a non-square matrix is not really useful

[willtebbutt]

@willtebbutt So are we going to implement kerneldiagmatrix and kerneldiagmatrix! for FBMKernel before merging or not?

Yes, but only the single-input ones.

[sharanry]

@theogf @willtebbutt I was confused because it always seems to return an array of ones. As it probably should. But I guess this is not the right place. I will discuss this over slack.

julia> kerneldiagmatrix(CosineKernel(), rand(3,3))
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

julia> kerneldiagmatrix(SqExponentialKernel(), rand(3,3))
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

julia> kerneldiagmatrix(ExponentialKernel(), rand(3,3))
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

[theogf]

Well it's because all these kernels are based on Euclidean distance and you will always get a distance of 0, hence all the ones! If you try with the LinearKernel, it won't be the case

[sharanry]

@theogf Thanks a lot for that! Now it makes sense :)

[sharanry]

@willtebbutt Please take a look. Implementing just _kernel specific to FBMKernel seems to suffice.

[willtebbutt]

LGTM!

[willtebbutt]

Actually, sorry, I've just thought: do we actually need to implement _kernel to make kerneldiagmatrix work? Looks like the _kernel implementations near here should do the trick.

[sharanry]

@willtebbutt Yeah but that uses kappa which we don't have.