矩阵相乘在3D变换中是被频繁用到的一种计算,但在矩阵相乘过程中用到了大量的乘法运算,而cpu中运算单元对于乘法的效率是比较低的,远低于加法运算, 所以,如果能找到一种用加法来替代乘法的方法实现矩阵相乘,将能大大提高我们程序的效率。我们的确有这种方法,这就是网上甚为流行的斯特拉森矩阵乘法,它是由v.斯特拉森在1969年提出的一个方法。
下面对其进行详细介绍.
一,推导
对于二阶矩阵
A = [a11 a12]
[a21 a22]
B = [b11 b12]
[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)来计算出最后结果。
A4 = [ma11 ma12]
[ma21 ma22]
B4 = [mb11 mb12]
[mb21 mb22]
其中
ma11 = [a11 a12]
[a21 a22]
ma12 = [a13 a14]
[a23 a24]
ma21 = [a31 a32]
[a41 a42]
ma22 = [a33 a34]
[a43 a44]
mb11 = [b11 b12]
[b21 b22]
mb12 = [b13 b14]
[b23 b24]
mb21 = [b31 b32]
[b41 b42]
mb22 = [b33 b34]
[b43 b44]
二,实现
typedef float Matrix44[ 4 ][ 4 ];
inline void Matrix22MulMatrix22(Matrix22 c, const Matrix22 & a, const Matrix22 & b)
{
float x1 = (a[ 0 ][ 0 ] + a[ 1 ][ 1 ]) * (b[ 0 ][ 0 ] + b[ 1 ][ 1 ]);
float x2 = (a[ 1 ][ 0 ] + a[ 1 ][ 1 ]) * b[ 0 ][ 0 ];
float x3 = a[ 0 ][ 0 ] * (b[ 0 ][ 1 ] - b[ 1 ][ 1 ]);
float x4 = a[ 1 ][ 1 ] * (b[ 1 ][ 0 ] - b[ 0 ][ 0 ]);
float x5 = (a[ 0 ][ 0 ] + a[ 0 ][ 1 ]) * b[ 1 ][ 1 ];
float x6 = (a[ 1 ][ 0 ] - a[ 0 ][ 0 ]) * (b[ 0 ][ 0 ] + b[ 0 ][ 1 ]);
float x7 = (a[ 0 ][ 1 ] - a[ 1 ][ 1 ]) * (b[ 1 ][ 0 ] + b[ 1 ][ 1 ]);
c[ 0 ][ 0 ] = x1 + x4 - x5 + x7;
c[ 0 ][ 1 ] = x3 + x5;
c[ 1 ][ 0 ] = x2 + x4;
c[ 1 ][ 1 ] = x1 + x3 - x2 + x6;
}
inline void Matrix44MulMatrix44(Matrix44 c, const Matrix44 & a, const Matrix44 & b)
{
Matrix22 x[ 7 ];
// (ma11 + ma22) * (mb11 + mb22)
Matrix22 a0 = {a[ 0 ][ 0 ] + a[ 2 ][ 2 ], a[ 0 ][ 1 ] + a[ 2 ][ 3 ], a[ 1 ][ 0 ] + a[ 3 ][ 2 ], a[ 1 ][ 1 ] + a[ 3 ][ 3 ]};
Matrix22 b0 = {b[ 0 ][ 0 ] + b[ 2 ][ 2 ], b[ 0 ][ 1 ] + b[ 2 ][ 3 ], b[ 1 ][ 0 ] + b[ 3 ][ 2 ], b[ 1 ][ 1 ] + b[ 3 ][ 3 ]};
Matrix22MulMatrix22(x[ 0 ], a0, b0);
// (ma21 + ma22) * mb11
Matrix22 a1 = {a[ 2 ][ 0 ] + a[ 2 ][ 2 ], a[ 2 ][ 1 ] + a[ 2 ][ 3 ], a[ 3 ][ 0 ] + a[ 3 ][ 2 ], a[ 3 ][ 1 ] + a[ 3 ][ 3 ]};
Matrix22 b1 = {b[ 0 ][ 0 ], b[ 0 ][ 1 ], b[ 1 ][ 0 ], b[ 1 ][ 1 ]};
Matrix22MulMatrix22(x[ 1 ], a1, b1);
// ma11 * (mb12 - mb22)
Matrix22 a2 = {a[ 0 ][ 0 ], a[ 0 ][ 1 ], a[ 1 ][ 0 ], a[ 1 ][ 1 ]};
Matrix22 b2 = {b[ 0 ][ 2 ] - b[ 2 ][ 2 ], b[ 0 ][ 3 ] - b[ 2 ][ 3 ], b[ 1 ][ 2 ] - b[ 3 ][ 2 ], b[ 1 ][ 3 ] - b[ 3 ][ 3 ]};
Matrix22MulMatrix22(x[ 2 ], a2, b2);
// ma22 * (mb21 - mb11)
Matrix22 a3 = {a[ 2 ][ 2 ], a[ 2 ][ 3 ], a[ 3 ][ 2 ], a[ 3 ][ 3 ]};
Matrix22 b3 = {b[ 2 ][ 0 ] - b[ 0 ][ 0 ], b[ 2 ][ 1 ] - b[ 0 ][ 1 ], b[ 3 ][ 0 ] - b[ 1 ][ 0 ], b[ 3 ][ 1 ] - b[ 1 ][ 1 ]};
Matrix22MulMatrix22(x[ 3 ], a3, b3);
// (ma11 + ma12) * mb22
Matrix22 a4 = {a[ 0 ][ 0 ] + a[ 0 ][ 2 ], a[ 0 ][ 1 ] + a[ 0 ][ 3 ], a[ 1 ][ 0 ] + a[ 1 ][ 2 ], a[ 1 ][ 1 ] + a[ 1 ][ 3 ]};
Matrix22 b4 = {b[ 2 ][ 2 ], b[ 2 ][ 3 ], b[ 3 ][ 2 ], b[ 3 ][ 3 ]};
Matrix22MulMatrix22(x[ 4 ], a4, b4);
// (ma21 - ma11) * (mb11 + mb12)
Matrix22 a5 = {a[ 2 ][ 0 ] - a[ 0 ][ 0 ], a[ 2 ][ 1 ] - a[ 0 ][ 1 ], a[ 3 ][ 0 ] - a[ 1 ][ 0 ], a[ 3 ][ 1 ] - a[ 1 ][ 1 ]};
Matrix22 b5 = {b[ 0 ][ 0 ] + b[ 0 ][ 2 ], b[ 0 ][ 1 ] + b[ 0 ][ 3 ], b[ 1 ][ 0 ] + b[ 1 ][ 2 ], b[ 1 ][ 1 ] + b[ 1 ][ 3 ]};
Matrix22MulMatrix22(x[ 5 ], a5, b5);
// (ma12 - ma22) * (mb21 + mb22)
Matrix22 a6 = {a[ 0 ][ 2 ] - a[ 2 ][ 2 ], a[ 0 ][ 3 ] - a[ 2 ][ 3 ], a[ 1 ][ 2 ] - a[ 3 ][ 2 ], a[ 1 ][ 3 ] - a[ 3 ][ 3 ]};
Matrix22 b6 = {b[ 2 ][ 0 ] + b[ 2 ][ 2 ], b[ 2 ][ 1 ] + b[ 2 ][ 3 ], b[ 3 ][ 0 ] + b[ 3 ][ 2 ], b[ 3 ][ 1 ] + b[ 3 ][ 3 ]};
Matrix22MulMatrix22(x[ 6 ], a6, b6);
// 第一块
c[ 0 ][ 0 ] = x[ 0 ][ 0 ][ 0 ] + x[ 3 ][ 0 ][ 0 ] - x[ 4 ][ 0 ][ 0 ] + x[ 6 ][ 0 ][ 0 ];
c[ 0 ][ 1 ] = x[ 0 ][ 0 ][ 1 ] + x[ 3 ][ 0 ][ 1 ] - x[ 4 ][ 0 ][ 1 ] + x[ 6 ][ 0 ][ 1 ];
c[ 1 ][ 0 ] = x[ 0 ][ 1 ][ 0 ] + x[ 3 ][ 1 ][ 0 ] - x[ 4 ][ 1 ][ 0 ] + x[ 6 ][ 1 ][ 0 ];
c[ 1 ][ 1 ] = x[ 0 ][ 1 ][ 1 ] + x[ 3 ][ 1 ][ 1 ] - x[ 4 ][ 1 ][ 1 ] + x[ 6 ][ 1 ][ 1 ];
// 第二块
c[ 0 ][ 2 ] = x[ 2 ][ 0 ][ 0 ] + x[ 4 ][ 0 ][ 0 ];
c[ 0 ][ 3 ] = x[ 2 ][ 0 ][ 1 ] + x[ 4 ][ 0 ][ 1 ];
c[ 1 ][ 2 ] = x[ 2 ][ 1 ][ 0 ] + x[ 4 ][ 1 ][ 0 ];
c[ 1 ][ 3 ] = x[ 2 ][ 1 ][ 1 ] + x[ 4 ][ 1 ][ 1 ];
// 第三块
c[ 2 ][ 0 ] = x[ 1 ][ 0 ][ 0 ] + x[ 3 ][ 0 ][ 0 ];
c[ 2 ][ 1 ] = x[ 1 ][ 0 ][ 1 ] + x[ 3 ][ 0 ][ 1 ];
c[ 3 ][ 0 ] = x[ 1 ][ 1 ][ 0 ] + x[ 3 ][ 1 ][ 0 ];
c[ 3 ][ 1 ] = x[ 1 ][ 1 ][ 1 ] + x[ 3 ][ 1 ][ 1 ];
// 第四块
c[ 2 ][ 2 ] = x[ 0 ][ 0 ][ 0 ] - x[ 1 ][ 0 ][ 0 ] + x[ 2 ][ 0 ][ 0 ] + x[ 5 ][ 0 ][ 0 ];
c[ 2 ][ 3 ] = x[ 0 ][ 0 ][ 1 ] - x[ 1 ][ 0 ][ 1 ] + x[ 2 ][ 0 ][ 1 ] + x[ 5 ][ 0 ][ 1 ];
c[ 3 ][ 2 ] = x[ 0 ][ 1 ][ 0 ] - x[ 1 ][ 1 ][ 0 ] + x[ 2 ][ 1 ][ 0 ] + x[ 5 ][ 1 ][ 0 ];
c[ 3 ][ 3 ] = x[ 0 ][ 1 ][ 1 ] - x[ 1 ][ 1 ][ 1 ] + x[ 2 ][ 1 ][ 1 ] + x[ 5 ][ 1 ][ 1 ];
}
三,分析
在标准的定义算法中我们需要进行n * n * n次乘法运算,新算法中我们需要进行7log2n次乘法,对于最常用的4阶矩阵:
原算法 新算法
加法次数 48 72(48次加法,24次减法)
乘法次数 64 49
需要额外空间 16 * sizeof(float) 28 * sizeof(float) (+2 * 4 * 7 * sizeof(float))
新算法要比原算法多了24次减法运算,少了15次乘法。但因为浮点乘法的运算速度要远远慢于加/减法运算,所以新算法的整体速度有所提高。
四,其他
这里列出了按通常公式计算矩阵乘法的函数,以作参考。感谢我的女朋友帮我完成了这两个函数:)值得一提的是我女朋友是学文科的,从不知道什么是矩阵,当然也没写过程序,但在我稍微指点了一下后,等我洗漱完回来,她已经写好了,经检查测试通过,把她高兴的...
{
c[ 0 ][ 0 ] = a[ 0 ][ 0 ] * b[ 0 ][ 0 ] + a[ 0 ][ 1 ] * b[ 1 ][ 0 ];
c[ 0 ][ 1 ] = a[ 0 ][ 0 ] * b[ 0 ][ 1 ] + a[ 0 ][ 1 ] * b[ 1 ][ 1 ];
c[ 1 ][ 0 ] = a[ 1 ][ 0 ] * b[ 0 ][ 0 ] + a[ 1 ][ 1 ] * b[ 1 ][ 0 ];
c[ 1 ][ 1 ] = a[ 1 ][ 0 ] * b[ 0 ][ 1 ] + a[ 1 ][ 1 ] * b[ 1 ][ 1 ];
}
inline void Matrix44MulMatrix44_(Matrix44 c, const Matrix44 & a, const Matrix44 & b)
{
c[ 0 ][ 0 ] = a[ 0 ][ 0 ] * b[ 0 ][ 0 ] + a[ 0 ][ 1 ] * b[ 1 ][ 0 ] + a[ 0 ][ 2 ] * b[ 2 ][ 0 ] + a[ 0 ][ 3 ] * b[ 3 ][ 0 ];
c[ 0 ][ 1 ] = a[ 0 ][ 0 ] * b[ 0 ][ 1 ] + a[ 0 ][ 1 ] * b[ 1 ][ 1 ] + a[ 0 ][ 2 ] * b[ 2 ][ 1 ] + a[ 0 ][ 3 ] * b[ 3 ][ 1 ];
c[ 0 ][ 2 ] = a[ 0 ][ 0 ] * b[ 0 ][ 2 ] + a[ 0 ][ 1 ] * b[ 1 ][ 2 ] + a[ 0 ][ 2 ] * b[ 2 ][ 2 ] + a[ 0 ][ 3 ] * b[ 3 ][ 2 ];
c[ 0 ][ 3 ] = a[ 0 ][ 0 ] * b[ 0 ][ 3 ] + a[ 0 ][ 1 ] * b[ 1 ][ 3 ] + a[ 0 ][ 2 ] * b[ 2 ][ 3 ] + a[ 0 ][ 3 ] * b[ 3 ][ 3 ];
c[ 1 ][ 0 ] = a[ 1 ][ 0 ] * b[ 0 ][ 0 ] + a[ 1 ][ 1 ] * b[ 1 ][ 0 ] + a[ 1 ][ 2 ] * b[ 2 ][ 0 ] + a[ 1 ][ 3 ] * b[ 3 ][ 0 ];
c[ 1 ][ 1 ] = a[ 1 ][ 0 ] * b[ 0 ][ 1 ] + a[ 1 ][ 1 ] * b[ 1 ][ 1 ] + a[ 1 ][ 2 ] * b[ 2 ][ 1 ] + a[ 1 ][ 3 ] * b[ 3 ][ 1 ];
c[ 1 ][ 2 ] = a[ 1 ][ 0 ] * b[ 0 ][ 2 ] + a[ 1 ][ 1 ] * b[ 1 ][ 2 ] + a[ 1 ][ 2 ] * b[ 2 ][ 2 ] + a[ 1 ][ 3 ] * b[ 3 ][ 2 ];
c[ 1 ][ 3 ] = a[ 1 ][ 0 ] * b[ 0 ][ 3 ] + a[ 1 ][ 1 ] * b[ 1 ][ 3 ] + a[ 1 ][ 2 ] * b[ 2 ][ 3 ] + a[ 1 ][ 3 ] * b[ 3 ][ 3 ];
c[ 2 ][ 0 ] = a[ 2 ][ 0 ] * b[ 0 ][ 0 ] + a[ 2 ][ 1 ] * b[ 1 ][ 0 ] + a[ 2 ][ 2 ] * b[ 2 ][ 0 ] + a[ 2 ][ 3 ] * b[ 3 ][ 0 ];
c[ 2 ][ 1 ] = a[ 2 ][ 0 ] * b[ 0 ][ 1 ] + a[ 2 ][ 1 ] * b[ 1 ][ 1 ] + a[ 2 ][ 2 ] * b[ 2 ][ 1 ] + a[ 2 ][ 3 ] * b[ 3 ][ 1 ];
c[ 2 ][ 2 ] = a[ 2 ][ 0 ] * b[ 0 ][ 2 ] + a[ 2 ][ 1 ] * b[ 1 ][ 2 ] + a[ 2 ][ 2 ] * b[ 2 ][ 2 ] + a[ 2 ][ 3 ] * b[ 3 ][ 2 ];
c[ 2 ][ 3 ] = a[ 2 ][ 0 ] * b[ 0 ][ 3 ] + a[ 2 ][ 1 ] * b[ 1 ][ 3 ] + a[ 2 ][ 2 ] * b[ 2 ][ 3 ] + a[ 2 ][ 3 ] * b[ 3 ][ 3 ];
c[ 3 ][ 0 ] = a[ 3 ][ 0 ] * b[ 0 ][ 0 ] + a[ 3 ][ 1 ] * b[ 1 ][ 0 ] + a[ 3 ][ 2 ] * b[ 2 ][ 0 ] + a[ 3 ][ 3 ] * b[ 3 ][ 0 ];
c[ 3 ][ 1 ] = a[ 3 ][ 0 ] * b[ 0 ][ 1 ] + a[ 3 ][ 1 ] * b[ 1 ][ 1 ] + a[ 3 ][ 2 ] * b[ 2 ][ 1 ] + a[ 3 ][ 3 ] * b[ 3 ][ 1 ];
c[ 3 ][ 2 ] = a[ 3 ][ 0 ] * b[ 0 ][ 2 ] + a[ 3 ][ 1 ] * b[ 1 ][ 2 ] + a[ 3 ][ 2 ] * b[ 2 ][ 2 ] + a[ 3 ][ 3 ] * b[ 3 ][ 2 ];
c[ 3 ][ 3 ] = a[ 3 ][ 0 ] * b[ 0 ][ 3 ] + a[ 3 ][ 1 ] * b[ 1 ][ 3 ] + a[ 3 ][ 2 ] * b[ 2 ][ 3 ] + a[ 3 ][ 3 ] * b[ 3 ][ 3 ];
}
当然, 这个用for循环写出来要简洁些,但是,这样更原汁原味:)