[eigen] request for comments / help

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]


I'm looking for a helpful soul who has a look at my first non-trivial Eigen function. It implements the QR-Cholesky method to solve the minimum-length least squares problem for possibly rank deficient matrices (or it applies the pseudo-inverse of a matrix to a vector)

I'm not too experienced in template programming.
My code is broken since it wouldn't work for matrices with a fixed number of columns and/or rows.

I don't understand the reasoning for this code snippet of ColPivHouseholder
namespace internal {

template<typename _MatrixType, typename Rhs>
struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
  : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs>

Do I need something similar for the return type of my 'solve' method below?

Please comment on errors, missed optimizations and style of programming.

Many thanks,

And here is my first try :

#include <Eigen/QR>

namespace Eigen {
template<typename _MatrixType> class QR_Cholesky : public ColPivHouseholderQR<_MatrixType>
  typedef _MatrixType MatrixType;
  typedef typename ColPivHouseholderQR<MatrixType>::Index Index;

QR_Cholesky() : ColPivHouseholderQR<MatrixType>(), my_isInitialized(false) {} QR_Cholesky(const MatrixType& matrix) : ColPivHouseholderQR<MatrixType>() {
  QR_Cholesky& compute(const MatrixType& matrix);
  MatrixType solve(const MatrixType& rhs) const;

  bool my_isInitialized;

  MatrixType V;
  LLT<MatrixType> CR;

template<typename MatrixType>
QR_Cholesky<MatrixType>& QR_Cholesky<MatrixType>::compute(const MatrixType& matrix) {
  Index q= this->dimensionOfKernel();
  Index Rank= this->rank();
  if ( q > 0 ) {
    V = this->matrixQR().topRightCorner(Rank,q);
this->matrixQR().topLeftCorner(Rank,Rank).template triangularView<Upper>().solveInPlace(V);
    CR.compute(V.transpose()*V + MatrixType::Identity(q,q));
  my_isInitialized= true;
  return *this;

template<typename MatrixType>
MatrixType QR_Cholesky<MatrixType>::solve(const MatrixType& rhs) const {
  eigen_assert(my_isInitialized && "QR_Cholesky is not initialized.");
  Index Rank= this->rank();
  MatrixType Sol(this->cols(),rhs.cols());
MatrixType y = (this->householderQ().inverse() * rhs).topLeftCorner(Rank,rhs.cols()); this->matrixQR().topLeftCorner(Rank,Rank).template triangularView<Upper>().solveInPlace(y);
  Index q = this->dimensionOfKernel();
  if ( q > 0 ) {
    MatrixType z = CR.solve(V.transpose()*y);
    Sol.topLeftCorner(Rank,rhs.cols())= y-V*z;
  } else
    Sol= y;
  return this->colsPermutation()*Sol;

}  // End of namespace Eigen

// T E S T

#include <iostream>
using std::cout;  using std::endl;
using Eigen::MatrixXd;  using Eigen::VectorXd;

int main() {
  MatrixXd A(3,4);
  VectorXd b(3);
  A << 1,2,3,4,

  b << 30,70,110;

  double eps=1E-8;
  Eigen::QR_Cholesky<MatrixXd> QRChol;
//   Eigen::QR_Cholesky<MatrixXd>  QRFac(A);
  VectorXd x = QRChol.solve(b);
  cout << "LSQ Solution: \n" << x << endl;

Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/