[eigen] Dense Matrix Decompositions/Solvers: Inherit from SolverBase

[Patrick Peltzer]

Good day,

as this is my first entry to this mailing list, I would like to shortly
introduce myself. My name is Patrick Peltzer and I'm a Computational
Engineering Science student at RWTH Aachen, Germany. I'm currently
working on my Bachelor thesis with the title "Efficient Adjoints of the
Linear Algebra Library
Eigen using dco/c++", in which I'm creating a fork of Eigen to optimize
its usage with the algorithmic differentiation software dco/c++. Amongst
other things, I'm implementing symbolic solver versions of the Eigen
build-in solvers. Symbolic solvers allow to reuse the already computed
matrix decomposition when doing the reverse part of algorithmic
differentiation, saving memory and speed by making the recording of the
decomposition calculation and the recalculation obsolete. If you are
interested, you can find more info here:
https://www.stce.rwth-aachen.de/research/publications/naumann2012dls

For these symbolic solvers, I've created a parent class which inherits
from SolverBase. The symbolic solver variants for FullPivLU and
PartialPivLU inherit from my parent class and they work like a charm.
However, I would like to include the other solvers listed here
https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html, e.g.
HouseholderQR. However, all these other decompositions do not inherit
from SolverBase. My questions is: Is there a reason why?

From the documentation I can see that all SolverBase subclasses must
have the same interface as shown here
https://eigen.tuxfamily.org/dox/classEigen_1_1SolverBase.html. Now, it
appears that only the transpose/adjoint solve implementation are missing
from the other decompositions. Is this the only reason why they do not
inherit from SolverBase? And if so, would it be desirable to implement them?

So far I have only looked at the HouseholderQR class, but changing it to
inherit from SolverBase appeared relatively easy by adding the trait
struct and the _solve_transposed implementation:


template<typename _MatrixType>
template<typename RhsType, typename DstType>
void HouseholderQR<_MatrixType>::_solve_impl_transposed(const RhsType
&rhs, DstType &dst) const
{
  const Index rank = (std::min)(rows(), cols());
  eigen_assert(rhs.rows() == rows());

  dst = householderSequence(m_qr.leftCols(rank),m_hCoeffs.head(rank)) *
m_qr.topLeftCorner(rank, rank).template
triangularView<Upper>().transpose().solve(rhs.topRows(rank));
}


I have only tested the above implementation for quadratic, real
matrices, so it might still need some changes. But it appears to me that
inheriting from SolverBase didn't break any other feature, as the qr
test still succeeds. So, can this be seen more like a feature request to
provide all the missing _solve_transposed implementations, or is there
anything else I'm missing out?


Best regards,

Patrick Peltzer

[Gael Guennebaud]

Hi Patrick,

sorry for the rather late reply.

On Tue, Feb 20, 2018 at 1:55 PM, Patrick Peltzer <

Post by Patrick Peltzer
as this is my first entry to this mailing list, I would like to shortly
introduce myself. My name is Patrick Peltzer and I'm a Computational
Engineering Science student at RWTH Aachen, Germany. I'm currently working
on my Bachelor thesis with the title "Efficient Adjoints of the Linear
Algebra Library
Eigen using dco/c++", in which I'm creating a fork of Eigen to optimize
its usage with the algorithmic differentiation software dco/c++. Amongst
other things, I'm implementing symbolic solver versions of the Eigen
build-in solvers. Symbolic solvers allow to reuse the already computed
matrix decomposition when doing the reverse part of algorithmic
differentiation, saving memory and speed by making the recording of the
decomposition calculation and the recalculation obsolete. If you are
interested, you can find more info here: https://www.stce.rwth-aachen.d
e/research/publications/naumann2012dls

This looks very interesting. Out of curiosity, have you considered
leveraging matrix differential rules, such as identifying expressions like
x = A.lu().solve(b) and instead of going through the whole factorization to
cary out the derivatives, simply use the rules for differentiating the
inverse of a matrix...??

Post by Patrick Peltzer
For these symbolic solvers, I've created a parent class which inherits
from SolverBase. The symbolic solver variants for FullPivLU and
PartialPivLU inherit from my parent class and they work like a charm.
However, I would like to include the other solvers listed here
https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html, e.g.
HouseholderQR. However, all these other decompositions do not inherit from
SolverBase. My questions is: Is there a reason why?
From the documentation I can see that all SolverBase subclasses must have
the same interface as shown here https://eigen.tuxfamily.org/do
x/classEigen_1_1SolverBase.html. Now, it appears that only the
transpose/adjoint solve implementation are missing from the other
decompositions. Is this the only reason why they do not inherit from
SolverBase? And if so, would it be desirable to implement them?

Yes, that's mostly a lack of time, and having all of them inherit
SolverBase is of course highly desirable.

Post by Patrick Peltzer
So far I have only looked at the HouseholderQR class, but changing it to
inherit from SolverBase appeared relatively easy by adding the trait struct
template<typename _MatrixType>
template<typename RhsType, typename DstType>
void HouseholderQR<_MatrixType>::_solve_impl_transposed(const RhsType
&rhs, DstType &dst) const
{
const Index rank = (std::min)(rows(), cols());
eigen_assert(rhs.rows() == rows());
dst = householderSequence(m_qr.leftCols(rank),m_hCoeffs.head(rank)) *
m_qr.topLeftCorner(rank, rank).template triangularView<Upper>().transp
ose().solve(rhs.topRows(rank));
}
I have only tested the above implementation for quadratic, real matrices,
so it might still need some changes. But it appears to me that inheriting
from SolverBase didn't break any other feature, as the qr test still
succeeds. So, can this be seen more like a feature request to provide all
the missing _solve_transposed implementations, or is there anything else
I'm missing out?

extending the respective unit tests to cover transpose/adjoint solving is
all what's missing to accept a pull-request.

gael

Post by Patrick Peltzer
Best regards,
Patrick Peltzer

[Patrick Peltzer]

Hi Gael,

thank you for your answer.

Post by Gael Guennebaud
This looks very interesting. Out of curiosity, have you considered
leveraging matrix differential rules, such as identifying expressions
like x = A.lu().solve(b) and instead of going through the whole
factorization to cary out the derivatives, simply use the rules for
differentiating the inverse of a matrix...??

This is exactly what the symbolic solvers are doing. Some words on how
the direct linear solvers (LU, QR, ...) can be differentiated:
The ordinary adjoint algorithmic differentiation approach would be to
save all necessary data in the forward run ( termed "taping") on an
external data structure (the "tape") so it can be accessed in the
reverse run when the adjoints are propagated. The forward run includes
the taping of the whole decomposition process, occupying memory in the
order of O(n^3) and accordingly requires the reverse passing of the
computation, causing an additional computational cost in the order of
O(n^3).
By exploiting the mathematical insight, it is possible to propagate the
adjoints in the reverse run without going through the decomposition
computation, given the finished decomposition was saved in the forward
run before. This reduces the computational complexity to O(n^2) in the
reverse section and the memory cost to O(n^2) as well.

I've implemented such symbolic variants for the FullPivLU, PartialPivLU
and HouseholderQR solvers so far. I'm currently collecting data and
things appear to look very promising (especially the additional memory
cost has been greatly reduced). Right now, the symbolic solvers need to
be explicitly instantiated, but I will try to find a nice way to
specialize the existing solvers for the dco/c++ data types in order to
automatically grant support for the hardcoded solver calls inside the
Eigen internals. Besides including the other dense solvers, other
symbolic variants of linear algebra operations are on my agenda for the
future.

Post by Gael Guennebaud
Yes, that's mostly a lack of time, and having all of them inherit
SolverBase is of course highly desirable.
...
extending the respective unit tests to cover transpose/adjoint solving
is all what's missing to accept a pull-request.
gael

Alright, great to know! You can also expect the other solvers to follow
in the future as well.

Best regards,
Patrick