dfm/george

Sometimes HODLR predictive samples are off

st-bender opened this issue · 17 comments

st-bender commented

Hello

I tried the model fitting tutorial notebook
https://george.readthedocs.io/en/latest/tutorials/model/
using the HODLR solver for the correlated noise fit by adding
solver=HODLRSolver to the gp initialisation in In[10].

The emcee fit works well and the numbers are close to the true
values and the BasicSolver results.
However, the predictive samples look like this:
hodlr_predictions

Is there any internal reason for this or is it a bug?
Changing tol does not help really, may it be related to the small sample size (50)?
But something like this happens also on larger data sets.
HODLR ist very nice to fit my large datasets in a resonable amount of
time but the strange predictive samples kept me from really trusting it.

Cheers.

dfm commented
st-bender commented
dfm commented

Very strange! I'll look into this this afternoon. Thanks for letting me know!

dfm commented

Sorry about the delay! I'm still looking into this. I can reproduce the bug when I run the notebook but I can't seem to reproduce it in isolation so something fishy is definitely going on. I'll keep on it and update when I figure it out.

st-bender commented

Hi, no problem.
For me this also happens when I run this script via normal python.
The only thing related to the solver would be solver.apply_inverse() in predict()
https://github.com/dfm/george/blob/master/george/gp.py#L502
Or maybe something related to kernel.get_value()?
I tried with return_cov=False and plot only the mean prediction and also cache=False in predict().
Both didn't help.
Cheers.

st-bender commented

Hi,
Thanks for looking into this. I found your test branch and gave it a go.
The tests don't give any error here too.
However, if I put compute() before set_parameter_vector(), the first assertion in https://github.com/dfm/george/blob/predict-bug/tests/test_solvers.py#L119 already fails, i.e. the alphas differ.
HTH, Stefan.

alexchgithub commented

Hi, I might have hit the same issue. Is it known what causes this?

dfm commented

I haven't had a chance to look into this in detail yet, but it's definitely a bug! I'll try to get to it next week (I'm just getting started back up after a cross-country move), but let me know if you come up with anything.

alexchgithub commented

I am not setup to debug c++, so not sure, but I think there is some caching in kernel classes with HODLR solver. If I re-create kernel, model=george.GP(kernel, solver = george.HODLRSolver), and then model.set_parameter_vector() it seems to work. If I do not recreate kernel seems not to. Also, emcee.EnsembleSampler seems to be wrong with HODLR solver, which would be related I guess, but it is harder to pinpont at the moment.

dfm commented
alexchgithub commented

For what it's worth, if I add self.kernel.get_value(self._x) in GP.compute, before self.solver.compute it seems to work. I think it works because it triggers creation of KernelInterface with new parameters in Kernel.kernel in kernels.py. Basic solver does that explicitely in BasicSolver.compute in basic.py. But hodlr solver is in c++, so I am not sure it is calling it and therefore KernelInterface creation is not triggered. I hope you find this useful and this is not a red herring.

dfm commented

That was exactly the clue I needed! Thanks, @alexchgithub. I've started working on a fix on the predict-bug branch. I'll post here again when I have a chance to implement some tests.

dfm commented

I've just pushed an update to the predict-bug branch that seems to fix this problem. There were indeed a few problems with the caching setup, but I'm really not entirely sure why this wasn't also a problem for the basic solver!

Take a look at that branch, see if it helps for your problems, and please let me know!

alexchgithub commented

Thanks for looking and fixing. I tried to locally build the branch but did not manage so far (I am on Windows). Will try again later. If you know some workaround or indeed if it does not work on Windows for some reason please let me know.

dfm commented

I'm not a Windows expert, but if you continue having issues with the Windows build, please open a separate issue and we'll try to figure it out.

st-bender commented

Hi!
Sorry for not answering lately but I saw you got some advise already and I also followed the updates.
I checked out the predict-bug branch and it really seems to solve the issue! Thanks for looking into it.

Some time ago you also wrote a test case for it, I found that it now works if you copy() the kernel object when creating the GP object. The following snippet now passes cleanly, after preparing Model, t, y, and yerr as in the model tutorial:

from copy import copy

params = dict([('amp', -0.9084520230247749),
               ('location', 0.014424970078361431),
               ('log_sigma2', -0.91748937072607739)])

kernel = kernels.ConstantKernel(log_constant=-2.826564897326464)
kernel *= kernels.Matern32Kernel(7.369719781508781)

gp0 = george.GP(copy(kernel), mean=Model(**params))
gp = george.GP(copy(kernel), mean=Model(**params), solver=george.HODLRSolver)

gp0.compute(t, yerr)
gp.compute(t, yerr)

p = np.array([-0.89532215,  0.04577885, -1.32793447, -1.32874009,  3.39837512])
gp0.set_parameter_vector(p)
gp.set_parameter_vector(p)

x = np.linspace(-5, 5, 500)

for _ in range(10):
    alpha0 = gp0.apply_inverse(y)
    alpha = gp.apply_inverse(y)
    assert np.allclose(alpha, alpha0)

for _ in range(10):
    mu0, cov0 = gp0.predict(y, x)
    mu, cov = gp.predict(y, x)
    assert np.allclose(mu, mu0)
dfm commented

Yes. The problem had something to do with the kernel thinking that it wasn't "dirty" because it had been used in another place. This shouldn't be a problem anymore and the performance hit was negligible. That's what I get for trying to be clever :-\

Thanks for tracking this down!