OK - I added a numpy implementation of the low rank term and I think it's right...

I've opened #8 to track the fact that we haven't implemented missing data yet, but I want to continue with this implementation without worrying about the details of that too much so I've removed that from the checklist.

I took a stab at the conditional mean matrices and it wasn't obvious how to make it work so I'm also going to punt on that for now.

I realized that there must be some way to convert Appendix C of the Celerite2D paper into the conditional_mean matrix format here (or update the algorithm a bit) so I'm going to add that back to the checklist. I have a draft of an initial tutorial and I'm working on an RV case study for exoplanet so we're getting close!

@tagordon @ericagol: In writing up the docs for this, I realized something that perhaps you guys have known all along. We've been artificially restricting ourselves with these two "cases" for the covariance. If we write: R = L * L^T where L has dimensions M, K where K <= M we capture the general R version and the alpha version from our paper. The total scaling of the method will be: O(N M J^2 K^2) which also covers both cases. So, in terms of implementation, for the general case we can take the Cholesky factorization of R to get L and then we don't need all the logic that we had before to deal with the two cases. I've now implemented this on this branch and it removed quite a bit of code and it supports other models where 1 < K < M! It's also interesting to think about what kinds of models correspond to low rank systems like this.

@dfm Nice; I hadn't realized this!

@dfm I guess I hadn’t thought it through that far! That works out nicely.