Cholesky decomposition rewrite: initial work #150

Andlon · 2017-02-11T21:04:01Z

I started some work on rewriting the Cholesky decomposition in the same style as LU. On the way there, I completely rewrote the algorithm using an algorithm from Golub and Van Loan, which describes the decomposition in terms of GAXPY operations. This meant that I can leverage your efficient dot implementation in the utils module. The preliminary benchmarks show a significant improvement (cholesky_nxn is the existing cholesky() implementation):

test linalg::cholesky::cholesky_100x100                                         ... bench:     201,593 ns/iter (+/- 1,126)
test linalg::cholesky::cholesky_500x500                                         ... bench:  21,765,711 ns/iter (+/- 208,204)
test linalg::cholesky::cholesky_decompose_unpack_100x100                        ... bench:      95,016 ns/iter (+/- 337)
test linalg::cholesky::cholesky_decompose_unpack_500x500                        ... bench:   9,607,424 ns/iter (+/- 535,354)

I just thought I'd put this up here, because I think some of this knowledge can also be applied to LU, although the situation is more complicated because of the additional work on U as well as the fact that one has to pivot.

EDIT: Note that this PR includes the entirety of #142! (will rebase after #142 eventually gets merged)

AtheMathmo · 2017-02-11T21:20:32Z

These initial numbers are very promising!

I should have some time to give #142 a final review tomorrow and hopefully this PR too.

Thanks!

Andlon · 2017-02-11T21:22:17Z

Take whatever time you need :) This PR is far from done, by the way. I've only wrapped the decomposition itself, but I still plan to provide e.g. solve, inverse and det, as in the case of LU.

AtheMathmo · 2017-02-17T20:15:41Z

#142 has been merged which I think means you need to rebase.

Benchmarks show significant performance gains.

Andlon · 2017-02-18T19:00:25Z

Rebasing took more of an effort than I expected, but at least I learned something new in the process. What's remaining is mostly:

An inverse() implementation
Deprecating .cholesky()
Write documentation for the new Cholesky type
Update documentation for the decomposition module so that the table of decompositions lists Cholesky as well.

Don't know when I actually have time to do this. I might find myself with almost no time or ample time in some cases. We'll see :)

AtheMathmo

The changes so far look good to me!

I think the only thing pending is the old cholesky function deprecation?

AtheMathmo · 2017-02-22T13:31:58Z

src/internal_utils.rs

@@ -0,0 +1,52 @@
+use matrix::BaseMatrixMut;


Maybe this module would be better placed in the matrix folder, so that it is src/matrix/utils.rs?

AtheMathmo · 2017-02-22T13:34:12Z

src/matrix/decomposition/cholesky.rs

+/// with the factor itself. Instead, we may wish to
+/// solve linear systems or compute the determinant or the
+/// inverse of a symmetric positive definite matrix.
+/// In this case, see the next subsections.///


Stray /// here

AtheMathmo · 2017-02-22T13:38:37Z

src/matrix/decomposition/cholesky.rs

+                        let k_row = a.row(k).raw_slice();
+                        dot(&k_row[0 .. j], &j_row[0 .. j])
+                    };
+                    a[[k, j]] = a[[k, j]] - kj_dot;


We might want to check whether using get_unchecked here improves vectorization (I doubt it).

Please see my other comment. The vast majority of the time should be spent in the kj_dot computation, so I think there's no vectorization to speak of for this single expression (?). However, dot allows for vectorization, and indeed this is where the speed gains demonstrated by the benchmarks come from :)

I had misunderstood what was going on when I made this comment :)

Can't blame you, this algorithm is not so intuitive because it doesn't follow what you would usually do by hand when doing Gaussian elimination, which basically corresponds to the algorithm based on outer products. However, according to Golub & Van Loan, the Gaxpy-rich version should only involve about half the memory traffic or so of the outer product version.

In the future we can however make it significantly more cache-efficient by using an algorithm based on (blocked) level 3 BLAS routines. This is a strong motivation for #155, because once we have some optimized and efficient BLAS-3 routines, we can significantly accelerate the various decompositions :)

AtheMathmo · 2017-02-22T13:44:42Z

src/matrix/decomposition/cholesky.rs

+        let mut a = matrix;
+
+        // Resolve each submatrix (j .. n, j .. n)
+        for j in 0 .. n {


Inside this loop I think we might be able to tidy things up in places by using a MatrixSliceMut and get_col_unchecked calls. In theory this should help with vectorization too, but I might be wrong as columnar slice iteration isn't the most efficient.

Another thought - here we are iterating over columns on the outside and rows on the inside. This means all of our inner loops have more cache misses (a[[k, j]] access pattern is worse than a[[j, k]] here). Is it possible to transpose the algorithm so that each inner operation is over column access?

These are just thoughts from someone who has spent 2 minutes reading the code. Please ignore these if not relevant.

Hmm, I just realized that the comment on line 128 is somewhat misleading. This is a left-looking algorithm based on GAXPY, instead of the usual right-looking algorithm based on outer products, which is probably how most people think of the LU decomposition. In this case it's not quite right to talk about the submatrix (j ..n, j .. n). I will remove this comment to avoid confusion.

If you look closely, you'll see that for each k, we compute the dot product of parts of the kth and jth row. This is where the majority of the time of the algorithm is spent. Since we're computing the dot products of two contiguous slices, this should be fairly cache friendly. Now, you'll also notice that when kj_dot has been computed, the last item in the kth row that was accessed in the dot product is the item with index j - 1, which means that the next item would be j, which corresponds exactly to a[[k, j]]. In other words, at the end of computing kj_dot, a[[k, j]] is probably already in cache due to prefetching.

My understanding is thus that this whole operation should be very cache friendly. It is of course possibly that I've misunderstood some parts :)

Also, I'm not entirely sure if I understand your comment about MatrixSliceMut. Could you perhaps elaborate?

I think the most sensible way to make the code cleaner is to rewrite rows 135 to 143 (the k loop) as a call to a gaxpy function. In fact, this is part of my motivation for #155, because this could be optimized individually and reused in several parts of the library. However, a difficulty here is that the "vector" y in the gaxpy expression y <- y + Ax is not a contiguous slice, but rather a column in a matrix.

Thanks for explaining - I agree that the current algorithm is better than a transposed version.

What I meant with the MatrixSliceMut comment is not applicable now that I understand what's going on. I was trying to make the case that if we slice from [j.. , j..] we can get slightly nicer code by doing things like:

*a.col_mut(k) -= kj_dot;

But given that we need to take the whole row of a at line 137-138 I think this would cause borrowing issues. And it only really saves us on that one line.

Ah, yes. That would be a good idea in that case. However, I think it's still possible doing something like that using split_at_mut(), but I think perhaps the code doesn't get clearer from this, as you'd have to spend additional lines of code in appeasing the borrow checker.

AtheMathmo · 2017-02-22T13:49:01Z

src/matrix/decomposition/cholesky.rs

-use libnum::Float;
+use libnum::{Zero, Float};
+
+/// Cholesky decomposition.


Really great documentation here!

Andlon · 2017-02-25T14:14:03Z

About the internal_utils module, moving it into matrix as utils is somewhat problematic, because it overlaps with e.g. use utils:: in matrix/mod.rs itself, and seems to create some potential confusion in the code.

I actually think utils should be made entirely private, in which case internal_utils could be merged with utils, but this would be a breaking change.

In any case, I'm not entirely sure what to do about it. Do you have any suggestions, @AtheMathmo?

I think I'm otherwise done with the rest! :)

AtheMathmo · 2017-03-07T21:34:46Z

Sorry for the delay with this.

I agree with what you said above. We can't really move the internal_utils module. I also think that making all of utils private might be the best move. It has been kept private thus far as I needed to use some of the functionality within rusty-machine. Now that the API is a little more developed this probably isn't necessary.

For now I think we should keep things as you have them in this PR - evaluating whether users might need the utils functions could be difficult.

I'll try to find time to give this PR one more look through either later today or tomorrow. From what I remember it should be ready to merge :)

Andlon · 2017-03-07T22:07:50Z

No worries!

I just realized I forgot to mention something I wanted to discuss, that I think impacts the merge.

I actually made the solve function return a Vector directly, because if the decomposition is successful, it should always be safe to solve the system. But this assumes that the decomposition and the triangular solve actually use the same criterion for "close to zero". I think they do this at the moment, but it seems a bit fragile. Perhaps to be on the same side we should make solve (and inverse) return Result? This is also more consistent with the LU decomposition API, which is probably a good idea.

AtheMathmo · 2017-03-07T22:13:02Z

I agree that using a Result is probably better. The consistency itself is a big enough motivating factor I think. It seems a good idea to have a Result heavy API with decomposition work anyway - lots of things can go wrong!

Andlon · 2017-03-07T22:16:24Z

All right, I'll take care of it as soon as I have time. Might not be before after next week. I'm in the last leg before thesis submission so I'm deliberately trying not to get too distracted by other endeavours at the moment! 😛

AtheMathmo · 2017-03-07T22:17:27Z

No worries. I've been swallowed by school work for the last couple weeks so you have my sympathy. Good luck with the submission!

Andlon · 2017-03-23T13:55:22Z

OK, I finally found the time to make the changes to solve and inverse so that they return Result. I suppose this PR should be complete now, barring any mistakes!

AtheMathmo

One small point.

Once we resolved that I am happy to merge :)

AtheMathmo · 2017-03-24T12:00:18Z

src/matrix/decomposition/cholesky.rs

+    }
+
+    /// Computes the determinant of the decomposed matrix.
+    pub fn det(&self) -> T {


What would happen if we tried to cholesky decompose an empty matrix here?

From what I can see - we would get an empty matrix within self and this function would return T::one(). I don't think that this is desired behaviour?

I see what you mean. Still, what would you have the determinant of an empty matrix A be? We could define it to be 0, but this would imply that it's singular, and since rank(A) == A.rows() == A.cols(), it is also invertible by the invertible matrix theorem, so we have a rather strange contradiction here. I think the best approach is to define the determinant to be any non-zero number, and so we might as well make it 1 for simplicity.

Oh, and due to https://en.wikipedia.org/wiki/Empty_product and the usual definition of determinant, one could argue that det(A) == 1 in any case :-)

Yes, I think this is a good argument.

Do you think it is worth adding a comment to state our choice here? It seems fairly niche but maybe it will save one person a headache in the future...

Yes, I agree! I will do this shortly.

Please could you add the same comment to the FullPivLU that we have just merged?

Andlon · 2017-03-24T12:15:52Z

I think perhaps you may want to delay merging this PR until #160 is in, so that @Prismatica doesn't have to make more tiny adjustments due to potential merge conflicts in decomposition documentation! He's already been forced to make so many small adjustments already, I'd rather I do this for this PR instead (if necessary. Perhaps git merges automatically).

Andlon · 2017-03-24T13:02:12Z

OK, I believe the outstanding issues should have been taken care of now!

EDIT: Note that @Prismatica had already added a note to the docs of FullPivLu::det about empty matrices, so I did not need to change that :-) I did update PartialPivLu::det though.

AtheMathmo · 2017-03-24T13:38:51Z

Awesome, thank you!

Andlon mentioned this pull request Feb 14, 2017

Rewrite forward_substitution and back_substitution #152

Merged

Andlon force-pushed the cholesky branch from 91b62c5 to 25046d6 Compare February 17, 2017 21:32

Andlon added 11 commits February 18, 2017 17:41

Add minimal internal_utils module

e1a1a0e

Basic implementation of Cholesky

5f04dcb

Rewrite lu unpack to use nullify from internal_utils

754de35

Add Cholesky benchmarks

f3807c7

Take raw_slice_mut() in nullify

38b08a3

Rewrite Cholesky in terms of dot()

29163f4

Benchmarks show significant performance gains.

Cholesky: validate diagonal entries

4520bd5

Implement det() for Cholesky

91b0751

Implement transpose_back_substitution

038f129

Implement Cholesky::solve

d243313

Add benchmarks for Cholesky::solve

8b398a1

Andlon force-pushed the cholesky branch from 25046d6 to 8b398a1 Compare February 18, 2017 18:53

Andlon mentioned this pull request Feb 19, 2017

Refactoring rulinalg in terms of a lower level BLAS-like interface #155

Open

AtheMathmo mentioned this pull request Feb 19, 2017

[WIP] Issue #152 use cholesky factorization AtheMathmo/rusty-machine#155

Open

Andlon added 5 commits February 20, 2017 12:04

Rewrite transpose_back_sub.. for better data access pattern

5b47cd5

Enable previously disabled test

5c8a566

Implement Cholesky::inverse()

2c45b7b

Documentation for Cholesky

5de2e50

Fix errors in Cholesky doctests

a4eb96c

AtheMathmo reviewed Feb 22, 2017

View reviewed changes

Andlon added 4 commits February 25, 2017 14:43

Remove stray ///

9e28cdd

Remove confusing comment

c222bb1

Add Cholesky to decomposition module overview

f387056

Deprecate .cholesky()

5af62fa

Add warning to Cholesky::decompose

04c2d77

Rewrite Cholesky solve/inverse to return Result

1916c71

AtheMathmo reviewed Mar 24, 2017

View reviewed changes

AtheMathmo mentioned this pull request Mar 24, 2017

Added LU decomposition with complete pivoting #160

Merged

Andlon added 2 commits March 24, 2017 13:44

Merge branch 'master' into cholesky

8edae5b

Explain determinant of empty matrix to Cholesky/PartialPivLu

8dc150b

Reorder Cholesky and FullPivLu in decomposition docs

9ec0d1c

AtheMathmo merged commit 91110a4 into AtheMathmo:master Mar 24, 2017

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Cholesky decomposition rewrite: initial work #150

Cholesky decomposition rewrite: initial work #150

Andlon commented Feb 11, 2017 •

edited

Loading

AtheMathmo commented Feb 11, 2017

Andlon commented Feb 11, 2017

AtheMathmo commented Feb 17, 2017

Andlon commented Feb 18, 2017

AtheMathmo left a comment

AtheMathmo Feb 22, 2017

AtheMathmo Feb 22, 2017

AtheMathmo Feb 22, 2017

Andlon Feb 22, 2017 •

edited

Loading

AtheMathmo Feb 22, 2017

Andlon Feb 23, 2017

AtheMathmo Feb 22, 2017

Andlon Feb 22, 2017 •

edited

Loading

AtheMathmo Feb 22, 2017

Andlon Feb 23, 2017

AtheMathmo Feb 22, 2017

Andlon commented Feb 25, 2017

AtheMathmo commented Mar 7, 2017

Andlon commented Mar 7, 2017

AtheMathmo commented Mar 7, 2017

Andlon commented Mar 7, 2017

AtheMathmo commented Mar 7, 2017

Andlon commented Mar 23, 2017

AtheMathmo left a comment

AtheMathmo Mar 24, 2017

Andlon Mar 24, 2017

Andlon Mar 24, 2017

AtheMathmo Mar 24, 2017

Andlon Mar 24, 2017

AtheMathmo Mar 24, 2017

Andlon Mar 24, 2017

Andlon commented Mar 24, 2017

Andlon commented Mar 24, 2017 •

edited

Loading

AtheMathmo commented Mar 24, 2017

Cholesky decomposition rewrite: initial work #150

Cholesky decomposition rewrite: initial work #150

Conversation

Andlon commented Feb 11, 2017 • edited Loading

AtheMathmo commented Feb 11, 2017

Andlon commented Feb 11, 2017

AtheMathmo commented Feb 17, 2017

Andlon commented Feb 18, 2017

AtheMathmo left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Andlon Feb 22, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Andlon Feb 22, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Andlon commented Feb 25, 2017

AtheMathmo commented Mar 7, 2017

Andlon commented Mar 7, 2017

AtheMathmo commented Mar 7, 2017

Andlon commented Mar 7, 2017

AtheMathmo commented Mar 7, 2017

Andlon commented Mar 23, 2017

AtheMathmo left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Andlon commented Mar 24, 2017

Andlon commented Mar 24, 2017 • edited Loading

AtheMathmo commented Mar 24, 2017

Andlon commented Feb 11, 2017 •

edited

Loading

Andlon Feb 22, 2017 •

edited

Loading

Andlon Feb 22, 2017 •

edited

Loading

Andlon commented Mar 24, 2017 •

edited

Loading