I'm trying to implement a periodic Gaussian Process:

GP = Table[Exp[-.1 Sin[\[Pi] (x - y)]^2], {x, 0., 3, .1}, {y, 0., 3, .1}];  


I'd like to be able to use this in the RandomVariate function, but can't because it is not positive definite.

Any help in pointing out anything I am doing wrong would be gratefully received.

Am using Mathematica, if that is relevant.

  • 2
    $\begingroup$ It looks like you may not have set up your matrix correctly. Look at GP//MatrixForm and verify that it has the values you expect. $\endgroup$ – bill s Jan 8 '15 at 12:44
  • $\begingroup$ GP = Table[ Exp[-.1 Sin[[Pi] (x - y)]^2], {x, 0., 3, .1}, {y, 0., 3, .1}] // Rationalize[#, 0.000001] &; GP // Det $\endgroup$ – chris May 8 '15 at 16:39

This is a question of numerical precision. Your matrix is symmetric, so it is positive definite if all eigenvalues are positive. Looking at the eigenvalues, however, shows that some of them are very close to zero. Check the Possible Issues section on PositiveDefiniteMatrixQ[]. Then you need to define your table as Table[...{0,3,1/10}] and use N[] with according precision. I just checked it on Alpha with

Eigenvalues[Table[Exp[-1/10 Sin[x-y]^2],{x,0,2,.1},{x,0,2,.1}]]


Eigenvalues[Table[Exp[-1/10 Sin[x-y]^2],{x,0,2,1/10},{x,0,2,1/10}]]

The first one giving the last 2 eigenvalues negative while the second one has all eigenvalues positive.

| improve this answer | |
  • $\begingroup$ Thanks Mikuszefski, that change to symbolic calculation and then taking N seems to do the trick, ie all eigenvalues positive. Sadly, Mathematica still doesn't like the covariance matrix (for use in RandomVariate on Multinormal Distribution), this time it complains about the Cholesky decomposition. I suspect I need to think of a slightly different way of doing this, given the extreme numerical instability issues. (Would upvote, but don't have enough reputation to do so at the moment) $\endgroup$ – Julian Francis Jan 9 '15 at 7:44
  • $\begingroup$ @JulianFrancis Surely you run into similar problems as the decoposition has similar requirements (Matrices need to be positive definite enough to overcome numerical roundoff). If this specific form of the matrix is not explicitly required, it is probably a good idea to choose one with somewhat bigger eigenvalues. $\endgroup$ – mikuszefski Jan 12 '15 at 17:24
  • $\begingroup$ reference.wolfram.com/language/howto/… This link will help you ;) $\endgroup$ – An old man in the sea. Aug 2 '17 at 20:16
  • $\begingroup$ @Anoldmaninthesea.I guess the link was meant for the OP? $\endgroup$ – mikuszefski Aug 3 '17 at 5:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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