14. Optimization

There are three main computational steps when doing a Gaussian process fit with lsqfitgp:

• Compute the prior covariance matrix using the kernel. This is \(O((n + m)^2)\) where n is the number of datapoints and m the number of additional points where the posterior is computed.

• Solve a linear system of size n. This is \(O(n^3)\).

• Take random samples from the posterior. This is \(O(m^3 + s)\) where s is the number of samples taken.

Since usually \(m \gg n\) because the plot is done on a finely spaced grid, the typical bottleneck is taking the samples, i.e., calling gvar.raniter(). This problem can be bypassed by plotting only the standard deviation band instead of taking samples, but it is less informative. To make gvar.raniter() faster, use its eps option: gvar.raniter(x, eps=1e-12). This forces it to use a Cholesky decomposition instead of a diagonalization.

In general the GP methods have options for doing everything without gvar, but don’t try to use all of them mindlessly before profiling the code to know where the bottleneck actually is. Python has the module cProfile for that, and in an IPython shell you can use %run -p. If you opt out of gvars, you can use lsqfitgp.raniter() to draw samples from an explicit mean vector and covariance matrix instead of gvar.raniter().

Once you have solved eventual gvar-related issues, if you have at least some hundreds of datapoints the next bottleneck is probably in GP.predfromdata(). Making it faster is quick: select a solver different from the default one when initializing the GP object, like GP(kernel, solver='chol'). This applies also when using empbayes_fit. And don’t forget to disable the positivity check: GP(kernel, solver='chol', checkpos=False).

Finally, if you have written a custom kernel, it may become a bottleneck. For example the letter counting kernel in A custom kernel: text classification was very slow. A quick way to get a 2x improvement is disabling the symmetry check in GP: GP(kernel, checksym=False).