CSC4120 W02 分治的应用

I. Applications of D&C

Convex Hull 凸包

给定平面上的 n 个点 S={(x_i , y_i) | i=1, 2…, n}

Convex Hull CH(S) 是包含所有这些点的最小凸多边形，最终可以用边界上的点按照顺时针顺序表示。

用橡皮筋套住平面的所有点，其形状就是凸包

现在假设这个平面没有三点共线 (方便起见)，如何找到平面的凸包？

——最直白的做法对任意两点 p_i , p_j，把它们连成一条直线，如果其余点都在其同侧，它就是凸包的一条边。但复杂度为 T(n) = $\binom{n}{2}\times (n-2)$ = Θ(n³)，which is unbearably slow。

“Geometric” D&C

自然的做法就是使用 D&C 的思路，将问题二分递归：

将所有点按 x 坐标大小二分 (按中值二分)，然后递归计算两边的凸包。但 Divide 完怎么 Combine 回来呢？

• 将两边 y 坐标最大/最小的点分别相连可以吗？想象左上图中，如果 B 区域在 y 轴负方向平移，a₄ 和 b₂ 的连线就会把 a₅ 留在包外，显然没这么简单。

• 但不难发现，合并时永远只需要画两条公共切线，⚠️ 且它们可以将所有其余任意连线线段夹在中间。

  假设 A (左凸包) 中 x 最大的点和 B (右凸包) 中 x 最小的点分别为 a₁ 和 b₁，那么直线 L = k, k∈(x_a1 , x_a2) 与平面中任意线段 a_ib_j 必然有交点。

  而且，因为两条公共切线在 y 方向上将所有线段夹在中间，使这个交点最高/低的 a_ib_j 就是我们要找的两条线！

• 我们仍需遍历 a_ib_j，但是简单的两点相连显然是 Θ(n₂) 的。但通过从 a₁b₁ 出发依次逆时针移动 a、顺时针移动 b，直到遇到交点 y 值的“峰”就找到了下切线。反之亦然。（如右上图）

• 在这种前提下，每个点最多“经过”常数次，Θ(n)

我们每次二分将问题减半，故总复杂度：T(n) = 2T(n/2) + Θ(n)，根据 Master Thm.，s=1，T(n) = Θ(nlogn)。

⚠️ Divide 不会天然让解题变快，如何便宜的完成 combine (merge) 在复杂问题中是关键。

Median Finding 中值问题

假设把数组 S 从小到大排序：s₍₁₎ , s₍₂₎ , … , s_(i) , … , s_(n)，将 s_(i) 称为 i-th order statistic 即 element of rank i。Lower median i = ⌊(n+1)/2⌋，upper median i = ⌈(n+1)/2⌉。

问题 Select(S, i) = find element of rank i in S（方便起见，假设数组无重复值）

——最直白的做法就是 sort(S) & return(S[1])，但在最简排序下也需要 Θ(nlogn)，但是其实不需要所有数的完整排序。

Partition + D&C

选择 Pivot x 将 S 分为两段：Partition(S, x) → A = {m | m ≤ x, m ∈ S}、{x}、B = {n | n ≥ x, n ∈ S} 三个无序数组，那么 x 就是第 |A|+1 小的元素。因为无需排序，T(n) = T(max{|A|, |B|}) + Θ(n) 复杂度。我们要找的 i-th 就在 A 或 B 中，平均来说，问题缩小约 1/2——

int Select(const vector<int>& S, int i, int x) {
    vector<int> B, A;
    for (size_t j = 1; j < S.size(); ++j) {
        if (S[j] < x) B.push_back(S[j]);
        else A.push_back(S[j]);
    }
		int k = static_cast<int>(A.size()) + 1;  // x 是 S 的 rank k 元素
		if (k == i) return x; 
    		else if (k > i) return Select(A, i);
    		else return Select(B, i - k);

但这套方法很吃我们选的 x 的选择：

• 极端差情况：每次都是剩余项中的最大/最小，T(n) = T(n-1) + Θ(n) = n+(n-1)+(n-2)+ … = Θ(n²)
• 非极端情况：每次都是剩余项的 10～90 百分位，max{|A|, |B|} ≤ 0.9n，T(n) ≤ T(.9𝑛) + Θ(n) = cn(1+0.9+0.9²+...) = Θ(n) 首项(根)主导

⚠️ 故 Selection 不要求 pivot 恰好是 median，不极端偏心即可保持线性复杂。

聪明地选 Pivot：Median of Medians

把所有元素 5 个一列，对每列求 Median M_i，x = Median(M_i)——中值的中值。

原理：至少有一半的列中值 ≥ x，其中每列又有至少 3 个 (超过一半) 元素 ≥ x，所以至少 n/4 个元素 ≥ x，vice versa 可知 n/4 < (|A|, |B|) < 3n/4 (对于足够大的 n > 200)。从而 x 一定处于 25～75 百分位，使得 T(n) ≤ T(.75𝑛) + Θ(n)。

MoM Partition + D&C

递归时，产生两个不同问题：1、分为 column-of-5 并找中值的中值从而确定 pivot；2、从 pivot 处分完后答案更高概率在较大 (≤3n/4) 侧，对其递归继续 partition——

$$ T(n) \le T(n/5)+T(3n/4)+{\frac{n}{5}\Theta(1)}+{\Theta(n)} $$

即：递归确定 MoM + 递归找答案 + 每列算 Median + 分组/Partition 等线性时间

所以： $$ \boxed{ T(n)\le T(n/5)+T(3n/4)+\Theta(n)=T(\frac15+\frac34)+\Theta(n)=20cn } $$

用数学归纳法(Induction)证明

⚠️ 证明复杂度基本都用数学归纳法，上周已经做过一次

猜 T(n) ≤ cn 对任意 n 有 T(n) ≤ T(n/5) + T(3n/4) + c₀n (n ≥ 200)：

假设对 200 ≤ k < n，T(k) ≤ ck 成立，然后证明对 k = n 成立——

那么 T(n) ≤ (1/5 + 3/4)cn +c₀n = 19cn/20 + c₀n，那么必须有 19cn/20 + c₀n < cn，即 c > 20c₀ 使得 T(n) < cn

Randomized Algorithms 随机算法

主线：Randomized Selection → Markov Inequality → Freivalds Algorithm → Monte Carlo / Las Vegas → Randomized Quicksort

刚才花了不少功夫，用 MoM 确定良好的 pivot，一定不能随机选吗？

我们的算法中 x 落在 25～50 百分位算“良好”，随机选就有 p = 1/2 的可能性“幸运”，那么得到一个“良好”pivot 的尝试次数 τ 的期望就是：E[τ]
=
p+(1-p)(1+E[τ]) → E[τ]
 = 1/p = 2—— $$ \bar{T}(n) = \bar{T}({3n\over 4}) + E[τ]\times cn=\bar{T}({3n\over 4})+2cn,\ thus\ \bar{T}(n)=\Theta(n) $$ “啊？既然随机选也是线性复杂度，为什么刚才要搞什么中值的中值？”

我知道你很急，但你先别急，但是代价是什么呢？虽然更简单高效，但是丢失了确定性。

Markov Inequality 马可夫不等式

对于一个非负随机变量 τ，它大于其期望 a 倍的概率存在上界：
P(τ > aE[τ]]) ≤ 
1/a （a > 1），即 P(x > a) ≤ E(τ)/a。

对于我们的问题，尝试超过 100 次的概率就是 P(τ > 100) ≤ 2/100 = 2%

Freivalds Algorithm

Randomization 不只是让算法跑得快，还能让一个原本昂贵的验证问题变得很便宜。

假设有一个 n x n 元素为 {-1, 0, 1} 的矩阵 D，我们想判断 D=0 也就是是不是所有元素都是 0。

李似不似傻，是不是 0 看一眼不就知道了，那如果 D 是一个计算结果，我想在不计算的前提下预测其为 0 的概率呢？比如等式成立的概率？

在 2-模运算下，1+1 mod 2 = 0，1+0 mod 2 = 1，0+0 mod 2 = 0，-1 ≡ 1 (mod 2)（-1 和 1 等价）

随机生成 r = [r₁ r₂ … r_n]^T 其中 P(r_i=0)=P(r_i=1) = 0.5，然后计算 D·r：

• 如果 D = 0 就一定有 Dr = 0

• 如果 D ≠ 0 那至少存在一行 d_i ≠ 0，d_ir mod 2 有至少 0.5 的概率不为 0

  P(Dr ≠ 0 | D ≠ 0) ≥ 0.5 ≡ P(Dr = 0 | D ≠ 0) ≤ 0.5，任意非零矩阵最多有50%概率能骗过这个随机测试

  更精确的：P(Dr = 0 | D ≠ 0) = 2^-rank_F2(D)（二分之一的“D 矩阵数域为 2 下线性独立的行数”次方，这可以用 rank-nullity 定理证明）

  ⚠️ 严格来说不能像课本上直接说有几行就是几次方！每增加一个线性独立的行约束，随机 (r) 能满足 (Dr=0) 的概率再减半。重复或线性相关的行不会增加新的约束 (乘 r 结果一样)，不会降低概率。

  原理：全 0 行 d_az·r mod 2 一定是 0；但非 0 行 1 ≡ -1 且无论非 0 元素个数，结果是 0 (没发现是 0) 还是 1 的概率都是 0.5。

从而，X = Y 可以转化为 Xr=Yr (mod2)，如相等前者一定等，不相等则至少 ≥ 0.5 可能不相等。

⚠️ 更重要的应用：AB = C (Θ(n³)) 可转为 A(Br) = Cr (mod2)，向量 x 矩阵只有 Θ(n²)。如果尝试 k 次，误判概率 ≤ 0.5^k，k 作为常数在复杂度中省略。

Monte Carlo vs Las Vegas

摩纳哥 VS 赌城，你们知道吗？就在刚刚…

相同 input，不同运行可能产生不同执行时间，甚至不同结果，我们把随机算法分成两类。

• Monte Carlo Algorithm：保证运行时间，答案可能错——固定尝试次数 τ，保证算法自身的多项式时间复杂度

• Las Vegas Algorithm：保证答案，运行时间随机——无法保证尝试次数，期望 E[T(n)] 为多项式时间复杂度

⚠️ 具有期望多项式运行时间的 Las Vegas 算法，可以通过设置 timeout 转化为 Monte Carlo 算法——

• 假设一个 Las Vegas 算法 E[T(n)] ≤ p(n)（p 为某多项式），任意常数 c，规定最多运行 cp(n) 步

  根据 Markov：P(T ≥ cp(n)) ≤ 1/c，这个新算法成功率 ≥ 1-1/c，固定时间 O(p(n))，显然被转为 Monte Carlo 算法。

Quicksort 快排序

“将中点作 pivot 左边放更小元素右边放更大元素，对左右区间递归；n 个元素放在二叉树中高 log₂n，平均/最优 O(nlogn)，如果每次选择的 pivot 都是最大/最小值，退化为 Selection Sort，最差 O(n²)。与有序性无关，但与 pivot 有关。”——CSC3060

Quicksort 和上文探讨的 Partition、Median Finding 的语境十分类似（只不过 pivot 的两侧都要递归）也有相似的“运气”问题……

• 极端差情况：每次都是剩余项中的最大/最小，T(n) = T(n-1) + T(0) + Θ(n) = n+(n-1)+(n-2)+ … = Θ(n²)

• 普通情况：每次都是剩余项的 10～90 百分位，max{|A|, |B|} ≤ 0.9n，T(n) ≤ T(.9𝑛) + T(.1n) + Θ(n) = Θ(nlogn)

• 最优情况：每次都是中值，完全平衡，T(n) = T(n/2) + T(n/2) + Θ(n) = Θ(nlogn)

平均：对于 T(n)，有 n 种 x 的可能，就有 $$ \bar{T}(n) = {(\bar{T}(0) + \bar{T}(n-1)) + (\bar{T}(1)+\bar{T}(n-1))+… + (\bar{T}(n-1) + \bar{T}(0)) \over n}+\Theta(n) $$ 即： $$ \bar{T}(n) = {1\over n}\sum_{k=0}^{n-1}2\bar{T}(k)+\Theta(n) = \Theta(n \log n)\ (Induction\ Prove) $$ 随机的 x 意味着程序无法被揣测并攻击 (即使使用最差的逆序数组)

II. Fast Fourier Transformation

多项式

任意一元 degree-d (d-次) 多项式 A(x)，都可以用两种形式表达：一种是 d+1 维系数向量 [a₀, a₁, … , a_d]；另一种是点值表示 A(x₀), A(x₁), …

次数不超过 d 的多项式，由 d+1 个系数唯一确定，也由 d+1 个不同点上的值唯一确定

两种表示之间的转换叫：

• 求值（evaluation）：系数 → 点值。
• 插值（interpolation）：点值 → 系数。

系数表示乘法较贵

对于两个 degree-d 多项式相乘： $$ A(x) = \sum_{k=0}^d a_kx^k,\ B(x)=\sum_{k=0}^d b_kx^k $$

$$ C(x)=A(x)B(x)=\sum_{i=0}^{2d}c_ix^i,\ where\ c_k=\sum_{j=0}^ka_jb_{k-j} $$

(a 和 b 系数和为 k，这些项相乘后求和作为 x^k 项系数)

对于 C_k→O(k)，直接计算系数 C(k)→O(d²)（要考虑 (d+1)² 对系数）

点值表示乘法便宜

系数→求点值→乘积的点值→差值求乘积的系数（evaluation, multiplication, interpolation）

1. 至少需要 2d+1 个不同点，通常选不小于 2d+1 的最小 2 的幂
2. 求值不能慢，如果在 n 个点上分别代入求值，总耗时仍是 O(n²)，没有获得加速
3. 逐点相乘只需要 O(n)
4. 解线性方程组的通用做法需要 O(n³)，但并非插值问题最短时间

求值：A(x) = M_n·a（中间这个矩阵就是 Vandermonde matrix 范德蒙德矩阵）

插值：a = M^-1_n(x)·A(ω)（高斯消元/矩阵求逆 下 O(n³)）

“计算一个多项式在很多点上的值”其实就是一个线性变换。

而只要 x₀,…,x_n-1 两两不同，Vandermonde 矩阵就是可逆的，因为 det V 
= Π_{0≤ i<j<n}(x_j-x_i) ≠ 0

线性系统

多项式乘法中的公式： $$ c_k=\sum_{j=0}^ka_jb_{k-j} $$ 叫作离散卷积 discrete convolution——一个函数在另一个函数上的加权叠加

例：考虑时间离散的线性-时不变 (Linear Time-Invariant: LTI) 系统。假设一个单位脉冲输入 δ(t)，会产生响应 b(t)。LTI 系统计算输出 c(t) 的过程，本身就是输入 δ(t) 与系统冲激响应 b(t) 的卷积。所以加快多项式乘法可以用于信号处理。

说到底如何快速计算 A(x_k) 的点值呢？

快速计算点值与 FFT

刚才我们说到高斯/求逆硬解插值需要恐怖的立方时间，比系数表示 (硬解卷积) 的 baseline 更差。FFT 故意选择一组极其特殊的 x_k：roots of unity 单位根——x_k =ω^k,
 ω=e^2πi/n

D&C with even/odd terms

将给定多项式 A(x) degree≤(n-1) 奇偶次项二分：A(x) = A_e(x²) + xA_o(x²)（A_e()、A_o() 为 degree≤(n/2-1) 的多项式）

虽然多项式系数数量减半，但我们还需要让每个子问题的求值点数量也减半到 n/2，那么 T(n) = 2T(n/2)+O(n) = O(nlogn)——

符合直觉来说：因为 x 和 -x 的平方相同，只需计算一次 A_e(x²) 和 xA_o(x²)，就得到两个答案……但是如何递归？

⚠️ 假设我们取点 ±1、±2、±3…，第一轮递归任务减半变成 1、4、9…可是它们之间无法继续形成正负对了！

所以我们需要的是：一组正负配对的点；平方去重后仍然可以正负配对；继续平方保持此结构直到剩下一点

那必然就需要平方为负数的数——

复数与极坐标

复平面上，负数 z = a + bi 投影于点 (a, b)，故可用极坐标表达 z = r(cosθ + isinθ) = re^iθ（r = |z| = (a² + b²)^1/2，θ ∈ [0, 2π) 表示 z 与实正半的逆时针夹角）

而且，极坐标点的乘法更容易：(r₁, θ₁) x (r₂, θ₂) = (r₁r₂, θ₁+θ₂)

• 考虑极坐标系单位圆 (r=1) 上的任意 z，zⁿ = (1, nθ)
• 复平面中的 1 是所有 (1,2πk) k∈N，故 (zⁿ = 1) ≡ ((1, nθ) = (1,2πk))，即 θ = 2πk/n，即 z_k = (1, 2πk/n) = e^2πi/n (0≤i≤n-1)

由此我们定义了 primitive root ω=e^2πi/n，n-次单位复根 x_k = ω_n^k = e^2πi·k/n (表示第 k 个复根) 就是单位圆上的所有 n 均匀等分点。比如 4 次单位根：1、i、-1、-i

(偶数) 单位复根的性质：

1. 正负配对：偶数次单位根在圆上两两中心对称—— $$ -ω = (1, {2πi·k\over n}+π) = (1, {2πi·(k+{n\over 2})\over n})\ \ while\ n∈even,\ -ω∈roots\ of\ unity $$ 换句话说，对于任意偶数 n，ω_n^n/2 = -1，故 ω_n^j+n/2 = -ω_n^j（所有点都有相反数点）

2. 平方减半：ω_n² = e^2πi/(n/2)=ω_n/2

   人话：每个 n/2 次单位根，都对应两个互为相反数的 n 次单位根——这两个 n 次单位根平方后得到同个 n/2 次单位根

   …… → {1,i,-1,-i} → {1,-1} → {1} （当 n 是 2 的幂这种结构就能一直递归下去）

快速傅里叶变换 FFT

对于 A(x) of degree≤n-1，且 n (总问题规模) 为 2 的幂，选择 n 个不同点 x_j =ω^j_n

FFT 想一次性算出所有 j-次单位根：A(ω^j_n)，j=0, … , n-1。换句话说：

• 输入：n 个系数
• 输出：在 n 个单位根上的 n 个点值。

奇偶二分：A(x) = A_e(x²) + xA_o(x²) → A(ω^j_n) = A_e(ω^j_n/2) + ω^j_nA_o(ω^j_n/2)，递归 A_e、A_o……递归完返回合并时，对于 A(ω^j_n) = E_j + ω^j_nO_j，那么奇偶二分的另一项无需单独算，就是：A(ω^j+n/2_n) = E_j - ω^j_nO_j。

• 向下递归时拆系数；向上返回时用正负配对合并点值
• 每一对输出只需要一次乘法、一次加法和一次减法——蝶形 butterfly 算法
• T(n) = 2T(n/2) + O(n) = O(nlogn)

FFT(a, ω) {
    n = a 的长度
    if (n = 1) return a
    E = FFT(a 的偶数下标元素, ω²)
    O = FFT(a 的奇数下标元素, ω²)
    for (j=0; j ≤ n/2 - 1; j++) {
        t = ω^j × O[j]
        y[j] = E[j] + t
        y[j + n/2] = E[j] - t
    }
    return y
}

Example of FFT

核心还是 A(x) = A_e(x²) + xA_o(x²)

例如 A(x) = 1 + 2x + 3x² + 4x³，输入系数 a = [1, 2, 3, 4]，4次单位负根为 1、i、-1、-i，要计算 A(单位负根)：

• 奇偶二分：A_o(x²) = 2+4x²；A_e(x²) = 1+3x²

• ⚠️ 合并返回：A_o(1) = 6，A_o(-1) = -2 → O = [6, -2]

  A_e(1) = 4，A_e(-1) = -2 → E = [4, -2]

  A(ω⁰₄) = E₀ + ω⁰₄O₀ = 4+1·6 = 10，那么 A(ω²₄) = E₀ - ω⁰₄O₀ = -2

  A(ω¹₄) = E₁ + ω¹₄O₁ = -2+i·(-2) = -2-2i，则 A(ω³₄) = E₁ - ω¹₄O₁ = -2+2i

• 得到 [10, -2-2i, -2, -2+2i]

离散傅里叶变换与点值表示的矩阵 (多项式) 乘法

求值：A(ω) = M_n(ω)·a = FFT[a, ω]

插值：a = M^-1_n(ω)·A(ω) = FFT[A, ω^-1]/n

也就是说，插值根本不需要真的做一般的 (O(n³) 矩阵求逆。对于单位根构成的这个特殊矩阵，只需要：ω → ω^-1 做一次 FFT，然后所有结果再除以 n 即可，这就是 Inverse FFT（IFFT）。

FFT 解多项式 (矩阵) 乘法

例如 A(x) = 1 + 2x；B(x) = 3 + 4x，乘完后最高系数为 2，至少需要 3 个点值，取最小 power-2，即 4

1. FFT 输入 a = [1, 2, 0, 0]；b = [3, 4, 0, 0]（补 0）和 ω = [1, i, -1, -i]

2. 正向 FFT 求值：A(ω) = [3, 1+2i, -1, 1-2i]；B(ω) = [7, 3+4i, -1, 3-4i]

3. 逐点相乘：C(ω) = [21, -5+10i, 1, -5-10i]

4. IFFT 输入 C(ω) 和 ω^-1 再做一次 FFT 然后除以 4：c = [3, 10, 8, 0]，即 3 + 10x + 8x²

   ωω^-1 = 1 取逆！不是求倒数！！！

FFT / IFFT 计算笔记

FFT：系数表示 → 点值表示。

对于系数 (a₀, a₁, ..., aₙ₋₁)，先写成：

A(x) = a₀ + a₁x + ... + aₙ₋₁xⁿ⁻¹

取 n 个 n 次单位根作为 x，分别代入 A(x)：

FFT(A) = (A(x₀), A(x₁), ..., A(xₙ₋₁))

例如 n=4 时（按本课约定）取值：

x = 1, i, -1, -i

IFFT：点值表示 → 系数表示。

计算方法几乎一样，但把 FFT 的取值点换成它们的共轭 / 逆，最后整体除以 n：

FFT  (n=4) 取：1,  i, -1, -i
IFFT (n=4) 取：1, -i, -1,  i
                         ↑
              每个单位根取逆（也等于取共轭）

最后所有结果 ÷ 4

所以一句话记忆：

FFT：代入 n 次单位根；IFFT：代入这些单位根的逆（共轭），再 ÷ n。

注意不是“取相反数”：i → -i 看起来是相反数，但 1 → 1、-1 → -1，准确说是 取逆 / 共轭。

例如要算 A(x) × B(x)：

① 补 0
选择足够大的 n（2 的幂），保证 n > 最终多项式次数

A coefficients → 补到 n 项
B coefficients → 补到 n 项

② 分别 FFT
FFT(A) = (A₀, A₁, ..., Aₙ₋₁)
FFT(B) = (B₀, B₁, ..., Bₙ₋₁)

③ 对应位置直接相乘（componentwise）
C = (A₀B₀, A₁B₁, ..., Aₙ₋₁Bₙ₋₁)

④ IFFT
coefficients of A(x)B(x) = IFFT(C)

Visualized FFT

• CSC4120