c - C中循环中的矩阵相乘

Question

我正在尝试使用 C 中的循环来乘以多个矩阵。我在 R 中获得预期答案，但无法在 C 中获得预期答案。我怀疑问题与 += 函数有关，该函数似乎在循环的第一次迭代后将乘积的值加倍。

我对 C 不太熟悉，无法用返回预期答案的函数替换 += 函数。

感谢您的建议。

首先，这是返回预期答案的 R 代码:

B0 = -0.40
B1 =  0.20

mycov1 = exp(B0 + -2 * B1) / (1 + exp(B0 + -2 * B1))
mycov2 = exp(B0 + -1 * B1) / (1 + exp(B0 + -1 * B1))
mycov3 = exp(B0 +  0 * B1) / (1 + exp(B0 +  0 * B1))
mycov4 = exp(B0 +  1 * B1) / (1 + exp(B0 +  1 * B1))

trans1 = matrix(c(1 - 0.25 - mycov1,     mycov1,   0.25 * 0.80,              0,
                                  0,   1 - 0.50,             0,    0.50 * 0.75,
                                  0,          0,             1,              0,
                                  0,          0,             0,              1), 
               nrow=4, ncol=4, byrow=TRUE)

trans2 = matrix(c(1 - 0.25 - mycov2,     mycov2,   0.25 * 0.80,              0,
                                  0,   1 - 0.50,             0,    0.50 * 0.75,
                                  0,          0,             1,              0,
                                  0,          0,             0,              1), 
               nrow=4, ncol=4, byrow=TRUE)

trans3 = matrix(c(1 - 0.25 - mycov3,     mycov3,   0.25 * 0.80,              0,
                                  0,   1 - 0.50,             0,    0.50 * 0.75,
                                  0,          0,             1,              0,
                                  0,          0,             0,              1), 
               nrow=4, ncol=4, byrow=TRUE)

trans4 = matrix(c(1 - 0.25 - mycov4,     mycov4,   0.25 * 0.80,              0,
                                  0,   1 - 0.50,             0,    0.50 * 0.75,
                                  0,          0,             1,              0,
                                  0,          0,             0,              1), 
               nrow=4, ncol=4, byrow=TRUE)

trans2b <- trans1  %*% trans2
trans3b <- trans2b %*% trans3
trans4b <- trans3b %*% trans4
trans4b

#
# This is the expected answer
#
#            [,1]      [,2]      [,3]      [,4]
# [1,] 0.01819965 0.1399834 0.3349504 0.3173467
# [2,] 0.00000000 0.0625000 0.0000000 0.7031250
# [3,] 0.00000000 0.0000000 1.0000000 0.0000000
# [4,] 0.00000000 0.0000000 0.0000000 1.0000000
#

这是我的C代码。 C 代码相当长，因为我不太了解 C 来提高效率:

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

char quit;

int main(){

int i, j, k, ii, jj, kk ;

double B0, B1, mycov ;

double trans[4][4]     = {0} ;

double prevtrans[4][4] = {{1,0,0,0},
                          {0,1,0,0},
                          {0,0,1,0},
                          {0,0,0,1}};

B0 = -0.40 ;
B1 =  0.20 ;

for (i=1; i <= 4; i++) {

          mycov = exp(B0 + B1 * (-2+i-1)) / (1 + exp(B0 + B1 * (-2+i-1))) ;

          trans[0][0] =     1 - 0.25 - mycov ;
          trans[0][1] =                mycov ;
          trans[0][2] =          0.25 * 0.80 ;
          trans[0][3] =                    0 ;

          trans[1][0] =                    0 ;
          trans[1][1] =             1 - 0.50 ;
          trans[1][2] =                    0 ;
          trans[1][3] =          0.50 * 0.75 ;

          trans[2][0] =                    0 ;
          trans[2][1] =                    0 ;
          trans[2][2] =                    1 ;
          trans[2][3] =                    0 ;

          trans[3][0] =                    0 ;
          trans[3][1] =                    0 ;
          trans[3][2] =                    0 ;
          trans[3][3] =                    1 ;

          for (ii=0; ii<4; ii++){

               for(jj=0; jj<4; jj++){

                    for(kk=0; kk<4; kk++){

                         trans[ii][jj] += trans[ii][kk] * prevtrans[kk][jj] ;

                    }
               }
          }

          prevtrans[0][0] =     trans[0][0] ;
          prevtrans[0][1] =     trans[0][1] ;
          prevtrans[0][2] =     trans[0][2] ;
          prevtrans[0][3] =     trans[0][3] ;
          prevtrans[1][0] =     trans[1][0] ;
          prevtrans[1][1] =     trans[1][1] ;
          prevtrans[1][2] =     trans[1][2] ;
          prevtrans[1][3] =     trans[1][3] ;
          prevtrans[2][0] =     trans[2][0] ;
          prevtrans[2][1] =     trans[2][1] ;
          prevtrans[2][2] =     trans[2][2] ;
          prevtrans[2][3] =     trans[2][3] ;
          prevtrans[3][0] =     trans[3][0] ;
          prevtrans[3][1] =     trans[3][1] ;
          prevtrans[3][2] =     trans[3][2] ;
          prevtrans[3][3] =     trans[3][3] ;

}

  printf("To close this program type 'quit' and hit the return key\n");
  printf("               \n");
  scanf("%d", &quit);

  return 0;

}

这是上述 C 代码返回的最终矩阵:

0.4821  3.5870  11.68  381.22
0       1       0      76.875
0       0       5      0
0       0       0      5

Answer 1

这一行

                     trans[ii][jj] += trans[ii][kk] * prevtrans[kk][jj] ;

不对。您正在就地修改 trans，同时仍使用它来计算结果矩阵。您需要另一个矩阵来临时存储乘法结果。然后使用:

      // Store the resultant matrix in temp.
      for (ii=0; ii<4; ii++){

           for(jj=0; jj<4; jj++){

                temp[ii][jj] = 0.0;

                for(kk=0; kk<4; kk++){

                     temp[ii][jj] += trans[ii][kk] * prevtrans[kk][jj] ;

                }
           }
      }

      // Transfer the data from temp to trans
      for (ii=0; ii<4; ii++){

           for(jj=0; jj<4; jj++){
                trans[ii][jj] = temp[ii][jj];
           }
      }