Recovering the full orthogonal Q matrix in QR decomposition

[Mikhail Katliar]

According to Blaze documentation, the qr() function returns matrix Q which is “a general m-by-min(m,n) matrix”. This does not match the definition of QR decomposition, which implies that the Q matrix must be orthogonal (and hence square): “QR decomposition, also known as a QR factorization, is a decomposition of a matrix A into a product A = QR of an orthogonal matrix Q and an upper triangular matrix R.” The same is stated in the documentation of SGEQRF LAPACK function, which is used to implement blaze::qr():

A is REAL array, dimension (LDA,N)
On entry, the M-by-N matrix A.
On exit, the elements on and above the diagonal of the array
contain the min(M,N)-by-N upper trapezoidal matrix R (R is
upper triangular if m >= n); the elements below the diagonal,
with the array TAU, represent the orthogonal matrix Q as a
product of min(m,n) elementary reflectors (see Further
Details).

There should be a way to obtain the full Q matrix from blaze::qr(). The org2r and orgr2 LAPACK functions can recover the full Q matrix from the result of geqrf().

Created from https://bitbucket.org/blaze-lib/blaze/issues/14#comment-52850773

[Mikhail Katliar]

My idea is to add a bool argument to blaze::qr() which specifies whether the full Q matrix or only its first min(m,n) columns should be returned:

blaze::DynamicMatrix<double> A(3, 2);
blaze::DynamicMatrix<double> Q1, R1, Q2, R2;

qr(A, Q1, R1, false);
assert(rows(Q1) == 3 && columns(Q1) == 2);
assert(rows(R1) == 2 && columns(R1) == 2);

qr(A, Q2, R2, true);
assert(rows(Q2) == 3 && columns(Q2) == 3);
assert(rows(R2) == 3 && columns(R2) == 2);
assert(Q2 * trans(Q2) == blaze::IdentityMatrix<double>(3));

[Mikhail Katliar]

Another option is to specify how many columns of Q need to be reconstructed:

blaze::DynamicMatrix<double> A(5, 2);
blaze::DynamicMatrix<double> Q, R, Q3, R3, Q5, R5;

// Old behavior
qr(A, Q, R);
assert(rows(Q1) == 5 && columns(Q1) == 2);
assert(rows(R1) == 2 && columns(R1) == 2);

// Recover 3 columns of Q
qr(A, Q3, R3, 3);
assert(rows(Q3) == 5 && columns(Q3) == 3);
assert(rows(R3) == 3 && columns(R3) == 2);
assert(trans(Q3) * Q3 == blaze::IdentityMatrix<double>(3));

// Recover all 5 columns of Q
qr(A, Q5, R5, 3);
assert(rows(Q5) == 5 && columns(Q5) == 5);
assert(rows(R5) == 5 && columns(R2) == 2);
assert(trans(Q5) * Q5 == blaze::IdentityMatrix<double>(5));

[Mikhail Katliar]

I have implemented a draft version of the qr() function which can reconstruct a specified number of columns of Q. Here is the pull request: https://bitbucket.org/blaze-lib/blaze/pull-requests/33.

[Klaus Iglberger]

Hi Mikhail!

Thanks a lot for creating the issue and already providing a pull request. I have added a couple of comments for the first draft.

Best regards,

Klaus!

[Klaus Iglberger]

Hi Mikhail!

With the last push we have added overloads for the LAPACK functions org2r(), orgr2(), org2l(), orgl2(), ung2r(), ungr2(), ung2l(), and ungl2(). Hope fully these help to finish pull request 33.

Best regards,

Klaus!