Adding Gibbskernel #374

Cyberface · 2021-09-29T07:04:57Z

Summary

This PR begins to implement the non-stationary Gibbs kernel what is being discussed in issue: #372

This is a modification of the standard Squared-Exponential kernel where the lengthscale is no longer constant but is a function of the input parameters.

In this PR I have simply added a function GibbsKernel to evaluate the gibbs kernel as well as added a test to check that when the lengthscale function just returns a constant value of 1.0 then it should be the same as the SqExponentialKernel().

The implementation of the gibbs kernel is probably going to change as I learn more about how your kernels are defined/implemented.

This is my first PR here and I am only just starting to learn julia so I'm sorry for some basic mistakes!
And I appriciate any comments to make this better thanks!

Proposed changes

Adds a new file "src/kernels/gibbskernel.jl"
include and export this in "src/KernelFunctions.jl"
add tests and include the test in "test/runtests.jl"

What alternatives have you considered?

A better implementation would be to use the kernel/transform infrastructure you have and I will work on this after this PR.

Breaking changes

None

just sketching out the code and doc string

willtebbutt

Thanks for opening this! It's a very nice well-sized PR.

The basic functionality all seems is present. The main thing that needs addressing is turning the GibbsKernel function into a struct (see comment below) -- other than that it's just formatting-related things and other minor details that need addressing.

src/kernels/gibbskernel.jl

test/kernels/gibbskernel.jl

src/kernels/gibbskernel.jl

Co-authored-by: willtebbutt <[email protected]>

Cyberface · 2021-09-29T17:26:31Z

I have updated the unit test - it was a little more confusing and complicated that I originally thought! Hopefully you can follow what I've done through the comments.

Let me know what you think, cheers!

test/kernels/gibbskernel.jl

theogf · 2021-09-29T17:33:33Z

test/kernels/gibbskernel.jl

+    k_gibbs = GibbsKernel(l_func)
+
+    # check they are equal
+    @test k_gibbs(x, y) == k(x, y)


for the test it would be great if you could be as faithful as possible to the latex formula, it makes it easier to see if it's doing the right thing!

I'm not really sure what you mean, could you elaborate or sketch out an idea of what you would like to see?

Thanks!

Something like

@test k_gibbs(x, y) == sqrt((2 * ell(x) * ell(y)) / (ell(x)^2 + ell(y)^2)) * exp( - norm(x - y)^2 / (ell(x)^2 + ell(y)^2))

n.b. probably better to use ≈ here rather than ==.

OK I see what you mean. I've added your suggestion but actually this test is failing and I'm not sure why!

I don't think I fully understand the implementation in https://github.com/Cyberface/KernelFunctions.jl/blob/gibbskernel/src/kernels/gibbskernel.jl#L37 so will take a look a bit later again.

Are you sure the using norm(x-y) is correct here? Calling this on

x = randn(2) y = randn(2) norm(x-y)

This returns a single number - shouldn't this be a 2x2 matrix?

Assuming the norm is from LinearAlgebra ?

Sorry for my confusion!

No no, this the L2 norm! So it's sqrt(sum((x_i - y_i)^2)), this is what we would expect I think

OK I get you now - I thought this was supposed to give the covariance matrix but that is the kernelmatrix function!

This is just the covariance between two points in n-dimensions.

Hopefully I got this correct :)

…ivial function for ell

added a broken test to seek help about why it's broken

theogf · 2021-09-30T11:42:12Z

src/kernels/gibbskernel.jl

+    l = hypot(lx, ly)
+    kernel = (sqrt(2 * lx * ly) / l) * with_lengthscale(SqExponentialKernel(), l)


Suggested change

l = hypot(lx, ly)

kernel = (sqrt(2 * lx * ly) / l) * with_lengthscale(SqExponentialKernel(), l)

l = hypot(lx, ly) / sqrt(2)

kernel = (sqrt(lx * ly) / l) * with_lengthscale(SqExponentialKernel(), l)

I think this should more correspond to the exact formula

Without promotion to Float64:

Suggested change

l = hypot(lx, ly)

kernel = (sqrt(2 * lx * ly) / l) * with_lengthscale(SqExponentialKernel(), l)

l = invsqrt2 * hypot(lx, ly)

kernel = (sqrt(lx * ly) / l) * with_lengthscale(SqExponentialKernel(), l)

Is this correct though because l is passed into with_lengthscale - doesn't the SqExponentialKernel already have the 1/2 in the exponential looking here

shouldn't it be

l = hypot(lx, ly) kernel = (sqrt(2 * lx * ly) / l) * with_lengthscale(SqExponentialKernel(), 1 / l^2)

Is this correct though because l is passed into with_lengthscale - doesn't the SqExponentialKernel already have the 1/2 in the exponential looking here

Yes but in the latex formula the factor 2 is not there, so this way we cancel it out!

I guess the 2 disappears because we already have l(x) + l(y) ?

Yeah the 2 would be there when l(x) = l(y) = constant.

I think I finally understand most of the things here but I am confused by with_lengthscale. Can you explain why we shouldn't square the argument to with_lengthscale? Looking at the with_lengthscale code it seems to just multiply by 1/l.

I guess the 2 disappears because we already have l(x) + l(y) ?

I'm not sure what 2 you refer to (and there is no l(x) + l(y) term anywhere 🤔). The 2 in the square root outside of the squared exponential kernel disappears because we already scale the square root of the denominator by 1/sqrt(2).

why we shouldn't square the argument to with_lengthscale

It will be squared by the squared exponential kernel. A squared exponential kernel with lengthscale l is of the form k(x, y) = exp(- (x - y)^2/(2*l^2)).

So the l(x) + l(y) term I was refering to is the denominator in the exponential i.e.
exp(-(x - y)^2 / (l(x)^2 + l(y)^2)

And when l(x) = l(y) = l = constant you get exp(-(x - y)^2 / (2l^2) ?

Also thanks for explaining to me about the with_lengthscale

Exactly, if they are constant you get exp(-(x-y)^2 / (2 * l^2)). Very explicitly, in general the implementation yields
$$k(x, y) = sqrt(l(x) * l(y)) / (sqrt(l(x)^2 + l(y)^2) / sqrt(2)) * exp(- (x - y)^2 / (2 * (sqrt(l(x)^2 + l(y)^2) / sqrt(2))^2)) = sqrt(2 * l(x) * l(y) / (l(x)^2 + l(y)^2)) * exp(- (x - y)^2 / (l(x)^2 + l(y)^2))$$
which is exactly the formula in the docstring (and the other issue).

theogf · 2021-09-30T11:46:50Z

test/kernels/gibbskernel.jl

+    ell(x) = exp(sum(sin, x))
+
+    # this is the gibbs lengthscale function.
+    # to check that we can recover the stationary SqExponentialKernel
+    # with a constant lengthscale we just set this function
+    # equal to a constant.
+    l_func(x) = ell(1)
+
+    # in order to compare to the Gibbs kernel we compute the
+    # equivalent lengthscale l(x)^2 + l(y)^2.
+    # See the denominator of the exponential term in the gibbs kernel.
+    lengthscale = hypot(l_func(x), l_func(y))
+
+    # create a SqExponentialKernel with a lengthscale
+    k = with_lengthscale(SqExponentialKernel(), lengthscale)
+
+    # create a gibbs kernel with our constant lengthscale function
+    k_gibbs = GibbsKernel(l_func)


Suggested change

ell(x) = exp(sum(sin, x))

# this is the gibbs lengthscale function.

# to check that we can recover the stationary SqExponentialKernel

# with a constant lengthscale we just set this function

# equal to a constant.

l_func(x) = ell(1)

# in order to compare to the Gibbs kernel we compute the

# equivalent lengthscale l(x)^2 + l(y)^2.

# See the denominator of the exponential term in the gibbs kernel.

lengthscale = hypot(l_func(x), l_func(y))

# create a SqExponentialKernel with a lengthscale

k = with_lengthscale(SqExponentialKernel(), lengthscale)

# create a gibbs kernel with our constant lengthscale function

k_gibbs = GibbsKernel(l_func)

# this is the gibbs lengthscale function.

ell(x) = exp(sum(sin, x))

# create a gibbs kernel with our specific lengthscale function

k_gibbs = GibbsKernel(ell)

That's more what I meant to use instead of a constant.

bugs ironed out hopefully with lots of help from the devs

Cyberface · 2021-10-01T09:27:33Z

OK I think this PR is basically done?

willtebbutt

I'm happy -- once @theogf is happy, we can merge :)

src/kernels/gibbskernel.jl

@willtebbutt

remove space thanks @willtebbutt Co-authored-by: willtebbutt <[email protected]>

theogf

Just a comment for the docstring. And could you also add the kernel in the docs (just add the kernel name in the list in docs/kernels.md)?

theogf · 2021-10-01T09:47:03Z

src/kernels/gibbskernel.jl

+kernel. The lengthscale parameter ``l`` becomes a function of
+position ``l(x)``.
+
+``l(x) = l`` then you recover the standard Squared-Exponential kernel


Suggested change

``l(x) = l`` then you recover the standard Squared-Exponential kernel

For a constant function``l(x) = c``, one recovers the standard Squared-Exponential kernel with lengthscale `0.5 c`

great thanks and that answers my other question about the square root / squaring about with_lengthscale

The lengthscale is just c, it seems: In this case we get exp(-(x-y)^2/(2*c^2)), both with our implementation and the formula it is based on which corresponds to a lengthscale of c.

@theogf

thanks @theogf Co-authored-by: Théo Galy-Fajou <[email protected]>

Cyberface · 2021-10-01T09:52:45Z

Just a comment for the docstring. And could you also add the kernel in the docs (just add the kernel name in the list in docs/kernels.md)?

Yeah sure thing. Shall I add a new section in base kernels called 'non-stationary kernels' ? Or shall I put it in 'Exponential Kernels' ?

willtebbutt · 2021-10-01T09:55:30Z

I'd just put it in alphabetical order in the base kernels list. Feels like the appropriate place to me (the base kernels aren't necessarily stationary)

Cyberface · 2021-10-01T10:31:39Z

great thanks - done!

src/kernels/gibbskernel.jl

test/kernels/gibbskernel.jl

@devmotion

thanks @devmotion Co-authored-by: David Widmann <[email protected]>

@devmotion

thanks @devmotion Co-authored-by: David Widmann <[email protected]>

Cyberface · 2021-10-01T11:55:55Z

The PR says there is a request to change something but I'm not sure what it's referring to.

devmotion · 2021-10-01T11:57:31Z

The PR says there is a request to change something but I'm not sure what it's referring to.

It was my request to fix the docstring.

test/kernels/gibbskernel.jl

devmotion

Some minor improvements to the docstring.

src/kernels/gibbskernel.jl

theogf · 2021-10-01T12:09:12Z

Great! All we need now is a version bump. In Project.toml can you change version = "0.10.20" to version = "0.10.21"

devmotion · 2021-10-01T12:46:21Z

Thanks a lot @Cyberface!

Cyberface · 2021-10-01T13:36:12Z

that's really great! Thanks so much for all your help with this! Next I'll work on getting the kernelmatrix working and eventually getting it so that the gibbs kernel length scale can be modelled with a separate GP and fit/optimise that!

Cyberface added 6 commits September 28, 2021 21:13

add gibbskernel.jl

d50f2b4

just sketching out the code and doc string

include gibbskernel.jl

2ca62f4

some modifications for gibbskernel and added a small test

2b83286

fix missing definition of ell function in test

e4d1414

fix gibbskernel test

f9f066a

add gibbskernel to test/runtests.jl

d6b53b9

willtebbutt reviewed Sep 29, 2021

View reviewed changes

Cyberface and others added 6 commits September 29, 2021 17:25

Update test/kernels/gibbskernel.jl

698c453

Co-authored-by: willtebbutt <[email protected]>

Update src/kernels/gibbskernel.jl

3a98028

Co-authored-by: willtebbutt <[email protected]>

Update src/kernels/gibbskernel.jl

0f757a3

Co-authored-by: willtebbutt <[email protected]>

implement changes based on PR discussion and update test

2e1df2d

add brackets back in

b5c6faa

fix unit test for gibbs kernel

19e6665

theogf reviewed Sep 29, 2021

View reviewed changes

Cyberface added 2 commits September 30, 2021 09:37

update gibbskernel test: implement theogf's idea about using a non-tr…

21de7d4

…ivial function for ell

ensure ell returns positive value in test

4ddafbe

added a broken test to seek help about why it's broken

theogf reviewed Sep 30, 2021

View reviewed changes

Cyberface added 2 commits September 30, 2021 21:16

gibbs kernel: fixed implementation and unittest

9fd1d14

bugs ironed out hopefully with lots of help from the devs

gibbs kernel: fix docstring

5531851

willtebbutt approved these changes Oct 1, 2021

View reviewed changes

src/kernels/gibbskernel.jl Show resolved Hide resolved

Update src/kernels/gibbskernel.jl

78faeac

remove space thanks @willtebbutt Co-authored-by: willtebbutt <[email protected]>

theogf reviewed Oct 1, 2021

View reviewed changes

Update src/kernels/gibbskernel.jl

ae8b3d4

thanks @theogf Co-authored-by: Théo Galy-Fajou <[email protected]>

add gibbskernel to docs

6b91a23

devmotion requested changes Oct 1, 2021

View reviewed changes

src/kernels/gibbskernel.jl Outdated Show resolved Hide resolved

devmotion reviewed Oct 1, 2021

View reviewed changes

test/kernels/gibbskernel.jl Outdated Show resolved Hide resolved

Cyberface and others added 4 commits October 1, 2021 12:34

Update src/kernels/gibbskernel.jl

a914957

thanks @devmotion Co-authored-by: David Widmann <[email protected]>

Update test/kernels/gibbskernel.jl

6eb4099

thanks @devmotion Co-authored-by: David Widmann <[email protected]>

gibbskernel test: formatting

79aa4f5

Merge branch 'master' into gibbskernel

0942447

devmotion reviewed Oct 1, 2021

View reviewed changes

test/kernels/gibbskernel.jl Outdated Show resolved Hide resolved

devmotion added 2 commits October 1, 2021 13:59

Fix format (hopefully)

de3e535

Add newline character

0caf2fd

devmotion reviewed Oct 1, 2021

View reviewed changes

Apply suggestions from code review

cfbca3d

devmotion approved these changes Oct 1, 2021

View reviewed changes

theogf added 2 commits October 1, 2021 14:23

Export invsqrt2

1fb7903

Version bump

f577700

devmotion merged commit 16f1149 into JuliaGaussianProcesses:master Oct 1, 2021

willtebbutt mentioned this pull request Oct 1, 2021

Implementing non-stationary gibbs kernel #372

Open

Cyberface deleted the gibbskernel branch October 1, 2021 13:36

		l = hypot(lx, ly)
		kernel = (sqrt(2 * lx * ly) / l) * with_lengthscale(SqExponentialKernel(), l)

	``l(x) = l`` then you recover the standard Squared-Exponential kernel
	For a constant function``l(x) = c``, one recovers the standard Squared-Exponential kernel with lengthscale `0.5 c`

Adding Gibbskernel #374

Adding Gibbskernel #374

Conversation

Cyberface commented Sep 29, 2021

willtebbutt left a comment

Choose a reason for hiding this comment

Cyberface commented Sep 29, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

theogf Sep 30, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Cyberface Sep 30, 2021 • edited Loading

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

theogf Sep 30, 2021 • edited Loading

Choose a reason for hiding this comment

Cyberface commented Oct 1, 2021

willtebbutt left a comment

Choose a reason for hiding this comment

theogf left a 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

Cyberface commented Oct 1, 2021

willtebbutt commented Oct 1, 2021

Cyberface commented Oct 1, 2021

Cyberface commented Oct 1, 2021

devmotion commented Oct 1, 2021

devmotion left a comment

Choose a reason for hiding this comment

theogf commented Oct 1, 2021

devmotion commented Oct 1, 2021

Cyberface commented Oct 1, 2021

theogf Sep 30, 2021 •

edited

Loading

Cyberface Sep 30, 2021 •

edited

Loading

theogf Sep 30, 2021 •

edited

Loading