矩阵运算是现代科学及工程计算的基石之一,而矩阵乘法则是其中最常见一种运算。对于二维矩阵A、B,如果A的列数等于B的行数,则矩阵A、B可乘,其结果矩阵C的行数等于A的行数,列数等于B的列数。

  形式化表达为:

e1

  其中,cij表示结果矩阵中第i行第j列的元素,其计算公式为:

e2

版权声明

本文可以在互联网上自由转载,但必须:注明出处(作者:海洋饼干叔叔)并包含指向本页面的链接。

本文不可以以纸质出版为目的进行改编、摘抄。

  下述示例中,一个3行2列的矩阵乘以一个2行3列的矩阵,其结果矩阵为3行3列。

m1

  对于结果矩阵的第i行第j列的元素,其值正好等于A阵的第i行和B阵的第j列的元素两两相乘,再求和。本例中,结果矩阵第2行第3列的元素值为49,它通过下列计算而得:

$$
4 \times 4+3 \times 11=49
$$


  请读者注意,上述表述中我们依从数学课本的习惯,行列从1开始计数。后续示例代码则按照程序设计的规则从0开始计数。

  下述代码中的multiMatrix()函数实现了不同行列数的二维矩阵的通用乘法。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
//Project - MulMatrix
#include <iostream>
#include <iomanip> //引入格式操作算子(manipulator),详见2.3节扩展阅读

bool multiMatrix(const int a[], int aRow, int aCol,
const int b[], int bRow, int bCol,
int c[], int cRow, int cCol){
if ((aCol!=bRow) || (aRow!=cRow) || (bCol!=cCol)
|| (aRow<=0) || (bRow<=0) || (bCol<=0))
return false;

for (int i=0;i<cRow;i++)
for (int j=0;j<cCol;j++){
c[i*cCol+j] = 0;
for (int k=0;k<aCol;k++)
c[i*cCol+j] += a[i*aCol+k] * b[k*bCol+j];
}
return true;
}

int main(){
using namespace std;
int a[] = { 5,7,4,3,6,1 };
int b[] = { -3,2,4,4,1,11 };
int c[9];
if (!multiMatrix(a,3,2,b,2,3,c,3,3)){
cout << "Illegal parameters for matrix multiply.";
return 0;
}

for (int i=0;i<3;i++){
for (int j=0;j<3;j++)
cout << setw(5) << c[i*3+j]; //set(w)设置输出宽度为5个字符
cout << "\n";
}
return 0;
}

上述代码的执行结果为:

1
2
3
 13   17   97
0 11 49
-14 13 35

  依本章稍早的讨论,当程序根据下标去访问二维数组的特定元素时,依赖于每行的元素个数,即列数来计算单个元素的存储位置。对于multiMatrix()函数,如果我们以二维数组的形式设计形参的话,就只能将其声明为下述形式:

1
2
3
bool multiMatrix(const int a[][2], int aRow,       //二维数组作为参数
const int b[][3], int bRow,
int c[][3], int cRow);

  该函数声明包含6个形参,其中a、b两个二维数组所代表的矩阵相乘的结果矩阵存储在二维数组c当中,aRow、bRow、cRow分别代表了三个矩阵的行数,矩阵的列数已包含在形参a、b、c的定义中。 请读者注意,上述声明中二维数组的行数被省略,这是允许的,因为程序依赖于列数来计算单个元素的存储位置,但不依赖于行数。

  显然,这种以二维数组来表示矩阵的方式使得multiMatrix()函数不够通用,其只能计算特定维度的矩阵的乘积。为了确保multiMatrix()函数的通用性,我们将二维数组“扁平化”,降维成一维数组。当我们在函数中按行号i和列号j访问矩阵的特定元素时,总是通过计算求得其在对应的一维数组中的下标来进行访问。

  不失一般性,对于二维数组C2,如果将其降维存储为一维数组C1,则有:

$$
C_{2}[i][j]=C_{1}[i \times \text { 列数 }+\mathrm{j}]
$$


🚩第5行:形参一维数组a代表被乘矩阵A,在程序中,其为一维数组,但读者需要将其理解为aRow行aCol列的二维数组。在函数执行过程中,矩阵A预期不应被修改,故a的元素类型设定为const int。请读者注意,const int a[]表示a为一个一维数组,其元素类型为const int,其元素个数不确定。当程序第26行的函数调用发生时,实参a向形参a的传递并非传递整个数组,实际被传递的是数组名,即数组首元素的地址。在函数内对形参数组a的访问,即是对第23行的实参数组a的访问。

🚩第6行:形参一维数组b代表乘数矩阵B,其含义同形参a。

🚩第7行:形参一维数组c代表结果矩阵C,其行数为cRow,列数为cCol。请读者注意,确保c数组有恰当的存储空间以容纳结果矩阵,是外部程序而不是本函数的职责。

🚩第8 ~ 10行:对参数的合法性进行检查,包括行列数是否大于0,A的列数是否等于B的行数,A的行数是否等于C的行数等。如不合法,返回false。

🚩第12 ~ 17行:按照计算公式逐一计算结果矩阵C中每一个元素cij 的值。第12 ~ 13行是一个双重循环,其内层循环即j循环的循环体用于计算cij 的值。第14行首先将cij 置0,第15 ~ 16行则依公式将矩阵A的第i行与矩阵B的第j列两两相乘,并累加至cij 。请读者留意第14 ~ 17行中一维数组a、b、c的下标,其是按“ 行号x列数+列号”的形式计算而得的,代表了元素在“假想的”的二维数组中的对应位置。

🚩第23行:3行2列矩阵A的扁平化一维数组。逻辑上,矩阵A的值如下:

m2

🚩第24行:2行3列矩阵B的扁平化一维数组。

🚩第25行:3行3列的结果矩阵C的扁平化一维数组,其包含9个元素。

🚩第26 ~ 29行:调用multiMatrix()函数计算矩阵A乘以矩阵B的积。如果计算失败,打印错误信息并返回。

🚩第31 ~ 35行:将结果矩阵打印出来。其中,setw(5)格式操作算子(manipulator)指定了整数输出的宽度为5个字符。第34行代码在每行数据输出完毕后添加一个换行符。同样地,第33行的下标[i*3+j]对应“假想”的二维数组下标[i][j],在该“假想”的二维数组中,列数为3。

  关于格式操作算子及通过cout进行复杂格式的输出,请参见2.3节扩展阅读部分。



扩展阅读📕


矩阵乘法的计算复杂性

  在上述示例中,结果矩阵包括m行n列共mn个元素,对于每个结果元素的计算,需要执行一个q轮的循环。因此,整个计算过程中第16行代码需要执行mnq次。在代码的第16行中,包括至少4次乘法,4次加法运算,如果我们将这些运算用“1”次来代替,则最终的估计值比实际值至少小8倍,但不会超过常数倍。如此简化之后,我们称一次矩阵乘法需要执行mnq次“基本运算”。进一步地,矩阵的行数列数之比不会超过常数倍,故我们用参数n来表征矩阵的行/列数,称一次矩阵乘法需要执行n3次“基本运算”。如此简化之后,估计值与实际值之比也不会超过常数倍。在《算法分析》这样的课程中,我们称矩阵乘法的计算复杂性为θ(n3) 。请读者注意,这里所谓的常数,是指任意小于无穷大的确定值。之所以在计算复杂性估计的过程中可以忽略常数倍的差异,理由之一是常数倍的计算量差异可以通过提高计算机速度等途径来弥补。