Linear systems using cholesky decomp

[Bonnevie]

While standard linear systems seem adequately supported, there seems to be no way to differentiate expressions involving the cholesky decomposition. I believe the decomposition itself might be supported, but the cho_factor, cho_solve and solve_triangular functions in scipy are not, and without them it's hard to take advantage of it.

[mattjj]

We have the solve_triangular grad, so you can use that with np.linalg.cholesky. With that, you can write cho_solve in terms of two solve_triangular calls (the LAPACK routine dpotrs does the same thing, calling dtrsm twice).

However, that might have significant overheads (I think I've benchmarked this kind of thing before). So we should probably look into wrapping cho_solve and cho_factor as you suggest.