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

Minor enhancements to custom LDLT implementations #489

Merged
merged 2 commits into from
Aug 2, 2024

Conversation

peddie
Copy link
Contributor

@peddie peddie commented Jun 21, 2024

This change is in support of PIC.

@peddie peddie requested a review from akleeman June 21, 2024 02:19
@peddie peddie force-pushed the peddie/ldlt-enhancements branch 2 times, most recently from 4fad1a3 to 6a141bb Compare June 21, 2024 02:30

template <typename Op, typename LhsType, typename RhsType>
inline Eigen::MatrixXd BlockDiagonalLDLT::solve(
const Eigen::CwiseBinaryOp<Op, LhsType, RhsType> &rhs_orig) const {
Copy link
Collaborator

Choose a reason for hiding this comment

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

What happens at the moment if you do:

Eigen::MatrixXd a, b;
block_diag.solve(a.cwiseProduct(b));

I would have guessed it'd find the solve(Matrix<_Scalar, _Rows, _Cols>) call here and do effectively the same thing as these new methods?

What if you change Matrix<...> to MatrixBase<...>?

Copy link
Contributor Author

@peddie peddie Jun 24, 2024

Choose a reason for hiding this comment

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

MatrixBase<Derived> doesn't cover CWiseBinaryOp for sure. But I think it could cover Block and SparseMatrix.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

There are a couple problems here. One is that we have to declare a type for output. So if we use MatrixBase<Derived>, then if you pass it a CWiseBinaryOp, it tries to declare one and assign to it. This is clearly no good as CWiseBinaryOps are never lvalues. But how else can we determine the return type so as to preserve sparsity? We don't care about that for rcond(), but we do care about it (I think) in PIC.

Copy link
Contributor Author

@peddie peddie Jun 24, 2024

Choose a reason for hiding this comment

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

The Right Thing to do here is to make BlockDiagonalLDLT cooperate with the Eigen expression system.

I think I'm going to just make it return MatrixXd for now. But also I was wrong above that MatrixBase covers Block and SparseMatrix.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

This change is up, which allows removing only the overload for CWiseBinaryOp.

const Eigen::CwiseBinaryOp<Op, LhsType, RhsType> &rhs_orig) const {
ALBATROSS_ASSERT(cols() == rhs_orig.rows());
Eigen::MatrixXd rhs{rhs_orig};
return solve(rhs);
Copy link
Collaborator

Choose a reason for hiding this comment

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

what happens if you do; solve(rhs_orig.eval())?

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 can test it, but I don't think forcing an expression type is enough here.

@peddie
Copy link
Contributor Author

peddie commented Jun 21, 2024

Thanks for the careful review. I wrote this code in mad prototyping mode and failed to go back through some of these small patches in detail. I'll rework this.

@peddie peddie force-pushed the peddie/ldlt-enhancements branch 3 times, most recently from 1572b11 to 1b9dca1 Compare June 25, 2024 06:36
@peddie
Copy link
Contributor Author

peddie commented Jun 25, 2024

This should be good to go if you are happy with the rewrite.

@peddie peddie requested a review from akleeman August 1, 2024 23:26
inline Eigen::Matrix<_Scalar, _Rows, _Cols> BlockDiagonalLDLT::solve(
const Eigen::Matrix<_Scalar, _Rows, _Cols> &rhs) const {
template <typename Derived>
inline Eigen::MatrixXd
Copy link
Collaborator

Choose a reason for hiding this comment

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

There must be a way to do this via Eigen types such that if the input is static sized the output will be too? Eigen::Matrix<_Scalar, _Rows, _Cols> would have allowed for that, but Eigen::MatrixXd will force all output to be dynamic doubles. I don't think it's a problem for us, and feel free to proceed without changing this, it just feels like there should be a way, Eigen::DenseBase<Derived>::PlainMatrix or something like that?

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 think the problem is not all expression types involved here derive from even DenseBase.

Copy link
Collaborator

Choose a reason for hiding this comment

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

MatrixBase inherits from DenseBase so it should right? https://gitlab.com/libeigen/eigen/-/blob/master/Eigen/src/Core/MatrixBase.h#L52

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes, but this gets used with things other than a MatrixBase.

Copy link
Collaborator

Choose a reason for hiding this comment

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

I guess I meant that if the input argument is a MatrixBase<Derived> then maybe you could get the PlainMatrix from it?

I doubt this would work but if Eigen::MatrixBase<Derived> works, then Eigen::DenseBase<Derived> should right?

template <typename Derived>
Eigen::DenseBase<Derived>::PlainMatrix BlockDiagonalLDLT::solve(const Eigen::MatrixBase<Derived> &rhs) const {

Anyway everytime I dig into Eigen I start to think I understand how to extend it, then I find something like Solve.h and get thoroughly lost. This isn't important enough to spend the time figuring 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.

Sorry . . . I think the long-term "right way" is to make our solve() members look just like Eigen's ones and only overload _solve_impl() or whatever it is that takes an RhsType reference.

@peddie peddie merged commit 877be3a into master Aug 2, 2024
9 checks passed
@peddie peddie deleted the peddie/ldlt-enhancements branch August 2, 2024 05:51
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants