Hello everyone,

I was porting a code from Matlab to C++ using Eigen and I have found some problems. When I perform a QR decomposition the R matrix I get from Eigen is very different from the R matrix I get from Matlab.

Code: Select all

HouseholderQR<MatrixXd> qr_left(W_left);
qr_left.matrixQR().block(0, 0, cols, cols).triangularView<Upper>();

I have also checked the magnitude of the residual of the reconstructed matrix Q*R = A1 compared to the original matrix A0, by using the piece of code below

Code: Select all

MatrixXd LR = qr_left.matrixQR().triangularView<Upper>();
MatrixXd LQ = qr_left.matrixQ();
MatrixXd A1 = LQ*LR;
double res = (A1 - A0).norm();

The residual I get from Matlab is around 4.9248e-07, while the one obtained from Eigen is something like 6.7e-6.

I thought both, Matlab and Eigen, where using LAPACK for estimating this decomposition. Does anybody know what can be the problem?
You can check this by yourself by computing QR to the matrices of these links matrix in matlab format (spaces) , matrix in Eigen format (commas). The matrix is 831x13

Best regards,
G. Ros

Eigen has it's own QR implementation. It slightly differs from LAPACK in the way Houseolder reflectors are computed thus explaining the sign differences on the diagonal of R. Here I get similar reconstruction errors:

Eigen ((Q*R-A).norm()):
4.79937e-07

Octave (LAPACK):
norm(Q * R - A) = 7.4467e-07
D = Q * R - A
sqrt(trace(D'*D)) = 8.0763e-07 (what Eigen's .norm() computes)

--Edit--

my test program:

Code: Select all

#include <iostream>
#include <Eigen/Dense>
#include <fstream>

using namespace Eigen;

int main(int argc, char **argv)
{
  MatrixXd A(831,13);
  std::fstream file(argv[1]);
  for(int i=0; i<A.rows(); ++i)
    for(int j=0; j<A.cols(); ++j) 
      file >> A(i,j);
   
  HouseholderQR<MatrixXd> qr(A);
  MatrixXd R = qr.matrixQR().triangularView<Upper>();
  MatrixXd Q = qr.householderQ();
  std::cout << R.topRows(13) << "\n";
  MatrixXd D = Q*R-A;
  std::cout << "\n" << (Q*R-A).norm() << "  " << sqrt((D.adjoint()*D).trace()) << "\n";
}


Thank you for your reply ggael. I don't get the same error magnitude with the code you provided. I think the problem might be the version of Eigen3 I use (3.0.5). I am going to check with a more recent version.

Best regards,
G. Ros

I am afraid that not even with the latest stable version of Eigen I get that. Did you use any special flag for compiling the code? Am I missing something? Are you using the hg version?

Best regards,
German Ros.

ggael wrote:Eigen has it's own QR implementation. It slightly differs from LAPACK in the way Houseolder reflectors are computed thus explaining the sign differences on the diagonal of R. Here I get similar reconstruction errors:

Eigen ((Q*R-A).norm()):
4.79937e-07

Octave (LAPACK):
norm(Q * R - A) = 7.4467e-07
D = Q * R - A
sqrt(trace(D'*D)) = 8.0763e-07 (what Eigen's .norm() computes)

--Edit--

my test program:

Code: Select all

#include <iostream>
#include <Eigen/Dense>
#include <fstream>

using namespace Eigen;

int main(int argc, char **argv)
{
  MatrixXd A(831,13);
  std::fstream file(argv[1]);
  for(int i=0; i<A.rows(); ++i)
    for(int j=0; j<A.cols(); ++j) 
      file >> A(i,j);
   
  HouseholderQR<MatrixXd> qr(A);
  MatrixXd R = qr.matrixQR().triangularView<Upper>();
  MatrixXd Q = qr.householderQ();
  std::cout << R.topRows(13) << "\n";
  MatrixXd D = Q*R-A;
  std::cout << "\n" << (Q*R-A).norm() << "  " << sqrt((D.adjoint()*D).trace()) << "\n";
}


I get 1.62442e-06 when disabling vectorization. Then Eigen version does not change anything for me.

If you're running a 32bits system, then compile with -msse2 to get SIMD speedups.

Thank you, it works now! However, there is something I don't understand. Enabling vectorization speeds up the program, ok, but why do I get a better result?

Best regards,
G. Ros

ggael wrote:I get 1.62442e-06 when disabling vectorization. Then Eigen version does not change anything for me.

If you're running a 32bits system, then compile with -msse2 to get SIMD speedups.