|
| 1 | +""" |
| 2 | + FBMKernel(; h::Real=0.5) |
| 3 | +
|
| 4 | +Fractional Brownian motion kernel with Hurst index h from (0,1) given by |
| 5 | +``` |
| 6 | + κ(x,y) = ( |x|²ʰ + |y|²ʰ - |x-y|²ʰ ) / 2 |
| 7 | +``` |
| 8 | +
|
| 9 | +For h=1/2, this is the Wiener Kernel, for h>1/2, the increments are |
| 10 | +positively correlated and for h<1/2 the increments are negatively correlated. |
| 11 | +""" |
| 12 | +struct FBMKernel{T<:Real} <: BaseKernel |
| 13 | + h::T |
| 14 | + function FBMKernel(;h::T=0.5) where {T<:Real} |
| 15 | + @assert h<=1.0 && h>=0.0 "FBMKernel: Given Hurst index h is invalid." |
| 16 | + return new{T}(h) |
| 17 | + end |
| 18 | +end |
| 19 | + |
| 20 | +_fbm(modX, modY, modXY, h) = (modX^h + modY^h - modXY^h)/2 |
| 21 | + |
| 22 | +function kernelmatrix(κ::FBMKernel, X::AbstractMatrix; obsdim::Int = defaultobs) |
| 23 | + @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))" |
| 24 | + modX = sum(abs2, X; dims = 3 - obsdim) |
| 25 | + modXX = pairwise(SqEuclidean(), X, dims = obsdim) |
| 26 | + return _fbm.(vec(modX), reshape(modX, 1, :), modXX, κ.h) |
| 27 | +end |
| 28 | + |
| 29 | +function kernelmatrix!(K::AbstractMatrix, κ::FBMKernel, X::AbstractMatrix; obsdim::Int = defaultobs) |
| 30 | + @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))" |
| 31 | + modX = sum(abs2, X; dims = 3 - obsdim) |
| 32 | + modXX = pairwise(SqEuclidean(), X, dims = obsdim) |
| 33 | + K .= _fbm.(vec(modX), reshape(modX, 1, :), modXX, κ.h) |
| 34 | + return K |
| 35 | +end |
| 36 | + |
| 37 | +function kernelmatrix( |
| 38 | + κ::FBMKernel, |
| 39 | + X::AbstractMatrix, |
| 40 | + Y::AbstractMatrix; |
| 41 | + obsdim::Int = defaultobs, |
| 42 | +) |
| 43 | + @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))" |
| 44 | + modX = sum(abs2, X, dims=3-obsdim) |
| 45 | + modY = sum(abs2, Y, dims=3-obsdim) |
| 46 | + modXY = pairwise(SqEuclidean(), X, Y,dims=obsdim) |
| 47 | + return _fbm.(vec(modX), reshape(modY, 1, :), modXY, κ.h) |
| 48 | +end |
| 49 | + |
| 50 | +function kernelmatrix!( |
| 51 | + K::AbstractMatrix, |
| 52 | + κ::FBMKernel, |
| 53 | + X::AbstractMatrix, |
| 54 | + Y::AbstractMatrix; |
| 55 | + obsdim::Int = defaultobs, |
| 56 | +) |
| 57 | + @assert obsdim ∈ [1,2] "obsdim should be 1 or 2 (see docs of kernelmatrix))" |
| 58 | + modX = sum(abs2, X, dims=3-obsdim) |
| 59 | + modY = sum(abs2, Y, dims=3-obsdim) |
| 60 | + modXY = pairwise(SqEuclidean(), X, Y,dims=obsdim) |
| 61 | + K .= _fbm.(vec(modX), reshape(modY, 1, :), modXY, κ.h) |
| 62 | + return K |
| 63 | +end |
| 64 | + |
| 65 | +function _kernel(κ::FBMKernel, x::Real, y::Real) |
| 66 | + _kernel(κ, [x], [y]) |
| 67 | +end |
| 68 | + |
| 69 | +## Apply kernel on two vectors ## |
| 70 | +function _kernel( |
| 71 | + κ::FBMKernel, |
| 72 | + x::AbstractVector, |
| 73 | + y::AbstractVector; |
| 74 | + obsdim::Int = defaultobs |
| 75 | + ) |
| 76 | + @assert length(x) == length(y) "x and y don't have the same dimension!" |
| 77 | + return κ(x,y) |
| 78 | +end |
| 79 | + |
| 80 | +#Syntactic Sugar |
| 81 | +function (κ::FBMKernel)(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}) |
| 82 | + modX = sum(abs2, x) |
| 83 | + modY = sum(abs2, y) |
| 84 | + modXY = sqeuclidean(x, y) |
| 85 | + return (modX^κ.h + modY^κ.h - modXY^κ.h)/2 |
| 86 | +end |
| 87 | + |
| 88 | +(κ::FBMKernel)(x::Real, y::Real) = (abs2(x)^κ.h + abs2(y)^κ.h - abs2(x-y)^κ.h)/2 |
| 89 | + |
| 90 | +function (κ::FBMKernel)(X::AbstractMatrix{<:Real}, Y::AbstractMatrix{<:Real}; obsdim::Integer=defaultobs) |
| 91 | + return kernelmatrix(κ, X, Y, obsdim=obsdim) |
| 92 | +end |
| 93 | + |
| 94 | +function (κ::FBMKernel)(X::AbstractMatrix{<:Real}; obsdim::Integer=defaultobs) |
| 95 | + return kernelmatrix(κ, X, obsdim=obsdim) |
| 96 | +end |
0 commit comments