矩阵相乘的快速算法

矩阵相乘的快速算法

算法介绍

矩阵相乘在进行3D变换的时候是经常用到的。在应用中常用矩阵相乘的定义算法对其进行计算。这个算法用到了大量的循环和相乘运算，这使得算法效率不高。而矩阵相乘的计算效率很大程度上的影响了整个程序的运行速度，所以对矩阵相乘算法进行一些改进是必要的。

这里要介绍的矩阵算法称为斯特拉森方法，它是由v.斯特拉森在1969年提出的一个方法。

我们先讨论二阶矩阵的计算方法。

对于二阶矩阵

a11 a12 b11 b12

A = a21 a22

B = b21 b22

先计算下面7个量(1)

x1 = (a11 + a22) * (b11 + b22);

x2 = (a21 + a22) * b11;

x3 = a11 * (b12 - b22);

x4 = a22 * (b21 - b11);

x5 = (a11 + a12) * b22;

x6 = (a21 - a11) * (b11 + b12);

x7 = (a12 - a22) * (b21 + b22);

再设C = AB。根据矩阵相乘的规则，C的各元素为(2)

c11 = a11 * b11 + a12 * b21

c12 = a11 * b12 + a12 * b22

c21 = a21 * b11 + a22 * b21

c22 = a21 * b12 + a22 * b22

比较(1)(2)，C的各元素可以表示为(3)

c11 = x1 + x4 - x5 + x7

c12 = x3 + x5

c21 = x2 + x4

c22 = x1 + x3 - x2 + x6

根据以上的方法，我们就可以计算4阶矩阵了，先将4阶矩阵A和B划分成四块2阶矩阵，分别利用公式计算它们的乘积，再使用(1)(3)来计算出最后结果。

ma11 ma12 mb11 mb12

A4 = ma21 ma22 B4 = mb21 mb22

其中

a11 a12 a13 a14 b11 b12 b13 b14

ma11 = a21 a22 ma12 = a23 a24 mb11 = b21 b22 mb12 = b23 b24

a31 a32 a33 a34 b31 b32 b33 b34

ma21 = a41 a42 ma22 = a43 a44 mb21 = b41 b42 mb22 = b43 b44

实现

// 计算2X2矩阵

void Multiply2X2(float& fOut_11, float& fOut_12, float& fOut_21, float& fOut_22,

float f1_11, float f1_12, float f1_21, float f1_22,

float f2_11, float f2_12, float f2_21, float f2_22)

{

const float x1((f1_11 + f1_22) * (f2_11 + f2_22));

const float x2((f1_21 + f1_22) * f2_11);

const float x3(f1_11 * (f2_12 - f2_22));

const float x4(f1_22 * (f2_21 - f2_11));

const float x5((f1_11 + f1_12) * f2_22);

const float x6((f1_21 - f1_11) * (f2_11 + f2_12));

const float x7((f1_12 - f1_22) * (f2_21 + f2_22));

fOut_11 = x1 + x4 - x5 + x7;

fOut_12 = x3 + x5;

fOut_21 = x2 + x4;

fOut_22 = x1 - x2 + x3 + x6;

}

// 计算4X4矩阵

void Multiply(CLAYMATRIX& mOut, const CLAYMATRIX& m1, const CLAYMATRIX& m2) {

float fTmp[7][4];

// (ma11 + ma22) * (mb11 + mb22)

Multiply2X2(fTmp[0][0], fTmp[0][1], fTmp[0][2], fTmp[0][3],

m1._11 + m1._33, m1._12 + m1._34, m1._21 + m1._43, m1._22 + m1._44,

m2._11 + m2._33, m2._12 + m2._34, m2._21 + m2._43, m2._22 + m2._44);

// (ma21 + ma22) * mb11

Multiply2X2(fTmp[1][0], fTmp[1][1], fTmp[1][2], fTmp[1][3],

m1._31 + m1._33, m1._32 + m1._34, m1._41 + m1._43, m1._42 + m1._44,

m2._11, m2._12, m2._21, m2._22);

// ma11 * (mb12 - mb22)

Multiply2X2(fTmp[2][0], fTmp[2][1], fTmp[2][2], fTmp[2][3],

m1._11, m1._12, m1._21, m1._22,

m2._13 - m2._33, m2._14 - m2._34, m2._23 - m2._43, m2._24 - m2._44);

// ma22 * (mb21 - mb11)

Multiply2X2(fTmp[3][0], fTmp[3][1], fTmp[3][2], fTmp[3][3],

m1._33, m1._34, m1._43, m1._44,

m2._31 - m2._11, m2._32 - m2._12, m2._41 - m2._21, m2._42 - m2._22);

// (ma11 + ma12) * mb22

Multiply2X2(fTmp[4][0], fTmp[4][1], fTmp[4][2], fTmp[4][3],

m1._11 + m1._13, m1._12 + m1._14, m1._21 + m1._23, m1._22 + m1._24,

m2._33, m2._34, m2._43, m2._44);

// (ma21 - ma11) * (mb11 + mb12)

Multiply2X2(fTmp[5][0], fTmp[5][1], fTmp[5][2], fTmp[5][3],

m1._31 - m1._11, m1._32 - m1._12, m1._41 - m1._21, m1._42 - m1._22,

m2._11 + m2._13, m2._12 + m2._14, m2._21 + m2._23, m2._22 + m2._24);

// (ma12 - ma22) * (mb21 + mb22)

Multiply2X2(fTmp[6][0], fTmp[6][1], fTmp[6][2], fTmp[6][3],

m1._13 - m1._33, m1._14 - m1._34, m1._23 - m1._43, m1._24 - m1._44,

m2._31 + m2._33, m2._32 + m2._34, m2._41 + m2._43, m2._42 + m2._44);

// 第一块

mOut._11 = fTmp[0][0] + fTmp[3][0] - fTmp[4][0] + fTmp[6][0];

mOut._12 = fTmp[0][1] + fTmp[3][1] - fTmp[4][1] + fTmp[6][1];

mOut._21 = fTmp[0][2] + fTmp[3][2] - fTmp[4][2] + fTmp[6][2];

mOut._22 = fTmp[0][3] + fTmp[3][3] - fTmp[4][3] + fTmp[6][3];

// 第二块

mOut._13 = fTmp[2][0] + fTmp[4][0];

mOut._14 = fTmp[2][1] + fTmp[4][1];

mOut._23 = fTmp[2][2] + fTmp[4][2];

mOut._24 = fTmp[2][3] + fTmp[4][3];

// 第三块

mOut._31 = fTmp[1][0] + fTmp[3][0];

mOut._32 = fTmp[1][1] + fTmp[3][1];

mOut._41 = fTmp[1][2] + fTmp[3][2];

mOut._42 = fTmp[1][3] + fTmp[3][3];

// 第四块

mOut._33 = fTmp[0][0] - fTmp[1][0] + fTmp[2][0] + fTmp[5][0];

mOut._34 = fTmp[0][1] - fTmp[1][1] + fTmp[2][1] + fTmp[5][1];

mOut._43 = fTmp[0][2] - fTmp[1][2] + fTmp[2][2] + fTmp[5][2];

mOut._44 = fTmp[0][3] - fTmp[1][3] + fTmp[2][3] + fTmp[5][3];

}

比较

在标准的定义算法中我们需要进行n * n * n次乘法运算，新算法中我们需要进行7log2n次乘法，对于最常用的4阶矩阵：原算法新算法

加法次数 48 72(48次加法，24次减法)

乘法次数 64 49

需要额外空间 16 * sizeof(float) 28 * sizeof(float)

新算法要比原算法多了24次减法运算，少了15次乘法。但因为浮点乘法的运算速度要远远慢于加/减法运算，所以新算法的整体速度有所提高。

一、两个矩阵相乘的经典算法:

若设Q=M*N其中，M是m1*n1矩阵，N是m2*n2矩阵。当n1=m2时有：

for (i=1;i

for ( j=1; j<=n2; ++j){

Q[i][j]=0;

for(k=1; k<=n1; ++k) Q[i][j]+=M[i][k]*N[k][j];

}

此算法的时间复杂度是O(m1*n1*n2)。

二、斯特拉森算法

斯特拉森方法，是由v.斯特拉森在1969年提出的一个方法。

我们先讨论二阶矩阵的计算方法。

对于二阶矩阵

a11 a12 b11 b12

A = a21 a22

B = b21 b22

先计算下面7个量(1)

x1 = (a11 + a22) * (b11 + b22);

x2 = (a21 + a22) * b11;

x3 = a11 * (b12 - b22);

x4 = a22 * (b21 - b11);

x5 = (a11 + a12) * b22;

x6 = (a21 - a11) * (b11 + b12);

x7 = (a12 - a22) * (b21 + b22);

再设C = AB。根据矩阵相乘的规则，C的各元素为(2)

c11 = a11 * b11 + a12 * b21

c12 = a11 * b12 + a12 * b22

c21 = a21 * b11 + a22 * b21

c22 = a21 * b12 + a22 * b22

比较(1)(2)，C的各元素可以表示为(3)

c11 = x1 + x4 - x5 + x7

c12 = x3 + x5

c21 = x2 + x4

c22 = x1 + x3 - x2 + x6

根据以上的方法，我们就可以计算4阶矩阵了，先将4阶矩阵A和B划分成四块2阶矩阵，分别利用公式计算它们的乘积，再使用(1)(3)来计算出最后结果。

ma11 ma12 mb11 mb12

A4 = ma21 ma22 B4 = mb21 mb22

其中

a11 a12 a13 a14 b11 b12 b13 b14

ma11 = a21 a22 ma12 = a23 a24 mb11 = b21 b22 mb12 = b23 b24

a31 a32 a33 a34 b31 b32 b33 b34

ma21 = a41 a42 ma22 = a43 a44 mb21 = b41 b42 mb22 = b43 b44

实现

// 计算2X2矩阵

void Multiply2X2(float& fOut_11, float& fOut_12, float& fOut_21, float& fOut_22, float f1_11, float f1_12, float f1_21, float f1_22,

float f2_11, float f2_12, float f2_21, float f2_22)

{

const float x1((f1_11 + f1_22) * (f2_11 + f2_22));

const float x2((f1_21 + f1_22) * f2_11);

const float x3(f1_11 * (f2_12 - f2_22));

const float x4(f1_22 * (f2_21 - f2_11));

const float x5((f1_11 + f1_12) * f2_22);

const float x6((f1_21 - f1_11) * (f2_11 + f2_12));

const float x7((f1_12 - f1_22) * (f2_21 + f2_22));

fOut_11 = x1 + x4 - x5 + x7;

fOut_12 = x3 + x5;

fOut_21 = x2 + x4;

fOut_22 = x1 - x2 + x3 + x6;

}

// 计算4X4矩阵

void Multiply(CLAYMATRIX& mOut, const CLAYMATRIX& m1, const CLAYMATRIX& m2)

{

float fTmp[7][4];

// (ma11 + ma22) * (mb11 + mb22)

Multiply2X2(fTmp[0][0], fTmp[0][1], fTmp[0][2], fTmp[0][3],

m1._11 + m1._33, m1._12 + m1._34, m1._21 + m1._43, m1._22 + m1._44,

m2._11 + m2._33, m2._12 + m2._34, m2._21 + m2._43, m2._22 + m2._44);

// (ma21 + ma22) * mb11

Multiply2X2(fTmp[1][0], fTmp[1][1], fTmp[1][2], fTmp[1][3],

m1._31 + m1._33, m1._32 + m1._34, m1._41 + m1._43, m1._42 + m1._44,

m2._11, m2._12, m2._21, m2._22);

// ma11 * (mb12 - mb22)

Multiply2X2(fTmp[2][0], fTmp[2][1], fTmp[2][2], fTmp[2][3],

m1._11, m1._12, m1._21, m1._22,

m2._13 - m2._33, m2._14 - m2._34, m2._23 - m2._43, m2._24 - m2._44);

// ma22 * (mb21 - mb11)

Multiply2X2(fTmp[3][0], fTmp[3][1], fTmp[3][2], fTmp[3][3],

m1._33, m1._34, m1._43, m1._44,

m2._31 - m2._11, m2._32 - m2._12, m2._41 - m2._21, m2._42 - m2._22);

// (ma11 + ma12) * mb22

Multiply2X2(fTmp[4][0], fTmp[4][1], fTmp[4][2], fTmp[4][3],

m1._11 + m1._13, m1._12 + m1._14, m1._21 + m1._23, m1._22 + m1._24,

m2._33, m2._34, m2._43, m2._44);

// (ma21 - ma11) * (mb11 + mb12)

Multiply2X2(fTmp[5][0], fTmp[5][1], fTmp[5][2], fTmp[5][3],

m1._31 - m1._11, m1._32 - m1._12, m1._41 - m1._21, m1._42 - m1._22,

m2._11 + m2._13, m2._12 + m2._14, m2._21 + m2._23, m2._22 + m2._24);

// (ma12 - ma22) * (mb21 + mb22)

Multiply2X2(fTmp[6][0], fTmp[6][1], fTmp[6][2], fTmp[6][3],

m1._13 - m1._33, m1._14 - m1._34, m1._23 - m1._43, m1._24 - m1._44,

m2._31 + m2._33, m2._32 + m2._34, m2._41 + m2._43, m2._42 + m2._44);

// 第一块

mOut._11 = fTmp[0][0] + fTmp[3][0] - fTmp[4][0] + fTmp[6][0];

mOut._12 = fTmp[0][1] + fTmp[3][1] - fTmp[4][1] + fTmp[6][1];

mOut._21 = fTmp[0][2] + fTmp[3][2] - fTmp[4][2] + fTmp[6][2];

mOut._22 = fTmp[0][3] + fTmp[3][3] - fTmp[4][3] + fTmp[6][3];

// 第二块

mOut._13 = fTmp[2][0] + fTmp[4][0];

mOut._14 = fTmp[2][1] + fTmp[4][1];

mOut._23 = fTmp[2][2] + fTmp[4][2];

mOut._24 = fTmp[2][3] + fTmp[4][3];

// 第三块

mOut._31 = fTmp[1][0] + fTmp[3][0];

mOut._32 = fTmp[1][1] + fTmp[3][1];

mOut._41 = fTmp[1][2] + fTmp[3][2];

mOut._42 = fTmp[1][3] + fTmp[3][3];

// 第四块

mOut._33 = fTmp[0][0] - fTmp[1][0] + fTmp[2][0] + fTmp[5][0];

mOut._34 = fTmp[0][1] - fTmp[1][1] + fTmp[2][1] + fTmp[5][1];

mOut._43 = fTmp[0][2] - fTmp[1][2] + fTmp[2][2] + fTmp[5][2];

mOut._44 = fTmp[0][3] - fTmp[1][3] + fTmp[2][3] + fTmp[5][3];

}

比较

在标准的定义算法中我们需要进行n * n * n次乘法运算，新算法中我们需要进行7log2n次乘法，对于最常用的4阶矩阵：原算法新算法

加法次数 48 72(48次加法，24次减法)

乘法次数 64 49

需要额外空间 16 * sizeof(float) 28 * sizeof(float)

新算法要比原算法多了24次减法运算，少了15次乘法。但因为浮点乘法的运算速度要远远慢于加/减法运算，所以新算法的整体速度有所提高。

三、Strassen矩阵乘法

矩阵乘法是线性代数中最常见的运算之一，它在数值计算中有广泛的应用。若A和B是2个n×n的矩阵，则它们的乘积C=AB同样是一个n×n的矩阵。A和B的乘积矩阵C中的元素C[i,j]定义为:

若依此定义来计算A和B的乘积矩阵C，则每计算C的一个元素C[i,j]，需要做n个乘法和n-1次加法。因此，求出矩阵C的n2个元素所需的计算时间为0(n3)。

60年代末，Strassen采用了类似于在大整数乘法中用过的分治技术，将计算2个n阶矩阵乘积所需的计算时间改进到O(n log7)=O(n2.18)。

首先，我们还是需要假设n是2的幂。将矩阵A，B和C中每一矩阵都分块成为4个大小相等的子矩阵，每个子矩阵都是n/2×n/2的方阵。由此可将方程C=AB重写为:

(1)

由此可得:

C11=A11B11+A12B21 (2)

C12=A11B12+A12B22 (3)

C21=A21B11+A22B21 (4)

C22=A21B12+A22B22 (5)

如果n=2，则2个2阶方阵的乘积可以直接用(2)-(3)式计算出来，共需8次乘法和4次加法。当子矩阵的阶大于2时，为求2个子矩阵的积，可以继续将子矩阵分块，直到子矩阵的阶降为2。这样，就产生了一个分治降阶的递归算法。依此算法，计算2个n阶方阵的乘积转化为计算8个n/2阶方阵的乘积和4个n/2阶方阵的加法。2个n/2×n/2矩阵的加法显然可以在c*n2/4时间内完成，这里c是一个常数。因此，上述分治法的计算时间耗费T(n)应该满足：

这个递归方程的解仍然是T(n)=O(n3)。因此，该方法并不比用原始定义直接计算更有效。究其原因，乃是由于式(2)-(5)并没有减少矩阵的乘法次数。而矩阵乘法耗费的时间要比矩阵加减法耗费的时间多得多。要想改进矩阵乘法的计算时间复杂性，必须减少子矩阵乘法运算的次数。按照上述分治法的思想可以看出，要想减少乘法运算次数，关键在于计算2个2阶方阵的乘积时，能否用少于8次的乘法运算。Strassen提出了一种新的算法来计算2个2阶方阵的乘积。他的算法只用了7次乘法运算，但增加了加、减法的运算次数。这7次乘法是:

M1=A11(B12-B22)

M2=(A11+A12)B22

M3=(A21+A22)B11

M4=A22(B21-B11)

M5=(A11+A22)(B11+B22)

M6=(A12-A22)(B21+B22)

M7=(A11-A21)(B11+B12)

做了这7次乘法后，再做若干次加、减法就可以得到: C11=M5+M4-M2+M6

C12=M1+M2

C21=M3+M4

C22=M5+M1-M3-M7

以上计算的正确性很容易验证。例如:

C22=M5+M1-M3-M7

=(A11+A22)(B11+B22)+A11(B12-B22)-(A21+A22)B11-(A11-A21)(B11+B12) =A11B11+A11B22+A22B11+A22B22+A11B12

-A11B22-A21B11-A22B11-A11B11-A11B12+A21B11+A21B12

=A21B12+A22B22

由(2)式便知其正确性。

至此，我们可以得到完整的Strassen算法如下：

procedure STRASSEN(n,A,B,C);

begin

if n=2 then MATRIX-MULTIPLY(A，B，C)

else begin

将矩阵A和B依(1)式分块;

STRASSEN(n/2,A11,B12-B22,M1);

STRASSEN(n/2,A11+A12,B22,M2);

STRASSEN(n/2,A21+A22,B11,M3);

STRASSEN(n/2,A22,B21-B11,M4);

STRASSEN(n/2,A11+A22,B11+B22,M5);

STRASSEN(n/2,A12-A22,B21+B22,M6);

STRASSEN(n/2,A11-A21,B11+B12,M7);

;

end;

end;

其中MATRIX-MULTIPLY(A，B，C)是按通常的矩阵乘法计算C=AB的子算法。

Strassen矩阵乘积分治算法中，用了7次对于n/2阶矩阵乘积的递归调用和18次n/2阶矩阵的加减运算。由此可知，该算法的所需的计算时间T(n)满足如下的递归方程:

按照解递归方程的套用公式法，其解为T(n)=O(n log7)≈O(n2.81)。由此可见，Strassen矩阵乘法的计算时间复杂性比普通矩阵乘法有阶的改进。

有人曾列举了计算2个2阶矩阵乘法的36种不同方法。但所有的方法都要做7次乘法。除非能找到一种计算2阶方阵乘积的算法，使乘法的计算次数少于7次，按上述思路才有可能进一步改进矩阵乘积的计算时间的上界。但

是Hopcroft和Kerr(197l)已经证明，计算2个2×2矩阵的乘积，7次乘法是必要的。因此，要想进一步改进矩阵乘法的时间复杂性，就不能再寄希望于计算2×2矩阵的乘法次数的减少。或许应当研究3×3或5×5矩阵的更好算法。在Strassen之后又有许多算法改进了矩阵乘法的计算时间复杂性。目前最好的计算时间上界是O(n2.367)。而目前所知道的矩阵乘法的最好下界仍是它的平凡下界Ω(n2)。因此到目前为止还无法确切知道矩阵乘法的时间复杂性。关于这一研究课题还有许多工作可做。