c - 这是 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。