c - 在 C 中使用 dgsev 时获取不正确的值

Question

我正在使用 C 中的 dgesv 库进行线性回归(带截距)并返回回归系数。我知道，给定下面代码中的 X 和 Y 值，我应该得到回归系数:

The regression coefficients: 26.851352 0.240842 0.119026

但是，当我运行代码时，我得到的回归系数为:

The regression coefficients: 10680696.608794, 112946.210120, -188705.310886

我相信我的矩阵乘法正在工作(调试时)。是否有什么遗漏导致我的回归系数值如此之大且不正确？

#include <stdio.h>
#define N 16
#define P 2
void dgesv_(int *M, int *NRHS, double *A, int *LDA, int *IPIV, double *B, int *LDB, int *INFO);

int main(){
/* longley dataset from R: Employed (Y) GNP.deflator and Population (X) */
double Y[N] = {60.323,61.122,60.171,61.187,63.221,63.639,64.989,63.761,66.019,67.857,68.169,66.513,68.655,69.564,69.331,70.551};
double X[N*(P+1)] = {1.0,83.0,107.608,1.0,88.5,108.632,1.0,88.2,109.773,1.0,89.5,110.929,1.0,96.2,112.075,1.0,98.1,113.27,1.0,99.0,115.094,1.0,100.0,116.219,1.0,101.2,117.388,1.0,104.6,118.734,1.0,108.4,120.445,1.0,110.8,121.95,1.0,112.6,123.366,1.0,114.2,125.368,1.0,115.7,127.852,1.0,116.9,130.081};

int i;
int j;
int k;

int n = P+1;
int nrhs = 1;
double XTX[(P+1)*(P+1)];
int lda = P+1;
int ipiv[P+1];
double XTY[P+1];
int ldb = P+1;
int info;

/* Matrix multiplication */
for (i=0; i<(P+1); i++){
    for (j=0; j<(P+1); j++){
        XTX[i*(P+1)+j]=0;
        for (k=0; k<N; k++){
            XTX[i*(P+1)+j] += X[(P+1)*k+i] * X[(P+1)*k+j];
        }
    }
}

/* Matrix multiplication */
for (i=0; i<(P+1); i++){
    XTY[i] = 0;
    for (j=0; j<N; j++){
        XTY[i] += X[(P+1)*j+1] * Y[j];
    }
}

dgesv_(&n, &nrhs, XTX, &lda, ipiv, XTY, &ldb, &info);

/* From docoumentation, the regression coefficients are returned in XTY */
printf("The regression coefficients: %f, %f, %f\n", XTY[0], XTY[1], XTY[2]);

return 0;
}

Answer 1

第二个矩阵乘法似乎不正确:

    XTY[i] += X[(P+1)*j+1] * Y[j];

为什么要添加1？不应该是:

    XTY[i] += X[(P+1)*j] * Y[j];