I use Lapacke. I am trying to do QR decomposition in C for complex data. For this I write the function (based on Haatschii code How to get the Q from the QR factorization output?):
// Q - input: matrix that we expand / output: Q matrix
// R - output: R matrix
// rows - input: number of rows of Q
// columns - input: number of columns of Q
// rows >= columns condition is always met
void QR(lapack_complex_double * Q, lapack_complex_double * R, size_t rows, size_t columns){
size_t i;
lapack_complex_double* tau = malloc(columns*sizeof(lapack_complex_double));
LAPACKE_zgeqrf(LAPACK_ROW_MAJOR, (int) rows, (int) columns, Q, (int) columns, tau); // returns the Q, R in a packed format
// Copy the upper triangular Matrix R (columns x columns).
for(i = 0; i < columns; ++i)
memcpy(R+i*columns+i, Q+i*columns+i, (columns-i)*sizeof(lapack_complex_double));
LAPACKE_zungqr(LAPACK_ROW_MAJOR, (int) rows, (int) columns, (int) columns, Q, (int) columns, tau); // returns the Q
free(tau);
}
It is worth noting that the author of Alireza also had problems with the function znugqr, but however he switched to the function zunmqr and it seems that happiness has come (LAPACK QR factorization). I believe that my problem is also related to LAPACKE_zungqr since the matrix R is the same as other methods and therefore LAPACKE_zgeqrf works successfully.
But in the end, comparing the similar result (QR decompostion) with Mathematica (QRDecomposition function) and Python (numpy.linalg.qr function), I see that the matrix Q is different, while the matrix R is the same.
input matrix, for simplicity 5×5:
1 + 3j 6 + 8j 11 + 13j 16 + 18j 21 + 23j
2 + 4j 7 + 9j 12 + 14j 17 + 19j 22 + 24j
3 + 5j 8 + 10j 13 + 15j 18 + 20j 23 + 25j
4 + 6j 9 + 11j 14 + 16j 19 + 21j 24 + 26j
5 + 7j 10 + 12j 15 + 17j 20 + 22j 25 + 27j
Output Q matrix (from my C code): first 3 columns:
-0.07254 - 0.21764j -0.61558 - 0.41039j 0.519770 - 0.06712j
-0.14509 - 0.29019j -0.35909 - 0.25649j -0.59817 + 0.211099j
-0.21764 - 0.36273j -0.10259 - 0.10259j 0.035755 - 0.18619j
-0.29019 - 0.43528j 0.153896 + 0.051298j -0.35605 + 0.007600j
-0.36273 - 0.50783j 0.410391 + 0.205195j 0.398709 + 0.034623j
last 2 columns:
-0.12316 - 0.06327j -0.11940 + 0.303152j
0.221078 + 0.491045j 0.084589 - 0.02148j
-0.06231 - 0.31146j 0.119483 - 0.80553j
-0.04594 - 0.59711j -0.01512 + 0.462905j
0.010343 + 0.480803j -0.06954 + 0.060958j
Output Q matrix (from Python code):
-0.11670 - 0.06185j -0.13105 + 0.301181j
0.223111 + 0.487988j 0.096454 - 0.02009j
-0.08117 - 0.30874j 0.138515 - 0.80184j
-0.04015 - 0.59906j -0.04217 + 0.459232j
0.014923 + 0.481676j -0.06174 + 0.061519j
(Here I am listing only the last 2 columns of this matrix. The first 3 columns are the same).
I calculated the maximum and average difference across the columns in these matrices. The conclusion is this: the first three columns differ at the level of 10^-11, and the difference between the last two is 10^-3, 10^-2, respectively (the difference is visible to the naked eye).
With an increase in the size of the matrix, an increase in the difference is observed and, as a rule, the first 2-3 columns coincide well.
Maybe someone can help me?
The input matrix A is 5x5 but it has rank 2. The first two columns are linearly independent, while the last three columns of A are linear combinations of the first two.
For a QR factorization, this means that there is not unicity for the last three columns of Q. An QR factorization implementation (e.g., LAPACK, numpy, etc. ) can return any three columns that are (1) mutually orthonormal and (2) in A^(perp), and that is a correct answer. There are many correct answers! The solution is not unique.
If you want to check the Q and R returned by LAPACK, (or any QR factorization implementation for the matter,) you can (1) compute Q'*Q and check that you get the 5x5 identity matrix, (please use BLAS HERK function to do so) and (2) compute Q'*A and check that you have the 5x5 upper triangular matrix R (as returned by QR). In your case you should see that the last 3 rows of R are all zeros, which indicates that the last three columns of A are linear combinations of the first two. To compute Q'*A, you can use BLAS GEMM.
I took your input A and your output Q from LAPACK and I checked for you that Q'*Q is identity and Q'*R is upper triangular with last three rows being fully zeros. So the Q output from LAPACK looks good to me. Yes, this Q is different from what is returned by another implementation and this is entirely possible.
In general we check the quality of a QR factorization by checking that: (1) || A - QR || / || A || is small, (2) || I - Q^T Q || is small, and R is upper triangular. (Take any norm that is easy to compute.)
Checking that two codes return the same output is not a good idea. Since there is not unicity of the output Q and R for two reasons: (1) when A is rank deficient, see the example that you give for example; (2) for any j, any column scaling with a complex number of modulus 1 for the j-th column of Q is possible as long as you scale the j-th row of R accordingly. This changes the Q and R factors but there are still valid Q and R factors.
In addition, assume that you force all diagonal elements of R to be real nonnegative, (by rescaling columns of Q and columns of R appropriately,) and that A is full rank, then you would expect unicity of Q and R. However checking that the pair (Q,R) is close to another one is related to the condition number of computing Q and R for a QR factorization and so you can see pairs that are far apart (forward error), even though they have a good backward error quality.