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}}]