Optimize! causes a DomainError

[arturgower]

A domain error is triggered when running

mZero = MeanZero()                             
kern = Mat(5/2,[0.0,0.0],0.0)     
gp = GP(X,Y,mZero,kern)
optimize!(gp; method=Optim.ConjugateGradient())

I'm using

GaussianProcesses             0.4.0+
Julia Version 0.5.1

I've struggled to find a minimal example to trigger the error, the best I could do is for the values of X and Y given at the end.

It seems to me that the error comes from line 107 in GP.jl
gp.mLL = -dot((gp.y - μ),gp.alpha)/2.0 - logdet(gp.cK)/2.0 - gp.nobsv*log(2π)/2.0

when returns NaN when gp.cK has a negative determinant.

Running this code on my laptop where

GaussianProcesses             0.4.0
Julia Version 0.5.1-pre+31

generates a different error when running optimize!(gp; method=Optim.ConjugateGradient())

ERROR: MethodError: no method matching set_params!(::GaussianProcesses.GP, ::Float64; noise=true, mean=true, kern=true)
Closest candidates are:
  set_params!(::GaussianProcesses.GP, ::Array{Float64,1}; noise, mean, kern) at /home/art/.julia/v0.5/GaussianProcesses/src/GP.jl:275
  set_params!{K<:GaussianProcesses.Kernel}(::GaussianProcesses.Masked{K<:GaussianProcesses.Kernel}, ::Any) at /home/art/.julia/v0.5/GaussianProcesses/src/kernels/masked_kernel.jl:55 got unsupported keyword arguments "noise", "mean", "kern"
 in #optimize!#19(::Bool, ::Bool, ::Bool, ::Optim.ConjugateGradient{Void,Optim.##29#31,LineSearches.#hagerzhang!}, ::Array{Any,1}, ::Function, ::GaussianProcesses.GP) at /home/art/.julia/v0.5/GaussianProcesses/src/optimize.jl:20
 in (::GaussianProcesses.#kw##optimize!)(::Array{Any,1}, ::GaussianProcesses.#optimize!, ::GaussianProcesses.GP) at ./<missing>:0

I'm using the data

Y = [-0.000160691 0.000561494 -0.000308228 4.14104e-5 0.000306943 6.24922e-5 -0.00013596 9.83276e-5 -0.000105637 -2.11221e-5 0.00373866 0.000200135 -0.000462546 -0.000230539 0.0003362 -0.000120488 0.000201228 0.000141567 6.60807e-5 -0.000240906 0.00527562 0.00112132 0.000880385 -0.000602714 -0.000203268 2.6165e-5 0.000257117 0.000272523 -0.000526565 -0.000142842 0.0113493 0.00218881 0.000720203 0.000591266 0.000606935 0.000629937 1.19301e-5 0.000336753 1.32784e-5 5.98963e-5 0.00307527 0.00158222 0.000883546 0.000434129 -0.000172114 0.000570647 -0.000293091 -0.000187017 0.000111851 0.00037517][:]

X =[0.01 0.2
0.01 0.4
0.01 0.6
0.01 0.8
0.01 1.0
0.01 1.2
0.01 1.4
0.01 1.6
0.01 1.8
0.01 2.0
0.02 0.2
0.02 0.4
0.02 0.6
0.02 0.8
0.02 1.0
0.02 1.2
0.02 1.4
0.02 1.6
0.02 1.8
0.02 2.0
0.03 0.2
0.03 0.4
0.03 0.6
0.03 0.8
0.03 1.0
0.03 1.2
0.03 1.4
0.03 1.6
0.03 1.8
0.03 2.0
0.04 0.2
0.04 0.4
0.04 0.6
0.04 0.8
0.04 1.0
0.04 1.2
0.04 1.4
0.04 1.6
0.04 1.8
0.04 2.0
0.05 0.2
0.05 0.4
0.05 0.6
0.05 0.8
0.05 1.0
0.05 1.2
0.05 1.4
0.05 1.6
0.05 1.8
0.05 2.0]'

[chris-nemeth]

This is interesting. It appears to be a numerical stability issue caused by setting the observation noise too low. Currently, if you don't specify the logNoise then it defaults to -1e8. We've never had problems with this before, but now the variance (exp(2*logNoise)) seems to be rounded to 0.0 and this is causing issues with Optim. For now, I've changed the default noise on the master branch to -5.0 and we'll look into this in more detail.

For the time being, if you run

gp = GP(X,Y,mZero,kern,-5.0)

everything should work. Let me know if you're still having problems.