分治法(1)：Karatsuba算法与大O表示法

本文将以大整数乘法的Karatsuba算法为例，介绍分治算法（Divide and Conquer），并介绍时间复杂度分析的大O表示法。

算法（algorithm） 一词来自于9世纪波斯数学家al-Khwarizmi。因为他撰写了一本《印度计算书》，其开篇为"Dixit Algorismi"，意为“al-Khwarizmi如是说”，而这个词逐渐演化到了今天我们熟知的algorithm一词。

那就让我们从一种数学运算——整数乘法开始学习算法。

Karatsuba算法

当我们计算简单的乘法，如$44 \times 34$时，用简单的竖式手算非常轻松（对于计算机更加轻松）。但当乘数的位数$n$不断增大，乘法的运算量也随之加大。那么有没有方法能减少大整数乘法中的运算量？

从大整数乘法开始

让我们先看看$n$位大整数乘法的运算量。如果我们以1位数的加法或乘法为一次运算，那么$n$位的大整数乘法需要多少次运算？由小学数学可知，大约有

• $n^2$次1位数乘法
• $n^2$次1位数加法（进位）
• $n$次$2n$位数加法

因此，我们总共要进行大约$5n^2$次一位数运算，可以写作$\bm{O(n^2)}$。

注意，我们在这里用大O表示法省略了常数项，是因为我们只关心$\bm{n}$很大的情况。

那么，我们是否有办法做得更好呢？

分治法的首次尝试

让我们尝试分治法！分治法，顾名思义，就是把原来的问题不断划分为规模更小的问题，直至最简单的情况。

算法设计

分治法的设计无非考虑两个问题：

1. 怎么分（如何将原问题划分为更小的子问题）
2. 怎么治（如何将子问题的解合并为原问题的解）

再这里，我们可以将乘数划分为高位和低位（怎么分），两两递归地计算乘积，再按权重相加（怎么治）。举个例子，我们用分治法计算$1234 \times 5678$：

$$\begin{aligned} 1234 \times 5678 =& \, (12 \times 100 + 34) \times (56 \times 100 + 78) \\ =& \, (12 \times 56) \cdot 10000 + (12 \times 78 + 34 \times 56) \cdot 100 + 34 \times 78 \end{aligned}$$

或者更一般地，

$$\begin{aligned} x \times y =& \, (a \cdot 10^{\frac{n}{2}} + b) \times (c \cdot 10^{\frac{n}{2}} + d) \\ =& \, ac \cdot 10^n + (ad + bc) \cdot 10^{\frac{n}{2}} + bd \end{aligned}$$

这样，我们就把一个$\bm{n}$位数乘法的问题，分解为了4个$\bm{\dfrac{n}{2}}$位数乘法的问题。而这4个子问题，将会按同样的方法继续递归调用，直到问题规模缩小到1位数乘法。

注意，我们在这里假定$n$是2的整数次幂。事实上，如果$n$不是2的整数幂，我们可以在高位补0。

算法分析

那么，我们的算法有没有比原来做得更好？也就是说，我们做了多少次1位数运算？在这里，我们用分治法中最经典的树形图进行分析。

这里，我们将每个大问题分出4个小问题，因此每一层的问题个数都是上一层的4倍。而由于我们每次能将问题规模缩小一半，即递归树的层数为$\log_2 n$，最终一共会分出$4^{\log_2 n}$个基本问题，即递归树有$4^{\log_2 n}$个叶结点。

在分治法中，我们的时间代价主要有两方面：

1. 分的代价：解决最终分出的所有最小规模的问题（即基本问题）
2. 治的代价：每次将小问题的解合并为大问题的解

在树形图中，“分的代价”即指解决最下面一层叶结点的代价，“治的代价”即指下层叶结点合并到上层父节点的代价。

由于解决每个基本问题的需要完成1次1位数乘法，我们一共需要完成

$$4^{\log_2 n} \cdot 1 = n^2$$

次1位数乘法！也就是说，我们需要做超过$\bm{n^2}$次1位数运算！

显然，我们的朴素分治算法并没有让计算量变得更小。这似乎也是符合直觉的，因为我们只是把要计算的乘法交换了计算顺序，但基本计算的总量并没有变。以$1234 \times 5678$为例，在几次递归之后，仍归结为每一位之间的两两相乘，即16个基本问题，并没有减少1位数运算的数量。

那么，有没有办法真正减少计算量呢？

Karatsuba算法

想要减少计算，我们必须看看我们是否做了冗余的计算，或者说有没有计算可以简化。

算法设计

前面已经推导过，

$$x \times y = ac \cdot 10^n + (ad + bc) \cdot 10^{\frac{n}{2}} + bd$$

也就是说，为了求出$xy$，我们需要知道$ac$、$ad + bc$、$bd$。而我们刚才计算了哪些呢？$ac$、$ad$、$bc$、$bd$。换句话说，其实我们不需要知道$\bm{ac}$和$\bm{bd}$分别是多少，只需要知道$\bm{ac + bd}$的值就够了。

那么有没有可能，可以不通过$\bm{ad}$和$\bm{bc}$来求得$\bm{ac + bd}$的值呢？相信你已经想到了，我们可以在算出$ac$和$bd$之后，再算出

$$z = (a + b)(c + d)$$

进而用

$$ad + bc = (a + b)(c + d) - ac - bd = z - ac - bd$$

来求得$ad + bc$，也就是

$$\begin{aligned} x \times y =& \, (a \cdot 10^{\frac{n}{2}} + b) \times (c \cdot 10^{\frac{n}{2}} + d) \\ =& \, ac \cdot 10^n + (ad + bc) \cdot 10^{\frac{n}{2}} + bd \\ =& \, ac \cdot 10^n + (z - ac - bd) \cdot 10^{\frac{n}{2}} + bd \end{aligned}$$

这就是Karatsuba算法。

原先要做4次乘法，现在只需要3次。利用Karatsuba算法，我们将一个$\bm{n}$位数乘法问题，分解为了3个$\bm{\dfrac{n}{2}}$位乘法问题，也就是每层分出的子问题数比之前少了1个。那么，这是否带来了时间复杂度的降低呢？

算法分析

我们先计算要解决基本问题的数量，即进行了多少次1位数乘法。与刚才类似，我们不难得到递归树的层数为$\log_2 n$，且每层的结点数为上一层的三倍，因此基本问题的数量（即叶结点个数）为

$$3^{\log_2 n} = n^{\log_2 3} \approx n^{1.6}$$

也就是说，我们只需要完成$\bm{n^{1.6}}$次1位数乘法！

可以验证，需要完成的1位数加法的数量也为$O(n^{1.6})$。

这样我们竟然得到了一个时间复杂度为$\bm{O(n^{1.6})}$的整数乘法算法！

事实上，我们可以做得更好！比如利用快速傅立叶变化（FFT），我们可以将两个$n$位整数相乘的时间复杂度降低到$O(n \log n)$。这一技巧我们将在后续小节介绍。

拓展到矩阵乘法（Strassen算法）

这一分治法的思想可以应用到矩阵乘法中。设有两个$n$阶方阵$X$和$Y$，要计算$Z = XY$。按常规计算方法，我们需要计算$Z$中的$n^2$个元素，每个元素要做$n$次乘法，因此时间复杂度为$\bm{O(n^3)}$。那么，我们可以用分治法做得更好吗？

不难想到，我们可以对X和Y分块。设

$$X = \begin{bmatrix} A & B \\ C & D \end{bmatrix}, \quad Y = \begin{bmatrix} E & F \\ G & H \end{bmatrix}$$

那么

$$Z = \begin{bmatrix} A & B \\ C & D \end{bmatrix} \begin{bmatrix} E & F \\ G & H \end{bmatrix} = \begin{bmatrix} AE+BG & AF+BH \\ CE+DG & CF+DH \end{bmatrix}$$

这样，我们将一个规模为$n$的矩阵乘法，分解为了8个规模为$\dfrac{n}{2}$的矩阵乘法。这能不能让计算次数减少呢？

相信你不会再次被骗，因为这和之前整数乘法中的朴素分治法很类似。事实上，我们仍需要做

$$8 ^ {\log_2 n} = n^3$$

次整数乘法，因此并不能降低运算量。

解决的思路也和之前类似。我们并不需要求出8个矩阵乘法的结果，只需要知道4个矩阵和即可，因此可以尝试将8次乘法减少为7次，从而让分解出的问题变少。

Strassen算法的具体计算不需要掌握，了解分治法的思想即可。

这个算法能让矩阵乘法的时间复杂度降低到$O(n^{2.86})$。

时间复杂度

为了定义时间复杂度，我们需要定义单位操作。例如，在整数乘法中，我们将1位数乘法定义为单位操作；在矩阵乘法中，我们将1次整数乘法定义为单位操作。那么在一般算法中，单位操作应该是什么呢？

计算模型

在这里，我们简单地介绍Word RAM 模型。我们设计算机的字长为$w$，那么可以认为内存中能存$2^w$个长度为$w$的数（因为内存地址的位数也为$w$）。我们定义对一个字进行一次操作为一次单位操作。那么，我们在算法分析中应该取多大的$w$？

以排序算法为例，假设我们的输入为$a_1, a_2, \ldots, a_n$。

• 地址空间限制：为了能存下所有的$n$个输入，我们要求$n \le 2^w$。
• 数据存储限制：为了能存下每一个输入的数字，我们要求$a_i \le 2^w$。

但同时，我们不希望$w$过大，因为这可能导致某些操作被过度优化。比如在大整数乘法中，如果取$w = n$，那么其时间复杂度可以被认为是$O(1)$，这显然不符合我们的期望。

因此，我们通常取

$$w = O(\max\{\log n, \log a_i\})$$

来确保时间复杂度的分析基于合理的字长假设。

大 O 表示法

为了了解一个算法的性能好坏，我们需要衡量一个算法的运行时间。那么如何衡量更加合理呢？

时间复杂度

我们知道，算法在一台计算机上的运行时间与算法本身、输入$\sigma(n)$的规模$n$有关，还与计算机处理器的性能有关。为了屏蔽处理器性能等因素对算法时间分析的影响，我们常用基本操作的数量来表征算法的运行时间$t(\sigma(n))$。

在算法分析中，我们常考虑最坏的输入情况，并定义时间复杂度

$$T(n) = \max_{\sigma(n)} t(\sigma(n))$$

即在输入规模为$\bm{n}$时需要的最多基本操作数量。

然而，单纯用$T(n)$来分析仍然有些复杂。例如，时间复杂度为$T_1(n) = 5n^2 + 3n + 1$ 和 $T_2(n) = n^2 + 2n + 100$ 的两个算法哪个性能更好？我们需要对$T(n)$作进一步化简。

大 O 表示法

我们在前面已经使用过大$O$表示法，比如我们可以将$T(n) = 5n^2 + 3n + 1$写作$T(n) = O(n^2)$。这表示当$\bm{n}$很大时，$\bm{T(n)}$最多和$\bm{n^2}$同阶，或者用更严格的语言表述，即当$n$足够大时，存在常数$C > 0$，使得$T(n) \le C n^2$。

注意，$T(n) = O(n^2)$可以等于$n^2 + 100n + 10000$，可以等于$1000n^2 + 10n$，也可以等于$n + 100$，因为大$O$表示法只规定了时间复杂度的上界。

我们可以对大$O$表示法进行严格定义：

• 大$O$表示法：$T(n) = O(g(n))$的定义为

$$\exist C, \, n_0, \, \text{s.t.} \, \forall n > n_0, \, T(n) \le C \cdot g(n)$$

• 大$\Omega$表示法：$T(n) = \Omega(g(n))$的定义为

$$\exist C, \, n_0, \, \text{s.t.} \, \forall n > n_0, \, T(n) \ge C \cdot g(n)$$

• 大$\Theta$表示法：$T(n) = \Theta(g(n))$的定义为

$$T(n) = O(g(n)) \quad \text{ and } \quad T(n) = \Omega(g(n))$$

可以看到，大$O$表示法刻画了$T(n)$的上界，大$\Omega$表示法给出了$T(n)$的下界，而大$\Theta$表示法给出了$T(n)$精确的阶，即$T(n)$与$g(n)$同阶。简单来说，一般我们对$T(n)$保留最高阶项即可得到大$O$表示法的$g(n)$。

其实，我们还有一种小$o$表示法：

• 小$o$表示法：$T(n) = o(g(n))$的定义为

$$\exist C, \, n_0, \, \text{s.t.} \, \forall n > n_0, \, T(n) < C \cdot g(n)$$

• 小$\omega$表示法：$T(n) = \omega(g(n))$的定义为

$$\exist C, \, n_0, \, \text{s.t.} \, \forall n > n_0, \, T(n) > C \cdot g(n)$$

不难发现，小$o$和小$\omega$表示法分别给出了$T(n)$的严格上界和严格下界。

以Karatsuba算法为例

我们可以以Karatsuba算法为例，练习时间复杂度的分析。设输入为n位整数的Karatsuba算法的时间复杂度为$T(n)$，那么我们可以其伪代码出发，逐行分析其时间复杂度。

由此，可以得到$T(n)$的递推公式

$$T(n) = 3T\left(\frac{n}{2}\right) + O(n)$$

这里，我们尝试通过逐步展开的方式求出$T(n)$。

$$\begin{aligned} T(n) =& \, 3T\left(\frac{n}{2}\right) + O(n) \\ =& \, 3^2 T\left(\frac{n}{2^2}\right) + O(1.5n + n) \\ =& \, 3^3 T\left(\frac{n}{2^3}\right) + O(1.5^2n + 1.5n + n) \\ =& \, \cdots \\ =& \, 3^{\log_2 n} T(1) + O(1.5^{\log_2 n - 1} + \cdots + 1.5^2n + 1.5n + n) \\ =& \, n^{\log_2 3} \cdot O(1) + O(n^{\log_2 3}) \\ =& O(n^{\log_2 3}) \end{aligned}$$

这和我们之前分析出的结果一致。

注：本文中所有图片均来自张宇昊老师的课程PPT。