矩阵方幂与矩阵乘法, 哪个更快?

写在前面

我们先提出一个问题.

Problem 1.1
有两个 n 行 n 列的矩阵 A, B. 定义两种运算, 第一种是计算 A 和 A 相乘的值, 第二种是计算 A 和 B 相乘的值. 求比较两种运算的复杂度.

即便我们不知道计算机是怎么实现矩阵乘法, 我们仍可以看到, 第一种运算是第二种运算的一个子问题, 因此复杂度不会高于后者. 

于是我们尝试计算一下矩阵的复杂度. 假定 A = [a[i][j]] 和 B = [b[i][j]] 的乘积是 C = [c[i][j]], 那么c[i][j]的值如下可求. 直观地讲, 计算出一个c[i][j]需要把 n 对 a[i][k] 和 b[k][j] 相乘, 因此需要至少 O(n) 的复杂度; 计算出整个 C 需要至少 O(n^3) 的复杂度.

$$ c_{ij} = \sum\limits_{k=1}^{n} a_{ik}b_{kj}$$

事实上我们可以采用一些巧妙的办法对矩阵乘法进行提速. 而传统运算中也有一些很经典的提速算法, 最典型的是快速幂. 我们将介绍一下运算操作的原理和一些提速的方法, 进而剖析一下矩阵复杂度的问题, 进而分析problem 1.1中蕴含的奥秘.

运算操作(arithmetic operation)

在前文中, 我们提到过, 补码可以忽略符号位直接进行竖式加法和减法, 结果不需要做任何转换. 我们继续这个话题. 计算机在实现加法的时候和竖式加法非常相似. 由于两个加数都是二进制数, 所以实现起来会省很多事. 对于整数 a 和 b 的加法, 我们只需要申请一个位 (bit) 的空间, 记录是否进位, 之后进行从低位到高位的按位加法就可以了. 那么, 从位运算的角度来讲, 因为一个数值为 n 的数有 log(n) 位, 所以两个数相加的复杂度是 O(log(n)).

计算实现整数的乘法, 要比竖式乘法简单很多. 假定我们区分被乘数和乘数(竖式第一行与第二行). 由于乘数每一位上的取值只有0和1, 所以我们只是把被乘数在乘数为1的位上抄写下来, 在最后把这些抄写下来的数加起来即可. 下面的竖式实现了二进制下 5 和 7 的乘法(实际的编码必然超过6位, 我们仅考虑后6位). 

\begin{array}{cccccccc} & 0 & 0 & 0 & 1 & 0 & 1 & (5)\\ * & 0 & 0 & 0 & 1 & 1 & 1 & (7)\\ - & - & - & - & - & - & - &\\ & & & & 1 & 0 & 1 &\\ & & & 1 & 0 & 1 & &\\ + & & 1 & 0 & 1 & & &\\ - & - & - & - & - & - & - & \\ & 1 & 0 & 0 & 0 & 1 & 1 & (35) \end{array}

因此乘法的本质只是 log(n) 次的移位和加法. 相比于加法的时间代价, 移位操作的代价可以忽略, 那么两个数相乘的复杂度是 O(log(n) * log(n)).

这里的移位操作相当于前文中提到的位运算中的左移. 

在位运算中, 左移一个数会在最右边添上一个新的"0". 这样一来用左移可以轻松实现一个数乘以 2 的操作, 从而在乘法操作中与竖式等价. 

而右移操作会在一个数的最左边添上一个"0"吗? 不是的. 负数的右移如果在左边添上一个0, 则会变成正数. 因此, 右移操作总是在最左边添上与符号位相同的一个位(把最左边的一位复制一遍). 这样的话, 右移操作就可以实现一个"除以 2 之后取下整"的操作了. 例如, -4(111...1100)和-5(111...1101)的右移结果都是-2(111..1110); 2(000...0010)和3(000...0011)的右移结果都是1(000...0001).

注意到竖式乘法的 log(n) 次加法操作仍然需要开一个位(bit)来储存借位的值. 因此, 在底层中实现乘以, 或者除以2的幂的时候, 往往用一次左移, 右移操作代替, 这样不仅化繁为简, 而且可以避免开额外的空间.

快速幂(exponentiation by squaring)

我们先讲幂(exponentiation), 再讲平方(squaring). 从算法优化的角度来讲, 后者甚至比前者难一些.

先假定两个数的乘法的时间代价是O(f(n)). 根据前文, f(n) = log(n) * log(n). 现在假定我们需要计算一个整数 a 的 m 次方的值. 朴素的算法是做 m-1 次乘法, 这个算法的复杂度是O(m*f(n)). 然而这个做法的优化空间很大.

考虑一次乘法之后, 我们得到了a^2. 之后我们再做一次乘法, 得到 a^3, 再做一次得到 a^4. 事实上, 从 a^2 经过一次乘法就可以得到 a^4, 即 a^2 * a^2 = a^4. 进一步说, 从 a 得到 a^16 只需要4次乘法, 从 a 得到 a^(2^k) 只需要 k 次乘法. 对于求 a 的 m 次方, 我们只需要做 log(m) 次乘法, 从而复杂度降到了 O(log(m) * f(n)).

那么如果 m 不是 2 的幂呢? 我们举一个例子, 现在 m 等于 20, 而我们在做 4 次乘法之后得到了 a^16. 注意到, 我们在做这 4 次乘法的时候, 可以把中间结果, 也就是 a^2, a^4, a^8 存起来. 这样的话, 我们只要做第 5 次乘法, a^16 * a^4, 就可以得到 a^20 了. 这个例子里面, 20 被拆成了两个 2 的幂的和, 16和4. 

因此我们需要把 m 拆成 2 的幂. 事实上, 由于 m 是二进制存储的, 我们可以直接从它的二进制位上读出它的 2 的幂的组合. 例如 19(000...010011) 有三个1, 对应16, 2, 1, 于是19 = 16+2+1. 我们在算 a 的 19 次方的时候, 只需要先用 4 次乘法分别得到a, a^2, a^4, a^8, a^16, 并且从右往左标在 19 的二进制位的下面; 再把二进制位为1的几个结果相乘起来即可.

\begin{array} {cccccccc} 0 & 1 & 0 & 0 & 1 & 1 & | & (19) \\ & a^{16} & a ^8 & a^4 & a^2 & a & | & a^{19} = a^{16} + a^2 + a \end{array}

快速幂的算法如下. 对于 a 的 m 次方, 先用一个 O(log(m)) 长度的列表, 在对 a 进行自平方的同时, 从右向左记录 a, a^2, a^4, a^8... 的值. 然后再用这张表对应 m 的二进制位, 把对应位为1的列表中的结果乘起来, 就得到了 a^m 的值. 在这个算法中, 我们进行了两轮乘法, 每轮最多做了 O(log(m)) 次, 因此总时间代价为 O(log(m) * f(n)).

在一类经典题目中, 经常会要求我们计算一个 a 的 b 次方模 m 的值. 比如"今天是星期天, 求123的45678次方天之后是星期几". 

这类问题我们应该对 a 的 b 次方进行快速幂运算, 而且在 m 不是大整数的情况下, 可以考虑每次都把中间结果对 m 取模. 这个算法叫做"快速幂取模", 适用于 a^b 是大整数, 甚至 a 或 b 是大整数的情况.

另外, 如果 m 也是大整数的话, 这时一次对 m 取模的时间复杂度很高, 需要在取模上进行优化.

平方(squaring)

我们可以看到, 快速幂有一半的时间都花在了平方上. 事实上, 我们需要指出, 对于一个整数 x , 对 x 的平方的运算是有很大优化空间的. (注意到 problem 1.1 关注了矩阵平方和矩阵乘法的复杂度比较) 我们在这里讲一个最基础的优化, 对后面的矩阵运算会有所启发.

根据前面的结论, 从位的角度来看, 对两个整数乘法的时间复杂度是O(log(n) * log(n)). 换言之, 如果两个整数的二进制码长度是 k , 那么对它们做乘法的时间复杂度就是O(k^2). 对它们做加法的时间代价却要小很多 (O(k)). 我们的思路是尽量减少使用乘法的代价, 尽可能地用加法来代替乘法操作. 

考虑平方和公式. 我们发现, 如果计算c = (a+b)^2, 只需要对 a 和 b 做共计三次乘法操作, 便得到 c 的值.

$$ c = (a+b)^2 = a^2 + 2ab + b^2 $$ $$ c = C_1 + 2C_2 + C_3 \begin{cases} C_1 = a^2 \\ C_2 = a*b \\ C_3 = b^2 \end{cases}$$

求c的值只需要一次乘法c = (a+b)*(a+b), 为什么要多此一举? 我们回到把整数 x 平方的例子. 假定 x 有 2k 位, 后 k 位所表示的值记做 b; x-b, 也就是前k位所表示的值乘以2^k, 记做 a. 注意 x 的有效位数是 2k, 对 x 直接平方的乘法时间代价是 O(4k^2), 而 a 和 b 的有效位数是 k , 对它们做三次乘法的时间代价是 O(3k^2). 当 k 充分大的时候, 加法的时间代价可以忽略.

注意, 对两个有效位个数是k的数做乘法, 尽管相乘之后的结果会有 2k 个有效位, 但是我们依然可以达到O(k^2)的时间代价. 读者可以通过重现竖式操作来明白这个问题.

这个方法在时间上的改进远不止3/4. 我们把求 c 的问题拆成了求 C1, C2, C3 三个子问题. 注意到求C1和C3也是在求整数的平方, 那么它们也可以做这样一个分割. 假定f(k)是对一个长度为k的整数求平方的算法的时间复杂度, 再假定C2这一类乘法不可优化, 我们便可以得到

$$ f(k) = 2f(\frac{k}{2}) + O(k^2) $$

根据Master Theorem(本质是数列的一类递推公式的一个总结性公式), 我们得到f(k) = 2*k^2. 解决一个长度为2k的整数平方, 利用这个基础的优化, 就可以省掉一半的乘法运算. 

那么, 以上的优化方法改进了快速幂的时间复杂度了吗? 答案是否定的. 快速幂有两个步骤, 第二个步骤的乘法不是平方乘法, 最坏的情况下也需要进行log(m) 次. 因此这个算法并没有对快速幂的第二步有任何的优化. 这里讲这个算法是为了引出Master Theorem和分治法的概念, 从而对矩阵乘法的时间复杂度分析有一定的启发.

事实上, 不仅仅是整数平方, 对整数乘法的优化也可以在log(n)级下进行. 这个上面的优化算法有很多. 举一个最有名的例子, Schönhage–Strassen Algorithm提供了一种基于快速傅立叶变换的算法来进行优化的方法. (感谢lxg的提示)

分治法(divide and conquer)

上面一个对整数平方的优化方法属于分治法. 

分治法是把母问题化归成子问题求解的方法. 由于子问题可以继续化归, 我们总是可以得到一个递推公式来计算它的时间代价. 根据Master Theorem的内容, 只要我们在化归的过程中产生额外的复杂度低于一个条件, 我们便总可以花更小的时间来解决这个问题. Master Theorem还可以直接求解分治法的时间复杂度.

Master Theorem 又译主定理, 是一类数列递推问题的总结性公式. 它处理的递推公式如下. 根据a, b 和c(n)的关系来判定f(n)的取值. 由于在很多算法问题上采用对母问题分而治之的方法求解, 因此得到的递推公式都属于Master Theorem定义的范畴. 

$$ f(n) = a f(\frac{n}{b}) + c(n) $$

分治法在算法领域内的应用非常广泛. 我们再这里考虑一个分块矩阵的例子. 假定我们有两个n*n的矩阵A 和 B 相乘, 为了简化分析过程, 我们假设n是2的幂(n = 2^k) . 它们相乘的结果是 C . 朴素的算法可以得到, 这个矩阵乘法需要O(n^3)的时间复杂度, 因为C中的每一个元素都需要n次乘法运算所得. 

Problem 0.0
假定两个2*2的矩阵相乘需要且只需要做 x 次乘法和 y 次加法, x和y是常数. 试分析是否可以用 x 和 y 表示两个 n*n 的矩阵乘法的时间复杂度.

解这个问题便需要对矩阵做分治法. 矩阵拥有分块的特性, 对于A*B = C, 如果我们把它们分别分块, 则得到这样的式子. 

$$ \begin{bmatrix} A_1 & A_2 \\ A_3 & A_4 \end{bmatrix} * \begin{bmatrix} B_1 & B_2 \\ B_3 & B_4 \end{bmatrix} = \begin{bmatrix} C_1 & C_2 \\ C_3 & C_4 \end{bmatrix}$$

假定我们得到了A1, A2, A3, A4和B1, B2, B3, B4, 我们是否可以通过 x 次乘法和 y 次加法来完成对 C1, C2, C3, C4的运算呢? 

如果我们可以通过 x 次乘法和 y 次加法来从A1, A2, A3, A4 和 B1, B2, B3, B4 来求出 C1, C2, C3, C4的话, 根据Master Theorem求解, 我们对整个C矩阵的求解使用的乘法个数和加法个数(n_Multiplications, n_Additions), 以及n*n矩阵相乘的算法复杂度Complexity, 就等于下式.

$$ (n_M, n_A) = ( x^k, \frac{y}{x-1}(x^k-1) ), k = \log _2 (n) $$ $$ Complexity = x^k = x^{\log _2 (n) } = n^{\log _2 (x)} $$

矩阵乘法(matrix multiplication)

我们现在需要回答刚刚提出的问题, 即能否通过 x 次乘法和 y 次加法来完成对C1, C2, C3, C4的运算. 如果可以的话, 我们就可以得到O(n^k)的复杂度, k = log(x). 一个很显然的结论, 就是 x 肯定不会超过 8, 因为朴素的算法允许 8 次乘法搞定2*2矩阵相乘. 如果我们可以得到 x 小于 8 的话, 这个复杂度是低于O(n^3)的. 

Strassen Algorithm 提供了一种 7 次乘法操作和 18 次加法操作来解决 2*2 的矩阵相乘方法. 按照这个方法, 我们可以用分治法求解矩阵相乘, 从而拿到一个不错的时间复杂度, O(n^k), k 等于 log_2(7), 约等于2.807. 

我们回到problem 1.1. 矩阵相乘可以在O(n^2.807)的时间内搞定, 那么矩阵平方会不会更简单一些呢? 从2*2的矩阵乘法可以看到, 对于2*2的矩阵平方, 我们是可以通过5次乘法来实现的. 

$$ \begin{bmatrix} a & b \\ c & d \end{bmatrix}^2= \begin{bmatrix} a^2+bc & ab+bd \\ ac+cd & bc+d^2 \end{bmatrix} $$ $$ = \begin{bmatrix} a*a & b*(a+d) \\ c*(a+d) & d*d \end{bmatrix} + \begin{bmatrix} b*c & 0 \\ 0 & b*c \end{bmatrix}$$

这是不是意味着我们可以用O(n^k), k = log_2(5)的时间来解决矩阵的平方问题呢? 是不是就意味着矩阵的平方可以算得比矩阵的乘法更快一些?

答案是否定的. 注意上面5次乘法的推导公式. ab+bd可以化成b*(a+d). 但是当我们拿着大的分块矩阵去做矩阵乘法的时候, 这个地方被完全改变了. AB+BD不等于B(A+D), 因为矩阵的乘法不满足交换律! 换而言之, 我们可以用5次乘法来实现2*2的矩阵平方, 却无法把它推广到矩阵平方的分治法上面.

我们回到Stassen Algorithm. 这个算法在推导过程中没有使用过乘法的交换律, 因而可以用在分治法的策略上面, 从而得到一个k = log_2(7)的算法. 

Answer 0.0
在用 x 次乘法和 y 次加法求解2*2矩阵乘法的过程中, 如果用到了矩阵的交换律, 那么我们不能用这个问题来对矩阵乘法进行分治; 否则我们可以用 x 和 y 来表示矩阵乘法的复杂度. Complexity = O(n^k), k = log_2(x). 一般地, 如果可以不用交换律地在z*z的矩阵上用x次乘法来完成运算, 那么矩阵乘法的复杂度表示如下, 其中对数函数的底数是z:

$$ Complexity = n^{\log _z (x)} $$

矩阵的乘法不满足交换律, 这是在推导矩阵乘法时需要注意的一个重要的地方.

矩阵乘法的复杂度下界是O(n^2). 在Strassen algorithm(1969)之后, 最有名的算法属于Coppersmith–Winograd algorithm(1990)算法, 这个算法把矩阵运算的复杂度优化到了O(n^2.376)上. 二十年间, CW算法一直是矩阵运算中最快的算法, 持续到2010年. 2014年, 有学者提出一种方法, 用O(n^2.3728639)的时间代价解决了矩阵乘法.

哪个更快(which one is faster)

本节会对Problem 1.1给出解答.

前文说过, 我们可以用 5 次乘法搞定 2*2 的矩阵平方, 却无法推广到更大的矩阵上. 事实上, 这个方法无法在提高矩阵方幂对矩阵乘法的复杂度. 我们倾向于矩阵方幂和矩阵乘法一样快.

考虑"平方(squaring)" 一节中提到的平方和公式. 如果我们把平方和公式中的数换成矩阵的话, 就会有这样一个公式.

$$ C = (A + B)^2 = A^2 + AB + BA + B^2 $$ $$ AB+BA = C^2 - A^2 - B^2 $$

由于矩阵乘法不满足交换律, 所以我们得到了AB+BA 而不是2AB. 那么AB+BA就等于三个矩阵平方进行加减得到的结果. 换而言之, 假如矩阵平方的时间代价是O(g(n)), 那么AB+BA的时间代价不高于O(g(n)). 

我们需要的是把O(g(n))和矩阵乘法的时间代价O(f(n))联系上, 现在我们已经证明AB+BA的时间代价不高于O(g(n)), 下一步只需要找到O(f(n))和AB+BA的关系. 首先, AB+BA的时间代价肯定也不高于O(f(n)), 因为我们总可以做两个矩阵乘法拿到这个式子. 

其次, 观察下面的构造. 对于任意一对矩阵 X, Y 相乘, 我们都可以构造另一对矩阵 A, B, 使得AB+BA足以求出XY的结果. 

$$A=\begin{bmatrix} X&0 \\ 0&0 \end{bmatrix}, B=\begin{bmatrix} 0&Y \\ 0&0 \end{bmatrix} $$ $$AB+BA = \begin{bmatrix} 0 & XY \\ 0 & 0 \end{bmatrix} + \begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix} = \begin{bmatrix} 0 & XY \\ 0 & 0 \end{bmatrix}$$

注意, A, B 相对于 X, Y来说, 把矩阵的大小扩大了常数倍. 由于时间代价也增长了一个常数, 因而不影响算法的复杂度.

从 X, Y 到A, B的构造时间复杂度也是O(n^2). 由于我们提到的矩阵平方, 矩阵乘法, 以及AB+BA的时间代价的下限都是O(n^2), 因此我们可以得出O(f(n))的时间代价不高于计算AB+BA的时间代价. 

这个重要的结论间接地说明了计算矩阵平方不会比计算矩阵乘法简单.

根据这个巧妙的构造, 我们已经可以对problem 1.1做出解答了.

Answer 1.1
假定矩阵乘法的时间代价是O(f(n)), 矩阵平方的时间代价是O(g(n)). 由于平方是乘法的子问题, 我们有g(n) 不高于 f(n).

另假设算出AB+BA的时间代价是O(h). 一方面, AB+BA可以由三次矩阵平方和加减法运算得到, 因此 h 不高于 g(n); 另一方面, 用上面的构造方法可以用AB+BA这个形式计算任意两个矩阵的乘法结果, 因此 f(n) 不高于 h. 两方面可得, f(n) 不高于 g(n).

综上所述, 矩阵乘法和矩阵平方的复杂度是相等的.

完.