Make GP-NARX optimization great again

[RJT1990]

Problem:

• Silly using default scipy L-BFGS for optimization since it uses numerical gradients and line-search which require multiple inversions of the covariance matrix. Need to replace with analytical gradients.

More generally:

• Need to profile the entire code and find out where the bottlenecks are, and fix them.

Other potential options:

• Multiple starting points for GP-NARX -> pick the one which ends up in the best region.