SmartPtr
本博客已搬至:http://www.cnblogs.com/baiyanhuang/
posts - 29,comments - 176,trackbacks - 0

By SmartPtr(http://www.cppblog.com/SmartPtr/)

  矩阵相乘在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 Matrix22[2][2];
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次乘法。但因为浮点乘法的运算速度要远远慢于加/减法运算,所以新算法的整体速度有所提高。

四,其他
这里列出了按通常公式计算矩阵乘法的函数,以作参考。感谢我的女朋友帮我完成了这两个函数:)值得一提的是我女朋友是学文科的,从不知道什么是矩阵,当然也没写过程序,但在我稍微指点了一下后,等我洗漱完回来,她已经写好了,经检查测试通过,把她高兴的... 

inline void Matrix22MulMatrix22_(Matrix22 c, const Matrix22& a, const Matrix22& b)
{
    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循环写出来要简洁些,但是,这样更原汁原味:)


posted on 2007-08-26 20:43 SmartPtr 阅读(5087) 评论(6)  编辑 收藏 引用

FeedBack:
# re: 矩阵快速乘法
2007-12-31 09:49 | kk
大哥,要是100阶的怎么办?  回复  更多评论
  
# re: 矩阵快速乘法
2008-05-04 16:21 | Seven
>在标准的定义算法中我们需要进行n * n * n次乘法运算,新算法中我们需要
>进行7log2n次乘法,对于最常用的4阶矩阵:

>新算法要比原算法多了24次减法运算,少了15次乘法。但因为浮点乘法的运算>速度要远远慢于加/减法运算,所以新算法的整体速度有所提高。
Hi 这是理论上的分析吧。。请问你有实际测试过这两种方法的实际执行效果吗? 因为编译器有自己的优化策略, 所以这样的改进不一定能够带来性能提高, 相反 我实际测试的结果倒是原来的乘法效率高。
请指点,谢谢!
  回复  更多评论
  
# re: 矩阵快速乘法
2012-07-12 21:40 | wx
@Seven
你可以将相同的理论应用到1000×1000的矩阵上测试,小矩阵的话误差会很大的  回复  更多评论
  
# re: 矩阵快速乘法
2013-12-20 20:05 | wu
@wx
要怎么推广到两个2^n*2^n的矩阵相乘?  回复  更多评论
  
# re: 矩阵快速乘法
2014-04-12 11:49 | yk
请问你是小学生吗,写的程序真幼稚  回复  更多评论
  
# re: 矩阵快速乘法
2015-09-08 18:17 | sdqxh
@yk
喷就不对了...  回复  更多评论
  

只有注册用户登录后才能发表评论。
【推荐】超50万行VC++源码: 大型组态工控、电力仿真CAD与GIS源码库
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理