算法重拾之路——strassen矩阵乘法

***************************************转载请注明出处：http://blog.csdn.net/lttree********************************************

第一章：分治与递归

STRASSEN矩阵乘法

算法描述：

        矩阵乘法是线性代数中最常见的问题之一，它在数值计算中有广泛的应用。设A 和 B 是两个n × n矩阵，它们的乘积AB同样是一个n×n矩阵。A和B的乘积矩阵C 中的元素cij定义为：

按照这个定义来看，计算A 与 B 矩阵乘法，至少需要 n 次乘法 与 n-1 次加法，所以可以知道，求矩阵C乘法的时间为O(n^3)

算法分析：

        strassen矩阵乘法采用类似于 大数乘法 中的分治技术，将计算2个n阶矩阵乘积所需的时间改进到O（n^log7) ≈ O（n^2.81）

这个算法有个前提条件，必须是 两个n×n阶矩阵相乘，而且n必须为2的幂。

这样我们每次都可以把大矩阵分割成四个小矩阵：

由上面图可知：

C11 = A11B11 + A12B21

C12 = A11B12 + A12B22

C21 = A21B11 + A22B21

C22 = A21B12 + A22B22

        如果 n=2 ，两个2×2阶矩阵乘法需要 8次乘法 和 4次加法。当子矩阵的阶数大于2时，为求两个子矩阵的乘积，可以继续将矩阵分块，直到子矩阵阶数降为2，这就是分治降阶的递归算法。

        按照这个算法，计算2个n阶矩阵的乘积转化为计算8个n/2阶矩阵的乘积和4个 n/2阶 矩阵的加法。 而 2个 n/2阶矩阵的加法显然可以再O（n^2）时间内完成，由此可知上述算法的时间耗费T（n）：

        T（n）   =         O（1）                                  当n=2

                        =         8T（n/2）+O（n^2）           当n>2

        Ok,解这个递归方程，发现 T（n） = O（n^3）。并没有减少！

        这是因为，这个方法并没有减少矩阵的乘法次数，要改进这个算法的时间，必须减少乘法的次数，Strassen就提出了 对于2阶矩阵的乘积方法，只用了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)

        做完这七次乘法，通过一些加减就可以得到，最后矩阵的值：

C11 = M5 + M4 - M2 + M6

C12 = M1 + M2

C21 = M3 + M4

C22 = M5 + M1 - M3 - M7

        这样做以后，它的时间T（n）为：

        T（n）    =        O（1）                                  当n=2

                        =         7T（n/2）+O（n^2）           当n>2

        解得 T（n） = O（n^log7） ≈ O（n^2.81）

程序代码：

<span style="font-family:Comic Sans MS;">const int N = 8;     //Define the size of the Matrix

// 输入矩阵
template<typename T>
void input(int n, T p[][N]) {

     for(int i=0; i<n; i++)
        for(int j=0; j<n; j++)
            cin>>p[i][j];
}

// 输出矩阵
template<typename T>
void output(int n, T C[][N]) {

    for(int i=0; i<n; i++) {
        for(int j=0; j<n; j++) {
            cout<<C[i][j]<<" ";
        }
        cout<<endl;
    }
}


// 普通的矩阵乘法
template<typename T>
void Matrix_Mul(T A[][N], T B[][N], T C[][N]) {
     for(int i=0; i<2; i++) {
        for(int j=0; j<2; j++) {
            C[i][j] = 0;
            for(int k=0; k<2; k++) {
               C[i][j] = C[i][j] + A[i][k]*B[k][j];
            }
        }
    }
}


// 矩阵加法  C=A+B
template <typename T>
void Matrix_Add(int n, T A[][N], T B[][N], T C[][N]) {
     for(int i=0; i<n; i++) {
        for(int j=0; j<n; j++) {
            C[i][j] = A[i][j] + B[i][j];
        }
    }
}


// 矩阵减法  C=A-B
template <typename T>
void Matrix_Sub(int n, T A[][N], T B[][N], T C[][N]) {
    for(int i=0; i<n; i++) {
        for(int j=0; j<n; j++) {
            C[i][j] = A[i][j] - B[i][j];
        }
    }
}

// strassen矩阵乘法
template <typename T>
void Strassen(int n, T A[][N], T B[][N], T C[][N]) {
    // 将第一个矩阵A分割成四个小矩阵
    T A11[N][N], A12[N][N], A21[N][N], A22[N][N];
    // 将第二个矩阵B分割成四个小矩阵
    T B11[N][N], B12[N][N], B21[N][N], B22[N][N];
    // M1~M7 先存起来，为后面最终C矩阵做准备
    T M1[N][N], M2[N][N], M3[N][N], M4[N][N], M5[N][N], M6[N][N], M7[N][N];
    // 四个小矩阵C  最后组合成大矩阵
    T C11[N][N], C12[N][N], C21[N][N], C22[N][N];
    // 中间变量
    T tempA[N][N], tempB[N][N];

    if(n == 2) {
        Matrix_Mul(A, B, C);
    }
    else {
        //将矩阵A和B分成阶数相同的四个子矩阵，即分治思想。
        for(int i=0; i<n/2; i++)
        {
            for(int j=0; j<n/2; j++)
            {
                A11[i][j] = A[i][j];
                A12[i][j] = A[i][j+n/2];
                A21[i][j] = A[i+n/2][j];
                A22[i][j] = A[i+n/2][j+n/2];

                B11[i][j] = B[i][j];
                B12[i][j] = B[i][j+n/2];
                B21[i][j] = B[i+n/2][j];
                B22[i][j] = B[i+n/2][j+n/2];
            }
        }

        /* ****** 分别计算 M1~M7 ****** */

        //Calculate M1 = (A0 + A3) × (B0 + B3)
        Matrix_Add(n/2, A11, A22, tempA);
        Matrix_Add(n/2, B11, B22, tempB);
        Strassen(n/2, AA, BB, M1);

        //Calculate M2 = (A2 + A3) × B0
        Matrix_Add(n/2, A21, A22, tempA);
        Strassen(n/2, tempA, B11, M2);

        //Calculate M3 = A0 × (B1 - B3)
        Matrix_Sub(n/2, B12, B22, tempB);
        Strassen(n/2, A11, tempB, M3);

        //Calculate M4 = A3 × (B2 - B0)
        Matrix_Sub(n/2, B21, B11, tempB);
        Strassen(n/2, A22, tempB, M4);

        //Calculate M5 = (A0 + A1) × B3
        Matrix_Add(n/2, A11, A12, tempA);
        Strassen(n/2, tempA, B22, M5);

        //Calculate M6 = (A2 - A0) × (B0 + B1)
        Matrix_Sub(n/2, A21, A11, tempA);
        Matrix_Add(n/2, B11, B12, tempB);
        Strassen(n/2, tempA, tempB, M6);

        //Calculate M7 = (A1 - A3) × (B2 + B3)
        Matrix_Sub(n/2, A12, A22, tempA);
        Matrix_Add(n/2, B21, B22, tempB);
        Strassen(n/2, tempA, tempB, M7);

        //Calculate C1 = M1 + M4 - M5 + M7
        Matrix_Add(n/2, M1, M4, tempA);
        Matrix_Sub(n/2, M7, M5, tempB);
        Matrix_Add(n/2, tempA, tempB, C11);

        //Calculate C2 = M3 + M5
        Matrix_Add(n/2, M3, M5, C12);

        //Calculate C3 = M2 + M4
        Matrix_Add(n/2, M2, M4, C21);

        //Calculate C4 = M1 - M2 + M3 + M6
        Matrix_Sub(n/2, M1, M2, tempA);
        Matrix_Add(n/2, M3, M6, tempB);
        Matrix_Add(n/2, tempA, tempB, C22);

        //将四个子矩阵组合回大矩阵 C[][N]
        for(int i=0; i<n/2; i++)
        {
            for(int j=0; j<n/2; j++)
            {
                C[i][j] = C11[i][j];
                C[i][j+n/2] = C12[i][j];
                C[i+n/2][j] = C21[i][j];
                C[i+n/2][j+n/2] = C22[i][j];
            }
        }
    }
}</span>

***************************************转载请注明出处：http://blog.csdn.net/lttree********************************************