initial working version of SolveCholesky and MatrixInverseCholesky

Hi

Here is a working implementation of a Cholesky based linear solver and a matrix inversion based on a Cholesky factor. This allow for effectively precomputing the solutions of linear systems (with psd matrix) without the need of ever storing a matrix inverse.

Currently, there is no gradient for the Cholesky inversion. While writing tests for my implementation, I discovered a possible bug in gradients of matrix inversions, where the matrix is computed from a Cholesky, see #2066 

I included a couple of unit tests. I am already using this code in a research project of mine and it works fine. It is much faster than always re-computing the linear solve, and storing the inverted matrix is not numerically stable in my example.

Let me know what you think, and what is needed to merge this.
I could do ```L = Cholesky(lower=True)(A)```

I dont get the comment on checking that nobody else uses L. The line was also there before. Could you elaborate a bit?