You do not give us any information about the sizes of the variables, which makes it difficult to test.
I notice that you always call kernelFunc() with (x, x, theta). If x is scalar or row vector then the result of the pdist2() call will be 0. If x is N x M for N > 1 then the result of the pdist2() will be N x N. The exp() will not change that, and multiplying by the scalar will not change that.
pinv() of N x N will be N x N . In order for pinv()*y to work, y must be N x P for some P, with the * giving an N x P result. The y' * before that would be * of a P x N, so that would be P x N * N * P, giving a P x P result. You multiply that by -0.5 and you subtract 0.5*log(abs(kernelFunc(x,x,theta))) where we have already determined that the kernelFunc returns an N x N result. P x P minus N x N will work under the circumstance that P is 1 or P is N; either way the subtraction going to return a N x N result .
Now, the objective function return value for fminsearch needs to be a scalar. For N x N to be scalar, then x would have had to have been N x M = 1 x M -- but in that case the pdist2() would return 0 making it essentially useless to make the pdist2() call.
This suggests that your formula is either fundamentally incorrect (but it does not look to be wrong), or else that you are asking for a matrix minimization, which fminsearch() cannot handle.
