4568 字

23 分钟

【算法设计与分析】分治法

2023-03-09

学习笔记/算法设计与分析

algorithm

/

分治法

前言#

这一章将介绍重要的算法设计策略——分治法（divide and conquer）。本章基于《Introduction to Algorithms》4th 第四章，包含分治法解决问题（以矩阵乘法为例）和分析分治法算法的运行时间——代入法、递归树法以及主方法。

矩阵乘法#

这一章举例是以N x N的两个矩阵A, B相乘，得到N x N矩阵C。

朴素解法#

在线性代数中我们学到，两个矩阵相乘，C的第 i 行第 j 列的数 $C_{ij}$ = A 的第 i 行 x B 的第 j 列（逐位相乘后逐位相加）。

抽象成数学语言，设 $C_{ijk}$ = A的第 i 行第 k 个数 x B 的第j列第 k 个数，于是有：

我们用伪代码很容易实现：

1
MATRIX-MULTIPLY(A, B, C, n)
2
    for i = 1 to n
3
        for j = 1 to n
4
            for k = 1 to n
5
                c_ij = c_ij + a_ik * b_kj

很明显，它有三层循环，时间复杂度为 $Θ(n^3)$ 。

我们可以考虑使用分治法改造它，使用分治法我们要做到：

分解(Divide)：分治法中的分，将原问题分解为若干子问题。
解决(Conquer)：分治法中的治，求出子问题的解。
合并(Combine)：将子问题的解合并为原问题的解。

分解#

我们可以把矩阵分为4个 $\frac{n}{2} \times \frac{n}{2}$ 的子矩阵（n是2的幂），分别为 $A_{11}, A_{12}, A_{21}, A_{22}$ 和 $B_{11}, B_{12}, B_{21}, B_{22}$ ，同时也会生成结果为4个n/2 X n/2的子矩阵： $C_{11}, C_{12}, C_{21}, C_{22}$

此时我们的公式就可以从 $C = A X B$ 改写为如下形式：

接下来我们要做的就是继续分割每个矩阵到最小的部分。每次都分成4个部分，直到无法再分。

解决#

我们已经分割成子问题去解决了，接下来我们只需要去解决子矩阵之间的乘法即可。

1
MATRIX-MULTIPLY-RECURSIVE(A, B, C, n)
2
    if n == 1
3
        c_11 = c_11 + a_11 · b_11
4
        return
5
    partition A, B, C into n/2 × n/2 submatrices A_11, A_12, A_21, A_22; B_11, B_12, B_21, B_22; C_11, C_12, C_21, C_22;
6
    MATRIX-MULTIPLY-RECURSIVE(A_11, B_11, C_11, n/2)
7
    MATRIX-MULTIPLY-RECURSIVE(A_11, B_12, C_12, n/2)
8
    MATRIX-MULTIPLY-RECURSIVE(A_21, B_11, C_21, n/2)
9
    MATRIX-MULTIPLY-RECURSIVE(A_21, B_12, C_22, n/2)
10
    MATRIX-MULTIPLY-RECURSIVE(A_12, B_21, C_11, n/2)
11
    MATRIX-MULTIPLY-RECURSIVE(A_12, B_22, C_12, n/2)
12
    MATRIX-MULTIPLY-RECURSIVE(A_22, B_21, C_21, n/2)
13
    MATRIX-MULTIPLY-RECURSIVE(A_22, B_22, C_22, n/2)

如之前图，我们执行了8次乘法运算得到了 $C$ 矩阵，运用了8个递归表达式。设运行时间为 $T(n)$ ，故8次递归的运行时间为 $8T(n/2)$ ，再加上其他语句，我们可以得到： $T(n) = 8T(n/2) + Θ(1)$

其时间复杂度为 $Θ(n^{\log 8}) = Θ(n^3)$ 但是这样的时间复杂度依旧很高，我们还可以进一步优化。

Strassen算法#

1969年，Volker Strassen发明了算法，将矩阵乘法的时间复杂度降低到了 $Θ(n^{\log 7})$ ，算法在运用分治法解决问题过程中，减少了递归树的茂密程度。Strassen算法将递归树的 8 个分枝减少为 7 个分枝，即只进行了7次递归即可解决。

从本质上来看，Strassen算法把原本计算时间长的8的乘法，转成了较多加减法的乘法，减少乘法的数量。这个分解思想来自于《On the Complexity of Matrix Multiplication》，这里只是介绍其过程，只需要知道这个算法将矩阵乘法的时间复杂度降低到了 $Θ(n^{\log 7})$ 。

分解#

和之前一样，我们先把矩阵分为4个 $n/2 \times n/2$ 的子矩阵，得到：

接下来，创建10个 $n/2 \times n/2$ 的矩阵 $S_1,S_2,\dots,S_{10}$ 保存上一步的子矩阵的和或差。定义为：

这 10 个矩阵分为两组，每组 5 个，都是 $A$ 的子矩阵或者 $B$ 的子矩阵内部之间相互加减。这个过程涉及矩阵加减法，时间复杂度为 $Θ(n^2)$ 。

解决#

接下来我们递归计算7个矩阵乘积 $P_1,P_2,\dots,P_7$ ，定义为：

最后可以直接计算出 $C_{11},C_{12},C_{21},C_{22}$ 。

Strassen算法用了 7 次矩阵乘法 18 次矩阵加减法。

无论如何，乘法运算的时间复杂度都高于加法运算的，所以通过增加常数个数的加法运算而减少乘法运算的策略都是非常划算的。

它的伪代码为：

1
STRASSEN(A, B, n)
2
    let C be a new n × n matrix
3
    if n == 1
4
        c[1, 1] = a[1, 1] · b[1, 1]
5
        return
6
    A[1, 1] = A[1 : n / 2][1 : n / 2]
7
    A[1, 2] = A[1 : n / 2][n / 2 + 1 : n]
8
    A[2, 1] = A[n / 2 + 1 : n][1 : n / 2]
9
    A[2, 2] = A[n / 2 + 1 : n][n / 2 + 1 : n]
10
    B[1, 1] = B[1 : n / 2][1 : n / 2]
11
    B[1, 2] = B[1 : n / 2][n / 2 + 1 : n]
12
    B[2, 1] = B[n / 2 + 1 : n][1 : n / 2]
13
    B[2, 2] = B[n / 2 + 1 : n][n / 2 + 1 : n]
14
    S[1] = B[1, 2] - B[2, 2]
15
    S[2] = A[1, 1] + A[1, 2]
16
    S[3] = A[2, 1] + A[2, 2]
17
    S[4] = B[2, 1] - B[1, 1]
18
    S[5] = A[1, 1] + A[2, 2]
19
    S[6] = B[1, 1] + B[2, 2]
20
    S[7] = A[1, 2] - A[2, 2]
21
    S[8] = B[2, 1] + B[2, 2]
22
    S[9] = A[1, 1] - A[2, 1]
23
    S[10] = B[1, 1] + B[1, 2]
24
    P[1] = STRASSEN(A[1, 1], S[1], n / 2)
25
    P[2] = STRASSEN(S[2], B[2, 2], n / 2)
26
    P[3] = STRASSEN(S[3], B[1, 1], n / 2)
27
    P[4] = STRASSEN(A[2, 2], S[4], n / 2)
28
    P[5] = STRASSEN(S[5], S[6], n / 2)
29
    P[6] = STRASSEN(S[7], S[8], n / 2)
30
    P[7] = STRASSEN(S[9], S[10], n / 2)
31
    C[1 : n / 2][1 : n / 2] = P[5] + P[4] - P[2] + P[6]
32
    C[1 : n / 2][n / 2 + 1 : n] = P[1] + P[2]
33
    C[n / 2 + 1 : n][1 : n / 2] = P[3] + P[4]
34
    C[n / 2 + 1 : n][n / 2 + 1 : n] = P[5] + P[1] - P[3] - P[7]
35
    return C

分治法分析#

求解分治法的运行时间没有这么直观，我们通过递归知道了分治法的时间表达式后，可以通过代入法、递归树法以及主方法来求解分治法的运行时间。

我们在分析中，存在递归式，一个递归式就是一个等式（方程式），它通过通常更小的函数参数值来描述一个函数。即是带有递归时间表达式。

代入法#

无渐进符号递归式#

代入法(substitution method) 分两步：

使用常数符号猜测解的形式。
使用数学归纳法求出解中常数，并证明解是正确的。

一般来说，我们的猜测有： $O(nlgn)$ 、 $O(n^{\lg n})$ 等等。接下来通过一个例子来解释代入法的过程。

例如e.g.1：

我们猜测 $T(n) = nlgn + n$ ①，因为 $2T(n/2)$ 通常是 $nlogn$ 的复杂度。进行数学归纳法的证明：

当 $n=1$ 时， $T(n) = 1lg1 +1 = 1$ ，符合 $T(n)$ 特征。
当 $n>1$ 时，当 $n=n/2$ 时， $T(n/2) = (n/2) lg (n/2) + n/2$ ② 把②带入原式 $T(n) = 2T(n/2) + n$ 中，得 $T(n) = n lg(n/2)+ 2n$ 。继续化简: $\begin{align*} T(n) &= n (\lg n - \lg 2) + 2n \\ \implies T(n) &= n (\lg n - 1) + 2n \\ \implies T(n) &= n \lg n + n \end{align*}$ 得到式子①，则证明成功。

此时 $T(n) = nlgn + n$ ，可以使用渐进符号表示为 $T(n) = Θ(nlgn)$ .

带有渐进符号的递归式#

这种是递归式中不包含渐进符号时我们做的，当然更加常见的是带有渐进符号的递归式，我们一般需要进行下面步骤：

用常数组合代替渐进符号
分别分析上下界，往往使用不同常数组合代替。

例如e.g.2： $T(n) = 2T(n/2) + Θ(n)$ .

首先我们如果想分析它的上界，即 $T(n) = 2T(n/2) + O(n)$ 。我们用 $cn$ 代替 $O(n)$ ，并且写成不等式的形式： $T(n) \le 2T(n/2) + cn$ ， $c$ 为正常数。

根据之前求到的，我们可以猜 $T(n) \le dn\ lgn$ ， $d$ 为其他正常数。和之前一样，对它进行数学归纳法证明：

当 $n=1$ 时， $T(n) \le d\ lg1= 0$ ，是个常数，符合T(n)特征。
当 $n>1$ 时，当 $n=n/2$ 时，

\begin{align*} T(n/2) &\le (dn/2) \lg(n/2) \\ \implies T(n) &\le 2T(n/2) + cn \\ &= dn \lg(n/2) + cn \\ &= dn (\lg n - 1) + cn \\ &= dn \lg n - dn + cn \\ &\le dn \lg n \end{align*}

（当且仅当 $-dn + cn \le 0 \implies c \le d$ 时成立，常数无所谓）

此时可以得到 $T(n) \le dn\ lgn$ ，写作 $T(n) = O(nlgn)$

同理我们还可以证明其下界，用 $T(n) \ge 2T(n/2) + cn$ 表示，猜测 $T(n) \ge dn\ lg n$ ，得到 $T(n) = Ω(nlgn)$ 。

因此，我们证明了 $T(n) = Θ(nlgn)$ 。

常数递归式减去低阶项#

当然，我们还可能遇到特殊的渐进符号，例如 $Θ(1)$ 。此时为了配平，右边减去低阶项，一阶项的低阶项为常数项。

例如e.g.3： $T(n) = 2T(n/2) + Θ(1)$ .

我们猜它： $T(n) = O(n)$ ，则此时写成 $T(n) \le cn - d$ ，而不是 $T(n) \le cn$ 。因为在验证中，若不等式写为 $T(n) \le cn$ ，则会出现下面过程：

\begin{align*} T(n) &\le 2(c(n/2)) + \Theta(1) \\ &= cn + \Theta(1) \\ &\not\le cn \end{align*}

此时这个不等式不一定成立，因为 $Θ(1)$ 不一定小于0，无法得到一个常数范围使其成立，证明不成立。我们应该使用 $T(n) \le cn - d$ ，作为配平。

\begin{align*} T(n) &\le cn - d, \text{ where } d \ge 0 \text{ is a constant:} \\ T(n) &\le 2(c(n/2) - d) + \Theta(1) \\ &= cn - 2d + \Theta(1) \\ &\le cn - d - (d - \Theta(1)) \\ &\le cn - d \end{align*}

此时我们只需要 $d\ -\ Θ(1)>0 \implies d>Θ(1)$ 即可成立。此时有一个范围使其成立，证明成立。

常见错误#

我们证明的时候可能会出现一些常见的问题，例如： $T(n) = 2T(n/2) + Θ(n)$ ，猜测 $T(n) = O(n)$

错误写法，直接带入渐进符号： $T(n) \le 2O(n/2) + Θ(n) = 2O(n) + Θ(n) = O(n)$ 因为此时 $O(n)$ 所隐藏的系数不断变化，也就是表达式不一定成立。
错误写法，不替换渐进符号： $T(n) \le 2c(n/2) + Θ(n) = cn + Θ(n) = O(n)$ 我们无法确定 $Θ(n)$ 所隐藏的系数一定大于 $c$ ，表达式不一定成立。因此我们需要使用常数组合表达渐进符号。

递归树法#

在递归树中，每个结点表示一个单一子问题的代价（需要时间），子问题对应某次递归函数的调用。我们将树中每层中的代价求和，得到每层代价，然后将所有层的代价求和，得到所有层次递归调用的总代价。

例如分析： $T(n) = 3T(n/4) + Θ(n^2)$ . 此时 $Θ(n^2)$ 是个二次匿名函数，我们可以用 $cn^2$ 表示，此时 $c$ 为常数且 $c>0$ 。故我们可以分析： $T(n) = 3T(n/4) + cn^2$ 。因为本质上递归中执行的为若干个 $cn^2$ （递归的终点是执行最后一次 $cn^2$ ），故我们通过分析 $cn^2$ 得到递归树。

因为我们递归有3个，故从 $T(n)$ 处可以分成三个分支：

$T(n/4)$ 又可以通过 $T(n) = 3T(n/4) + cn^2$ 得到， $T(n/4) = 3T(n/16) + cn^2$ . 故每个 $T(n/4)$ 又可以有三个分支：

同理我们可以继续分 $T(n/16)$ ，直到最后 n 分到 1 为止，此时递归树为：

我们发现到了最后的时候，全都是 $Θ(1)$ 的运行时间。我们分析每一层，第0层运行时间为 $cn^2$ ，第1层加起来为 $(3/16)cn^2$ ，第2层则为 $(3/16)^2cn^2$ 。此时它们与层数呈等比数列：第i层运行时间为 $(3/16)^icn^2$ 。

每个子结点规模为父母结点的1/4，最深处我们定义为 $d$ ，此时最深处子节点规模为1， $n/(4)^d = 1 \implies d = \log_4 n$ 。在树中，高度与深度的数值一样，故高度为 $h = \log_4 n$ 。

我们知道了高度，就可以求高度为 $h_i$ 处的结点的位置，那一层的结点数为 $3^{h_i}$ 。例如第0层，结点数为 $3^0 = 1$ ，第一层为 $3^1 = 3$ ，第二层为 $3^2 = 9$ . 在最后一层，结点数为 $3^{\log_4 n}$ 。

此时根据数学知识得，我们有一个指数公式： $a^{\log_b c} = c^{\log_b a}$ 。

所以可以得到 $3^{\log_4 n} = n^{\log_4 3}$ 。最后一层全都是 $Θ(1)$ 常数，故最后一层运行时间为： $Θ(1) \times n^{\log_4 3} =$ $Θ(n^{\log_4 3})$ 。

我们把所有层的运行时间加起来，除了最后一层外，其他层的和是一个等比数列求和。

此时我们可以求它的上下界。例如求其上界，继续对它化简得：

这里运用了放缩，使用一个无穷范围以消除指数。这样我们就求到了它的上界。

当然有些时候我们可能遇到树叶结点参差不齐的情况，例如： $T(n) = T(n/3) + T(2n/3) + Θ(n)$ .

这棵递归树是不平衡的，很明显一个父结点的两个子结点的代价不一样，左子结点的数据规模是父结点的 1/3，右子结点的数据规模是父结点的 2/3，从根结点到不同的叶结点的路径长度也不一样。

同样的，我们改为 $T(n) = T(n/3) + T(2n/3) + cn$ 进行分析。

我们发现，从根节点开始的路径在2/3的分支上更长，因为我们假设均为1/3分，或均为2/3分，得到节点数为 $2^{\log_3 n}$ 和 $2^{\log_{3/2} n}$ ，此时根据数学知识 $2^{\log_{3/2} n}$ 更加大，故此时 $2^{\log_{3/2} n}$ 的分支更长。

注： $(2/3)^i n=1 \implies i = \log_{3/2} n$ ，完全二叉树的情况下，节点数为 $2^i$

故我们得到这样的递归树：

此时最深的深度为 $\log_{3/2} n$ 。因为有叶子节点的存在，导致每一层的代价和小于 $cn$ （没到叶子节点的每层代价和为 $cn$ ，到了叶子节点的小于 $cn$ ）。此时我们计算总的代价和，可以用 $O(cn\log_{3/2}n) = O(n\lg n)$ 表示（因为假设最多的情况为完全满二叉树，高度为 $\log_{3/2} n$ ，实际上达不到，作为一个粗略的上界）。

因为是一个粗略的上界，我们需要使用代入进行验证。即证明： $T(n)\le dn\ lgn$ ， $d>0$ 是个常数。

\begin{align*} T(n) &\le T(n/3) + T(2n/3) + \Theta(n) \\ &\le d(n/3) \lg(n/3) + d(2n/3) \lg(2n/3) + cn \\ &= d(n/3) (\lg n - \lg 3) + d(2n/3) (\lg n - \lg(3/2)) + cn \\ &= dn \lg n - (d(n/3)\lg 3 + d(2n/3) \lg(3/2)) + cn \\ &= dn \lg n - dn((2\lg 3)/3 + (2\lg 3/2)/3) + cn \\ &\le dn \lg n \text{ 当且仅当 } c/(\lg 3 - (2/3)) \le d \end{align*}

此时得到的一个粗略上界即可，无需对递归树的每层代价进行精确计算。

同理我们还可以对它进行下界的估计，发现 $T(n)\ge dn\ lgn$ 在 $d>0$ 时也成立，故 $T(n) = Θ(nlogn)$ 。

Master method#

主方法主要来自于递归树法的总结，适用于形如 $T(n) = aT(n/b) + f(n)$ 的递归式。此时 $a ≥ 1$ , $b > 1$ ，把 $f(n)$ 称为驱动函数（driving function），它是一个渐进非负函数。

同时，我们的递归式还可以推导为： $T(n) = a'\ T(⌊n/b⌋) + a''\ T(⌈n/b⌉) + f(n)$ ， $a'$ , $a'' ≥ 0$ and $a' + a'' = a$ .

则分为下面三种情况：

如果存在 $ε>0$ ， $f(n) = O(n^{\log_b a-ε})$ ，则 $T(n) = Θ(n^{\log_b a})$ ；
如果存在 $k≥0$ ， $f(n) = Θ(n^{\log_b a} \lg^k n)$ ，则 $T(n) = Θ(n^{\log_b a} \lg^{k+1} n)$ ；
如果存在 $ε>0$ ， $f(n) = Ω(n^{\log_b a+ε})$ ，且对于某个常数 $c<1$ 和所有足够大的 $n$ 有 $af(n/b) ≤ cf(n)$ 则 $T(n) = Θ(f(n)).$

主方法的含义是： $f(n)$ 与 $n^{\log_b a}$ 两个函数中较大的那个决定递归式的解。其中f(n)通常是使用常数c的组合表示，在递归树法中分析代价所用，而 $n^{\log_b a}$ 常常指递归树中最后一层节点的代价（ $Θ(n^{\log_b a})$ ）。

当 $f(n) < n^{\log_b a}$ 时，递归式的解 $T(n) = Θ(n^{\log_b a})$ ；这种情况下代价主要来源于叶子节点。
当 $f(n) ≈ n^{\log_b a}$ 时，递归式的解乘上一个对数因子， $T(n) = Θ(n^{\log_b a} \lg^{k+1} n) = Θ(f(n) \lg^{k+1} n)$ ；这种情况下每一层的代价基本相等为 $n^{\log_b a} \lg^k$ ，而其高度一般为 $Θ(\lg n)$ 。
当 $f(n) > n^{\log_b a}$ 时，递归式的解 $T(n) = Θ(f(n))$ 。这种情况下代价主要来源于根节点。

此时的大小比较时多项式的渐进大小比较，而不是常数项与低次项的比较。实际上比的是两个函数渐近界的大小，即哪个函数增长得快，它们至少要相差 $n^ε$ 个因子。

下面分情况分析如何使用主方法。

eg.1： $T(n) = 2T(n/4) + 1$ .

分析可知： $f(n) = 1$ ， $n^{\log_b a} = n^{\log_4 2} = n^{1/2}$ 。很明显， $f(n) < n^{\log_b a}$ ，可以使用主方法中的情况1，得到 $T(n) = Θ(n^{1/2})$

eg.2： $T(n) = 2T(n/4) + \sqrt{n}\lg^2 n$

分析可知： $f(n) = \sqrt{n}\lg^2 n$ ， $n^{\log_b a} = n^{\log_4 2} = n^{1/2}$ 。此时在增长率上， $f(n) = Θ(n^{1/2} \lg^2 n) = Θ(n^{\log_b a} \lg^2 n)$ ，它们增长率相当，适用于情况2，得到 $T(n) = Θ(n^{1/2} \lg^3 n)$ 。

eg.3： $T(n) = 2T(n/4) + n^2$ .

分析可知： $f(n) = n^2$ ， $n^{\log_b a} = n^{\log_4 2} = n^{1/2}$ 。很明显， $f(n) > n^{\log_b a}$ ，且此时满足 $2f(n/4) = (1/8)n^2 ≤ cf(n)$ ，此时 $0<c<1$ 。可以使用主方法中的情况3，得到 $T(n) = Θ(n^2)$ 。

主方法并不能覆盖所有的递归式，它会出现空隙无法使用。

当 $f(n)$ 与 $n^{\log_b a}$ 两个函数无法比较时。
当 $f(n) = o(n^{\log_b a})$ 时，情况一和情况二之间会有一个间隙， $n^{\log_b a}$ 多项式的增长速率不大于 $f(n)$ 的多项式的增长速率。
当 $f(n) = ω(n^{\log_b a})$ 时，情况二和情况三之间会有一个间隙， $n^{\log_b a}$ 多项式的增长速率不小于 $f(n)$ 的多项式的增长速率。

解决这类递归式可以尝试使用其他方法，或Akra-Bazzi 递归式，这里不在正文处涉及。

后记#

本篇主要记录了如何使用分治法——以矩阵乘法为例，以及分析分治法的运行时间的三种方法——代入法、递归树法以及主方法。

References#

基于algorithm2e的LaTeX算法流程使用 (含中英文模板)_latex 算法流程-CSDN博客. https://blog.csdn.net
LaTeX算法排版 - huapyuan - 博客园. https://www.cnblogs.com
LaTeX公式引用技巧-CSDN博客. https://blog.csdn.net
Linear format equations using UnicodeMath and LaTeX in Word - Microsoft Support. https://support.microsoft.com
Welts. https://welts.xyz
LeTeX 写算法流程【algorithm,algorithmic】【algorithm, algorithmicx, algpseudocode】_latex写算法流程-CSDN博客. https://blog.csdn.net
LaTeX算法排版技巧. https://blog.csdn.net
Markdown 插入 LaTex 数学公式-腾讯云开发者社区-腾讯云. https://cloud.tencent.cn
Latex-算法伪代码 - 知乎. https://zhuanlan.zhihu.com
LaTeX 中的算法与伪代码实现-CSDN博客. https://blog.csdn.net