您的位置：首页 > 其它

算法导论-矩阵乘法-strassen算法

2017-08-23 11:23 633 查看


目录                                                                                               


     1、矩阵相乘的朴素算法


     2、矩阵相乘的strassen算法


     3、完整测试代码c++


     4、性能分析


     5、参考资料


内容                                                                                                


     1、矩阵相乘的朴素算法 T(n) = Θ(n3)                                                    

         朴素矩阵相乘算法，思想明了，编程实现简单。时间复杂度是Θ(n^3)。伪码如下

1 for i ← 1 to n
2     do for j ← 1 to n
3         do c[i][j] ← 0
4             for k ← 1 to n
5                 do c[i][j] ← c[i][j] + a[i][k]⋅ b[k][j]



     2、矩阵相乘的strassen算法 T(n)=Θ(nlog7) =Θ (n2.81)                       

       矩阵乘法中采用分治法，第一感觉上应该能够有效的提高算法的效率。如下图所示分治法方案，以及对该算法的效率分析。有图可知，算法效率是Θ(n^3)。算法效率并没有提高。下面介绍下矩阵分治法思想：



             
鉴于上面的分治法方案无法有效提高算法的效率，要想提高算法效率，由主定理方法可知必须想办法将2中递归式中的系数8减少。Strassen提出了一种将系数减少到7的分治法方案，如下图所示。


                             效率分析如下：



                  伪码如下：

1 Strassen (N,MatrixA,MatrixB,MatrixResult)
2
3     //splitting input Matrixes, into 4 submatrices each.
4             for i  <-  0  to  N/2
5                 for j  <-  0  to  N/2
6                     A11[i][j]  <-  MatrixA[i][j];                   //a矩阵块
7                     A12[i][j]  <-  MatrixA[i][j + N / 2];           //b矩阵块
8                     A21[i][j]  <-  MatrixA[i + N / 2][j];           //c矩阵块
9                     A22[i][j]  <-  MatrixA[i + N / 2][j + N / 2];//d矩阵块
10
11                     B11[i][j]  <-  MatrixB[i][j];                    //e 矩阵块
12                     B12[i][j]  <-  MatrixB[i][j + N / 2];            //f 矩阵块
13                     B21[i][j]  <-  MatrixB[i + N / 2][j];            //g 矩阵块
14                     B22[i][j]  <-  MatrixB[i + N / 2][j + N / 2];    //h矩阵块
15             //here we calculate M1..M7 matrices .
17             //递归求M1
18             HalfSize  <-  N/2
19             AResult  <-  A11+A22
20             BResult  <-  B11+B22
21             Strassen( HalfSize, AResult, BResult, M1 );   //M1=(A11+A22)*(B11+B22)          p5=(a+d)*(e+h)
22             //递归求M2
23             AResult  <-  A21+A22
24             Strassen(HalfSize, AResult, B11, M2);          //M2=(A21+A22)B11                 p3=(c+d)*e
25             //递归求M3
26             BResult  <-  B12 - B22
27             Strassen(HalfSize, A11, BResult, M3);         //M3=A11(B12-B22)                  p1=a*(f-h)
28             //递归求M4
29             BResult  <-  B21 - B11
30             Strassen(HalfSize, A22, BResult, M4);         //M4=A22(B21-B11)                  p4=d*(g-e)
31             //递归求M5
32             AResult  <-  A11+A12
33             Strassen(HalfSize, AResult, B22, M5);         //M5=(A11+A12)B22                  p2=(a+b)*h
34             //递归求M6
35             AResult  <-  A21-A11
36             BResult  <-  B11+B12
37             Strassen( HalfSize, AResult, BResult, M6);     //M6=(A21-A11)(B11+B12)          p7=(c-a)(e+f)
38             //递归求M7
39             AResult  <-  A12-A22
40             BResult  <-  B21+B22
41             Strassen(HalfSize, AResult, BResult, M7);      //M7=(A12-A22)(B21+B22)          p6=(b-d)*(g+h)
42
43             //计算结果子矩阵
44             C11  <-  M1 + M4 - M5 + M7;
45
46             C12  <-  M3 + M5;
47
48             C21  <-  M2 + M4;
49
50             C22  <-  M1 + M3 - M2 + M6;
51             //at this point , we have calculated the c11..c22 matrices, and now we are going to
52             //put them together and make a unit matrix which would describe our resulting Matrix.
53             for i  <-  0  to  N/2
54                 for j  <-  0  to  N/2
55                     MatrixResult[i][j]                  <-  C11[i][j];
56                     MatrixResult[i][j + N / 2]          <-  C12[i][j];
57                     MatrixResult[i + N / 2][j]          <-  C21[i][j];
58                     MatrixResult[i + N / 2][j + N / 2]  <-  C22[i][j];



     3、完成测试代码                                                                                     

    Strassen.h


 Strassen.h

 

 

 

   Strassen.cpp 

1 #include <iostream>
2 #include <ctime>
3 #include <Windows.h>
4 using namespace std;
5 #include "Strassen.h"
6
7 int main()
8 {
9     Strassen_class<int> stra;//定义Strassen_class类对象
10     int MatrixSize = 0;
11
12     int** MatrixA;    //存放矩阵A
13     int** MatrixB;    //存放矩阵B
14     int** MatrixC;    //存放结果矩阵
15
16     clock_t startTime_For_Normal_Multipilication ;
17     clock_t endTime_For_Normal_Multipilication ;
18
19     clock_t startTime_For_Strassen ;
20     clock_t endTime_For_Strassen ;
21     srand(time(0));
22
23     cout<<"\n请输入矩阵大小(必须是2的幂指数值(例如:32,64,512,..): ";
24     cin>>MatrixSize;
25
26     int N = MatrixSize;//for readiblity.
27
28     //申请内存
29     MatrixA = new int *[MatrixSize];
30     MatrixB = new int *[MatrixSize];
31     MatrixC = new int *[MatrixSize];
32
33     for (int i = 0; i < MatrixSize; i++)
34     {
35         MatrixA[i] = new int [MatrixSize];
36         MatrixB[i] = new int [MatrixSize];
37         MatrixC[i] = new int [MatrixSize];
38     }
39
40     stra.FillMatrix(MatrixA,MatrixB,MatrixSize);  //矩阵赋值
41
42   //*******************conventional multiplication test
43         cout<<"朴素矩阵算法开始时钟:  "<< (startTime_For_Normal_Multipilication = clock());
44
45         stra.MUL(MatrixA,MatrixB,MatrixC,MatrixSize);//朴素矩阵相乘算法 T(n) = O(n^3)
46
47         cout<<"\n朴素矩阵算法结束时钟: "<< (endTime_For_Normal_Multipilication = clock());
48
49         cout<<"\n矩阵运算结果... \n";
50         stra.PrintMatrix(MatrixC,MatrixSize);
51
52   //*******************Strassen multiplication test
53         cout<<"\nStrassen算法开始时钟: "<< (startTime_For_Strassen = clock());
54
55         stra.Strassen( N, MatrixA, MatrixB, MatrixC ); //strassen矩阵相乘算法
56
57         cout<<"\nStrassen算法结束时钟: "<<(endTime_For_Strassen = clock());
58
59
60     cout<<"\n矩阵运算结果... \n";
61     stra.PrintMatrix(MatrixC,MatrixSize);
62
63     cout<<"矩阵大小 "<<MatrixSize;
64     cout<<"\n朴素矩阵算法: "<<(endTime_For_Normal_Multipilication - startTime_For_Normal_Multipilication)<<" Clocks.."<<(endTime_For_Normal_Multipilication - startTime_For_Normal_Multipilication)/CLOCKS_PER_SEC<<" Sec";
65     cout<<"\nStrassen算法:"<<(endTime_For_Strassen - startTime_For_Strassen)<<" Clocks.."<<(endTime_For_Strassen - startTime_For_Strassen)/CLOCKS_PER_SEC<<" Sec\n";
66     system("Pause");
67     return 0;
68
69 }


 

            输出：


     


       4、性能分析                                                                                                                 

        
矩阵大小朴素矩阵算法（秒）Strassen算法（秒）
320.0030.003
640.0040.004
1280.0210.071
2560.090.854
5120.7826.408
10248.90852.391


  可以发现：可以看到使用Strassen算法时，耗时不但没有减少，反而剧烈增多，在n=512时计算时间就无法忍受，效果没有朴素矩阵算法好。网上查阅资料，现罗列如下：

  1）采用Strassen算法作递归运算，需要创建大量的动态二维数组，其中分配堆内存空间将占用大量计算时间，从而掩盖了Strassen算法的优势

  2）于是对Strassen算法做出改进，设定一个界限。当n<界限时，使用普通法计算矩阵，而不继续分治递归。需要合理设置界限，不同环境（硬件配置）下界限不同

  3）矩阵乘法一般意义上还是选择的是朴素的方法，只有当矩阵变稠密，而且矩阵的阶数很大时，才会考虑使用Strassen算法。

分析原因：（网上总结的说法）

http://blog.csdn.net/handawnc/article/details/7987107

仔细研究后发现，采用Strassen算法作递归运算，需要创建大量的动态二维数组，其中分配堆内存空间将占用大量计算时间，从而掩盖了Strassen算法的优势。于是对Strassen算法做出改进，设定一个界限。当n<界限时，使用普通法计算矩阵，而不继续分治递归。

改进后算法优势明显，就算时间大幅下降。之后，针对不同大小的界限进行试验。在初步试验中发现，当数据规模小于1000时，下界S法的差别不大，规模大于1000以后，n取值越大，消耗时间下降。最优的界限值在32～128之间。

因为计算机每次运算时的系统环境不同（CPU占用、内存占用等），所以计算出的时间会有一定浮动。虽然这样，试验结果已经能得出结论Strassen算法比常规法优势明显。使用下界法改进后，在分治效率和动态分配内存间取舍，针对不同的数据规模稍加试验可以得到一个最优的界限。

 

http://www.cppblog.com/sosi/archive/2010/08/30/125259.html

时间复杂度就马上降下来了。。但是不要过于乐观。

从实用的观点看，Strassen算法通常不是矩阵乘法所选择的方法：

1 在Strassen算法的运行时间中，隐含的常数因子比简单的O(n^3)方法常数因子大

2 当矩阵是稀疏的时候，为稀疏矩阵设计的算法更快

3 Strassen算法不像简单方法那样子具有数值稳定性

4 在递归层次中生成的子矩阵要消耗空间。

所以矩阵乘法一般意义上还是选择的是朴素的方法，只有当矩阵变稠密，而且矩阵的阶数>20左右，才会考虑使用Strassen算法。


       5、参考资料                                                                                              

        【1】http://blog.csdn.net/xyd0512/article/details/8220506

 

        【2】http://blog.csdn.net/zhuangxiaobin/article/details/36476769

        【3】http://blog.csdn.net/handawnc/article/details/7987107

        【4】http://www.xuebuyuan.com/552410.html

        【5】http://blog.csdn.net/chenhq1991/article/details/7599824

 
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 点击这里给我发消息
标签： 
相关文章推荐
新的分享
章节导航