矩阵乘法(一):基本运算
矩阵,是线性代数中的基本概念之一。一个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;
}
推荐阅读
-
矩阵乘法(一):基本运算
-
矩阵乘法(二):利用矩阵快速幂运算完成递推
-
矩阵乘法(七):其它一些典型应用
-
python人脸识别项目之基础学习(三):矩阵的基本运算 + 张量的阶和形态
-
矩阵乘法(一):基本运算
-
20.7.17 笔记算数运算符 复合运算符重载 比较运算重载 多态 设计原则 类的单一职责 依赖倒置 组合复用原则 里氏替换 迪米特法则 矩阵转置原理
-
矩阵乘法(二):利用矩阵快速幂运算完成递推
-
在python中的实现矩阵的乘法 、矩阵的转置、学习numpy中的矩阵运算方法
-
FPGA学习(第5节)-看电路图写出Verilog代码(乘法运算+自加一+模块实例化等)
-
Go语言基础(一):变量与常量、基本数据类型、运算符