qr factorization in LAPACK / Orthonormalization

[ganders]

Dear all,

I need to do an orthonormalization for a set of vectors in 3D space. I.e., I have 3 vectors or a 3x3 matrix as input respectively. Since it may happen that sometimes 3 vectors lay almost on a line, I thought Gram - Schmidt is to unstable ( read in http://www.wikipedia.org/wiki/QR_decomposition). Than I read that QR factorization is a different approach for this orthonormalization. I found in this forum a C++ template (thanks a lot to the author) and modified it. Here it ( lapack_forum_example.cpp) is:

#include <iostream>

using namespace std;

extern "C" void dgeqrf_(int *m, int *n, double da[], int *lda, double dtau[], double dwork[], int *ldwork, int *info);

int main(int argv, char** argc)
{
int n = 3;
int m = 3;
int lda = 3;
int info;
double work[3];
double tau[3];
double A[9];

A[0] = 12.0;
A[1] = 6.0;
A[2] = -4.0;
A[3] = -51.0;
A[4] = 167.0;
A[5] = 24.0;
A[6] = 4.0;
A[7] = -68.0;
A[8] = -41.0;

cout << "input A:" << endl;
for(int i=0; i < m; i++)
{
for(int j=0; j < n; j++)
{
cout << A[i*n +j] << endl;
}
cout << ""<< endl;
}

cout << "------------------------------------------------------------\n"<< endl;

dgeqrf_(&m,&n,A,&lda,tau,work,&n,&info);

if (info == 0)
{
cout << "output A:"<< endl;
for(int i=0; i < m; i++)
{
for(int j=0; j < n; j++)
{
cout << A[i*n+j] << endl;
}

cout << "" << endl;
}
}
else cout << "Problem!!! info = " << info << "\n";

cout << "------------------------------------------------------------\n"<< endl;

return 0;
}

Compilation was done successfully via:

g++ -o lapack_forum_example lapack_forum_example.cpp -llapack

Obtained output is:

input A:
12
6
-4

-51
167
24

4
-68
-41

------------------------------------------------------------

output A:
-14
0.230769
-0.153846

-21
-175
0.0555556

14
70
-35

------------------------------------------------------------

Ifaik that "output A" corresponds to the matrix R of the QR decomposition. I cross-checked this output of R on wikipedia. And the values of the lower triangle are all 0 - fine.
Even taking rounding errors for the non - zero values of the output shown above, the difference might be significant/problematic for further matrix multiplications. I also cross-checked with octave - same output as above.
Why are the values in the lower triangle of R not (closer to) 0?
Is there something to improve in this source code?
Alternatively: Is there another / a better LAPACK - based way to obtain Orthonormalization?

Thanks a lot in advance,

Gerd

[Julien Langou]

Hello
with input A:

Code: Select all

input = [
    12   -51     4
     6   167   -68
    -4    24   -41
]


The LAPACK output is:

Code: Select all

output = [
  -14.0000  -21.0000     14.0000
    0.2308 -175.0000   70.0000
   -0.1538    0.0556  -35.0000
]


The way you need to interpret this is that the upper triangle matrix R is
in the upper triangle so you get:

Code: Select all

R = [
  -14.0000   -21.0000    14.0000
    0       -175.0000    70.0000
    0          0        -35.0000
]


The three digits in the lower part are put there on purpose.

Code: Select all

 = [
     *         *          *     
    0.2308     *          *     
   -0.1538    0.0556      *
]


They are useful if you need the Q-factor for example.

[ganders]

Hello Julien,

thanks a lot for your fast and already very helpful answer. I still have a question to your reply:
You mentioned, that these three lower triangle values are useful if I "need the Q-factor for example". Indeed, actually I am looking for the Q matrix instead of the R matrix , since Q is the output of the Gram Schmidt Orthonormalization I want to apply - via QR factorization.
Thus, I would simply do a suitable matrix multiplacation now, like AR* = Q (or similiar) to obtain Q. Are these three lower triangle values useful for that or do you mean something different (another operation)?

Thanks a lot in advance,

Gerd

[Julien Langou]

You can compute Q via A/R (where the / operation is a call to DTRSM) but
this is not a stable way to obtain Q. For ill-conditionned matrix A, you will
get a Q with columns that are not at all orthogonal.

There is a routine to compute the Q factor in LAPACK. After having call
DGEQRF, you need to call DORGQR. It uses the three values that I
previously mentionned.

Julien.