Hi all,

I'm trying to replace Octave in a custom simulation by Eigen. In total, it was very smooth, but I got totally stuck with the QR decomposition.

For the determination of some residual, I need both Q and R matrices, but I did not make it by matrixQ() from FullPivHouseholderQR.
Also, I was not able to determine the data format returned by matrixQR(). The only hint is that it is compatible to Lapack (even though I found a discussion, whether this is really the case...). But also getting the format Lapack uses, was out of reach for me (sorry for my incompetence...).

So, my question is, can anyone tell me, how to extract Q from the return of matrixQR?

Here is my test program, with which i tried to figure out, what is happening:

Code: Select all

#include <iostream>
#include <eigen3/Eigen/Dense>
using namespace Eigen;
int main()
{
    MatrixXd A = MatrixXd::Random(3,3);

    FullPivHouseholderQR<MatrixXd> qr = A.fullPivHouseholderQr();
    MatrixXd R = qr.matrixQR().triangularView<Upper>();
    MatrixXd Q = qr.matrixQ();

    // If I am not totally wrong, I should get Q also by
    MatrixXd q = A*R.inverse();
    // but inverting R is not acceptable in the simulations

    std::cout << qr.matrixQR() << std::endl << std::endl;
    std::cout << R << std::endl << std::endl;
    std::cout << Q << std::endl << std::endl;
    std::cout << q << std::endl << std::endl;

    std::cout << A << std::endl << std::endl;
 
    std::cout << Q*R << std::endl << std::endl;
    std::cout << q*R << std::endl << std::endl;
 
      return 0;
}



The output is

Code: Select all

-1.18321 0.0932189 -0.434249
 0.297473 -0.905219  0.606159
-0.301468  0.283037  0.194758

 -1.18321 0.0932189 -0.434249
        0 -0.905219  0.606159
        0         0  0.194758

-0.504459 -0.803563 -0.315923
-0.695816  0.161697  0.699782
 0.511235 -0.572836  0.640701

-0.575026 -0.718592 -0.737725
 0.178527 -0.891113   5.92602
-0.478528  0.618954  -5.27544

 0.680375   0.59688 -0.329554
-0.211234  0.823295  0.536459
 0.566198 -0.604897 -0.444451

  0.59688  0.680375 -0.329554
 0.823295 -0.211234  0.536459
-0.604897  0.566198 -0.444451

 0.680375   0.59688 -0.329554
-0.211234  0.823295  0.536459
 0.566198 -0.604897 -0.444451



Any help?!?!?!
Thanks a lot!
Hecke

matrixQ() returns a proxy object representing the matrix Q as a combination of a permutation and q householder sequence, but you can simply do:

MatrixXd Q = qr.matrixQ();

to get a standard matrix.BTW, do you really need full pivoting?

hi Ggael,

thanks a lot for your reply.

I tried

Code: Select all

MatrixXd Q = qr.matrixQ();

as you can see in my example code.
But unfortunately, the result does not seem to be the Matrix that fits into A = Q*R which should be the same as Q = A*R.inverse() which is matrix q (which is not equal to the output of MatrixXd Q = qr.matrixQ(); ) in my example. What did I get wrong?

Full pivoting... hm. I am dealing with a quite complex flow solver which was written by someone else, and I just jumped into it. So I cannot really say. Maybe we can drop the accuracy of the full pivoting and go to the simple HouseholderQr(), but there is no function qr.matrixQ() available at all, right?

thanks again
Hecke

OH, as it was Friday, I did not see that in my example the second to last matrix (Q*R) is A with the first and second column switched...

So, maybe Q uses a different sorting of the columns according to the pivoting, correct? How can I adjust for this?

thanks
Hecke

the factorisation is A * P = Q * R;

you can get P via the colsPermutation() function. You can apply it to a matrix as usual, and apply the inverse using the ... .inverse() method!

Hi Ggael,

thanks again for your help. In the meantime I dug a bit deeper in the code, and saw that the QR decomposition was primarily used to solve a linear equation system by hand, so I can use the solve() function!

Code: Select all

 y = A.householderQr().solve(b);

This seems to do the job, but now I have a different problem, presumably not with Eigen...

thanks a lot and keep up the good work!
Hecke