欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页  >  IT编程

矩阵乘法(一):基本运算

程序员文章站 2023-01-04 13:10:59
矩阵,是线性代数中的基本概念之一。一个m×n的矩阵就是m×n个数排成m行n列的一个数阵。在计算机中,一个矩阵实际上就是一个二维数组。因此,可以将矩阵定义为一个结构体: struct Matrix { int mat[110][110]; // 存储矩阵中各元素 int row,col; // 矩阵的 ......

      矩阵,是线性代数中的基本概念之一。一个m×n的矩阵就是m×n个数排成m行n列的一个数阵。在计算机中,一个矩阵实际上就是一个二维数组。因此,可以将矩阵定义为一个结构体:

      struct matrix
      {
            int  mat[110][110]; // 存储矩阵中各元素
            int  row,col; // 矩阵的大小,row行,col列
      };

     矩阵相乘是矩阵的一种基本运算。

     设a为m×n矩阵,b为n×k矩阵,则它们的乘积ab(有时记做a·b)是一个m×k矩阵。

     其乘积矩阵a·b的第i行第j列的元素为第一个矩阵a第i行上的n个数与第二个矩阵b第j列上的n个数对应相乘后所得的n个乘积之和。即:

矩阵乘法(一):基本运算

 

      需要注意的是:只有当矩阵a的列数与矩阵b的行数相等时,矩阵a×b才有意义。因此,矩阵相乘不满足交换律。设a是3×4矩阵,b是4×5矩阵,a与b相乘后,a·b是3×5矩阵;但b·a根本就无法运算。

      矩阵乘法满足结合律。

【例1】矩阵的乘法。
      输入矩阵a和矩阵b的数据,输出新的矩阵c=a*b。

      例如,样例输入

4 3
1 2 3
4 5 6
7 8 9
10 11 12
3 5
7 8 9 10 11
4 5 6 7 8
1 2 3 4 5
样例输出
18 24 30 36 42
54 69 84 99 114
90 114 138 162 186
126 159 192 225 258

      (1)编程思路。

       按照矩阵乘法的定义,用一个三重循环完成运算。

      (2)源程序。

#include <stdio.h>
#include <string.h>
struct matrix
{
      int mat[110][110]; // 存储矩阵中各元素
      int row,col; // 矩阵的大小,row行,col列
};
matrix matmul(matrix a ,matrix b)      // 矩阵a*b
{
      matrix c;
      c.row=a.row;
      c.col=b.col;
      memset(c.mat,0,sizeof(c.mat));
      int i,j,k;
      for (i = 0; i<=a.row ; i++)
         for (j=0 ;j<b.col; j++)
              for (k = 0 ;k<a.col;k++)
                   c.mat[i][j] += a.mat[i][k] * b.mat[k][j];
      return c;
}
int main()
{
      int i,j,x,y;
      matrix a,b,c;
      scanf("%d%d",&x,&y);
      a.row=x;
      a.col=y;
      for (i=0;i<x;i++)
          for (j=0;j<y;j++)
              scanf("%d" ,&a.mat[i][j]);
      scanf("%d%d",&x,&y);
      b.row=x;
      b.col=y;
      for (i=0;i<x;i++)
           for (j=0;j<y;j++)
               scanf("%d" ,&b.mat[i][j]);
      c=matmul(a,b);
      for (i = 0 ;i <c.row;i++)
      {
           for (j=0;j<c.col;j++)
               printf("%5d" ,c.mat[i][j]);
           printf("\n");
      }
      return 0;
}

      在实际应用中,我们经常会用到矩阵的幂运算。

      n个矩阵a相乘称为a的n次方,或称a^n为矩阵a的n次幂。

      求矩阵a的n次方通常采用快速幂运算。下面我们来探讨快速幂运算的思路。

      由于矩阵乘法具有结合律,因此 a^4 = a * a * a * a = (a*a) * (a*a) = a^2 * a^2。由此可以得到这样的结论:当n为偶数时,a^n = a^(n/2) * a^(n/2);当n为奇数时,a^n = a^(n/2) * a^(n/2) * a (其中n/2取整)。这样,我们可以采用一种类似于二分的思想快速求得矩阵的幂。

       例如,a^9 = a*a*a*a*a*a*a*a*a     (一个一个乘,要乘9次)

                        = a*(a*a)*(a*a)*(a*a)*(a*a)

                        = a*(a^2)^4  

                        = a*((a^2)^2)^2   (a平方后,再平方,再平方,再乘上剩下的一个a,要乘4次)

       设c=a^k,c初始化为一个单位矩阵,即c矩阵中除了对角线的元素为1外,其余全部元素为0。

       c.mat[i][i]=1 ,          c,mat[i][j]=0  (i!=j)。

       任何一个矩阵乘以单位矩阵就是它本身。即可以把单位矩阵等价为整数1。因此,矩阵快速幂的算法描述为:

       while (k!=0)

       {

               if  (k%2==1)  c=c*a;     // c=c*a,表示矩阵c与a相乘,结果送c

               a=a*a;

               b=b/2;

       }

       为加深理解,以c=a^9模拟手算一下。

       k=9,  k!=0                       c=c*a (运算结果 c=a)   a=a*a  (运算结果 a=a^2)

       k=k/2  k=4!=0  4%2==0                                              a=a*a   (运算结果a=a^4)

       k=k/2 k=2!=0  2%2==0                                               a =a*a   (运算结果a=a^8) 

       k=2/2 k=1!=0  1%2==1   c=c*a (运算结果 c=a*a^8=a^9)   a=a*a  (运算结果 a=a^16)

       k=1/2  k=0  算法结束。

       可以看出,上述手算过程正好和9的二进制数表示1001相契合。

 

【例2】矩阵快速幂。

      给定n*n的矩阵a,求a^k。

      输入格式:
      第1行, n,k

      第2至n+1行,每行n个数,第i+1行第j个数表示矩阵第i行第j列的元素

      输出格式:
      共n行,每行n个数,第i行第j个数表示矩阵第i行第j列的元素,每个元素模10^5+7

      (1)编程思路。

      因为矩阵的幂运算参与运算的矩阵一定是n*n方阵。因此,在下面的程序中我们将结构体定义简化,去掉表示矩阵行列的变量row和col。

      另外,矩阵乘法运算后,所得结果通常会很大,所以一般采用64位整数表示。同时一般会在计算过程中不断取模,避免高精度运算。

      (2)源程序。

#include <stdio.h>
#include <string.h>
#define modnum 100007
struct matrix
{
      __int64 mat[101][101]; // 存储矩阵中各元素
};
matrix matmul(matrix a ,matrix b,int n)
{
      matrix c;
      memset(c.mat,0,sizeof(c.mat));
      int i,j,k;
      for (i = 1; i<=n ; i++)
          for (j=1 ;j<=n ; j++)
             for (k = 1 ;k<=n ;k++)
             {
                   c.mat[i][j]=(c.mat[i][j]+a.mat[i][k] * b.mat[k][j]) % modnum;
             }
      return c;
}
matrix quickmatpow(matrix a ,int n,int b)   // n阶矩阵a快速b次幂
{
      matrix c;
      memset(c.mat ,0 ,sizeof(c.mat));
      int i;
      for (i = 1 ;i <= n ;i++)
            c.mat[i][i] = 1;
      while (b!=0)
      {
           if (b & 1)
                c = matmul(c ,a ,n);      // c=c*a;
           a = matmul(a ,a ,n);          // a=a*a
           b /= 2;
      }
      return c;
}
int main()
{
      int i,j,n,k;
      matrix a;
      scanf("%d%d",&n,&k);
      for (i=1;i<=n;i++)
           for (j=1;j<=n;j++)
               scanf("%i64d",&a.mat[i][j]);
      a=quickmatpow(a,n,k);
      for (i = 1 ;i <=n;i++)
      {
            for (j=1;j<=n;j++)
                 printf("%i64d ",a.mat[i][j]);
            printf("\n");
      }
      return 0;
}