《算法导论》——矩阵乘法的Strassen算法

﹏ヽ暗。殇╰゛Y 2021-11-23 02:48 621阅读 0赞

前言：

很多朋友看到我写的《算法导论》系列，可能会觉得云里雾里，不知所云。这里我再次说明，本系列博文时配合《算法导论》一书，给出该书涉及的算法的c++实现。请结合《算法导论》一书阅读该系列博文。我这里有该书的电子版，有需要的朋友可以留言。

正题：

今天讨论的算法是矩阵乘法的Strassen算法，该算法的精髓在于减少n/2矩阵\*n/2矩阵的次数。首先，作一些写该算法的基础工作：

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

/*
        * 矩阵的加法运算
        */
        void Add(int** matrixA, int** matrixB, int** matrixResult,int length)
        {
            for(int i = 0; i < length; i++) {
                for(int j = 0; j < length; j++) {
                    matrixResult[i][j] = matrixA[i][j] + matrixB[i][j];
                }
            }
        }
    
        /*
        * 矩阵的减法运算
        */
        void Sub(int** matrixA, int** matrixB, int** matrixResult,int length)
        {
            for(int i = 0; i < length; i++) {
                for(int j = 0; j < length; j++) {
                    matrixResult[i][j] = matrixA[i][j] - matrixB[i][j];
                }
            }
        }
    
        /*
        * 矩阵乘法
        */
        void Mul(int** matrixA, int** matrixB, int** matrixResult){
            for(int i = 0; i < 2; ++i) {
                for(int j = 0; j < 2; ++j) {
                    matrixResult[i][j] = 0;
                    for(int k = 0; k < 2; ++k) {
                        matrixResult[i][j] += matrixA[i][k] * matrixB[k][j];
                    }
                }
            }
        }

接着进入核心部分：

void Strassen(int** matrixA, int** matrixB, int** matrixResult,int length)
        {
            int halfLength=length/2;
            int** a11=new int*[halfLength];
            int** a12=new int*[halfLength];
            int** a21=new int*[halfLength];
            int** a22=new int*[halfLength];
    
            int** b11=new int*[halfLength];
            int** b12=new int*[halfLength];
            int** b21=new int*[halfLength];
            int** b22=new int*[halfLength];
    
            int** s1=new int*[halfLength];
            int** s2=new int*[halfLength];
            int** s3=new int*[halfLength];
            int** s4=new int*[halfLength];
            int** s5=new int*[halfLength];
            int** s6=new int*[halfLength];
            int** s7=new int*[halfLength];
    
            int** matrixResult11=new int*[halfLength];
            int** matrixResult12=new int*[halfLength];
            int** matrixResult21=new int*[halfLength];
            int** matrixResult22=new int*[halfLength];        
    
            int** temp=new int*[halfLength];
            int** temp1=new int*[halfLength];
            if(halfLength==1){
                Mul(matrixA, matrixB, matrixResult);
            }else{
                //首先将矩阵A，B 分为4块
                for(int i = 0; i < halfLength; i++) {    
                    a11[i]=new int[halfLength];
                    a12[i]=new int[halfLength];
                    a21[i]=new int[halfLength];
                    a22[i]=new int[halfLength];
    
                    b11[i]=new int[halfLength];
                    b12[i]=new int[halfLength];
                    b21[i]=new int[halfLength];
                    b22[i]=new int[halfLength];
    
                    s1[i]=new int[halfLength];
                    s2[i]=new int[halfLength];
                    s3[i]=new int[halfLength];
                    s4[i]=new int[halfLength];
                    s5[i]=new int[halfLength];
                    s6[i]=new int[halfLength];
                    s7[i]=new int[halfLength];
    
                    matrixResult11[i]=new int[halfLength];
                    matrixResult12[i]=new int[halfLength];
                    matrixResult21[i]=new int[halfLength];
                    matrixResult22[i]=new int[halfLength];
    
                    temp[i]=new int[halfLength];
                    temp1[i]=new int[halfLength];
                    for(int j = 0; j < halfLength; j++) {
                        a11[i][j]=matrixA[i][j];
                        a12[i][j]=matrixA[i][j+halfLength];
                        a21[i][j]=matrixA[i+halfLength][j];
                        a22[i][j]=matrixA[i+halfLength][j+halfLength];
                        b11[i][j]=matrixB[i][j];
                        b12[i][j]=matrixB[i][j+halfLength];
                        b21[i][j]=matrixB[i+halfLength][j];
                        b22[i][j]=matrixB[i+halfLength][j+halfLength];
                    }                
                }
    
                //计算s1
                Sub(b12, b22, temp,halfLength);
                Strassen(a11, temp, s1,halfLength);
                //计算s2
                Add(a11, a12, temp,halfLength);
                Strassen(temp, b22, s2,halfLength);
                //计算s3
                Add(a21, a22, temp,halfLength);
                Strassen(temp, b11, s3,halfLength);
                //计算s4
                Sub(b21, b11, temp,halfLength);
                Strassen(a22, temp, s4,halfLength);
                //计算s5
                Add(a11, a22, temp1,halfLength);
                Add(b11, b22, temp,halfLength);
                Strassen(temp1, temp, s5,halfLength);
                //计算s6
                Sub(a12, a22, temp1,halfLength);
                Add(b21, b22, temp,halfLength);
                Strassen(temp1, temp, s6,halfLength);
                //计算s7
                Sub(a11, a21, temp1,halfLength);
                Add(b11, b12, temp,halfLength);
                Strassen(temp1, temp, s7,halfLength);
    
                //计算matrixResult11
                Add(s5, s4, temp1,halfLength);
                Sub(temp1, s2, temp,halfLength);
                Add(temp, s6, matrixResult11,halfLength);
                //计算matrixResult12
                Add(s1, s2, matrixResult12,halfLength);
                //计算matrixResult21
                Add(s3, s4, matrixResult21,halfLength);
                //计算matrixResult22
                Add(s5, s1, temp1,halfLength);
                Sub(temp1, s3, temp,halfLength);
                Sub(temp, s7, matrixResult22,halfLength);
    
                //结果送回matrixResult中
                for(int i = 0; i < halfLength; i++) {
                    for(int j = 0; j < halfLength; j++) {
                        matrixResult[i][j]=matrixResult11[i][j];
                        matrixResult[i][j+halfLength]=matrixResult12[i][j];
                        matrixResult[i+halfLength][j]=matrixResult21[i][j];
                        matrixResult[i+halfLength][j+halfLength]=matrixResult22[i][j];
                    }
                    delete(a11[i]);
                    delete(a12[i]);
                    delete(a21[i]);
                    delete(a22[i]);
    
                    delete(b11[i]);
                    delete(b12[i]);
                    delete(b21[i]);
                    delete(b22[i]);
    
                    delete(s1[i]);
                    delete(s2[i]);
                    delete(s3[i]);
                    delete(s4[i]);
                    delete(s5[i]);
                    delete(s6[i]);
                    delete(s7[i]);
    
                    delete(matrixResult11[i]);
                    delete(matrixResult12[i]);
                    delete(matrixResult21[i]);
                    delete(matrixResult22[i]);
    
                    delete(temp[i]);
                    delete(temp1[i]);
                }
                delete(a11);
                delete(a12);
                delete(a21);
                delete(a22);
    
                delete(b11);
                delete(b12);
                delete(b21);
                delete(b22);
    
                delete(s1);
                delete(s2);
                delete(s3);
                delete(s4);
                delete(s5);
                delete(s6);
                delete(s7);
    
                delete(matrixResult11);
                delete(matrixResult12);
                delete(matrixResult21);
                delete(matrixResult22);
    
                delete(temp);
                delete(temp1);
            }
        }

该算法看着或许有些冗长，几乎一半都在进行动态指针的初始化和删除。利用该算法计算矩阵乘的时间复杂度为θ（n^lg7）。

测试一下吧：

#include "stdafx.h"
    #include <iostream>
    #include "SquareMatrix.h"
    
    using namespace std;
    using namespace dksl;
    
    
    //STRASSEN矩阵乘法算法
    
    
    const int N=8; //常量N用来定义矩阵的大小
    int _tmain(int argc, _TCHAR* argv[])
    {
        int **a=new int*[4];
        int **b=new int*[4];
        int **c=new int*[4];
        for(int i=0;i<4;i++)
            {
                a[i]=new int[4];
                b[i]=new int[4];
                c[i]=new int[4];
            for(int j=0;j<4;j++)
            {
                a[i][j]=1;
                b[i][j]=2;
            }
        }
        Strassen(a,b,c,4);
        for(int i=0;i<4;i++)
        {
            for(int j=0;j<4;j++)
                cout<<c[i][j]<<" ";
            cout<<endl;
        }
        system("PAUSE");
        return 0;
    }

![22224016-33b507f52d7c48dbba8911ea4d27c886.png][]

转载于:https://www.cnblogs.com/DKSL/p/3150338.html

[ContractedBlock.gif]: https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif
[ExpandedBlockStart.gif]: /images/20211122/cc0149546b4b492f9d4614affe78b6d5.png
[22224016-33b507f52d7c48dbba8911ea4d27c886.png]: /images/20211122/9c273343249e4f8c89ba515a411e7b34.png