Is this a lapack problem or maybe a bug in my code?
https://stackoverflow.com/questions/5852598
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianQuestion
I am writing c code to get get the QR decomposition for a matrix A using lapack lib in R. My result is different from the one i get using R language command.
Is this a lapack thing or maybe a bug in my code ?
for matrix (row major) :
1, 19.5, 43.1, 29.1,
1, 24.7, 49.8, 28.2,
1, 30.7, 51.9, 37.0,
1, 29.8, 54.3, 31.1,
1, 19.1, 42.2, 30.9,
1, 25.6, 53.9, 23.7,
1, 31.4, 58.5, 27.6,
1, 27.9, 52.1, 30.6,
1, 22.1, 49.9, 23.2,
1, 25.5, 53.5, 24.8,
1, 31.1, 56.6, 30.0,
1, 30.4, 56.7, 28.3,
1, 18.7, 46.5, 23.0,
1, 19.7, 44.2, 28.6,
1, 14.6, 42.7, 21.3,
1, 29.5, 54.4, 30.1,
1, 27.7, 55.3, 25.7,
1, 30.2, 58.6, 24.6,
1, 22.7, 48.2, 27.1,
1, 25.2, 51.0, 27.5
r result is :
V1 V2 V3 V4
[1,] -4.4721360 -113.16740034 -2.288392e+02 -123.52039508
[2,] 0.2236068 -21.89587861 -2.107945e+01 -7.27753395
[3,] 0.2236068 0.29484219 8.733781e+00 -14.04825478
[4,] 0.2236068 0.25373857 7.566965e-02 -1.55436071
[5,] 0.2236068 -0.23493787 2.999600e-01 0.26555995
[6,] 0.2236068 0.06192165 -3.343037e-01 0.12188660
[7,] 0.2236068 0.32681168 -2.315941e-01 0.40765540
[8,] 0.2236068 0.16696425 1.213823e-01 -0.18580207
[9,] 0.2236068 -0.09792578 -2.561224e-01 -0.16369010
[10,] 0.2236068 0.05735458 -2.993562e-01 0.52588892
[11,] 0.2236068 0.31311047 -4.660317e-02 0.29409317
[12,] 0.2236068 0.28114099 -1.340150e-01 0.16746961
[13,] 0.2236068 -0.25320615 -2.357881e-01 0.26072358
[14,] 0.2236068 -0.20753545 1.360745e-01 0.18135493
[15,] 0.2236068 -0.44045600 -2.456167e-01 0.15393180
[16,] 0.2236068 0.24003736 3.166468e-02 -0.02119950
[17,] 0.2236068 0.15783011 -2.667146e-01 0.32042553
[18,] 0.2236068 0.27200685 -3.732646e-01 0.05926317
[19,] 0.2236068 -0.07052337 3.634501e-03 -0.20518296
[20,] 0.2236068 0.04365337 -4.566657e-02 -0.03457051
my result:
-4.4721, -113.1674, -228.8392, -123.5204,
0.1827, -21.8959, -21.0794, -7.2775,
0.1827, 0.2888, 8.7338, -14.0483,
0.1827, 0.2486, 0.0523, -1.5544,
0.1827, -0.2301, 0.2071, 0.2316,
0.1827, 0.0607, -0.2309, 0.1063,
0.1827, 0.3201, -0.1599, 0.3555,
0.1827, 0.1636, 0.0838, -0.1620,
0.1827, -0.0959, -0.1769, -0.1427,
0.1827, 0.0562, -0.2067, 0.4586,
0.1827, 0.3067, -0.0322, 0.2565,
0.1827, 0.2754, -0.0925, 0.1460,
0.1827, -0.2480, -0.1628, 0.2274,
0.1827, -0.2033, 0.0940, 0.1581,
0.1827, -0.4315, -0.1696, 0.1342,
0.1827, 0.2351, 0.0219, -0.0185,
0.1827, 0.1546, -0.1842, 0.2794,
0.1827, 0.2665, -0.2578, 0.0517,
0.1827, -0.0691, 0.0025, -0.1789,
0.1827, 0.0428, -0.0315, -0.0301,
#include <stdio.h>
#include <R.h>
#include <R_ext/BLAS.h>
#include <R_ext/Lapack.h>
int min(int x, int y) {
if (x <= y)
return x;
else
return y;
}
int max(int x, int y) {
if (x >= y)
return x;
else
return y;
}
void transpose(int* nrow, int* ncol, double* a) {
int i, j, index, k = 0;
double* atransp = malloc(*nrow**ncol * sizeof (double));
//compute transpose
for (i = 0; i<*ncol; i++) {
index = i;
atransp[k] = a[index];
k++;
for (j = 0; j<*nrow - 1; j++) {
index += *ncol;
atransp[k] = a[index];
k++;
}
}
//copy transpose in array a
for (i = 0; i<*nrow**ncol; i++)
a[i] = atransp[i];
//free memory
free(atransp);
}
void getQR(int* rowX, int* colX, double* X, double* Tau) {
const int m = *rowX;
const int n = *colX;
double* a = X;
const int lda = max(1, m);
double* tau = malloc(min(m, n) * sizeof (double));
const int lwork = max(1, n);
double* work = malloc(max(1, lwork) * sizeof (double));
int info;
F77_NAME(dgeqrf)(&m, &n, a, &lda, tau, work, &lwork, &info);
printf("\n dgeqrf() ended with : %d\n",info);
copyTo(min(m, n), tau, Tau);
free(work);
free(tau);
}
int main() {
int rX = 20, cX = 4;
double X[] = {
1, 19.5, 43.1, 29.1,
1, 24.7, 49.8, 28.2,
1, 30.7, 51.9, 37.0,
1, 29.8, 54.3, 31.1,
1, 19.1, 42.2, 30.9,
1, 25.6, 53.9, 23.7,
1, 31.4, 58.5, 27.6,
1, 27.9, 52.1, 30.6,
1, 22.1, 49.9, 23.2,
1, 25.5, 53.5, 24.8,
1, 31.1, 56.6, 30.0,
1, 30.4, 56.7, 28.3,
1, 18.7, 46.5, 23.0,
1, 19.7, 44.2, 28.6,
1, 14.6, 42.7, 21.3,
1, 29.5, 54.4, 30.1,
1, 27.7, 55.3, 25.7,
1, 30.2, 58.6, 24.6,
1, 22.7, 48.2, 27.1,
1, 25.2, 51.0, 27.5
};
//column major
transpose(&rX,&cX,X);
//tau is needed to extract Q later
double* tau = malloc(min(rX, cX) * sizeof (double));
double* QR = malloc(rX*cX*sizeof(double));
copyTo(rX*cX,X,QR);
getQR(&rX, &cX, QR, tau);
//printmat(cX,rX,QR,"QR");
return 0;
}
Solution
First of all, R uses the LINPACK routine DQRDC2 as a default. You can use the qr() command in R with the LAPACK=TRUE option to use a LAPACK routine :
> QR <- qr(Mat,LAPACK=TRUE)
> QR
$qr
V3 V4 V2 V1
[1,] -229.9739116 -123.04447843 -114.61600066 -4.45006998
[2,] 0.1823682 -19.23736803 -1.81750327 -0.25177355
[3,] 0.1900584 0.41052447 -12.08962672 0.36451274
[4,] 0.1988473 0.04298840 0.18492760 0.02485389
[5,] 0.1545369 0.37519893 -0.14576039 0.48992952
[6,] 0.1973825 -0.32149888 -0.01277324 -0.02631498
[7,] 0.2142277 -0.25359559 0.19391630 0.15368409
[8,] 0.1907908 0.07984478 0.13051129 -0.21692785
[9,] 0.1827344 -0.23371181 -0.11707414 -0.19513972
[10,] 0.1959177 -0.25431801 -0.01531797 0.38150533
[11,] 0.2072699 -0.07795264 0.21041668 0.11596815
[12,] 0.2076361 -0.16711575 0.17604756 -0.02208373
[13,] 0.1702836 -0.14766629 -0.23298857 0.30390902
[14,] 0.1618609 0.20180495 -0.14729488 0.33577876
[15,] 0.1563679 -0.12647957 -0.37138938 0.29164143
[16,] 0.1992135 -0.01062557 0.17041958 -0.12582222
[17,] 0.2025093 -0.25954266 0.06531149 0.14217017
[18,] 0.2145939 -0.40877853 0.13742162 -0.20783552
[19,] 0.1765090 0.01244890 -0.06067115 -0.15702073
[20,] 0.1867626 -0.04646282 0.01506042 -0.06657167
However, you should be aware that the function qr() in this case uses the LAPACK routine DGEQP3. Contrary to the DGEQRF routine you're using, DGEQP3 calculates the QR decomposition of a matrix with column pivoting.
So pretty normal you get different results, as you're not using the same methods. You should remember that a QR decomposition is not a unique solution. To know whether your QR decomposition is correct, you can simply check if the Q and R matrices fulfill the requirements. eg in R :
> Q <- qr.Q(QR)
> round( t(Q) %*% Q , 10 )
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 0 1 0 0
[3,] 0 0 1 0
[4,] 0 0 0 1
> all.equal(Q %*% qr.R(QR),Mat)
[1] TRUE
See also ?qr.
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
