矩阵乘法的不准确性- C

Question

我正在用 C 语言做一个 DS 项目。我一直在与矩阵乘法函数作斗争：

void mat_mult(double **mat1, double **mat2, double **res, int n) {
    int i, j, k;
    
    for (i = 0; i < n; i++) {
        
        for (j = 0; j < n; j++) {
            assert(res[i][j] == 0);
            for (k = 0; k < n; k++)
                res[i][j] += mat1[i][k] * mat2[k][j];;
        }
    }
}

我要相乘的矩阵是：

0.00000000, 0.00007291, 0.00000000 
0.00007291, 0.00000000, 0.00000001 
0.00000000, 0.00000001, 0.00000000

和：

117.11258762, 0.00000000, 0.00000000 
0.00000000, 117.10123871, 0.00000000 
0.00000000, 0.00000000, 8087.59220061

结果是：

0.00000000, 0.00853872, 0.00000007 
0.00853790, 0.00000000, 0.00000172 
0.00000467, 0.00011897, 0.00000000

而不是：

0.00000000 0.00853868 0.00000000
0.00853785 0.00000000 0.00000117
0.00000000 0.00008088 0.00000000

因此，误差很小，但可能很重要。有没有人知道如何解决这个问题或在尝试解决这个问题时要寻找什么？

***编辑-生成代码***

int main(){
int i,j;
int dim=3;
int n=3;
double* a;
double** wam_res, **ddg_res, **lnorm_res, **test_res;
double** vecs;
vecs= (double**)calloc(n, sizeof(double*));
assert(vecs!=NULL);

for(i=0; i<n; i++)
{
    a= (double*)calloc(dim,sizeof(double));
    assert(a!=NULL);
    for(j=0; j<dim; j++)
    {
        a[j]= rand() %50;
    }
    vecs[i]= a;
    
}
wam_res= wam(vecs,dim,n);
printf("%s","wam:\n");
print_matrix(wam_res, n);
ddg_res= ddg(wam_res, n);
printf("%s","ddg:\n");
print_sym_mat(ddg_res, n); 

lnorm_res= lnorm(ddg_res, wam_res, n);
printf("%s","lnorm:\n"); }


double** wam(double** vecs, int dim, int n) {
 /*Calculating wam for the given set of vectors- n vectors of dimension dim */
double** res;
int i,j, curr_index=0;
double norm;
res= (double**)calloc(n, sizeof(double*));
for(i=0; i< n; i++)
{
    res[i]= (double*)calloc(n, sizeof(double));

    /*For vector i, calculating wij for all 0<=j<n*/
    for(j=0; j<n; j++)
    {
        if(i!=j)

        {
            
            norm= l2_norm2vec(*(vecs+i),*(vecs+j),dim); /* function for l2 norm */
            
            res[i][j]= exp((-1*norm/2));
            
            curr_index++;
        }
        else
        {
            /* no self loops */
            res[i][j]=0;
        }
    }
}
return res; }

double** ddg(double** wam, int n) {
int i, j; 
double** res;
double d=0; 
res= (double**)calloc(n*n, sizeof(double*)); /* a n*n array with just 0s */
for(i=0; i<n; i++)
{
    res[i]= (double*)calloc(n, sizeof(double));
    for(j=0; j<n; j++)
    {
        d+=wam[i][j];
    }
    res[i][i]= 1/sqrt(d);
    d=0;

}
return res; }

double** lnorm(double** ddg, double** wam,  int n) {
double** res, **tmp;
int i;
res= (double**)calloc(n, sizeof(double));
tmp= (double**)calloc(n, sizeof(double));
for(i=0; i<n; i++)
{
    res[i]= (double*)calloc(n, sizeof(double));
    tmp[i]= (double*)calloc(n, sizeof(double));
}

mat_mult(ddg,wam, res, n);
printf("%s", "step 1 results: \n");
print_sym_mat(res,n); /* errors happened here */
mat_mult(copyMatrix(res,n), ddg,tmp,n);  /*copyMatrix- simple loop to allocate an identical matrix */
res= copyMatrix(tmp,n);
printf("%s", "step 2 results: \n");
print_sym_mat(res,n); /* errors happened here too */
mat_subI(res,n); 
printf("%s", "step 3 results: \n");
print_sym_mat(res,n);
return res; }

调用矩阵乘法的主要函数是lnorm，出现在代码末尾。它接收不同图形表示的两个向量计算

Answer 1

不准确不应该归咎于此：第二个矩阵是对角矩阵，因此矩阵乘法涉及结果矩阵的每个元素的单次乘法。您的代码仅输出 8 位小数，而值具有更多数据。

这里是编译和输出发布结果的代码的修改版本，显示了问题：矩阵中的值四舍五入到小数点后 8 位，这对于某些元素来说是一个显着差异：

#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

double **wam(double **vecs, int dim, int n);
double **ddg(double **wam, int n);
double **lnorm(double **ddg, double **wam, int n);

double l2_norm2vec(double *vec1, double *vec2, int n) {
    double res = 0;
    for (int i = 0; i < n; i++) {
        double d = vec2[i] - vec1[i];
        res += d * d;
    }
    return sqrt(res);
}

#define print_sym_mat print_matrix

double **copyMatrix(double **res, int n) {
    double **mat = malloc(sizeof(*mat) * n);
    for (int i = 0; i < n; i++) {
        mat[i] = malloc(sizeof(*mat[i]) * n);
        for (int j = 0; j < n; j++) {
            mat[i][j] = res[i][j];
        }
    }
    return mat;
}

double **mat_subI(double **mat, int n) {
    for (int i = 0; i < n; i++) {
        mat[i][i] -= 1;
    }
    return mat;
}

void mat_mult(double **mat1, double **mat2, double **res, int n) {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            res[i][j] = 0;
            for (int k = 0; k < n; k++)
                res[i][j] += mat1[i][k] * mat2[k][j];
        }
    }
}

void print_matrix(double **mat, int n) {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            printf("%19.15f,", mat[i][j]);
        }
        printf("\n");
    }
    printf("\n");
}

int main() {
    int i, j;
    int dim = 3;
    int n = 3;
    double *a;
    double **wam_res, **ddg_res, **lnorm_res;
    double **vecs;
    vecs = (double**)calloc(n, sizeof(double*));
    assert(vecs != NULL);

    for (i = 0; i < n; i++) {
        a = (double*)calloc(dim, sizeof(double));
        assert(a != NULL);
        for (j = 0; j < dim; j++) {
            a[j] = rand() % 50;
        }
        vecs[i] = a;
    }
    wam_res = wam(vecs, dim, n);
    printf("%s", "wam:\n");
    print_matrix(wam_res, n);
    ddg_res = ddg(wam_res, n);
    printf("%s", "ddg:\n");
    print_sym_mat(ddg_res, n);

    lnorm_res = lnorm(ddg_res, wam_res, n);
    printf("%s", "lnorm:\n");
    print_matrix(lnorm_res, n);
    return 0;
}

double **wam(double **vecs, int dim, int n) {
    /* Calculating wam for the given set of vectors- n vectors of dimension dim */
    double **res;
    int i, j, curr_index = 0;
    double norm;
    res = (double **)calloc(n, sizeof(double*));
    for (i = 0; i < n; i++) {
        res[i] = (double*)calloc(n, sizeof(double));

        /*For vector i, calculating wij for all 0<=j<n*/
        for (j = 0; j < n; j++) {
            if (i != j) {
                norm = l2_norm2vec(*(vecs+i), *(vecs+j), dim); /* function for l2 norm */
                res[i][j] = exp((-1 * norm / 2));
                curr_index++;
            } else {
                /* no self loops */
                res[i][j] = 0;
            }
        }
    }
    return res;
}

double **ddg(double **wam, int n) {
    int i, j;
    double **res;
    double d = 0;
    res = (double**)calloc(n*n, sizeof(double*)); /* a n*n array with just 0s */
    for (i = 0; i < n; i++) {
        res[i] = (double*)calloc(n, sizeof(double));
        for (j = 0; j < n; j++)  {
            d += wam[i][j];
        }
        res[i][i] = 1 / sqrt(d);
        d = 0;
    }
    return res;
}

double **lnorm(double **ddg, double **wam, int n) {
    double **res, **tmp;
    int i;
    res = (double**)calloc(n, sizeof(double));
    tmp = (double**)calloc(n, sizeof(double));
    for (i = 0; i < n; i++) {
        res[i] = (double*)calloc(n, sizeof(double));
        tmp[i] = (double*)calloc(n, sizeof(double));
    }

    mat_mult(ddg,wam, res, n);
    printf("%s", "step 1 results: \n");
    print_sym_mat(res, n); /* errors happened here */
    mat_mult(copyMatrix(res, n), ddg, tmp, n);  /*copyMatrix- simple loop to allocate an identical matrix */
    res = copyMatrix(tmp, n);
    printf("%s", "step 2 results: \n");
    print_sym_mat(res, n); /* errors happened here too */
    mat_subI(res, n);
    printf("%s", "step 3 results: \n");
    print_sym_mat(res, n);
    return res;
}

带 8 位小数的输出：

wam:
  0.00000000,  0.00007291,  0.00000000,
  0.00007291,  0.00000000,  0.00000001,
  0.00000000,  0.00000001,  0.00000000,

ddg:
117.11258762,  0.00000000,  0.00000000,
  0.00000000,117.10123871,  0.00000000,
  0.00000000,  0.00000000,8087.59220061,

step 1 results:
  0.00000000,  0.00853872,  0.00000007,
  0.00853790,  0.00000000,  0.00000172,
  0.00000467,  0.00011897,  0.00000000,

step 2 results:
  0.00000000,  0.99989517,  0.00054713,
  0.99989517,  0.00000000,  0.01393205,
  0.00054713,  0.01393205,  0.00000000,

step 3 results:
 -1.00000000,  0.99989517,  0.00054713,
  0.99989517, -1.00000000,  0.01393205,
  0.00054713,  0.01393205, -1.00000000,

lnorm:
 -1.00000000,  0.99989517,  0.00054713,
  0.99989517, -1.00000000,  0.01393205,
  0.00054713,  0.01393205, -1.00000000,

更多的地方，看起来wam中的系数与0.00000001中的系数有很大不同：

wam:
  0.000000000000000,  0.000072910387337,  0.000000000577653,
  0.000072910387337,  0.000000000000000,  0.000000014710728,
  0.000000000577653,  0.000000014710728,  0.000000000000000,

ddg:
117.112587619035523,  0.000000000000000,  0.000000000000000,
  0.000000000000000,117.101238706028283,  0.000000000000000,
  0.000000000000000,  0.000000000000000,8087.592200608456551,

step 1 results:
  0.000000000000000,  0.008538724125301,  0.000000067650464,
  0.008537896671658,  0.000000000000000,  0.000001722644500,
  0.000004671823707,  0.000118974371006,  0.000000000000000,

step 2 results:
  0.000000000000000,  0.999895172041842,  0.000547129363250,
  0.999895172041842,  0.000000000000000,  0.013932046219111,
  0.000547129363250,  0.013932046219111,  0.000000000000000,

step 3 results:
 -1.000000000000000,  0.999895172041842,  0.000547129363250,
  0.999895172041842, -1.000000000000000,  0.013932046219111,
  0.000547129363250,  0.013932046219111, -1.000000000000000,

lnorm:
 -1.000000000000000,  0.999895172041842,  0.000547129363250,
  0.999895172041842, -1.000000000000000,  0.013932046219111,
  0.000547129363250,  0.013932046219111, -1.000000000000000,

矩阵乘法的不准确性- C

Inaccuracies in matrix multiplication- C

c

matrix

double