《数据结构与算法分析：C语言描述》复习——第十章“算法设计技巧”——Strassen矩阵乘法

2014.07.06 18:52

简介：

　　给定两个大小相同的方阵A和B，我们要计算AXB。方阵的大小是2的整次幂，比如2^k。对于这种特殊大小的方阵乘法，Strassen算法能够带来一定程度的加速，矩阵越大，加速效果越明显。

描述：

　　这个例子其实也是分治法的典型算法，通过矩阵分块进行乘法，然后合并结果。

　　不过我觉得这个例子，更突出的特点其实是“数学”二字。我想作者应该是想告诉我们，数学和计算机的关联，就在于用数学工具来思考问题可以把某些你想不清楚的东西一步一步推导清楚。不论是代数、几何、微积分（数学系的如果路过请勿鄙视，我们作为码农学好这些已经算是恪守本分了）、数值分析、科学计算等等，都不可能让你觉得很实用，但只要某一个瞬间能从这些交过学费的课程里得到些灵感，那就不算白学了。

　　根据矩阵乘法的基本定义，两个nXn的方阵相乘需要O(n^3)级别的算术运算。而Strassen算法的时间复杂度大概O(n^2.81)，这样看起来不很显著的数量级差异需要较大的矩阵才能显示出性能的差异。因此Strassen算法的应用并不十分广泛，但思想的确是很有用的（至少Volker Strassen教授告诉我们学好代数是非常重要的）。初学这个算法时的感受，就和学习斐波那契数的对数级算法一样——数学要好好学，因为你学的那点数学实在很浅。

　　矩阵的分块我们在高中的代数里就学过，如果一个方阵的尺寸恰好是2的整次幂，那么我们可以一直进行四等分，直到尺寸为1。

　　于是矩阵A可以等分为四块儿：

　　　　┏━━┳━━┓

　　　　┃A_1,1┃A_1,2┃

　　　　┣━━╋━━┫

　　　　┃A_2,1┃A_2,2┃

　　　　┗━━┻━━┛

　　矩阵B也一样：

　　　　┏━━┳━━┓

　　　　┃B_1,1┃B_1,2┃

　　　　┣━━╋━━┫

　　　　┃B_2,1┃B_2,2┃

　　　　┗━━┻━━┛

 　　本来2X2X2=8，如果对于这些子矩阵需要进行8次乘法，那就没有效率上的提升。但Herr Strassen设计了一种神奇的计算方法，使用7次子矩阵乘法，请看下面的公式：

　　　　1. M₁=(A₁₂-A₂₂)(B₂₁+B₂₂)

　　　　2. M₂=(A₁₁+A₂₂)(B₁₁+B₂₂)

　　　　3. M₃=(A₁₁-A₂₁)(B₁₁+B₁₂)

　　　　4. M₄=(A₁₁+A₁₂)B₂₂

　　　　5. M₅=A₁₁(B₁₂-B₂₂)

　　　　6. M₆=A₂₂(B₂₁-B₁₁)

　　　　7. M₇=(A₂₁+A₂₂)B₁₁

　　这7个小矩阵M₁~M₇都是n/2大小的方阵。通过下面四个公式又可以将它们组合成最终结果C：

　　　　1. C11=M1+M2-M4+M6

　　　　2. C12=M4+M5

　　　　3. C21=M6+M7

　　　　4. C22=M2-M3+M5-M7

　　　　┏━━┳━━┓

　　　　┃C_1,1┃C_1,2┃

　　　　┣━━╋━━┫

　　　　┃C_2,1┃C_2,2┃

　　　　┗━━┻━━┛

　　七个式子是乘法，属于子问题的分解；四个式子是加减法，复杂度是严格O(n^2)的。两者的复杂度在下面会进行推导。

　　接下来开启高中模式，用数列知识来证明O(n^2.81)的复杂度（高考之前的我们应该是对数列问题最敏感的）：

　　　　T(n)=7T(n/2)+O(n^2)，前半部分表示7次子矩阵乘法，后半部分是4次加减法。

　　　　接下来把7T(n/2)这部分逐层拆成49T(n/4)、343T(n/8)...

　　　　T(n)=7(7T(n/4)+O(n^2/4))+O(n^2)

　　　　T(n)=49T(n/4)+O(n^2/4)+O(n^2)

　　　　...

　　　　T(n)=7^k * T(n/(2^k))+O((n/(2^k))^2)+...+O((n/2)^2)+O(n^2)

　　　　此处的k等于多少呢？等于log₂n。为什么呢？因为2^k=n，这样拆了k次以后就不能再拆了，而k不一定要是整数。

　　　　再来看看后面的等比数列部分O(n^2)+O(n^2/4)+...+O(n^2/(4^k))。这部分的和S=O(n^2) * (4/3-1/(3 * 4^k))

　　　　由于2^k=n，于是4^k=n^2。那么S=O(n^2) * (4/3-1/(3 * n^2))=O(n^2)-O(1)=O(n^2)。

　　　　所以后面的等比数列总体的时间复杂度是O(n^2)级别的。

　　　　再来看看前面的部分7^k * T(n/(2^k))

　　　　T(n/(2^k))=O(1)，表示这个公式展开了k次以后就不能再展开了。

　　　　7^k=7^(log₂n)=7^(log₇n/log₇2)=7^(log₇n * log₂7)=(7^log₇n)^log₂7=n^log₂7

　　　　又因为log₂7≈2.81，于是O(n^2.81)+O(n^2)=O(n^2.81)。

　　由此看来Strassen算法的O(n^2.81)的时间复杂度是由于化8次子问题为7次才得到的。

　　不过Srassen算法的局限性也是不少的：矩阵的尺寸限制的太严格，而且处理子问题的并行能力不太强。因此这个算法更多是出现在教材和ACM题库里，实际应用则没那么多。

实现：

// My implementation for Strassen Algorithm. Matrix size is limited to powers of 2 only.
#include <iostream>
#include <vector>
using namespace std;

void getSubmatrix(const vector<vector<int> > &m, vector<vector<int> > &sub, int part)
{
    int n = (int)m.size();
    int i, j;
    int n1 = n / 2;
    
    int top = part / 2 * n1;
    int left = part % 2 * n1;
    for (i = 0; i < n1; ++i) {
        for (j = 0; j < n1; ++j) {
            sub[i][j] = m[top + i][left + j];
        }
    }
}

void addMatrix(const vector<vector<int> > &a, const vector<vector<int> > &b, 
    vector<vector<int> > &c, int n)
{
    int i, j;
    
    for (i = 0; i < n; ++i) {
        for (j = 0; j < n; ++j) {
            c[i][j] = a[i][j] + b[i][j];
        }
    }
}

void subtractMatrix(const vector<vector<int> > &a, const vector<vector<int> > &b, 
    vector<vector<int> > &c, int n)
{
    int i, j;
    
    for (i = 0; i < n; ++i) {
        for (j = 0; j < n; ++j) {
            c[i][j] = a[i][j] - b[i][j];
        }
    }
}

void setSubmatrix(const vector<vector<int> > &sub, vector<vector<int> > &m, int part)
{
    int n = (int)m.size();
    int i, j;
    int n1 = n / 2;
    
    int top = part / 2 * n1;
    int left = part % 2 * n1;
    for (i = 0; i < n1; ++i) {
        for (j = 0; j < n1; ++j) {
            m[top + i][left + j] = sub[i][j];
        }
    }
}

void matrixMultiplicationRecursive(const vector<vector<int> > &a, 
    const vector<vector<int> > &b, vector<vector<int> > &c, int n)
{
    if (n == 1) {
        c[0][0] = a[0][0] * b[0][0];
        return;
    }
    
    int i;
    int n1 = n / 2;
    
    vector<vector<int> > aa[4];
    vector<vector<int> > bb[4];
    
    for (i = 0; i < 4; ++i) {
        aa[i].resize(n1, vector<int>(n1));
        bb[i].resize(n1, vector<int>(n1));
    }
    
    for (i = 0; i < 4; ++i) {
        getSubmatrix(a, aa[i], i);
        getSubmatrix(b, bb[i], i);
    }
    
    vector<vector<int> > x, y;
    vector<vector<int> > m[7];
    
    x.resize(n1, vector<int>(n1));
    y.resize(n1, vector<int>(n1));
    for (i = 0; i < 7; ++i) {
        m[i].resize(n1, vector<int>(n1));
    }
    
    subtractMatrix(aa[1], aa[3], x, n1);
    addMatrix(bb[2], bb[3], y, n1);
    matrixMultiplicationRecursive(x, y, m[0], n1);
    
    addMatrix(aa[0], aa[3], x, n1);
    addMatrix(bb[0], bb[3], y, n1);
    matrixMultiplicationRecursive(x, y, m[1], n1);
    
    subtractMatrix(aa[0], aa[2], x, n1);
    addMatrix(bb[0], bb[1], y, n1);
    matrixMultiplicationRecursive(x, y, m[2], n1);
    
    addMatrix(aa[0], aa[1], x, n1);
    matrixMultiplicationRecursive(x, bb[3], m[3], n1);
    
    subtractMatrix(bb[1], bb[3], y, n1);
    matrixMultiplicationRecursive(aa[0], y, m[4], n1);
    
    subtractMatrix(bb[2], bb[0], y, n1);
    matrixMultiplicationRecursive(aa[3], y, m[5], n1);
    
    addMatrix(aa[2], aa[3], x, n1);
    matrixMultiplicationRecursive(x, bb[0], m[6], n1);
    
    addMatrix(m[0], m[1], x, n1);
    subtractMatrix(x, m[3], x, n1);
    addMatrix(x, m[5], x, n1);
    setSubmatrix(x, c, 0);
    
    addMatrix(m[3], m[4], x, n1);
    setSubmatrix(x, c, 1);
    
    addMatrix(m[5], m[6], x, n1);
    setSubmatrix(x, c, 2);
    
    subtractMatrix(m[1], m[2], x, n1);
    addMatrix(x, m[4], x, n1);
    subtractMatrix(x, m[6], x, n1);
    setSubmatrix(x, c, 3);

    for (i = 0; i < 4; ++i) {
        aa[i].clear();
        bb[i].clear();
    }
    for (i = 0; i < 7; ++i) {
        m[i].clear();
    }
    x.clear();
    y.clear();
}

void matrixMultiplication(const vector<vector<int> > &a, 
    const vector<vector<int> > &b, vector<vector<int> > &c)
{
    int n = (int)a.size();
    
    matrixMultiplicationRecursive(a, b, c, n);
}

int main()
{
    int n;
    int i, j;
    vector<vector<int> > a, b, c;
    
    while (cin >> n && n > 0) {
        a.resize(n, vector<int>(n));
        b.resize(n, vector<int>(n));
        c.resize(n, vector<int>(n));
        
        for (i = 0; i < n; ++i) {
            for (j = 0; j < n; ++j) {
                cin >> a[i][j];
            }
        }

        for (i = 0; i < n; ++i) {
            for (j = 0; j < n; ++j) {
                cin >> b[i][j];
            }
        }
        
        matrixMultiplication(a, b, c);

        for (i = 0; i < n; ++i) {
            for (j = 0; j < n; ++j) {
                cout << c[i][j] << ' ';
            }
            cout << endl;
        }
        cout << endl;

        a.clear();
        b.clear();
        c.clear();
    }
    
    return 0;
}