Hello,
I'm comparing our rather naive LA implementation with Eigen.

We want to solve a linear system like A x = b, with x being a matrix Z and b being the unity matrix I. Afterward we perform some matrix matrix multiplication
Dimensions are n = 9600, m = 20.

With Eigen we do like that:

Code: Select all

  MatrixXd J(n, n);
  MatrixXd V(n, m);
  MatrixXd W(n, m);
  // fill with values
  MatrixXd Z = V.householderQr().solve(Eigen::MatrixXd::Identity(n, n));
  MatrixXd Jn = J + (W - J * V) * Z;


The householderQr took 3600ms.

The multiplication took 1970ms which is much faster then with our own implementation which 5040ms took.

But the QR decomposition took just 50m in hour implementation. We used another methods, but I wonder how can I optimize Eigen further, especially by exploiting the unity matrix at th rhs.

That's how we did it the QR:

QR:

Code: Select all

  modifiedGramSchmidt(V, Q, R);

  for (int i = 0; i < Q.rows(); i++) {
    for (int j = 0; j < Q.cols(); j++) {
      iterVec.append( Q(i,j) );  // copy i-ith row (we use Q^T = Q^-1) to iterVec
    }
    backSubstitution(R, iterVec, tmpVec);
    Z.append(tmpVec);
    iterVec.clear();
  }

So we exploit that the right hand side is the unity matrix here and it runs blazingly fast. How can I optimize Eigen towards that? I looked into triangularView but was not able to make it work.

Thanks for any idea!

You should extract the thinQ first:

Code: Select all

HouseholderQR<MatrixXd> qr(V);
MatrixXd thinQ = MatrixXd::Identity(n, m); 
thinQ.applyOnTheLeft(qr.householderQ().transpose());


and then you can either compute Z as:

Code: Select all

MatrixXd Z = qr.matrixQR().topLeftCorner(m,m).triangularView<Upper>().solve(thinQ.transpose());


or better not explicitly compute the pseudo inverse of V but rather apply solve to the right:

Code: Select all

MatrixXd tmp = qr.matrixQR().topLeftCorner(m,m).triangularView<Upper>().solve<OnTheRight>(W-J*V);
MatrixXd Jn = J + tmp * thinQ.transpose();