Massive Algorithms: Matrix Multiply

Matrix Multiply - 矩阵乘法

两个矩阵 A，B：A 的列数等于 B 的行数，则A、B可以相乘。即，如果 A = （a_ij）是一个m * n的矩阵，B =（b_jk）是一个 n * p 的矩阵，则它们的乘积 C = AB 是 m * p 矩阵 C = （c_ik）。

学过计算机图形学，会发现，不管是二维图形的平移，旋转，缩放，三维图形的取景变换，投影变换等都是通过矩阵乘法来实现，例如，二维点P（x，y）平移（tx，ty）后得到 P’(x’，y’)，可以通过矩阵计算：

$\begin{pmatrix} x'\\ y'\\ 1 \end{pmatrix}= \begin{pmatrix} 1 & 0 & tx\\ 0 & 1 & ty\\ 0 & 0 & 1 \end{pmatrix}\begin{pmatrix} x\\ y\\ 1 \end{pmatrix}$

求矩阵相乘的通用公式为：

$AB_{ij} = \sum_{r = 1}^{n}a_{ir}b_{rj} =a_{i1}b_{1j} + a_{i2}b_{2j}+……+a_{in}+b_{nj}$

还有一个常见用法，求斐波那契数列：

$\inline \begin{pmatrix} 1 & 1\\ 1& 0 \end{pmatrix} * \begin{pmatrix} Fi\\ Fi-1 \end{pmatrix} = \begin{pmatrix} Fi + Fi-1\\ Fi \end{pmatrix} = \begin{pmatrix} Fi+1\\ Fi \end{pmatrix}$

$\inline \begin{vmatrix} 1 &1 \\ 1 & 0 \end{vmatrix} * \begin{vmatrix} 8\\ 5 \end{vmatrix}=\begin{vmatrix} 13\\ 8 \end{vmatrix}$

所以当求 Fn 时可以用求幂来快速求出，而求幂也是建立在矩阵乘法的基础上：

$\begin{vmatrix} 1 &1 \\ 1&0 \end{vmatrix} * \begin{vmatrix} 1 &1 \\ 1&0 \end{vmatrix} * \begin{vmatrix} 1 &1 \\ 1&0 \end{vmatrix} * ……* \begin{vmatrix} 1 &1 \\ 1&0 \end{vmatrix} * \begin{vmatrix} F1\\ F0 \end{vmatrix} = \begin{vmatrix} Fn\\ Fn-1 \end{vmatrix}$

$\inline A = \begin{vmatrix} F1\\ F0 \end{vmatrix}B =\begin{vmatrix} 1 &1 \\ 1& 0 \end{vmatrix}C = \begin{vmatrix} Fn\\ Fn-1 \end{vmatrix}$

$C = B^{n-1}*A$

矩阵相乘有两种方法，普通的矩阵乘法（Matrix Multiply）和Strassen算法。

最普通的矩阵乘法

Matrix operator * (Matrix a, Matrix b) {

    Matrix ret; 

    ret.init(a.n, b.m);

    for (int i = 0; i < a.n; i++) {

        for (int k = 0; k < a.m; k++) if (a.mat[i][k]) {

            for (int j = 0; j < b.m; j++) if (b.mat[k][j]) {

                ret.mat[i][j] = ret.mat[i][j] + a.mat[i][k] * b.mat[k][j];  

                if (ret.mat[i][j] >= mod) {

                    ret.mat[i][j] %= mod;

                }//if

            }//for(j)

        }//for(k)

    }//for(i)

    return ret;

}//乘法

矩阵加法只需简单的将两个矩阵的元素相加：

Matrix operator + (Matrix a, Matrix b) {

    Matrix ret; 

    ret.init(a.n, a.m);

    for (int i = 0; i < a.n; i++) {

        for (int j = 0; j < a.m; j++) {

            ret.mat[i][j] = a.mat[i][j] + b.mat[i][j];

            if (ret.mat[i][j] >= mod) {

                ret.mat[i][j] %= mod;

}

}

}

    return ret;

}//加法

矩阵求幂

Matrix operator ^ (Matrix a, int b) {

    Matrix ret = a,  tmp = a;

    ret.init_e();

    for ( ; b; b >>= 1) {

        if (b & 1) {

            ret = ret * tmp;

}

        tmp = tmp * tmp;

}

    return ret;

}//

幂求和，即求S = A + A² + A³ +… + A^k ，同样的二分思想，但是利用递归，可以很快求出：

//递归幂求和

//用二分法求矩阵和,递归实现  

Matrix Power_Sum1(Matrix a, int b) {

    Matrix ret = a;

    ret.init_e();

    if (b == 1) {

        return a;

    } else if (b & 1) {

        return (a ^ b) + Power_Sum1(a, b - 1);

    } else {

        return Power_Sum1(a, b >> 1) * ((a ^ (b >> 1)) + ret);

}

}

//非递归幂求和

Matrix Power_Sum2(Matrix a, int b) {

    int k = 0 ,ss[32];

    Matrix tp1, tp2, ret;

    tp1 = tp2 = ret = a;

    ret.init_e();

    while (b) { 

        ss[k++] = b & 1;

        b >>= 1;

}

    for (int i = k - 2; i >= 0; i--) {

        tp1 = tp1 * (tp2 + ret);

        tp2 = tp2 * tp2;

        if (ss[i]) {

            tp2 = tp2 * a;

            tp1 = tp1 + tp2;

}

}

    return tp1;

}

二、Strassen算法

Strassen算法核心思想是分治，是一种递归算法，运行时间为O（n^lg7） = O（n^2.81），当 n 很大时，优化很明显，在普通的矩阵乘法中，C = A * B，按照：

$C[i,j] = \sum_{k = 1}^{n}A[i, k] * B[k, j]$

每计算一个元素C[i，j]，需要做 n 个乘法和 n – 1 次加法。因此，求出矩阵 C 的 n²个元素所需的计算时间为0（n^3）。Strassen算法的分治体现在：假设 n 是 2 的幂，将将矩阵A，B和C中每一矩阵都分块成为 4 个大小相等的子矩阵，每个子矩阵都是n / 2 × n / 2的方阵。由此可将方程C = AB重写为:

$\begin{vmatrix} C_{0} &C_{1}\\ C_{2} &C_{3} \end{vmatrix} = \begin{vmatrix} A_{0} &A_{1}\\ A_{2} &A_{3} \end{vmatrix} \times \begin{vmatrix} B_{0} &B_{1}\\ B_{2} &B_{3} \end{vmatrix}$

由此可得

$C_{0} = A_{0} * B_{0} + A_{1} * B_{2}\\~~~~~ C_{1} = A_{0} * B_{1} + A_{1} * B_{3}\\~~~~~ C_{2} = A_{2} * B_{0} + A_{3} * B_{2}\\~~~~~ C_{3} = A_{2} * B_{1} + A_{3} * B_{3}$

可以看出，进行了8次乘法，4次加法，当子矩阵的阶大于2时，为求2个子矩阵的积，可以继续将子矩阵分块，直到子矩阵的阶降为2，利用这个简单的分治策略，最后可以得出：T（n） = 8T（n / 2) + O（n²），但是这个式子的解任然为T（n）= O（n³），和普通的方法效率一样，没有任何提高，原因是上边的四个式子并没有减少矩阵乘法的次数（乘法极其耗费时间，学过底层二进制计算的，必然了解，而加减操作非常轻松），所以改进算法的关键是，减少乘法次数。

Strassen算法的高效之处，就在于，它成功的减少了乘法次数，将8次乘法，减少到7次。不要小看这减少的一次，每递归计算一次，效率就可以提高1 / 8，比如一个大的矩阵递归5次后，（7 / 8）⁵ = 0.5129，效率提升一倍。不过，这只是理论值，实际情况中，有隐含开销，并不是最常用算法，《算法导论》中给出四条理由：1）隐含的常数因子比简单的O（n3）方法中的常数因子要大。2）矩阵是稀疏矩阵时，为稀疏矩阵设计的方法更快。还有两点已经被缓解，可以不考虑。所以，稠密矩阵上的快速矩阵乘法实现一般采用Strassen算法。

M1 = (A0 + A3) × (B0 + B3)

M2 = (A2 + A3) × B0

M3 = A0 × (B1 – B3)

M4 = A3 × (B2 – B0)

M5 = (A0 + A1) × B3

M6 = (A2 – A0) × (B0 + B1)

M7 = (A1 – A3) × (B2 + B3)

C0 = M1 + M4 – M5 + M7

C1 = M3 + M5

C2 = M2 + M4

C3 = M1 – M2 + M3 + M6

求解M1，……，M7总共7次乘法，其他都是加法和减法，比如将C0扩展开后，最后结果是，C0 = A0 * B0 + A1 * B2，《算法导论》里有一句奇怪的话：“现在我们还不清楚Strassen当时是如何找出算法正常运行的关键——子矩阵乘积”，一次乘法的消失过程真的这么吊诡？

Please read full article from Matrix Multiply - 矩阵乘法