02 Divide-and-Conquer

乘法的优化

暴力的乘法复杂度是$O(n^2)$的。

如果两个$n$位数$x,y$相乘（我们保证$n$是2的幂），那么我们按位把$x$分成$x_L,x_R$两半，$y$分成$y_L,y_R$，使得$x=x_L \cdot 2^{n/2}+x_R,y=y_L \cdot 2^{n/2}+y_R$。那么$xy=x_Ly_L \cdot 2^n + (x_Ly_R+x_Ry_L)\cdot 2^{n/2}+x_Ry_R$。加法和乘以2的幂次都是$O(n)$的，而余下的乘法可以递归的完成。设两个$n$位数的乘法所需时间是$T(n)$。那么有递推关系$T(n)=4T(n/2)+O(n)$。可以画出递归树解得$T(n)=O(n^2)$。

利用Guass's Trick，我们可以进行优化。我们可以计算出$x_ly_L,x_Ry_R$，再计算出$(x_L+x_R)(y_L+y_R)$，用第三个数减去前两个数就可以得到交叉项。这样四次乘法就被优化到了三次。此时$T(n)=3T(n/2)+O(n)$。这时解变成了$T(n)=O(n \cdot \left(\dfrac{3}{2}\right)^{\log n})=T(3^{\log_2 n})=T(n^{\log_2 3})\approx T(n^{1.59})$

主定理

从上述例子中可以发现，分治算法的复杂度分析一般式为：$T(n)=aT(n/b)+O(n^d)$。

其中假设$n$是$b$的幂。如果不是，我们可以把它“补”成$b$的幂。

我们依然画出递归树。得到（$k=\log_b n$）

$T(n)=O\left(n^d+a \cdot \dfrac{n^d}{b^d}+a^2 \cdot \dfrac{n^d}{(b^2)^d}+\cdots+a^k \cdot \dfrac{n^d}{(b^k)^d}\right)$

这是等比数列，公比为$\dfrac{a}{b^d}$。因此需要分类讨论。

当$a>b^d$时，$T(n)=O(n^d \cdot \dfrac{\frac{a^k}{b^{dk}}-1}{\frac{a}{b^d}-1})=O(n^d \cdot \dfrac{a^{\log_bn}}{b^{d\log_bn}})=O(a^{\log_b n})=O(n^{\log_b a})$

当$a=b^d$时，$T(n)=O(n^d \cdot k)=O(n^d \cdot \log_b n)=O(n^d \log n)$

当$a<b^d$时，$T(n)=O(n^d \cdot \dfrac{1-\frac{a^k}{b^{dk}}}{1-\frac{a}{b^d}})=O(n^d)$

排序

归并排序满足$T(n)=2T(n/2)+O(n)$，解得$T(n)=O(n \log n)$。我们可以证明基于比较的排序算法不可能比这更优。

考虑给一个全排列$\sigma$排序，第一次比较$\sigma_{i_1},\sigma_{j_1}$两个元素。$\sigma$分为两类，一类满足$\sigma_{i_1}<\sigma_{j_1}$，一类满足$\sigma_{i_1}>\sigma_{j_1}$；在每一类中，又可以分为$\sigma_{i_2}<\sigma_{j_2}$和$\sigma_{i_2}>\sigma_{j_2}$。如果把这种分类看成一颗二叉树的话，它的叶节点就有$n!$个。排序算法的优劣就在于怎么选择每一次的$i_k,j_k$，因为这决定了二叉树的形态。为了让排序算法最优秀，一定要让最坏的也就是深度最大的叶节点的深度尽量小。因此这应当尽量是一颗平衡树（满二叉树），此时最大深度为$\log_2 (n!)$。根据Stirling's Formula，$n! \sim \sqrt{2n\pi}\left(\dfrac{n}{e}\right)^n$，因此$O(\log n!)=O(n \log n)$。可见任何排序算法的比较次数都不可能小于$O(n \log n)$这个量级。

Medians

求中位数的复杂度是可以优于排序的复杂度的。我们每次随机选择一个数$v$，把小于$v$的和大于$v$的挑出来，转化成在新的数组里找第$k$小的问题。我们不停随机这个$k$，直到挑到$v$的排名在$1/4$到$3/4$之间。由于概率为$1/2$，随机次数的期望是$1/2+1/4+\cdots=2$。因此可以写出期望复杂度$T(n)=T(3/4n)+O(n)$，解得$T(n)=O(n)$。

在理论计算机导论中，还学过一个并不随机的算法Median of Medians，复杂度是$T(n)=T(\dfrac{n}{5})+T(\dfrac{7n}{10})+O(n)$，解得$T(n)=O(n)$。

矩阵乘法的优化

传统的矩阵乘法$c_{ij} = \sum\limits_{k=1}^{n}a_{ik}b_{kj}$复杂度为$O(n^3)$（认为整数乘法是$O(1)$的）。我们可以一个一个算，也可以从分块矩阵

$X Y=\left[\begin{array}{ll} A & B \\ C & D \end{array}\right]\left[\begin{array}{ll} E & F \\ G & H \end{array}\right]=\left[\begin{array}{ll} A E+B G & A F+B H \\ C E+D G & C F+D H \end{array}\right]$

的角度写出分治的复杂度$T(n) = 8T(n/2)+O(n^2)$，得出$T(n) = O(n^{\log_2 8}) = O(n^3)$。

Strassen's trick指出$X Y=\left[\begin{array}{cc} P_{5}+P_{4}-P_{2}+P_{6} & P_{1}+P_{2} \\ P_{3}+P_{4} & P_{1}+P_{5}-P_{3}-P_{7} \end{array}\right]$，其中

$\begin{array}{l} P_{1}=A(F-H) \quad P_{5}=(A+D)(E+H) \\ P_{2}=(A+B) H \quad P_{6}=(B-D)(G+H) \\ P_{3}=(C+D) E \quad P_{7}=(A-C)(E+F) \\ P_{4}=D(G-E) \\ \end{array}$

于是复杂度优化为$T(n)=7 T(n / 2)+O\left(n^{2}\right)=O\left(n^{\log _{2} 7}\right) \approx O\left(n^{2.81}\right)$

快速傅里叶变换(FFT)

快速傅里叶变换把多项式乘法的复杂度优化到了$O(n \log n)$。

基于乘法分配律，多项式乘法本质上完成的工作是一个称为系数的“卷积”的运算：$c_i = \sum\limits_{k=0}^{i}a_kb_{i-k}$。基于这个表达式，似乎最快的复杂度就是$O(n^2)$了。事实上我们发现，多项式不仅有“系数表达式”这一种表示方法。只要给定$n+1$个点（自变量与函数值），就可以唯一确定一个$n$次多项式。因此这些点也可以唯一对应一个多项式，称为“点值表达式”。点值表达式的乘法是容易的，只需要$O(n)$就可以完成。因此，关键是要找到系数表达式与点值表达式之间的“变换”方法。

点值表达式中的点的选取是随意的。因此算法的效率很可能取决于点的选取策略。我们发现当我们正负正负成对的选取时，会发生有趣的事。假设我们选择$\pm x_0,\pm x_1,\cdots,\pm x_{n/2}$，那么假设有多项式

$A(x) = 3+4x+6x^2+2x^3+x^4+10x^5$

那么

$A(x_i) = 3+4x_i+6x_i^2+2x_i^3+x_i^4+10x_i^5$

$A(-x_i) = 3-4x_i+6x_i^2-2x_i^3+x_i^4-10x_i^5$

正是我们选取了正负成对的点这个事实决定了得到的多项式只在奇数次系数上出现符号不同，而在其余方面两个结果是完全相同的——冗余信息对算法的优化是至关重要的！如果我们用分治的观点，把奇数次项和偶数次项看作两个子多项式$A_o(x)=4+2x+10x^2$和$A_e(x)=3+6x+x^2$，那么就会发现

$A(x) = xA_o(x^2)+A_e(x^2),A(-x)=-xA_o(x^2)+A_e(x^2)$

这样只要我们计算出了$A_o,A_e$的点值表达式，就可以得到$A$的点值表达式。可惜这样是有问题的。我们所需要的$A_o,A_e$的$n/2$个点值都是关于$x^2$的，因此符号都是正的，没法满足原来“正负成对”这一特殊的取点规则了！

而当我们引入“复数”后，我们就可以顺利解决这个问题。我们希望找到$n$个正负成对的数（依然假设$n$是$2$的幂）；将每个数平方之后，得到$n/2$个数依然正负成对；再平方，得到$n/4$个正负成对的数……直到最后变成唯一的一个数。我们就令这个数是1，那么最初的$n$个数$x_i$必须满足$(x_i)^{2^{\log_2 n}}=x_i^n=1$。根据代数基本定理，这$n$个数只能是单位根$\omega_i$。我们可以直接想象出单位根是满足我们所需要的性质的，一个单位根对应着一个角度$\dfrac{2\pi}{n}\cdot k$。那么在$[0,n)$中，起初$0,\cdots,n-1$都有，然后只剩下偶数，然后只剩下$4$的倍数，然后只剩下$8$的倍数……最后只剩下0。无论在那一轮，单位根都是均匀对称分布的，即始终是正负成对的。

这样我们就能通过分治把系数表达式转换成点值表达式了。当$n = 1$时点值表达式就是对应的系数；对于$n>1$，$T(n) = 2T(n/2)+O(n)=O(n \log n)$。

我们把这个变换（FFT）写作矩阵的形式：

$\begin{bmatrix}A(\omega_0) \\ A(\omega_1) \\ \vdots \\ A(\omega_{n-1})\end{bmatrix}=\begin{bmatrix}1 & \omega_0 & \omega_0^2 & \cdots & \omega_{0}^{n-1}\\1 & \omega_1 & \omega_1^2 & \cdots & \omega_{1}^{n-1}\\& \vdots & & \ddots & \vdots\\1 & \omega_{n-1} & \omega_{n-1}^2 & \cdots & \omega_{n-1}^{n-1}\end{bmatrix}\begin{bmatrix}a_0 \\ a_1 \\ \vdots \\ a_{n-1}\end{bmatrix}$

所以FFT其实干了这么一件事：对于一个多项式（一个客观实体），它有两种表达方式。这两种表达方式都可以表示为一个$n$维向量。这两个向量之间存在一个线性映射，而普通的计算矩阵乘法需要消耗$O(n^2)$的时间，而FFT用分治的方法只用$O(n \log n)$就完成了这个矩阵乘法。

中间那个矩阵有着特殊的结构（它被称为“范德蒙德矩阵”）——容易验证在这里它的列向量是两两正交的，因此它一定是可逆的。它的逆矩阵也容易求出，只要给每个数都取倒数（-1次方），和原来的矩阵相乘就会得到单位矩阵的$n$倍。

对于$A=Ma$，要完成点值表达式到系数表达式的转换只需要$a=M^{-1}A$。而我们知道

$M^{-1}=\dfrac{1}{n}\begin{bmatrix}1 & \omega_0^{-1} & \omega_0^{-2} & \cdots & \omega_{0}^{1-n}\\1 & \omega_1^{-1} & \omega_1^{-2} & \cdots & \omega_{1}^{1-n}\\& \vdots & & \ddots & \vdots\\1 & \omega_{n-1}^{-1} & \omega_{n-1}^{-2} & \cdots & \omega_{n-1}^{1-n}\end{bmatrix}$。而$\omega_k^{-1} = \omega_{n-k}$。所以我们只需要把原本的第$i$个单位根换成第$n-i$个，把$A$“当作”系数表达式再做一次FFT，把得到的结果除以$n$，就完成了逆变换。

总复杂度$O(n \log n)$。