[eigen] request for comments / help

[Helmut Jarausch <jarausch@xxxxxxxxxxxxxxxxxxx>]

Hi,

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>
{
  EIGEN_MAKE_SOLVE_HELPERS(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,
Helmut.

And here is my first try :


#include <Eigen/QR>


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

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

protected:
  bool my_isInitialized;

  MatrixType V;
  LLT<MatrixType> CR;
};

template<typename MatrixType>
QR_Cholesky<MatrixType>& QR_Cholesky<MatrixType>::compute(const  
MatrixType& matrix) {
  ColPivHouseholderQR<MatrixType>::compute(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.");
  eigen_assert(this->rows()==rhs.rows());
  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;
    Sol.bottomLeftCorner(q,rhs.cols())=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,
       5,6,7,8,
       9,10,11,12;

  b << 30,70,110;

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