Re: [eigen] QR factorization question

[Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>]

On Wed, 26 Aug 2009, Eric Chu wrote:

> great! that solved my problem.
> i have another issue though: A is a skinny matrix (m > n), and the matrix Q
> is square (m x m).
>
> how do i retrieve the appropriate Q matrix that goes with my (now) upper
> triangular R? (i.e. A = [Q1 Q2] [R1; 0]--or the "economy size decomposition"
> in Matlab)

To get "economy size decomposition", you just have to select the right 
submatrices. The standard QR decomposition writes A = QR where Q is m x m 
and R is m x n, but the last (m-n) rows of R are zero. The economy-size 
decomposition throws away the last (m-n) rows of R and the last (m-n) 
columns of Q, so you get A = QR where Q is m x n and R = n x n.

Here's some code to get you going:

  const int m = 3;
  const int n = 2;
  MatrixXd A(m,n);
  A << 1,6,  3,8,  7,3;
  cout << "A = " << endl << A << endl;
  HouseholderQR<MatrixXd> qrOfA(A);
  MatrixXd Q = qrOfA.matrixQ().block(0,0,m,n);
  cout << "Q = " << endl << Q << endl;
  MatrixXd R = qrOfA.matrixQR().block(0,0,n,n).triangularView<UpperTriangular>();
  cout << "R = " << endl << R << endl;

Output:

  A =
  1 6
  3 8
  7 3
  Q =
  -0.130189 -0.637408
  -0.390567 -0.671066
  -0.911322  0.378658
  R =
  -7.68115 -6.63963
         0   -8.057

I believe that this is rather inefficient, in that Eigen stores the full 
QR decomposition (Q is m x m), while you need only the economy-size one.

Matlab commands for comparison:

  >> A = [1,6;3,8;7,3]
  A =
       1     6
       3     8
       7     3
  >> [Q,R] = qr(A,0)
  Q =
     -0.1302   -0.6374
     -0.3906   -0.6711
     -0.9113    0.3787
  R =
     -7.6811   -6.6396
           0   -8.0570

Hope this helps.

Jitse