cholesky_decompose!

[Rob Trangucci]

Hey guys,

So I finally figured out the cholesky derivatives for our autodiff using Marcus' pattern for f: matrix -> matrix. I implemented both the matrix algebra version (in cholesky_decompose_eigen.hpp) and the scalar version (cholesky_decompose.hpp). The former algo is outlined here: http://mathoverflow.net/questions/150427/the-derivative-of-the-cholesky-factor

and the latter is about halfway through this paper: https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf

Frustratingly enough, the scalar version is about as fast wall-time-wise (I tested for 10x10 to 500x500) as the naive autodiff through Eigen's llt. The upside is that it uses about 1/10th of the memory, so I'd probably be able to run a 1000x1000 GP on my machine (albeit very slowly!), whereas I'd max out my RAM with our current version. 

The cholesky_decompose_eigen.hpp version is about 5x slower for 10x10 and likely an order of magnitude slower for bigger matrices because you end up doing a bunch of operations that are on order of dim(input matrix) ^ 3. I haven't tried any of the sparse solvers (you'll see that a bunch of the operations in chain() are sparse), so maybe that would provide a speedup. Right now I'm doing the lower triangular view of a dense matrix and solving in place. Maybe QR would be better? I'm not sure right now. My hunch is that this version has more headroom for improvement, so even if we could get it to be 10% faster than Eigen, it'd be great. The memory usage for this version is also about 1/10th of the naive AD through Eigen's cholesky.

Any thoughts on how we can speed up the Eigen version? 

Rob

[Ben Goodrich]

On Saturday, July 25, 2015 at 8:14:48 PM UTC-4, robert.trangucci wrote:

So I finally figured out the cholesky derivatives for our autodiff using Marcus' pattern for f: matrix -> matrix. I implemented both the matrix algebra version (in cholesky_decompose_eigen.hpp) and the scalar version (cholesky_decompose.hpp). The former algo is outlined here: http://mathoverflow.net/questions/150427/the-derivative-of-the-cholesky-factor

and the latter is about halfway through this paper: https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf

This is awesome. If nothing else, now that we (and by we I mean you) have 1 piece of Giles' paper implemented, it would probably be "easy" to do those other matrix decompositions and to do fwd mode and at least get a big RAM improvement. In particular, if we could do section 2.3.5 on the matrix exponential, we could utilize Eigen's matrix exponential function (which doesn't do custom scalar types currently). I remember forgetting some tricks about commutation and duplication matrices and their inverses that would probably help the other implementation. Also, as discussed on

https://github.com/stan-dev/stan/pull/1405

there are different ways of approaching operations with sparse matrices.

Ben

[Ben Goodrich]

On Saturday, July 25, 2015 at 8:14:48 PM UTC-4, robert.trangucci wrote:

The cholesky_decompose_eigen.hpp version is about 5x slower for 10x10 and likely an order of magnitude slower for bigger matrices because you end up doing a bunch of operations that are on order of dim(input matrix) ^ 3. I haven't tried any of the sparse solvers (you'll see that a bunch of the operations in chain() are sparse), so maybe that would provide a speedup. Right now I'm doing the lower triangular view of a dense matrix and solving in place. Maybe QR would be better? I'm not sure right now. My hunch is that this version has more headroom for improvement, so even if we could get it to be 10% faster than Eigen, it'd be great. The memory usage for this version is also about 1/10th of the naive AD through Eigen's cholesky.

Attached is an R file with a make_d1() function based on the accepted answer

http://mathoverflow.net/questions/150427/the-derivative-of-the-cholesky-factor

that does the equivalent of cholesky_decompose_eigen.hpp without a communation matrix, elimination matrix, or Kronecker product. If you have a lower-triangular matrix already, then a lower-triangular solve with the identity matrix is the most efficient way to get its inverse, but we might get some speed improvements from using a sparse solver, particularly if you statically allocate it so that we just have to update its non-zero elements every time the log-posterior is called.

Ben

[Bob Carpenter]

We saw the same thing w.r.t. matrix multiply in that
naively calling Eigen was about the same speed or even
a bit faster than our custom approach, but the custom
approach uses far less memory.

My suspicion is that Eigen's maintaining memory
locality better than we are and that we're spending all
the time in our algorithm pulling from memory to CPU.
This is reinforced by the fact that CppAD, which is normally
slow compared to the other autodiff algorithms is very fast
in this case relative to everything but Stan, and CppAD is
the only other package with the traits defined to tell Eigen
the sizes of autodiff variables so that it can sort out blocking
for memory locality efficiently.

I believe Marcus already had the Giles-style matrix derivatives
for left and right division.

If I understand, there are three systems in play:

1. naive autodiff of Eigen's LLT

2. analytic scalar diff (cholesky_decompose.hpp)

3. analytic matrix diff (cholesky_decompose_eigen.hpp)

* 1 and 2 have the same speed, but 3 is 5x slower.

* 2 and 3 use the same amount of memory, but 1 uses 10x as much.

For an N x N matrix, the amount of memory to store the expression
graph for approach 1 grows as O(N^3). Can it be kept
to O(N^2) for 2 and 3? Peak memory usage is going to be determined
by the LLT algorithm --- is that O(N^2) or O(N^3)? Of course, time's
going to hit O(N^3) whatever we do.

And keep in mind that a huge advantage of doing all of our ops
as doubles is that we can speed things up with GPUs, etc. I
think Eigen's planing around with CUDA implementations.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <cholesky_decompose_eigen.hpp><cholesky_decompose.hpp>

[Robert Trangucci]

This is great! I'll implement this tonight and see how much of a boost we get. 

Rob

<cholesky_decompose_eigen.R>

[Ben Goodrich]

On Sunday, July 26, 2015 at 8:28:11 PM UTC-4, robert.trangucci wrote:

This is great! I'll implement this tonight and see how much of a boost we get. 

We may need a hand-optimized forwardsove to make this competitive, but with either approach you can get a little bit more speed using the .coeffRef() method for assignment and the coeff() method for extraction, neither of which check bounds.

http://eigen.tuxfamily.org/dox/group__QuickRefPage.html

Ben

[Rob Trangucci]

I'll also look into this. Marcus had a bunch of great suggestions for the scalar version that may speed it up. 

To Bob: I had my accounting off earlier on the memory usage for the expression graph. Both 2 and 3 use O(N^2) and 1 is O(N^3). 

Rob

On Jul 26, 2015, at 8:40 PM, Ben Goodrich <goodri...@gmail.com> wrote:

[Robert Trangucci]

FYI, using coeffRef throughout the chain method in version 2 (the scalar version) yields about a 40% speed up over our current implementation on 50 and 500 dimensional matrices.

Rob

[Ben Goodrich]

On Sunday, July 26, 2015 at 10:31:06 PM UTC-4, robert.trangucci wrote:

FYI, using coeffRef throughout the chain method in version 2 (the scalar version) yields about a 40% speed up over our current implementation on 50 and 500 dimensional matrices.

Do you get another 40% if you use the .coeff() method for extraction?

[Rob Trangucci]

I tried this last night, it does impart a small speedup, but so much of the time is spent in the chain() method that it's not much of a speedup compared to all the boosts we saw from the O(N^3) loop.

That being said, I'm going to try to do away with two of the O(N^2) loops in chain() later this afternoon. Will let you know how it goes...

Rob

[Ben Goodrich]

On Monday, July 27, 2015 at 3:10:36 PM UTC-4, robert.trangucci wrote:

I tried this last night, it does impart a small speedup, but so much of the time is spent in the chain() method that it's not much of a speedup compared to all the boosts we saw from the O(N^3) loop.

I was more referring to assignment and extraction inside the chain() method for Giles' version.

Ben

[Marcus Brubaker]

Just FYI, most of my comments would apply to both scalar and Eigen versions.

Cheers,

Marcus