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.

Sign up for GitHub

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
Merged

Minor enhancements to custom LDLT implementations #489

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
2b48438
Sparse matrix support in block-diagonal LDLT
peddie Jun 5, 2024
e0ba89c
Condition number estimation for custom LDLTs
peddie Jun 5, 2024
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
15 changes: 12 additions & 3 deletions include/albatross/src/eigen/serializable_ldlt.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,10 @@ namespace Eigen {
// See LDLT.h in Eigen for a detailed description of the decomposition
class SerializableLDLT : public LDLT<MatrixXd, Lower> {
public:
using RealScalar = double;
using Scalar = double;
using MatrixType = MatrixXd;

SerializableLDLT() : LDLT<MatrixXd, Lower>(){};

SerializableLDLT(const MatrixXd &x) : LDLT<MatrixXd, Lower>(x.ldlt()){};
Expand All @@ -43,6 +47,11 @@ class SerializableLDLT : public LDLT<MatrixXd, Lower> {

bool is_initialized() const { return this->m_isInitialized; }

double l1_norm() const {
ALBATROSS_ASSERT(is_initialized() && "Must initialize first!");
return this->m_l1_norm;
}

/*
* Computes the inverse of the square root of the diagonal, D^{-1/2}
*/
Expand Down Expand Up @@ -88,10 +97,10 @@ class SerializableLDLT : public LDLT<MatrixXd, Lower> {
* Computes the product of the square root of A with rhs,
* D^{-1/2} L^-1 P rhs
*/
template <class Rhs>
Eigen::MatrixXd sqrt_solve(const MatrixBase<Rhs> &rhs) const {
template <class Rhs> Eigen::MatrixXd sqrt_solve(const Rhs &rhs) const {
return diagonal_sqrt_inverse() *
this->matrixL().solve(this->transpositionsP() * rhs);
this->matrixL().solve(this->transpositionsP() *
Eigen::MatrixXd(rhs));
peddie marked this conversation as resolved.
Show resolved Hide resolved
}

Eigen::MatrixXd sqrt_transpose() const {
Expand Down
103 changes: 89 additions & 14 deletions include/albatross/src/linalg/block_diagonal.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,30 +22,49 @@ struct BlockDiagonalLDLT;
struct BlockDiagonal;

struct BlockDiagonalLDLT {
using RealScalar = Eigen::SerializableLDLT::RealScalar;
using Scalar = Eigen::SerializableLDLT::Scalar;
using MatrixType = Eigen::SerializableLDLT::MatrixType;
std::vector<Eigen::SerializableLDLT> blocks;

template <class _Scalar, int _Rows, int _Cols>
Eigen::Matrix<_Scalar, _Rows, _Cols>
solve(const Eigen::Matrix<_Scalar, _Rows, _Cols> &rhs) const;
template <typename Derived>
inline Eigen::MatrixXd solve(const Eigen::MatrixBase<Derived> &rhs) const;

template <class _Scalar, int _Rows, int _Cols>
Eigen::Matrix<_Scalar, _Rows, _Cols>
sqrt_solve(const Eigen::Matrix<_Scalar, _Rows, _Cols> &rhs) const;
template <typename Derived>
inline Eigen::MatrixXd
sqrt_solve(const Eigen::MatrixBase<Derived> &rhs) const;

template <class _Scalar, int _Options, typename _StorageIndex>
Eigen::Matrix<_Scalar, Eigen::Dynamic, Eigen::Dynamic>
solve(const Eigen::SparseMatrix<_Scalar, _Options, _StorageIndex> &rhs) const;

template <class _Scalar, int _Options, typename _StorageIndex>
Eigen::Matrix<_Scalar, Eigen::Dynamic, Eigen::Dynamic> sqrt_solve(
const Eigen::SparseMatrix<_Scalar, _Options, _StorageIndex> &rhs) const;

template <class _Scalar, int _Rows, int _Cols>
Eigen::Matrix<_Scalar, _Rows, _Cols>
solve(const Eigen::Matrix<_Scalar, _Rows, _Cols> &rhs,
ThreadPool *pool) const;

template <typename XprType, int BlockRows, int BlockCols, bool InnerPanel>
inline Eigen::MatrixXd
solve(const Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel> &rhs_orig)
const;

template <class _Scalar, int _Rows, int _Cols>
Eigen::Matrix<_Scalar, _Rows, _Cols>
sqrt_solve(const Eigen::Matrix<_Scalar, _Rows, _Cols> &rhs,
ThreadPool *pool) const;

const BlockDiagonalLDLT &adjoint() const;

std::map<size_t, Eigen::Index> block_to_row_map() const;

double log_determinant() const;

double rcond() const;

Eigen::Index rows() const;

Eigen::Index cols() const;
Expand Down Expand Up @@ -74,12 +93,12 @@ struct BlockDiagonal {
/*
* BlockDiagonalLDLT
*/
template <class _Scalar, int _Rows, int _Cols>
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.

BlockDiagonalLDLT::solve(const Eigen::MatrixBase<Derived> &rhs) const {
ALBATROSS_ASSERT(cols() == rhs.rows());
Eigen::Index i = 0;
Eigen::Matrix<_Scalar, _Rows, _Cols> output(rows(), rhs.cols());
Eigen::MatrixXd output(rows(), rhs.cols());
for (const auto &b : blocks) {
const auto rhs_chunk = rhs.block(i, 0, b.rows(), rhs.cols());
output.block(i, 0, b.rows(), rhs.cols()) = b.solve(rhs_chunk);
Expand All @@ -88,12 +107,53 @@ inline Eigen::Matrix<_Scalar, _Rows, _Cols> BlockDiagonalLDLT::solve(
return output;
}

template <class _Scalar, int _Rows, int _Cols>
inline Eigen::Matrix<_Scalar, _Rows, _Cols> BlockDiagonalLDLT::sqrt_solve(
const Eigen::Matrix<_Scalar, _Rows, _Cols> &rhs) const {
template <typename Derived>
inline Eigen::MatrixXd
BlockDiagonalLDLT::sqrt_solve(const Eigen::MatrixBase<Derived> &rhs) const {
ALBATROSS_ASSERT(cols() == rhs.rows());
Eigen::Index i = 0;
Eigen::Matrix<_Scalar, _Rows, _Cols> output(rows(), rhs.cols());
Eigen::MatrixXd output(rows(), rhs.cols());
for (const auto &b : blocks) {
const auto rhs_chunk = rhs.block(i, 0, b.rows(), rhs.cols());
output.block(i, 0, b.rows(), rhs.cols()) = b.sqrt_solve(rhs_chunk);
i += b.rows();
}
return output;
}

template <class _Scalar, int _Options, typename _StorageIndex>
inline Eigen::Matrix<_Scalar, Eigen::Dynamic, Eigen::Dynamic>
BlockDiagonalLDLT::solve(
const Eigen::SparseMatrix<_Scalar, _Options, _StorageIndex> &rhs) const {
ALBATROSS_ASSERT(cols() == rhs.rows());
Eigen::Index i = 0;
Eigen::Matrix<_Scalar, Eigen::Dynamic, Eigen::Dynamic> output(rows(),
rhs.cols());
for (const auto &b : blocks) {
const auto rhs_chunk = rhs.block(i, 0, b.rows(), rhs.cols());
output.block(i, 0, b.rows(), rhs.cols()) = b.solve(rhs_chunk);
i += b.rows();
}
return output;
}

template <typename XprType, int BlockRows, int BlockCols, bool InnerPanel>
inline Eigen::MatrixXd BlockDiagonalLDLT::solve(
const Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel> &rhs_orig)
const {
ALBATROSS_ASSERT(cols() == rhs_orig.rows());
Eigen::MatrixXd rhs{rhs_orig};
return solve(rhs);
}

template <class _Scalar, int _Options, typename _StorageIndex>
inline Eigen::Matrix<_Scalar, Eigen::Dynamic, Eigen::Dynamic>
BlockDiagonalLDLT::sqrt_solve(
const Eigen::SparseMatrix<_Scalar, _Options, _StorageIndex> &rhs) const {
ALBATROSS_ASSERT(cols() == rhs.rows());
Eigen::Index i = 0;
Eigen::Matrix<_Scalar, Eigen::Dynamic, Eigen::Dynamic> output(rows(),
rhs.cols());
for (const auto &b : blocks) {
const auto rhs_chunk = rhs.block(i, 0, b.rows(), rhs.cols());
output.block(i, 0, b.rows(), rhs.cols()) = b.sqrt_solve(rhs_chunk);
Expand Down Expand Up @@ -157,6 +217,17 @@ inline double BlockDiagonalLDLT::log_determinant() const {
return output;
}

inline double BlockDiagonalLDLT::rcond() const {
// L1 induced norm is just the maximum absolute column sum.
// Therefore the L1 induced norm of a block-diagonal matrix is the
// maximum of the L1 induced norms of the individual blocks.
double l1_norm = -INFINITY;
peddie marked this conversation as resolved.
Show resolved Hide resolved
for (const auto &b : blocks) {
l1_norm = std::max(l1_norm, b.l1_norm());
}
return Eigen::internal::rcond_estimate_helper(l1_norm, *this);
}

inline Eigen::Index BlockDiagonalLDLT::rows() const {
Eigen::Index n = 0;
for (const auto &b : blocks) {
Expand All @@ -173,6 +244,10 @@ inline Eigen::Index BlockDiagonalLDLT::cols() const {
return n;
}

inline const BlockDiagonalLDLT &BlockDiagonalLDLT::adjoint() const {
return *this;
}

/*
* Block Diagonal
*/
Expand Down
Loading