Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Updates LLT to inplace decomposition per eigen 3.3 doc #549

Merged
merged 2 commits into from
May 7, 2017
Merged
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Next Next commit
Updates LLT to inplace decomposition per eigen 3.3 doc
  • Loading branch information
rtrangucci committed May 4, 2017
commit 5cf225f0e3a50195136e4ac6348e326d78ee27be
6 changes: 1 addition & 5 deletions stan/math/rev/mat/fun/cholesky_decompose.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -233,8 +233,6 @@ namespace stan {
* Internally calls llt rather than using cholesky_decompose in order to
* use selfadjointView<Lower> optimization.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this doc broken now? You've removed the selfadjointView<Lower> bit of the code.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed the doc. Thanks for the review.

*
* TODO(rtrangucci): Use Eigen 3.3 inplace Cholesky when possible
*
* Note chainable stack varis are created below in Matrix<var, -1, -1>
*
* @param A Matrix
Expand All @@ -246,10 +244,8 @@ namespace stan {
check_symmetric("cholesky_decompose", "A", A);

Eigen::Matrix<double, -1, -1> L_A(value_of_rec(A));
Eigen::LLT<Eigen::MatrixXd> L_factor
= L_A.selfadjointView<Eigen::Lower>().llt();
Eigen::LLT<Eigen::Ref<Eigen::MatrixXd> > L_factor(L_A);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a reason we're no longer using a selfadjointView of L_A? Have the input requirements changed now that we use all of it?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good question! I'm not sure if I can do a Ref to a triangular view of a matrix. I'll see if I can figure it out.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I used the selfadjointViewEigen::Lower() to cut down on the copy when L_factor was instantiated, but I'm fairly certain that LLT defaults to a lower triangular view. In any case, the inplace decomposition obviates the need for any selfadjointView() instantiation here. I've added an explicit Eigen::Lower template argument, but that's the default for the class.

check_pos_definite("cholesky_decompose", "m", L_factor);
L_A = L_factor.matrixL();

// Memory allocated in arena.
// cholesky_scalar gradient faster for small matrices compared to
Expand Down