Re: [eigen] request for comments / help

[Helmut Jarausch <jarausch@xxxxxxxxxxxxxxxxxxx>]

On 06/11/2012 04:14:14 PM, Gael Guennebaud wrote:
> Hi,
>
> I guess the problem is on the type of V and thus CR which have to be
> of dynamic size:
>
> enum {
>   MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
> };
>
> typedef  
> Matrix<Scalar,Dynamic,Dynamic,MaxColsAtCompileTime,MaxColsAtCompileTime>
> MatrixVType;
>
> MatrixVType V;
> LLT<MatrixVType> CR;
>
>
> Same for the matrix y...

Thanks Gael,

I've attached the corrected version of QR_Cholesky.
Is it impudent to suggest adding this to Eigen?

I think it's the standard way of solving the minimum length least  
squares problem.

Helmut.

#include <Eigen/QR>

namespace Eigen {
template<typename _MatrixType> class QR_Cholesky : public ColPivHouseholderQR<_MatrixType>
{
public:
  typedef _MatrixType MatrixType;
  typedef typename ColPivHouseholderQR<MatrixType>::Index Index;
  typedef typename ColPivHouseholderQR<MatrixType>::Scalar Scalar;
  enum { 
    Options = MatrixType::Options,
    MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
    MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
  };
  typedef Matrix<Scalar,Dynamic,Dynamic,Options,MaxColsAtCompileTime,MaxColsAtCompileTime> MatrixVType;
  typedef Matrix<Scalar,Dynamic,Dynamic,Options,MaxColsAtCompileTime,MaxRowsAtCompileTime> MatrixRType;

  QR_Cholesky() : ColPivHouseholderQR<MatrixType>(), my_isInitialized(false) {}
  QR_Cholesky(const MatrixType& matrix) : ColPivHouseholderQR<MatrixType>() {
    compute(matrix);
  }
  QR_Cholesky& compute(const MatrixType& matrix);  
  MatrixRType solve(const MatrixType& rhs) const;
  
protected:
  bool my_isInitialized;
  
  MatrixVType V;
  LLT<MatrixVType> 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>
typename QR_Cholesky<MatrixType>::MatrixRType 
  QR_Cholesky<MatrixType>::solve(const MatrixType& rhs) const {
  eigen_assert(my_isInitialized && "QR_Cholesky is not initialized.");
  eigen_assert(this->rows()==rhs.rows());
  Index q =   this->dimensionOfKernel();
  Index Rank= this->rank();
  MatrixRType Sol(this->cols(),rhs.cols());
  MatrixRType y = (this->householderQ().inverse() * rhs).topLeftCorner(Rank,rhs.cols());
  this->matrixQR().topLeftCorner(Rank,Rank).template triangularView<Upper>().solveInPlace(y);
  if ( q > 0 ) {
    MatrixRType 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

#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;
}