4
$\begingroup$

I am trying to implement a Gaussian Process and have problems to maximize the loglikelihood because the function consists of a determinant of a matrix that is dependent of 2 variables and the product of the inverse of the matrix with 2 vectors. The full problem I am working on has 6 instead of 2 variables but I already fail with the 2. Do I have to use a different representation of the loglikehood? Maybe work with the partial derivatives of the loglikelihood with respect to the hyperparamters {h0,h1}? I think there are implementations of a GP in Matlab but not for Mathematica.

 gpdata = {{-1.5, -1.8}, {-1., -1.2}, {-0.75, -0.4}, {-0.4, 
0.1}, {-0.25, 0.5}, {0., 0.8}};

kernelfunction[i_, j_, h0_, h1_] := 
 h0*h0*Exp[-(gpdata[[i, 1]] - gpdata[[j, 1]])^2/(2*h1^2)] + 
 KroneckerDelta[i, j]*0.09;

covariancematrix[h0_, h1_] = 
ParallelTable[kernelfunction[i, j, h0, h1], {i, 1, 6}, {j, 1, 6}];

loglikelihood[h0_, h1_] := -0.5*
gpdata[[All, 2]].Inverse[covariancematrix[h0, h1]].gpdata[[All, 2]] - 
0.5*Log[Det[covariancematrix[h0, h1]]] - 3*Log[2*Pi];

FindMaximum[loglikelihood[a, b], {{a, 1}, {b, 1.1}}]
Improve this question
$\endgroup$
1
  • $\begingroup$ have you succeed with your Gaussian Process code? For now, I'm first reading to understand what GP exactly is. I have interest on it because I want to fit a set of {xi,yi} data using GP, since I heard it works really well. The way I've been doing is to use a sum of polynomials to 'N' specific sections of my set of data points. Them I join the fitted polynomial curves in their intersection points to give me the final Fit to the {xi,yi}.The reason why I need to select ranges is that the sum of polynomials do not work well for the complete graph {xi,yi}. It is fine for just one plot of {xi,yi} bu $\endgroup$
    – user27755
    Apr 14, 2015 at 21:50

1 Answer 1

Reset to default
3
$\begingroup$

FindMaximum[loglikelihood[a, b], {{a, 1}, {b, 1.1}}]

FindMaximum::lstol:The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient increase in the function. You may need more than MachinePrecision digits of working precision to meet these tolerances. More…

{-4.42065, {a -> 1.45198, b -> 1.10692}}

The maximum of the likelihood function is poorly defined; that is, there are many values of the parameters that will give almost the same likelihood. Such cases often work better with two starting points.

FindMaximum[loglikelihood[a, b], {{a, 1, 2}, {b, 1.1, 1.2}}] //InputForm

{-4.420652200472999, {a -> 1.451981780418416, b -> 1.1069143572795543}}

FindMaximum[loglikelihood[a, b], {{a, 1.4, 1.5}, {b, 1.10, 1.11}}] //InputForm

{-4.420652200453491, {a -> 1.4519808900841398, b -> 1.1069162467332068}}

EDIT - I use v5.2 on a 1.6 GHz Mac PowerPC G5. I reran the analysis, with Timing. The results were the same as above; the result with the error message took 7.07 sec, the good results took 3.06 sec and 3.17 sec. The only alteration to the code was to change ParallelTable to Table.

I also ran highly altered code that minimizes -2*loglikelihood without the constant and solves for H0 = Log[h0^2]. (I looked for but failed to find analytic solutions for the inverse and determinant of the covariance matrix.)

{x,y} = Transpose@gpdata;
cov[H0_,h1_] = Exp[H0 - .5 Outer[Subtract,x/h1,x/h1]^2] + .09 IdentityMatrix@Length@x;
f[H0_?NumericQ, h1_?NumericQ] := y.LinearSolve[cov[H0,h1],y] + Log@Det@cov[H0,h1]

Timing@FindMinimum[f[A, b], {{A, .5, 1}, {b, 1., 1.2}}] //InputForm

{0.06*Second, {-2.185957997477882, {A -> 0.7458602818030642, b -> 1.1069163090371128}}}

Sqrt@Exp@A /. %[[2,2]] //InputForm

1.4519829020731116

EDIT 2 - This is on the order of twice as fast as my previous code and gives the same results.

{x,y} = Transpose@gpdata; xx = -.5 Outer[Subtract,x,x]^2; i9 = .09 IdentityMatrix@Length@x;
f2[H0_?NumericQ, h1_?NumericQ] := y.LinearSolve[#,y] + Log@Det@# &[Exp[H0 + xx/h1^2] + i9]
Improve this answer
$\endgroup$
3
  • $\begingroup$ Dear Ray Koopman, thank you very much for your help. How much time does it need to evaluate? My program is already running for 10 minutes on a Dual Core 2,4 GHz ,4GB Ram with SSD. I use Mathematica 9. $\endgroup$
    – neverender
    Nov 17, 2013 at 1:46
  • $\begingroup$ I also get different values for the loglikelihood when I run the loglikelihoodfunction with the parameters you got from the maximization. 'loglikelihood[1.451981780418416, 1.1069143572795543]=-9.93425'. Ok so you took out the constant log(2*Pi) but it is still not evaluating. Have you made other changes in the code? $\endgroup$
    – neverender
    Nov 17, 2013 at 5:05
  • $\begingroup$ Your altered code works fine on Mathematica 9, but my old code does not evaluate at all under Mathematica 9, even if I use Table instead of ParallelTable. Very strange. I will try to test your solution now on my full problem with the 6 variables. Thank you very much for your help! $\endgroup$
    – neverender
    Nov 17, 2013 at 9:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.