# FindMaximum of loglikelihood in Gaussian Process

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}}]
• 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
– user27755
Apr 14, 2015 at 21:50

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]
• 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. Nov 17, 2013 at 1:46
• 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? Nov 17, 2013 at 5:05
• 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! Nov 17, 2013 at 9:32