QR Decomposition results in eigen library differs from Matlab

Question

I am posed with a problem whereby I am trying to decompose a rank deficient matrix into its Q-R components. This is the matrix we are interested in:

1 -2 4 -8 1 -1 1 -1 1 0 0 0 1 1 1 1 1 2 4 8 Let's call it mymatrix.

When I call the qr(mymatrix, 0) function in MATLAB, here is the R-matrix that gets returned to me:

-2.2361 -0.0000 -4.4721 -0.0000 0 3.1623 0.0000 10.7517 0 0 3.7417 0.0000 0 0 0 3.7947

Now, in eigen c++, I tried to do the same but, alas, I am overwhelmed with an altogether different matrix:

11.4018 -8.88178e-16 -2.22045e-16 2.982 0 -5.83095 -1.71499 3.33067e-16 0 0 -1.43486 9.41987e-17 0 0 0 -1.05247

Here's the function that returns my QR matrix in eigen:

  #include <Eigen/LU>
  #include <Eigen/Dense>
  #include <Eigen/QR>
  #include <cmath>

  using namespace Eigen;
  using namespace std;

  int main(){
  //inter is initialized as mymatrix
  ColPivHouseholderQR<MatrixXd> qr(inter.rows(), inter.cols());
  qr.compute(mymatrix);

  FullPivLU<MatrixXd>lu_decomp(inter);
  int Rank = lu_decomp.rank() ;

  cout <<"\nInter's Rank is: " << Rank << endl;

  if (Rank == inter.rows() + 1 || Rank == inter.cols() + 1 )
  {
  MatrixXd R = qr.matrixR().template triangularView<Upper>();            //retrieve the R - Matrix
  cout << "\nR Matrix: \n" <<  R << endl;
  }
 else          //For rank deficient matrices
 { 
   FullPivLU<MatrixXd>lu_decomp(inter);
   MatrixXd R = qr.matrixR().topLeftCorner(lu_decomp.rank(),        lu_decomp.rank()).template triangularView<Upper>();
   cout <<"\n R with rank deficiency: \n" << R << endl; 
  }
   return 0;
  }

What could I possibly be missing?

Answer 1

The documentation is rather sparse for Eigen's QR documentation. You could retrieve the Q matrix from the householder QR module as outlined here. After retrieving the QR matrix, extracting the R matrix is simply a matter of extracting the upper triangular matrix of matrixQR.

Suppose we have a full rank matrix, we could retrieve the R matrix as follows:

  /*Construct the QR factorization from a given matrix (mymatrix).*/
  HouseholderQR<MatrixXd> qr(mymatrix);
  qr.compute(mymatrix);

  MatrixXd R = qr.matrixQR().template  triangularView<Upper>();

If the original matrix is rank deficient, this would do:

FullPivLU<MatrixXd>lu_decomp(mymatrix);
int Rank = lu_decomp.rank();                 //retrieve rank of matrix
MatrixXd R = qr.matrixQR().topLeftCorner(Rank, Rank).template triangularView<Upper>();

I hope this would help someone down the lane.

Answer 2

Note that the different underlying matrix packages (LAPACK and Eigen) most likely implement a different type of QR algorithm.

Most notable differences are the choice of the pivoting strategy and transformation type (Housholder transformation / Givens rotation / Gram schmidt).

Answer 3

Although this question is a bit old, if I may add something here. The ColPivHouseholderQR should provide a decomposition of matrix A into matrices P, Q and R such that AP = QR, were P is a permutation matrix. (from Eigen documentation). Hence you should not expect to get the same result as Matlab. In order to achieve that you need to use HouseholderQR, as mentioned in in Calorified's answer. Nevertheless, you can get the permutation matrix using

  P = qr.colsPermutation();

then

   A = Q*R*P.inverse();

Using ColPivHouseholderQR the Q and P matrices should be:

      -0.702   -0.686   0.123   -0.0877    0.12 
      -0.0877  -0.171  -0.492    0.702    -0.478 
    Q= 0         0     -0.697     0        0.717 
       0.0877  -0.171  -0.492   -0.702    -0.478 
       0.702   -0.686   0.123    0.0877    0.12 

       0 0 1 0        
    P= 0 0 0 1      
       0 1 0 0 
       1 0 0 0