topic The Lapack routines ?geqrf() in Intel® oneAPI Math Kernel Library

Problems with dgeqrf and dorgqr

Paulo_G_1 — Thu, 28 Apr 2016 13:45:33 GMT

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

There is no problem. Just as

mecej4 — Thu, 28 Apr 2016 14:11:20 GMT

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.

yeap, I know that the given

Paulo_G_1 — Thu, 28 Apr 2016 15:21:00 GMT

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

The Lapack routines ?geqrf()

mecej4 — Thu, 28 Apr 2016 17:28:34 GMT

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.

I see, I did that already,

Paulo_G_1 — Thu, 28 Apr 2016 17:37:51 GMT

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.