Strassen 矩阵相乘算法

偶尔在算法课本上面看到矩阵相乘的算法，联想到自己曾经在蓝桥杯系统上曾经做过一道矩阵相乘的题目，当时用的是普通的矩阵相乘的方法，效率极低，勉强通过编译。所以决定研究一下Strassen矩阵相乘算法，由于本人比较懒，所以就从网上找了一些相关的资料供大家参考；

下面内容均转自 https://i.cnblogs.com/EditPosts.aspx?opt=1 ,如需转载请注明出处，https://www.cnblogs.com/zhuchenglin/p/6555495.html

题目描述

    请编程实现矩阵乘法，并考虑当矩阵规模较大时的优化方法。

思路分析

    根据wikipedia上的介绍：两个矩阵的乘法仅当第一个矩阵B的列数和另一个矩阵A的行数相等时才能定义。如A是m×n矩阵和B是n×p矩阵，它们的乘积AB是一个m×p矩阵，它的一个元素其中 1 ≤ i ≤ m, 1 ≤ j ≤ p。

    

    值得一提的是，矩阵乘法满足结合律和分配率，但并不满足交换律，如下图所示的这个例子，两个矩阵交换相乘后，结果变了：

     下面咱们来具体解决这个矩阵相乘的问题。

解法一、暴力解法

    其实，通过前面的分析，我们已经很明显的看出，两个具有相同维数的矩阵相乘，其复杂度为O（n^3），参考代码如下：

//矩阵乘法，3个for循环搞定
void Mul(int** matrixA, int** matrixB, int** matrixC)
{
    for(int i = 0; i < 2; ++i)
    {
        for(int j = 0; j < 2; ++j)
        {
            matrixC[i][j] = 0;
            for(int k = 0; k < 2; ++k)
            {
                matrixC[i][j] += matrixA[i][k] * matrixB[k][j];
            }
        }
    }
}

解法二、Strassen算法

    在解法一中，我们用了3个for循环搞定矩阵乘法，但当两个矩阵的维度变得很大时，O（n^3）的时间复杂度将会变得很大，于是，我们需要找到一种更优的解法。

    一般说来，当数据量一大时，我们往往会把大的数据分割成小的数据，各个分别处理。遵此思路，如果丢给我们一个很大的两个矩阵呢，是否可以考虑分治的方法循序渐进处理各个小矩阵的相乘，因为我们知道一个矩阵是可以分成更多小的矩阵的。

    如下图，当给定一个两个二维矩阵A B时：

    这两个矩阵A B相乘时，我们发现在相乘的过程中，有8次乘法运算，4次加法运算：

    矩阵乘法的复杂度主要就是体现在相乘上，而多一两次的加法并不会让复杂度上升太多。故此，我们思考，是否可以让矩阵乘法的运算过程中乘法的运算次数减少，从而达到降低矩阵乘法的复杂度呢？答案是肯定的。

    1969年，德国的一位数学家Strassen证明O（N^3）的解法并不是矩阵乘法的最优算法，他做了一系列工作使得最终的时间复杂度降低到了O(n^2.80)。

    他是怎么做到的呢？还是用上文A B两个矩阵相乘的例子，他定义了7个变量：

    如此，Strassen算法的流程如下：

• 两个矩阵A B相乘时，将A, B, C分成相等大小的方块矩阵：

；

• 可以看出C是这么得来的：

• 现在定义7个新矩阵（读者可以思考下，这7个新矩阵是如何想到的）：

• 而最后的结果矩阵C 可以通过组合上述7个新矩阵得到：

    表面上看，Strassen算法仅仅比通用矩阵相乘算法好一点，因为通用矩阵相乘算法时间复杂度是，而Strassen算法复杂度只是。但随着n的变大，比如当n >> 100时，Strassen算法是比通用矩阵相乘算法变得更有效率。

具体实现的伪代码如下：

Strassen (N,MatrixA,MatrixB,MatrixResult)
          
    //splitting input Matrixes, into 4 submatrices each.
            for i  <-  0  to  N/2
                for j  <-  0  to  N/2
                    A11[i][j]  <-  MatrixA[i][j];                   //a矩阵块
                    A12[i][j]  <-  MatrixA[i][j + N / 2];           //b矩阵块
                    A21[i][j]  <-  MatrixA[i + N / 2][j];           //c矩阵块
                    A22[i][j]  <-  MatrixA[i + N / 2][j + N / 2];//d矩阵块
                                
                    B11[i][j]  <-  MatrixB[i][j];                    //e 矩阵块
                    B12[i][j]  <-  MatrixB[i][j + N / 2];            //f 矩阵块
                    B21[i][j]  <-  MatrixB[i + N / 2][j];            //g 矩阵块
                    B22[i][j]  <-  MatrixB[i + N / 2][j + N / 2];    //h矩阵块
            //here we calculate M1..M7 matrices .                                                                                                                       
            //递归求M1
            HalfSize  <-  N/2    
            AResult  <-  A11+A22
            BResult  <-  B11+B22                                                                     
            Strassen( HalfSize, AResult, BResult, M1 );   //M1=(A11+A22)*(B11+B22)          p5=(a+d)*(e+h)    
            //递归求M2
            AResult  <-  A21+A22    
            Strassen(HalfSize, AResult, B11, M2);          //M2=(A21+A22)B11                 p3=(c+d)*e
            //递归求M3
            BResult  <-  B12 - B22   
            Strassen(HalfSize, A11, BResult, M3);         //M3=A11(B12-B22)                  p1=a*(f-h)
            //递归求M4
            BResult  <-  B21 - B11  
            Strassen(HalfSize, A22, BResult, M4);         //M4=A22(B21-B11)                  p4=d*(g-e)
            //递归求M5
            AResult  <-  A11+A12    
            Strassen(HalfSize, AResult, B22, M5);         //M5=(A11+A12)B22                  p2=(a+b)*h
            //递归求M6
            AResult  <-  A21-A11
            BResult  <-  B11+B12      
            Strassen( HalfSize, AResult, BResult, M6);     //M6=(A21-A11)(B11+B12)          p7=(c-a)(e+f)
            //递归求M7
            AResult  <-  A12-A22
            BResult  <-  B21+B22      
            Strassen(HalfSize, AResult, BResult, M7);      //M7=(A12-A22)(B21+B22)          p6=(b-d)*(g+h)

            //计算结果子矩阵
            C11  <-  M1 + M4 - M5 + M7;

            C12  <-  M3 + M5;

            C21  <-  M2 + M4;

            C22  <-  M1 + M3 - M2 + M6;
            //at this point , we have calculated the c11..c22 matrices, and now we are going to
            //put them together and make a unit matrix which would describe our resulting Matrix.
            for i  <-  0  to  N/2
                for j  <-  0  to  N/2
                    MatrixResult[i][j]                  <-  C11[i][j];
                    MatrixResult[i][j + N / 2]          <-  C12[i][j];
                    MatrixResult[i + N / 2][j]          <-  C21[i][j];
                    MatrixResult[i + N / 2][j + N / 2]  <-  C22[i][j];