Re: [eigen] QR factorization question |

[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]

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