這是一個lapack問題，或者可能是我的代碼中的錯誤？

Question

我正在寫c代碼來獲得使用R中的lapack庫獲得矩陣A的QR分解。 我的結果與使用R語言命令得到的結果不同。

這是一個lapack的東西，或者可能是我的代碼中的錯誤？

矩陣（行主要）：

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結果是：

   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 

我的結果：

-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; 
} 

Answer 1

首先，R使用LINPACK常規DQRDC2爲默認。在R中可以使用qr()命令與LAPACK=TRUE選項使用LAPACK例行：

> 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 

但是，你應該知道，在這種情況下，功能qr()使用LAPACK例行DGEQP3。與您正在使用的DGEQRF例程相反，DGEQP3通過列旋轉來計算矩陣的QR分解。

非常正常，你會得到不同的結果，因爲你沒有使用相同的方法。你應該記住QR分解是而不是一個獨特的解決方案。要知道您的QR分解是否正確，您可以簡單地檢查Q和R矩陣是否滿足要求。例如在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 

另請參見?qr。