Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Paulo_G_1
Beginner
92 Views

Problems with dgeqrf and dorgqr

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

 

0 Kudos
4 Replies
mecej4
Black Belt
92 Views

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.

Paulo_G_1
Beginner
92 Views

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

mecej4
Black Belt
92 Views

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.

Paulo_G_1
Beginner
92 Views

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.


 

Reply