- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Hello,

I'm having some unexpected results with the function LAPACKE_dgeqrf. Apparently I'm unable to get the appropriate QR decomposition at some cases, I'm rather obtaining a QR decomposition with some unexpected vector orientations for the orthogonal matrix Q.

Here is a MWE of the problem:

#include <stdio.h> #include <stdlib.h> #include "mkl.h" #define N 2 int main() { double *x = (double *) malloc( sizeof(double) * N * N ); double *tau = (double *) malloc( sizeof(double) * N ); int i, j; /* Pathological example */ x[0] = 4.0, x[1] = 1.0, x[2] = 3.0, x[3] = 1.0; printf("\n INITIAL MATRIX\n\n"); for (i = 0; i < N; i++) { for (j = 0; j < N; j++) { printf(" %3.2lf\t", x[i*N+j]); } printf("\n"); } LAPACKE_dgeqrf ( LAPACK_ROW_MAJOR, N, N, x, N, tau); printf("\n R MATRIX\n\n"); for (i = 0; i < N; i++) { for (j = 0; j < N; j++) { if ( j >= i ){ printf(" %3.2lf\t", x[i*N+j]); }else{ printf(" %3.2lf\t", 0.0); } } printf("\n"); } LAPACKE_dorgqr ( LAPACK_ROW_MAJOR, N, N, N, x, N, tau); printf("\n Q MATRIX\n\n"); for (i = 0; i < N; i++) { for (j = 0; j < N; j++) { printf(" %3.2lf\t", x[i*N+j]); } printf("\n"); } printf("\n"); return 0; }

With this example, the output I get is:

INITIAL MATRIX

4.00 1.00

3.00 1.00

R MATRIX

-5.00 -1.40

0.00 0.20

Q MATRIX

-0.80 -0.60

-0.60 0.80

However, the expected QR decomposition would be:

R MATRIX

5.00 1.40

0.00 0.20

Q MATRIX

0.80 -0.60

0.60 0.80

I have found this problem with other Initial matrices as well.

Thanks in advance,

Paulo

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

There is no problem. Just as (-2) X 3 and 2 X (-3) are both acceptable factorizations of -6, some columns of Q and the corresponding rows of R may have their signs flipped.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Hello mecej4, thanks for the reply,

yeap, I know that the given factorization is acceptable. My point (and I probably should have mentioned that explicitly in the description of the problem) is that it is not the expected factorization obtained typically by the gram-schimidt process. It may seem irrelevant, but in the particular application I'm interested it is very important that the directions of the orthonormalized column vectors of Q are preserved, so as the diagonal elements of R are positive.

So, in other terms, is it possible to force the library to obtain the expected QR decomposition by GS?

Thank you

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

The Lapack routines ?geqrf() do not use Gram-Schmidt or Modified Gram-Schmidt. In fact, after calling ?geqrf() the input matrix has been overwritten by the Householder reflectors that were produced by the factorization.

In other words, Q is not stored in the usual matrix convention, but as a sequence of reflectors from which, if desired, one can calculate the usual representation of Q by calling ?orgqr(). However, in many algorithms one does not want Q explicitly, but wishes to obtain the product of Q and another matrix, using ?ormqr().

If you really wish to obtain Q explicitly and insist on a convention (e.g., all diagonal elements of R should be positive, as you specified), it is easy to flip the signs of the corresponding columns of Q and rows of R to suit.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

I see, I did that already, but I was hoping it could be done by the library, so I wouldn't need the extra loop for flipping the signs

Unfortunately, I do need the Q matrix explicitly and also need its columns to be aligned so all diagonal elements of R are positive.

Anyway, thank you very much for your help.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page