Problems with dgeqrf and dorgqr

[Paulo_G_1]

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

[mecej4]

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]

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]

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]

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.