算法设计与分析

动态规划

计算机学院    张腾

tengzhang@hust.edu.cn

课程大纲

斐波那契数

$0, ~ 1, ~ 1, ~ 2, ~ 3, ~ 5, ~ 8, ~ 13, \ldots$

$$ \begin{align*} \quad F(n) = \begin{cases} n, & n \le 1 \\ F(n-1) + F(n-2), & n > 1 \end{cases} \end{align*} $$

分治法实现

def fibo_rec(n):
    if n <= 1:
        return n
    else:  # 递归求解两个子问题
        return fibo_rec(n - 1) + fibo_rec(n - 2)

效率如何？

分治法 时间复杂度

设$T(n)$为求$F(n)$需调用 fibo_rec() 的次数

• 边界：$T(0) = T(1) = 1$
• $T(n) = T(n-1) + T(n-2) + 1$

$$ \begin{align*} \quad \underbrace{T(n) + 1}_{G(n)} & = \underbrace{T(n-1) + 1}_{G(n-1)} + \underbrace{T(n-2) + 1}_{G(n-2)}, ~ \ldots, ~ G(0) = G(1) = 2 \\[10pt] & \Longrightarrow G(n) = 2 F(n+1) \\[5pt] & \Longrightarrow T(n) = 2 F(n+1) - 1 \end{align*} $$

$F(n) = \frac{1}{\sqrt{5}} ( (\frac{1 + \sqrt{5}}{2})^n - (\frac{1 - \sqrt{5}}{2})^n )$，$T(n)$亦为指数增长

分治法 时间复杂度

问题 F(n-k) 的个数为 F(k+1)

备忘录

保持分治递归的形式，每个子问题第一次求解后将结果存下来，每碰到一个子问题，先检查是否已求解过，若否则求解保存

def fibo_dp_memoized(n):
    F = [-float("inf")] * (n + 1)  # 备忘录初始化为负无穷
    return fibo_dp_memoized_aux(n, F)


def fibo_dp_memoized_aux(n, F):
    if n <= 1:
        return n
    else:
        if F[n] >= 0:
            return F[n]
        else:
            F[n] = fibo_dp_memoized_aux(n - 1, F) + fibo_dp_memoized_aux(n - 2, F)
            return F[n]

备忘录英文为Memoization，注意i前面没有r，由英国人工智能先驱唐纳德·米奇于1967年提出，但其实米奇本人只提出了memo function的概念，并没有用这个词，而更早的是香农在他1950年制造的机械鼠Theusus中就使用该词了

备忘录

• F[ ] 的填写顺序为 F[2]、F[3]、F[4]、……、F[n]
• 计算 F[2]、F[3]、F[4]、……、F[n] 只需做一次加法
• F[2]、F[3]、F[4]、……、F[n-2] 会被查询一次
• 时间复杂度为$O(n)$

动态规划

既然 F[ ] 的填写顺序为 F[2]、F[3]、F[4]、……、F[n]，为何不直接用循环顺序填写？

def fibo_dp_iter(n):
    F = [0] * (n + 1)
    F[1] = 1
    for i in range(2, n + 1):
        F[i] = F[i - 1] + F[i - 2]
    return F[n]

只有一重循环，时间复杂度为$O(n)$

动态规划 有选择地记

有时也未必需要记录所有子问题的结果，节省存储空间

def fibo_dp_iter2(n):
    prev, curr = 0, 1
    for i in range(n - 1):
        next = prev + curr
        prev, curr = curr, next
    return curr

动态规划 轶闻

动态规划 (dynamic programming, DP) 由美国数学家理查德·贝尔曼于 50 年代中期正式提出

programming 意为规划、计划、调度，不是编程的意思

• 线性规划 (linear programming, LP)
• 电视节目单 (television programming)
• 学位课程计划 (degree program)

1940 年左右

• Pierre Massé 使用动态规划优化法国水电站的运行，发表于 1944 年
• 冯·诺伊曼使用动态规划研究完美信息博弈，发表于 1944 年
• 图灵使用动态规划破译密码，80 年代中期解密

动态规划 轶闻

贝尔曼提出动态规划用于优化美国空军的训练和后勤保障

当时的国防部长是通用电气前首席执行官威尔逊，他非常讨厌数学研究，因此贝尔曼在起名上想方设法隐藏动态规划的数学味

We had a very interesting gentleman in Washington named Wilson. He was secretary of Defense, and he actually had a pathological fear and hatred of the word "research". I'm not using the term lightly; I'm using it precisely. His face would suffuse, he would turn red, and he would get violent if people used the term "research" in his presence. You can imagine how he felt, then, about the term "mathematical" ...... I felt I had to do something to shield Wilson and the Air Force from the fact that I was really doing mathematics inside the RAND Corporation. What title, what name, could I choose?
                              --- Bellman's autobiography, 1984

事实上贝尔曼首次用动态规划这个词是在 1952 年，早于威尔逊去五角大楼任职数月

还能更快吗

根据递推关系有

$$ \begin{align*} \quad & \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix} \begin{bmatrix} F(n-1) \\ F(n) \end{bmatrix} = \begin{bmatrix} F(n) \\ F(n-1) + F(n) \end{bmatrix} = \begin{bmatrix} F(n) \\ F(n+1) \end{bmatrix} \\[4pt] & \qquad \Longrightarrow \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix}^n \begin{bmatrix} F(0) \\ F(1) \end{bmatrix} = \begin{bmatrix} F(n) \\ F(n+1) \end{bmatrix} \end{align*} $$

利用矩阵乘法的结合律

$$ \begin{align*} \quad & \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix}^2 = \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix} \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix}, \quad \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix}^4 = \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix}^2 \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix}^2, \quad \ldots \end{align*} $$

时间复杂度为$O(\lg n)$？

还能更快吗 实现

根据数学归纳法易证

$$ \begin{align*} \quad & F(n) = F(m) F(n-m-1) + F(m+1) F(n-m), \quad \forall n,m \in \Nbb \\[5pt] & \qquad \Longrightarrow \begin{cases} F(2n-1) = F(n-1)^2 + F(n)^2 \\ F(2n) = F(n) (2 \cdot F(n-1) + F(n)) \end{cases} \end{align*} $$

def fibo_rec_faster(n):
    if n == 1:
        return 0, 1
    m = int(n / 2)
    prev_, curr_ = fibo_rec_faster(m)
    prev = prev_**2 + curr_**2
    curr = curr_ * (2 * prev_ + curr_)
    next = prev + curr
    if n & 1:
        return curr, next
    else:
        return prev, curr

还能更快吗 总结

$F(n) = \frac{1}{\sqrt{5}} ( (\frac{1 + \sqrt{5}}{2})^n - (\frac{1 - \sqrt{5}}{2})^n )$，对于充分大的$n$

• $\lg F(n) \approx n \lg (\frac{1 + \sqrt{5}}{2}) \approx 2n /3$
• $\log_{10} F(n) \approx n \log_{10} (\frac{1 + \sqrt{5}}{2}) \approx n /5$

无论二进制还是十进制，打印一个斐波那契数就要$\Theta(n)$

• 在 fibo_rec_faster() 中，每次递归，需要计算若干次两个相邻斐波那契数的乘积，目前最快的两个$n$位数的乘法时间开销为$O(n \lg n)$，因此$T(n) = T(n/2) + O(n \lg n)$，根据主定理情形 3，$T(n) = \Theta(n \lg n)$
• 前面的一重循环中有两个$n$位数相加，加法时间开销为$\Theta(n)$，实际时间复杂度为$\Theta(n^2)$

动态规划未必一定比分治递归快，当分治递归产生的子问题个数更少且没有重复子问题时，动态规划就没有优势了

动态规划

动态规划是一种聪明的、不求解重复子问题的递归方式

通过保存部分子问题的解，以空间换时间，动态规划是一种时空权衡 (time-memory trade-off)

动态规划常用来求解最优化问题

• 最大收益的钢条切割方案
• 最小开销的矩阵连乘方案
• 最长公共子序列
• 最少编辑次数(编辑距离)
• 最长递增子序列
• 最优二叉搜索树
• 最大连续子数组

动态规划 一般步骤

第一阶段：递归 (recursive) 表示

1. 用清晰的语言描述问题
2. 将问题的解用规模更小的子问题的解表示出来

第二阶段：递推 (recurrent) 求解

1. 确定保存子问题结果的数据结构，通常为多维表格
2. 除边界情况外，确定子问题间的依赖关系，这是一个偏序
3. 根据上一步的依赖关系确定子问题的求解顺序，偏序 => 线序
4. 若除最优值外还需最优解本身，在第 1 步的多维表格里维护一些额外信息

钢条切割

假设切割钢条的工序本身没有成本

• 输入：长度$n$的钢条，价格表$p[]$
• 输出：最优分解$n = n_1 + \cdots + n_k$使$r_n = p[n_1] + \cdots + p[n_k]$最大化

如何将$n$的最优分解用$m ~ (< n)$的最优分解表示出来？

钢条切割 最优子结构

设最优分解为$n = n_1 + \cdots + n_k$

任取$\Ical \subseteq \{1, \ldots, k\}$，记$\overline{\Ical} = \{1, \ldots, k\} \setminus \Ical$，设$m = \sum_{i \in \Ical} n_i$，则

• $m = \sum_{i \in \Ical} n_i$是$m$的最优分解
• $n - m = \sum_{i \in \overline{\Ical}} n_i$是$n-m$的最优分解

反设$\sum_{j \in \Jcal} n'_j = m$是$m$的更优分解，则$\sum_{j \in \Jcal} n'_j + \sum_{i \in \overline{\Ical}} n_i$是$n$的更优分解

最优解的每个子部分都是子问题的最优解，称为最优子结构性

钢条切割 分治

双子问题：

• 分：$n = m + (n-m)$
• 治：分别求$m$、$n-m$的最优分解对应的收益$r_m$、$r_{n-m}$

分解为$m$、$n-m$两部分未必是最优的，需遍历所有分解取最大

$$ \begin{align*} \quad r_n = \begin{cases} \max \{ r_1 + r_{n-1}, ~ r_2 + r_{n-2}, ~ \ldots, ~ r_{\lfloor n/2 \rfloor} + r_{\lceil n/2 \rceil}, ~ p[n] \}, & n \le \text{len}(p) \\ \max \{ r_1 + r_{n-1}, ~ r_2 + r_{n-2}, ~ \ldots, ~ r_{\lfloor n/2 \rfloor} + r_{\lceil n/2 \rceil} \}, & n > \text{len}(p) \end{cases} \end{align*} $$

从$(1,n-1)$到$(\lfloor n/2 \rfloor, \lceil n/2 \rceil)$即可遍历所有分解，在$n$不超过$p[]$的长度时还需考虑$p[n]$，它对应的是$r_0 + r_n$，但不能以这样的形式写，否则会出现无限递归

钢条切割 分治

单子问题：将最优分解中的某个数视为一部分，其余为另一部分

• 分：$n = i + (n-i)$
• 治：只对$n-i$求最优分解对应的收益$r_{n-i}$

分解为$i$、$n-i$两部分未必是最优的，需遍历所有分解取最大

$$ \begin{align*} \quad r_n = \begin{cases} \max \{ p[1] + r_{n-1}, ~ \ldots, ~ p[n-1] + r_1, ~ p[n] \}, & n \le \text{len}(p) \\ \max \{ p[1] + r_{n-1}, ~ \ldots, ~ p[\text{len}(p)] + r_{n - \text{len}(p)} \}, & n > \text{len}(p) \end{cases} \end{align*} $$

单子问题产生的子问题更少，效率更优

钢条切割 分治 实现

def cut_rod_rec2(n, p):  # 分治 双子问题
    if n <= 1:
        return p[n]
    else:
        if n < len(p):
            v = p[n]
        else:
            v = -float("inf")
        for i in range(1, math.ceil(n / 2)):  # i = 1, 2, ...
            v = max(v, cut_rod_rec2(i, p) + cut_rod_rec2(n - i, p))
    return v

def cut_rod_rec(n, p):  # 分治 单子问题
    if n <= 1:
        return p[n]
    else:
        v = -float("inf")
        for i in range(1, min(n + 1, len(p))):  # i = 1, 2, ...
            v = max(v, p[i] + cut_rod_rec(n - i, p))
    return v

钢条切割 分治 复杂度

令$T(n)$表示求 cut_rod_rec(n,p) 的递归调用次数，则

$$ \begin{align*} \quad T(n) & = \begin{cases} 1, & n = 1 \\ 1 + \sum_{i=1}^{n-1} T(i), & 1 < n \le \text{len}(p) \\ 1 + \sum_{i=n-\text{len}(p)}^{n-1} T(i), & n > \text{len}(p) \end{cases} \end{align*} $$

$T(n)$的增长速度快于斐波那契数列，至少是指数增长速度

钢条切割 动态规划

带备忘录的自顶向下法 (top-down with memoization)

def cut_rod_dp_memoized(n, p):
    r = [-float("inf")] * (n + 1)  # 备忘录初始化为负无穷
    return cut_rod_dp_memoized_aux(n, r, p)


def cut_rod_dp_memoized_aux(n, r, p):
    if r[n] >= 0:  # 查备忘录 若之前已计算过就直接用
        return r[n]
    if n <= 1:
        r[n] = p[n]
    else:
        for i in range(1, min(n + 1, len(p))):  # i = 1, 2, ...
            r[n] = max(r[n], p[i] + cut_rod_dp_memoized_aux(n - i, r, p))
    return r[n]

每个子问题求解需执行第 12 ~ 13 行的一重 for 循环

规模为$1,2,\ldots,n$的子问题各求解一次，时间复杂度$\Theta(n^2)$

钢条切割 动态规划

自底向上法 (bottom-up method)

• 子问题按规模从小到大的顺序依次求解
• 当求解某个子问题时，它所依赖的更小子问题都已求解完毕

def cut_rod_dp_bottom_up(n, p):
    r = [0] * (n + 1)
    for j in range(1, n + 1):  # 依次求解 r[1], r[2], ...
        v = -float("inf")
        for i in range(1, min(j + 1, len(p))):  # 求解 r[j] 时遍历 i = 1, 2, ...
            v = max(v, p[i] + r[j - i])  # 此时 r[j-1], r[j-2], ... 均已求好
        r[j] = v
    return r[n]

二重 for 循环时间复杂度$\Theta(n^2)$，一维表格空间复杂度$\Theta(n)$

自顶向下 vs. 自底向上

两者有相同的渐进时间复杂度$\Theta(n^2)$

自顶向下

• 无需考虑求子问题的顺序，递归时碰到就求
• 实现容易，在分治递归的基础上加上备忘录查询/保存即可
• 存在递归的开销

自底向上

• 需自己事先指定求子问题的顺序
• 循环实现，没有递归的开销

子问题图

动态规划子问题之间的依赖关系，可用子问题图描述

• 子问题图是一个有向图，每个顶点唯一地对应一个子问题
• 子问题 x 的最优解依赖子问题 y，则有一条从 x 到 y 的有向边
• 递归树中相同子问题的结点合并为一，相关的边从父结点指向子结点

子问题图

自底向上处理子问题的顺序：

• 对任意子问题，先求解它指向的子问题
• 当它指向的所有子问题都求解完才会求解它

算法运行时间等于所有子问题求解时间之和

子问题的求解时间与对应顶点的出度成正比

动态规划的运行时间一般与顶点和边的数量至少呈线性关系

子集和数

输入：$n$个正整数的集合$W = \{ w_1, w_2, \ldots, w_n \}$和正整数$M$

输出：$W$中的和数等于$M$的所有子集

例 1：

• $n=4$，$(w_1, w_2, w_3, w_4) = (7, 11, 13, 24)$，$M = 31$
• 子集：$(7, 11, 13)$、$(7, 24)$

例 2：

• $n=6$，$(w_1, w_2, w_3, w_4, w_5, w_6) = (5, 10, 12, 13, 15, 18)$，$M = 30$
• 子集：$(5, 10, 15)$、$(5, 12, 13)$、$(12,18)$

子集和数 分治

子问题：$\{ w_i, w_2, \ldots, w_n \}$的子集和为$m$

令$\text{ss}[i, m]$表示子问题是否有解，$\text{ss}[1, M]$即为所求目标

• 若$w[i]$不是和为$m$的加数之一，则$\text{ss}[i, m] = \text{ss}[i+1, m]$
• 若$w[i]$是和为$m$的加数之一，则$\text{ss}[i, m] = \text{ss}[i+1, m-w[i]]$

综上，$\text{ss}[i, m] = \text{ss}[i+1, m] \vee \text{ss}[i+1, m-w[i]]$

$i$从 1 开始每增大 1，产生两个新的子问题，子问题总数为$2^n$

子集和数 动态规划

子问题：$\{ w_i, w_2, \ldots, w_n \}$的子集和为$m$

令$\text{ss}[i, m]$表示子问题是否有解，$\text{ss}[1, M]$即为所求目标

$\text{ss}[i, m] = \text{ss}[i+1, m] \vee \text{ss}[i+1, m-w[i]]$

边界：

• 若$m = 0$，$\text{ss}[i, m]$恒为真
• 若$m < 0$或$i>n$，$\text{ss}[i, m]$恒为假

优化：

• 若$m < w[i]$，则$w[i]$必然不是加数之一，$\text{ss}[i, m] = \text{ss}[i+1, m]$

子集和数 动态规划

子问题：$\{ w_i, w_2, \ldots, w_n \}$的子集和为$m$

令$\text{ss}[i, m]$表示子问题是否有解，$\text{ss}[1, M]$即为所求目标

$$ \begin{align*} \quad \text{ss}[i, m] = \begin{cases} \text{True}, & m = 0 \\ \text{False}, & m < 0 ~ \text{or} ~ i > n \\ \text{ss}[i+1, m], & m < w_i \\ \text{ss}[i+1, m] \vee \text{ss}[i+1, m-w_i], & \ow \end{cases} \end{align*} $$

• 表$\text{ss}$的大小为$n \times (M+1)$，元素为布尔值
• 表第$0$列全部为真
• 第$i$行依赖第$i+1$行的元素，填表时从最后一行往前填

子集和数 实现

令$\text{ss}[i, m]$表示子问题是否有解

$$ \begin{align*} \quad \text{ss}[i, m] = \begin{cases} \text{True}, & m = 0 \\ \text{False}, & m < 0 ~ \text{or} ~ i > n \\ \text{ss}[i+1, m], & m < w_i \\ \text{ss}[i+1, m] \vee \text{ss}[i+1, m-w_i], & \ow \end{cases} \end{align*} $$

def subset_sum_dp(w, M):
    n = len(w)
    ss = np.zeros((n, M + 1), dtype="bool")
    ss[:, 0] = True  # 第0列全为真
    if w[-1] <= M:
        ss[-1, w[-1]] = True  # 填最后一行
    for i in range(n - 2, 0, -1):  # 从下往上填
        for m in range(1, M + 1):  # 每行从左往右填
            if m < w[i]:  # 如果目标和小于w[i]
                ss[i, m] = ss[i + 1, m]
            else:
                ss[i, m] = ss[i + 1, m] or ss[i + 1, m - w[i]]
    return ss

子集和数 打印解

def print_sol(ss, w, i, m, M, list):
    if w[i] == m:  # 如果当前目标和等于w[i] 则找到一个解
        for j in list:
            print("%d + " % j, end="")
        print("%d = %d" % (w[i], M))
    else:
        if ss[i + 1, m - w[i]]:  # 如果w[i]是加数之一
            list.append(w[i])  # 入栈
            print_sol(ss, w, i + 1, m - w[i], M, list)
            list.pop()  # 出栈
        if ss[i + 1, m]:
            print_sol(ss, w, i + 1, m, M, list)


w, M = [0, 5, 10, 12, 13, 15, 18], 30
ss = subset_sum_dp(w, M)
print(ss[1, M])
print_sol(ss, w, 1, M, M, [])
------------------------------------------
True
5 + 10 + 15 = 30
5 + 12 + 13 = 30
12 + 18 = 30

矩阵连乘

设$\Av \in \Rbb^{p \times q}$、$\Bv \in \Rbb^{q \times r}$，则$\Cv = \Av \Bv \in \Rbb^{p \times r}$

计算$c_{ij} = \sum_{1 \le k \le q} a_{ik} b_{kj}$需做$q$次标量乘法

计算$\Cv$需$pqr$次标量乘法

设$\Av_1 \in \Rbb^{10 \times 100}$、$\Av_2 \in \Rbb^{100 \times 5}$、$\Av_3 \in \Rbb^{5 \times 50}$

• 计算$(\Av_1 \Av_2) \Av_3$共$10 \times 100 \times 5 + 10 \times 5 \times 50 = 7,500$次标量乘法
• 计算$\Av_1 (\Av_2 \Av_3)$共$100 \times 5 \times 50 + 10 \times 100 \times 50 = 75,000$次标量乘法

矩阵连乘时，怎么两两相乘才能使总的代价最小呢？

矩阵连乘

输入：$n$个矩阵，$\Av_1 \in \Rbb^{p_0 \times p_1}, \Av_2 \in \Rbb^{p_1 \times p_2}, \ldots, \Av_n \in \Rbb^{p_{n-1} \times p_n}$
输出：计算代价最小的连乘顺序对应的加括号方案

令$P(n)$表示$n$个矩阵连乘时总的括号化方案的数量

$$ \begin{align*} \quad P(n) = \begin{cases} 1, & n = 1 \\ \sum_{k=1}^{n-1} P(k) P(n-k), & n > 1 \end{cases} \end{align*} $$

可以证明$P(n) = \binom{2n}{n} / (n+1) = \Omega(4^n / n^{3/2}) $，穷举法不可行！

$P(n)$称为卡特兰数 (Catalan number)

1. 给定进栈顺序的$n$个元素的不同出栈序列个数
2. $n+1$个叶子结点构成的不同(国际)满二叉树的个数
3. $2n$个高矮不同的人排成两队，每队$n$个人，每排都是从低到高排，且第二排的第$i$个人比第一排中第$i$个人高，总的排队方式

矩阵连乘 最优子结构

子问题：$\Av_i \cdots \Av_j$的最优加括号方案

用上划线表示加括号方案$\overline{\Av_i \cdots \Av_j}$，$\star$表示最优

假设$\Av_i \cdots \Av_j$的最优加括号方案在$\Av_k$和$\Av_{k+1}$之间分开：

$$ \begin{align*} \quad \overline{\Av_i \cdots \Av_j}^\star = (\overline{\Av_i \cdots \Av_k}) (\overline{\Av_{k+1} \cdots \Av_j}) \end{align*} $$

则必然有

• $\overline{\Av_i \cdots \Av_k} = \overline{\Av_i \cdots \Av_k}^\star$，否则$(\overline{\Av_i \cdots \Av_k}^\star)(\overline{\Av_{k+1} \cdots \Av_j})$优于$\overline{\Av_i \cdots \Av_j}^\star$
• 同理可得$\overline{\Av_{k+1} \cdots \Av_j} = \overline{\Av_{k+1} \cdots \Av_j}^\star$

问题：如何确定最优的分开位置$k$呢？

矩阵连乘 动态规划

$\Av_i \in \Rbb^{p_{i-1} \times p_i}$、$\Av_k \in \Rbb^{p_{k-1} \times p_k}$、$\Av_{k+1} \in \Rbb^{p_k \times p_{k+1}}$、$\Av_j \in \Rbb^{p_{j-1} \times p_j}$

$\Av_i \cdots \Av_k \in \Rbb^{p_{i-1} \times p_k}$、$\Av_{k+1} \cdots \Av_j \in \Rbb^{p_k \times p_j}$

令$m[i,j]$表示计算$\Av_i \cdots \Av_j$所需的最少标量乘法次数

$$ \begin{align*} \quad m[i,j] = \begin{cases} 0, & n = 1 \\ \min_{i \le k < j} \{ m[i,k] + m[k+1,j] + p_{i-1} p_k p_j \}, & n > 1 \end{cases} \end{align*} $$

采用自底向上法，子问题$\Av_i \cdots \Av_j$的矩阵个数为$l = 2 \rightarrow n$

• $l = 2$：$\Av_1 \Av_2$、$\Av_2 \Av_3$、……、$\Av_{n-1} \Av_n$
• $l = 3$：$\Av_1 \Av_2 \Av_3$、$\Av_2 \Av_3 \Av_4$、……、$\Av_{n-2} \Av_{n-1} \Av_n$
• $\qquad \vdots$
• $l = n-1$：$\Av_1 \Av_2 \cdots \Av_{n-1}$、$\Av_2 \Av_3 \cdots \Av_n$
• $l = n$：$\Av_1 \Av_2 \cdots \Av_{n-1} \Av_n$

矩阵连乘 动态规划

$m[i,j] = \begin{cases} 0, & n = 1 \\ \min_{i \le k < j} \{ m[i,k] + m[k+1,j] + p_{i-1} p_k p_j \}, & n > 1 \end{cases}$

$s[i,j]$表示$\Av_i \cdots \Av_j$的最优分开位置，用来重构最优解

def matrix_chain(p, n):
    m, s = np.zeros((n + 1, n + 1)), np.zeros((n + 1, n + 1))
    for l in range(2, n + 1):           # 子问题长度 l = 2 -> n
        for i in range(1, n - l + 2):   # 从第i个矩阵开始
            j = i + l - 1               # 到第j个矩阵结束
            m[i, j] = float("inf")      # 初始化为无穷大
            for k in range(i, j):       # 遍历最优分开位置
                cost = m[i, k] + m[k + 1, j] + p[i - 1] * p[k] * p[j]
                if cost < m[i, j]:
                    m[i, j] = cost
                    s[i, j] = k
    return m.astype(int), s.astype(int)

三重 for 循环时间复杂度$\Theta(n^3)$，二维表格空间复杂度$\Theta(n^2)$

矩阵连乘 例子

• $m[i,j]$表示计算$\Av_i \cdots \Av_j$所需的最少标量乘法次数
• $s[i,j]$表示$\Av_i \cdots \Av_j$的最优分开位置，用来重构最优解

采用自底向上法，子问题$\Av_i \cdots \Av_j$的矩阵个数为$l = 1 \rightarrow 6$

$$ \begin{align*} \quad l=1: m[1,1] & = m[2,2] = \cdots = m[6,6] = 0, ~ s[i,i] = null \\[5pt] l=2: m[1,2] & = 30 \times 35 \times 15 = 15750, ~ s[1,2] = 1 \\ m[2,3] & = 35 \times 15 \times 5 = 2625, ~ s[2,3] = 2 \\ m[3,4] & = 15 \times 5 \times 10 = 750, ~ s[3,4] = 3 \\ m[4,5] & = 5 \times 10 \times 20 = 1000, ~ s[4,5] = 4 \\ m[5,6] & = 10 \times 20 \times 25 = 5000, ~ s[5,6] = 5 \end{align*} $$

矩阵连乘 例子 填表

递推关系$m[i,j] = \min_{i \le k < j} \{ m[i,k] + m[k+1,j] + p_{i-1} p_k p_j \}$

$$ \begin{align*} \quad m[1,1] & = m[2,2] = \cdots = m[6,6] = 0 \\ m[1,2] & = 15750, ~ m[2,3] = 2625, ~ m[3,4] = 750 \\ m[4,5] & = 1000, ~ m[5,6] = 5000 \\[5pt] m[1,3] & = \min \{m[1,1] + m[2,3] + p_0 p_1 p_3, ~ m[1,2] + m[3,3] + p_0 p_2 p_3 \} \\ & = \min \{ 2625 + 30 \times 35 \times 5, ~ 15750 + 30 \times 15 \times 5 \} \\ & = \min \{ \class{blue}{2625 + 5250}, ~ 15750 + 2250 \} = 7875, ~ s[1,3] = 1 \\ m[2,4] & = \min \{m[2,2] + m[3,4] + p_1 p_2 p_4, ~ m[2,3] + m[4,4] + p_1 p_3 p_4 \} \\ & = \min \{ 750 + 35 \times 15 \times 10, ~ 2625 + 35 \times 5 \times 10 \} \\ & = \min \{ 750 + 5250, ~ \class{blue}{2625 + 1750} \} = 4375, ~ s[2,4] = 3 \end{align*} $$

矩阵连乘 例子 填表

递推关系$m[i,j] = \min_{i \le k < j} \{ m[i,k] + m[k+1,j] + p_{i-1} p_k p_j \}$

$$ \begin{align*} \quad m[1,1] & = m[2,2] = \cdots = m[6,6] = 0 \\ m[1,2] & = 15750, ~ m[2,3] = 2625, ~ m[3,4] = 750 \\ m[4,5] & = 1000, ~ m[5,6] = 5000 \\ m[1,3] & = 7875, ~ m[2,4] = 4375, ~ m[3,5] = 2500, ~ m[4,6] = 3500 \\[5pt] m[1,4] & = \min \{m[1,1] + m[2,4] + p_0 p_1 p_4, ~ m[1,2] + m[3,4] + p_0 p_2 p_4, \\ & \qquad \qquad m[1,3] + m[4,4] + p_0 p_3 p_4 \} \\ & = \min \{ 4375 + 30 \times 35 \times 10, ~ 15750 + 750 + 30 \times 15 \times 10, \\ & \qquad \qquad 7875 + 30 \times 5 \times 10 \} \\ & = \min \{ 14875, ~ 21000, ~ \class{blue}{9375} \} = 9375, ~ s[1,4] = 3 \end{align*} $$

矩阵连乘 例子 填表

递推关系$m[i,j] = \min_{i \le k < j} \{ m[i,k] + m[k+1,j] + p_{i-1} p_k p_j \}$

$m[1,6] = 15,125$即为$\Av_1 \cdots \Av_6$的最优加括号方案的乘法次数

矩阵连乘 构造最优解

$s[i,j]$表示$\Av_i \cdots \Av_j$的最优分开位置

• $s[1,6]=3$，故$\Av_1 \cdots \Av_6$在$\Av_3$后分开
• $s[1,3]=1$，故$\Av_1 \cdots \Av_3$在$\Av_1$后分开
• $s[4,6]=5$，故$\Av_4 \cdots \Av_6$在$\Av_5$后分开

最优方案：$(\Av_1 (\Av_2 \Av_3)) ((\Av_4 \Av_5) \Av_6)$

def print_sol(s, i, j):
    if i == j:
        print("A%d" % (i), end='')
    else:
        print("(", end='')
        print_sol(s, i, s[i, j])
        print_sol(s, s[i, j] + 1, j)
        print(")", end='')

矩阵连乘 实现

p = [30, 35, 15, 5, 10, 20, 25]
n = len(p) - 1  # 矩阵个数
m, s = matrix_chain(p, n)
print(m[1, n])
print(m[1:, 1:])
print(s[1:, 1:])
print_sol(s, 1, n)
------------------------------------------
15125
[[    0 15750  7875  9375 11875 15125]
 [    0     0  2625  4375  7125 10500]
 [    0     0     0   750  2500  5375]
 [    0     0     0     0  1000  3500]
 [    0     0     0     0     0  5000]
 [    0     0     0     0     0     0]]
[[0 1 1 3 3 3]
 [0 0 2 3 3 3]
 [0 0 0 3 3 3]
 [0 0 0 0 4 5]
 [0 0 0 0 0 5]

小结

第一阶段：递归 (recursive) 表示

1. 用清晰的语言描述问题
2. 将问题的解用规模更小的子问题的解表示出来

第二阶段：递推 (recurrent) 求解

1. 确定保存子问题结果的数据结构，通常为多维表格
2. 除边界情况外，确定子问题间的依赖关系，这是一个偏序
3. 根据上一步的依赖关系确定子问题的求解顺序，偏序 => 线序
4. 若除最优值外还需最优解本身，在第 1 步的多维表格里维护一些额外信息

对最优化问题，递归表示中第二步就是在判断最优子结构性

• 最优钢条切割方案的某一部分也是最优的
• 最优矩阵连乘方案的某一部分也是最优的

最优子结构性成立

有向图上的无权最短路径具有最优子结构性

设$u \ne v$，$u$到$v$的最短路径为$p$，记为$u \overset{p}{\rightsquigarrow} v$

路径$p$上的中间结点$w$将$p$分为$p_1$、$p_2$两部分，$u \overset{p_1}{\rightsquigarrow} w \overset{p_2}{\rightsquigarrow} v$

• $p_1$是从$u$到$w$的最短路径，否则设$p'_1$更短，则$u \overset{p'_1}{\rightsquigarrow} w \overset{p_2}{\rightsquigarrow} v$也更短
• $p_2$是从$w$到$v$的最短路径，否则设$p'_2$更短，则$u \overset{p_1}{\rightsquigarrow} w \overset{p'_2}{\rightsquigarrow} v$也更短

最优子结构性均可用上述反证法来证，书上称剪切-粘贴技术

1. 假设子问题的解不是其自身的最优解，最优解另有其解
2. 将该解剪掉，并将最优解粘进来，从而得到原问题的一个更优解
3. 这与最初的解是原问题的最优解矛盾

最优子结构性不成立

有向图上的无权最长简单(无环)路径不具有最优子结构性

• ①=>②=>③ 是 ① 到 ③ 的一条最长路径
• ①=>② 不是 ① 到 ② 的最长路径 (①=>④=>③=>②)
• ②=>③ 不是 ② 到 ③ 的最长路径 (②=>①=>④=>③)
• ①=>② 的最长路径连上 ②=>③ 的最长路径不是简单路径

区别的根源：子问题无关性，同一个原问题的子问题相互独立

最短路径的两个子问题相互独立，$p_1$、$p_2$除$w$外没有公共点，否则设公共点为$x$，则$u \rightsquigarrow x \rightsquigarrow w \rightsquigarrow x \rightsquigarrow v$可优化为$u \rightsquigarrow x \rightsquigarrow v$

最长路径的两个子问题不独立，子问题 ① 到 ② 的最长路径经过的点，在子问题 ② 到 ③ 的最长路径中不能再经过，否则构成环，换言之子问题 ① 到 ② 的求解影响了子问题 ② 到 ③ 的求解

LCS 最优子结构

序列$X = \langle x_1, x_2, \ldots, x_m \rangle$，序列$Y = \langle y_1, y_2, \ldots, y_n \rangle$

记$X$的前$i$个元素构成的前缀子序列为$X_i = \langle x_1, x_2, \ldots, x_i \rangle$

设序列$Z = \langle z_1, z_2, \ldots, z_k \rangle$是$X$和$Y$的任意 LCS

可否考虑序列的第一个元素和后缀？

LCS 最优子结构

① 若$x_m = y_n$，则$z_k = x_m = y_n$且$Z_{k-1}$是$X_{m-1}$和$Y_{n-1}$的 LCS

反设$z_k \ne x_m = y_n$，将$x_m$接到$Z$的末尾可以得到一个长为$k+1$的公共子序列，这与$Z$的最优性矛盾

已证$z_k = x_m = y_n$，故$Z_{k-1}$是$X_{m-1}$和$Y_{n-1}$的公共子序列，若其不是最长，则另有一个长度大于$k-1$的公共子序列，将$x_m$接到其末尾可以得到$X$和$Y$的一个长度大于$k$的公共子序列，矛盾

LCS 最优子结构

② 若$x_m \ne y_n$，则$z_k \ne x_m$蕴含$Z$是$X_{m-1}$和$Y$的一个 LCS

若$z_k \ne x_m$，$Z$显然是$X_{m-1}$和$Y$的公共子序列，若其不是最长，则另有一个长度大于$k$的子序列，该子序列也是$X$和$Y$的子序列，这与$Z$的最优性矛盾

③ 若$x_m \ne y_n$，则$z_k \ne y_n$蕴含$Z$是$X$和$Y_{n-1}$的一个 LCS

LCS 递推关系

• 若$x_m = y_n$，则$z_k = x_m = y_n$且$Z_{k-1}$是$X_{m-1}$和$Y_{n-1}$的 LCS
• 若$x_m \ne y_n$，则$z_k \ne x_m$蕴含$Z$是$X_{m-1}$和$Y$的一个 LCS
• 若$x_m \ne y_n$，则$z_k \ne y_n$蕴含$Z$是$X$和$Y_{n-1}$的一个 LCS

记$c[i,j]$为$X_i$和$Y_j$的 LCS 的长度：

$$ \begin{align*} \quad c[i,j] = \begin{cases} 0, & i = 0 ~ \text{或} ~ j = 0 \\ c[i-1,j-1]+1, & i,j > 0 ~ \text{且} ~ x_i = y_j \\ \max \{ c[i,j-1], ~ c[i-1,j] \}, & i,j > 0 ~ \text{且} ~ x_i \ne y_j \end{cases} \end{align*} $$

LCS 实现

$c[i,j]$为$X_i$和$Y_j$的 LCS 的长度，$b[i,j]$记录$c[i,j]$的计算情况

$$ \begin{align*} \quad c[i,j] = \begin{cases} 0, & i = 0 ~ \text{或} ~ j = 0 \\ c[i-1,j-1]+1, & i,j > 0 ~ \text{且} ~ x_i = y_j \\ \max \{ c[i,j-1], ~ c[i-1,j] \}, & i,j > 0 ~ \text{且} ~ x_i \ne y_j \end{cases} \end{align*} $$

def lcs(X, Y):
    m, n = len(X), len(Y)
    c, b = np.zeros((m + 1, n + 1), dtype="int"), np.empty((m, n), dtype="str")
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            if X[i - 1] == Y[j - 1]:  # 如果当前两个子串的最后一个字符相同
                c[i, j] = c[i - 1, j - 1] + 1
                b[i - 1, j - 1] = "↖"
            elif c[i - 1, j] >= c[i, j - 1]:
                c[i, j] = c[i - 1, j]
                b[i - 1, j - 1] = "↑"
            else:
                c[i, j] = c[i, j - 1]
                b[i - 1, j - 1] = "←"
    return c, b

LCS 例子

$X = \langle A, B, C, B, D, A, B \rangle$，长度 7
$Y = \langle B, D, C, A, B, A \rangle$，长度 6

$c[i,j]$为$X_i$和$Y_j$的 LCS 的长度
$b[i,j]$记录$c[i,j]$的计算情况

for i in range(1, m + 1):
    for j in range(1, n + 1):
        if X[i - 1] == Y[j - 1]:
            c[i, j] = c[i - 1, j - 1] + 1
            b[i - 1, j - 1] = '↖'
        elif c[i - 1, j] >= c[i, j - 1]:
            c[i, j] = c[i - 1, j]
            b[i - 1, j - 1] = '↑'
        else:
            c[i, j] = c[i, j - 1]
            b[i - 1, j - 1] = '←'

二重 for 循环时间复杂度$\Theta(n^2)$，二维表格空间复杂度$\Theta(n^2)$

LCS 例子

$X = \langle A, B, C, B, D, A, B \rangle$，长度 7
$Y = \langle B, D, C, A, B, A \rangle$，长度 6

$c[i,j]$为$X_i$和$Y_j$的 LCS 的长度
$b[i,j]$记录$c[i,j]$的计算情况

打印最长公共子序列

def restore_lcs(i, j, LCS):
    if i == -1 or j == -1:
        return
    if b[i, j] == "↖":
        restore_lcs(i - 1, j - 1, LCS)
        LCS.append(X[i])
    elif b[i, j] == "↑":
        restore_lcs(i - 1, j, LCS)
    else:
        restore_lcs(i, j - 1, LCS)

编辑距离

将一个序列转换成另一个所需的最少的插入/删除/替换次数

$X = \langle A, B, C, B, D, A, B \rangle$、$Y = \langle B, D, C, A, B, A \rangle$

$3$次删除、$2$次插入

$2$次删除、$2$次替换、$1$次插入

LCS 的长度与不带替换的编辑距离是相关的，从前一个例子可以看出

编辑距离

$X$、$Y$各插入一些空格变成等长，然后以某种方式对齐，产生的最少的不相同的列数就是编辑距离

最优子结构性：将$X$、$Y$按编辑距离的方式对齐后，去掉最后一列，剩余列的对齐方式依然是前缀的编辑距离

记$c[i,j]$为$X_i$和$Y_j$的编辑距离

最后一列分 3 种情况：

• 删除：$x_m$对空格，$c[i,j] = c[i-1,j] + 1$
• 插入：空格对$y_n$，$c[i,j] = c[i,j-1] + 1$
• 替换：$x_m$对$y_n$，$c[i,j] = c[i-1,j-1] + \Ibb(x_m \neq y_n)$，其中$\Ibb(\cdot)$为指示函数

编辑距离 递推关系

记$c[i,j]$为$X_i$和$Y_j$的编辑距离

最后一列分 3 种情况：

• 删除：$x_m$对空格，$c[i,j] = c[i-1,j] + 1$
• 插入：空格对$y_n$，$c[i,j] = c[i,j-1] + 1$
• 替换：$x_m$对$y_n$，$c[i,j] = c[i-1,j-1] + \Ibb(x_m \neq y_n)$，其中$\Ibb(\cdot)$为指示函数

3 种情况取最小即为最优

$$ \begin{align*} \quad c[i,j] = \begin{cases} i, & j = 0 \\ j, & i = 0 \\ \min \left\{ \begin{array}{l} c[i-1,j] + 1, \\ c[i,j-1] + 1, \\ c[i-1,j-1] + \Ibb(x_m \neq y_n), \end{array} \right\} & \ow \end{cases} \end{align*} $$

编辑距离 实现

def edit_dis(X, Y):
    m, n = len(X), len(Y)
    c, b = np.empty((m + 1, n + 1), dtype='int'), np.empty((m, n), dtype='str')
    for i in range(m + 1):
        c[i, 0] = i
    for j in range(n + 1):
        c[0, j] = j
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            a = [c[i - 1, j - 1] + int(X[i - 1] != Y[j - 1]), c[i, j - 1] + 1, c[i - 1, j] + 1]
            index = np.argmin(a)
            c[i, j] = a[index]
            if index == 0:
                b[i - 1, j - 1] = '↖'
            elif index == 1:
                b[i - 1, j - 1] = '←'
            else:
                b[i - 1, j - 1] = '↑'
    return c, b

编辑距离 实现

def restore_sol(i, j, b, l):
    if i == -1:
        for k in range(j + 1):
            l.append([' ', Y[j]])
        return
    if j == -1:
        for k in range(i + 1):
            l.append([X[i], ' '])
        return
    if b[i, j] == '↖':
        restore_sol(i - 1, j - 1, b, l)
        l.append([X[i], Y[j]])
    elif b[i, j] == '↑':
        restore_sol(i - 1, j, b, l)
        l.append([X[i], ' '])
    else:
        restore_sol(i, j - 1, b, l)
        l.append([' ', Y[j]])


X = ['A', 'B', 'C', 'B', 'D', 'A', 'B']
Y = ['B', 'D', 'C', 'A', 'B', 'A']
c, b = edit_dis(X, Y)
l = []
restore_sol(len(X) - 1, len(Y) - 1, b, l)
l = np.transpose(l)

print(c)
print(b)
print(c[len(X) - 1, len(Y) - 1])
print(l)
----------------------------------------
[[0 1 2 3 4 5 6]
 [1 1 2 3 3 4 5]
 [2 1 2 3 4 3 4]
 [3 2 2 2 3 4 4]
 [4 3 3 3 3 3 4]
 [5 4 3 4 4 4 4]
 [6 5 4 4 4 5 4]
 [7 6 5 5 5 4 5]]
[['↖' '↖' '↖' '↖' '←' '↖']
 ['↖' '↖' '↖' '↖' '↖' '←']
 ['↑' '↖' '↖' '←' '←' '↖']
 ['↖' '↖' '↖' '↖' '↖' '←']
 ['↑' '↖' '↖' '↖' '↖' '↖']
 ['↑' '↑' '↖' '↖' '↖' '↖']
 ['↖' '↑' '↖' '↖' '↖' '←']]
5
[['A' 'B' 'C' 'B' 'D' 'A' 'B' ' ']
 [' ' 'B' ' ' 'D' 'C' 'A' 'B' 'A']]

最长递增子序列

输入：元素各不相同的序列$X$

输出：$X$的最长递增子序列 (longest increasing subsequence, LIS)

最优子结构性：设以$X[i]$结尾的 LIS 为$X[i_1], \ldots, X[i_k], X[i]$

则$X[i_1], \ldots, X[i_k]$必然是以$X[i_k]$结尾的 LIS

记$d[i]$为以$X[i]$结尾的 LIS 的长度，显然所有$d[i]$的最小值是 1

$d[0] = 1$，以$X[0]$结尾的 LIS 就是$X[0]$

$d[1]$根据$X[1]$是否可以接在$X[0]$后面分两种情况：

• 若$X[0] < X[1]$，则$d[1] = 2$，此时 LIS 就是$X[0, 1]$
• 若$X[0] > X[1]$，则$d[1] = 1$，此时 LIS 就是$X[1]$

LIS 递推关系式

$X[i]$接在以$X[0], X[1], \ldots, X[i-1]$结尾的哪个 LIS 后面？

在所有可接的 LIS 后面选一个最长的，若无则为$1$：

$$ \begin{align*} \quad d[i] = \begin{cases} 1, & i = 0 \\ \max_{j < i, X[j] < X[i]} \{ 1, d[j] + 1 \}, & i > 0 \\ \end{cases} \end{align*} $$

def lis_dp(X, n):
    d, b = [1] * n, [-1] * n
    for i in range(1, n):
        for j in range(i):
            if X[j] < X[i] and d[i] < d[j] + 1:  # X[i]可接在X[j]后面
                d[i] = d[j] + 1
                b[i] = j  # 更新前一个元素的索引
    return d, b

二重 for 循环时间复杂度$\Theta(n^2)$，一维表格空间复杂度$\Theta(n)$

LIS 例子

$$ \begin{align*} \quad d[i] = \begin{cases} 1, & i = 0 \\ \max_{j < i, X[j] < X[i]} \{ 1, d[j] + 1 \}, & i > 0 \\ \end{cases} \end{align*} $$

def restore_lis_dp(X, b, index, LIS):
    if index == -1:
        return
    restore_lis_dp(X, b, b[index], LIS)
    LIS.append(X[index])

LIS 动态规划 改进

维护表$e[]$，初始为空

$e[j]$始终记录当前长度为$j$的递增子序列的末位元素最小值

以$X = [2, 8, 4, 9, 1, 6, 7, 3, 0, 5]$为例，最终$e = [0, 3, 5, 7]$

$e$的长度为 4 代表 LIS 的长度为$4$

$e[-1] = 7$代表 LIS 末位元素为$7$，$[2, 4, 6, 7]$是一个 LIS

$e$的元素是单调递增的，这是提升速度的关键

LIS 动态规划 改进

$e[j]$记录长度为$j$的递增子序列的末位元素最小值，初始为空

LIS 动态规划 改进

在任意时刻$e$都是严格单调递增的，若$i<j$且$e[i] > e[j]$，则有

$$ \begin{align*} \quad \text{LIS}_i:~ \underbrace{\square ~ \square ~ \cdots ~ \square ~ e[i]}_{\text{len}~=~i}, \quad \text{LIS}_j:~ \underbrace{\square ~ \square ~ \square ~ \square ~ \cdots ~ \square ~ \square ~ e[j]}_{\text{len}~=~j} \end{align*} $$

去掉$\text{LIS}_j$前$j-i$个元素可得长度$i$、末位元素$< e[i]$的递增序列

在遍历$X$时，对于$X[i]$

• 若$X[i] > e[-1]$，其可接在当前 LIS 后产生更长的 LIS，$e = [e, [X[i]]$
• 若$e[j-1] < X[i] < e[j]$，其可接在长度为$j-1$的递增子序列后产生长度为$j$末位元素更小的递增子序列，$e[j] = X[i]$，确定$j$可用二分查找

二分查找时间复杂度为$O(\lg n)$，总时间复杂度为$O(n \lg n)$

LIS 动态规划 改进

def lis_dp_plus(X, n):
    e = []
    b = [-1] * n  # b[i]记录i在LIS中的前一个元素的值
    for i in range(n):
        j = bisect_left(e, X[i])  # 在e中对X[i]进行二分查找
        if j == len(e):  # X[i] > e[-1]
            if j > 0:
                b[i] = e[-1]
            e.append(X[i])  # 将其接在e后面表示找到了更长的递增子序列
        else:  # e[j-1] < X[i] < e[j]
            if j > 0:
                b[i] = e[j - 1]
            e[j] = X[i]  # 将e[j]改为X[i] 可以改进现有长度为j的递增子序列
    return e, b


def restore_lis_dp_plus(X, n, e, b):
    dict = {key: value for key, value in zip(X, range(n))}  # 倒排索引字典 X[i]: i
    pre = e[-1]  # e的最后一个元素是LIS的最后一个元素
    LIS = []
    restore_lis_dp_plus_aux(dict, b, pre, LIS)  # 从最后一个元素开始 向前将LIS构造出来
    return LIS


def restore_lis_dp_plus_aux(dict, b, pre, LIS):
    if pre < 0:
        return
    restore_lis_dp_plus_aux(dict, b, b[dict[pre]], LIS)
    LIS.append(pre)

语言翻译

对每个单词，在字典里找该词的翻译，有的能找到，有的找不到

Daenerys Stormborn of House Targaryen, the First of Her Name, Queen of the Andals, the Rhoynar and the First Men, Queen of Meereen, Protector of the Realm, Lady Regnant of the Seven Kingdoms, Mother of Dragons, Khaleesi of the Great Grass Sea, the Unburnt, Breaker of Chains.

二叉搜索树

方法：以单词作关键字建二叉搜索树

目标：尽快找到单词，期望搜索时间最少

• 第$i$层的关键字比较$i$次可以找到
• 左子树的所有元素比根结点小
• 右子树的所有元素比根结点大
• 左、右子树也是二叉搜索树

本页15个结点的二叉树根据字典序而来

二叉搜索树

频繁被查询的单词靠近根可减少期望比较次数

$n$个关键字产生$n+1$个区间$d_i = (k_i, k_{i+1})$

• 不在词典中的词会落入其中之一的区间
• $n+1$个区间对应二叉搜索树的叶节点的子结点

这些区间称为伪关键字

最优二叉搜索树

输入：

• $n$个排好序的关键字$k_1 < k_2 < \cdots < k_n$及其概率$\Pbb [x = k_i] = p_i$
• $n+1$个伪关键字$d_0, d_1, \ldots, d_n$及其概率$\Pbb [x \in d_i] = q_i$

输出：以$k_1, k_2, \ldots, k_n$为内部结点、$d_0, d_1, \ldots, d_n$为叶子结点的最优二叉搜索树 (optimal binary search tree, OBST)，期望搜索时间最少

$$ \begin{align*} \quad \Ebb & [\text{搜索时间}] = \sum_{i=1}^n (\dep (k_i) + 1) \cdot p_i + \sum_{i=0}^n (\dep (d_i) + 1) \cdot q_i \\ & \overset{\sum_{i=1}^n p_i + \sum_{i=0}^n q_i = 1}{=} 1 + \sum_{i=1}^n \dep (k_i) \cdot p_i + \sum_{i=0}^n \dep (d_i) \cdot q_i \end{align*} $$

OBST 例子

有$n = 5$个关键字，$6$个伪关键字

OBST 例子

OBST 最优子结构性

如果$\Tcal$是关于关键字$k_1, \ldots, k_n$和伪关键字$d_0, \ldots, d_n$的 OBST，则$\Tcal$中仅包含$k_i, \ldots, k_j$和$d_{i-1}, \ldots, d_j$的子树$\Tcal'$必然也是 OBST

反设$\Tcal'$不是最优，则存在更优的$\Tcal''$，将$\Tcal$中的$\Tcal'$替换成$\Tcal''$可以得到比$\Tcal$更优的 OBST

子问题：仅包含$k_i, \ldots, k_j$和$d_{i-1}, \ldots, d_j$的 OBST

• 边界情况：$j=i-1$，只有一个伪关键字
• 非边界情况：设最优根结点是$k_r$
• 左子树是仅包含$k_i, \ldots, k_{r-1}$的 OBST
• 右子树是仅包含$k_{r+1}, \ldots, k_j$的 OBST

OBST 递推关系

子问题：仅包含$k_i, \ldots, k_j$和$d_{i-1}, \ldots, d_j$的 OBST

设$e[i,j]$是子问题 OBST 的期望搜索代价

• 边界情况：只有一个伪关键字$q_{i-1}$，$e[i,j] = q_{i-1}$
• 非边界情况：设最优根结点是$k_r$

$$ \begin{align*} \quad e[i,j] & = \overbrace{p_r}^{\text{根}} + \overbrace{e[i,r-1] + w (i,r-1)}^{\text{左子树}} + \overbrace{e[r+1,j] + w (r+1,j)}^{\text{右子树}} \\ & = e[i,r-1] + e[r+1,j] + w(i,j) \end{align*} $$

其中左、右子树分别$+ w (i,r-1)$、$+w (r+1,j)$是因为它们成为$k_r$的子结点时，所有结点深度加$1$

$w(i,j) = \sum_{l=i}^j p_l + \sum_{l=i-1}^j q_l$是结点概率和

OBST 实现

设$root[i,j]$是仅包含$k_i, \ldots, k_j$和$d_{i-1}, \ldots, d_j$的 OBST 的根节点

$$ \begin{align*} \quad e[i,j] = \begin{cases} q_{i-1}, & j = i - 1 \\ \min_{i \le r \le j} ~ \{ e[i,r-1] + e[r+1,j] + w(i,j) \}, & i \le j \end{cases} \end{align*} $$

def optimal_bst(p, q):
    n = len(p)
    e = np.full((n + 2, n + 1), float("inf"))
    w = np.zeros((n + 2, n + 1))
    root = np.zeros((n + 1, n + 1), dtype="int")
    for i in range(1, n + 2):
        e[i, i - 1] = w[i, i - 1] = q[i - 1]  # 边界情况
    for l in range(1, n + 1):  # 子问题 l = 1 -> n
        for i in range(1, n - l + 2):  # 从第i个关键字
            j = i + l - 1  # 到第j个关键字
            w[i, j] = w[i, j - 1] + p[j - 1] + q[j]  # 填写w表
            for r in range(i, j + 1):  # 遍历寻找最优根节点
                t = e[i, r - 1] + e[r + 1, j] + w[i, j]
                if t < e[i, j]:
                    e[i, j] = t  # 更新期望搜索代价
                    root[i, j] = r  # 更新最优根节点
    return e, w, root

三重 for 循环时间复杂度$\Theta(n^3)$，二维表格空间复杂度$\Theta(n^2)$

OBST 例子 填表

$$ \begin{align*} w(i,j) & = q_{i-1} + p_i + q_i + \cdots + p_j + q_j \\[5px] w(1,0) & = q_0 = 0.05, ~ w(2,1) = q_1 = 0.1, ~ w(3,2) = q_2 = 0.05 \\ w(4,3) & = q_3 = 0.05, ~ w(5,4) = q_4 = 0.05, ~ w(6,5) = q_5 = 0.1 \\[5px] w(1,1) & = q_0 + p_1 + q_1 = 0.05 + 0.15 + 0.1 = 0.3, ~ \ldots, ~ w(5,5) = 0.35 \\[5px] w(1,2) & = q_0 + p_1 + q_1 + p_2 + q_2 = 0.45, ~ \ldots, ~ w(4,5) = 0.5 \\[5px] w(1,3) & = q_0 + p_1 + q_1 + p_2 + q_2 + p_3 + q_3 = 0.55, ~ \ldots, ~ w(3,5) = 0.6 \\[5px] w(1,4) & = q_0 + p_1 + q_1 + p_2 + q_2 + p_3 + q_3 + p_4 + q_4 = 0.7, ~ w(2,5) = 0.8 \\[5px] w(1,5) & = 1 \end{align*} $$

OBST 例子 填表

$$ \begin{align*} e[i,j] & = \begin{cases} q_{i-1}, & j = i - 1 \\ \min_{i \le r \le j} ~ \{ e[i,r-1] + e[r+1,j] + w(i,j) \}, & i \le j \end{cases} \\[5px] e[1,0] & = q_0 = 0.05, ~ e[2,1] = q_1 = 0.1, ~ e[3,2] = q_3 = 0.05 \\ e[4,3] & = q_4 = 0.05, ~ e[5,4] = q_5 = 0.05, ~ e[6,5] = q_6 = 0.1 \\[5px] e[1,1] & = e[1,0] + e[2,1] + w(1,1) = 0.05 + 0.1 + 0.3 = 0.45, ~ root[1,1] = 1 \\ e[2,2] & = e[2,1] + e[3,2] + w(2,2) = 0.1 + 0.05 + 0.25 = 0.4, ~ root[2,2] = 2 \\ e[3,3] & = e[3,2] + e[4,3] + w(3,3) = 0.05 + 0.05 + 0.15 = 0.25, ~ root[3,3] = 3 \\ e[4,4] & = e[4,3] + e[5,4] + w(4,4) = 0.05 + 0.05 + 0.2 = 0.3, ~ root[4,4] = 4 \\ e[5,5] & = e[5,4] + e[6,5] + w(5,5) = 0.05 + 0.1 + 0.35 = 0.5, ~ root[5,5] = 5 \\[5px] e[1,2] & = \min \{ e[1,0] + e[2,2], ~ e[1,1] + e[3,2]\} + w(1,2) \end{align*} $$

OBST 例子 填表

$$ \begin{align*} e[i,j] & = \begin{cases} q_{i-1}, & j = i - 1 \\ \min_{i \le r \le j} ~ \{ e[i,r-1] + e[r+1,j] + w(i,j) \}, & i \le j \end{cases} \\[5px] e[1,2] & = \min \{ e[1,0] + e[2,2], ~ e[1,1] + e[3,2]\} + w(1,2) \\ & = \min \{ \class{blue}{0.05 + 0.4}, ~ 0.45 + 0.05 \} + 0.45 = 0.9, ~ root[1,2] = 1 \\ e[2,3] & = \min \{ e[2,1] + e[3,3], ~ e[2,2] + e[4,3]\} + w(2,3) \\ e[3,4] & = \min \{ e[3,2] + e[4,4], ~ e[3,3] + e[5,4]\} + w(3,4) \\ e[4,5] & = \min \{ e[4,3] + e[5,5], ~ e[4,4] + e[6,5]\} + w(4,5) \\[5px] e[1,3] & = \min \{ e[1,0] + e[2,3], ~ e[1,1] + e[3,3], ~ e[1,2] + e[4,3] \} + w(1,3) \\ e[2,4] & = \min \{ e[2,1] + e[3,4], ~ e[2,2] + e[4,4], ~ e[2,3] + e[5,4] \} + w(2,4) \\ & \qquad \vdots \end{align*} $$

OBST 例子 填表

再看最大子数组

跨越中点的最大子数组涉及求以$A[mid]$作结尾的最大子数组

左子问题递归涉及求以$A[mid/2]$作结尾的最大子数组

对任意$i \le mid/2$，求和$A[i] + \cdots + A[mid/2]$会被重复计算

分治法求最大子数组的时间复杂度$\Theta(n \lg n)$

再看最大子数组

子问题$i$：求以$A[i]$作结尾的最大子数组$A[j, \ldots, i]$的和

将这$n$个子问题求解完，其中和最大的就是原问题的最大子数组

最优子结构性：设$A[j, \ldots, i-1]$是以$A[i]$作结尾的最大子数组，其中$j < i$，则$A[j, \ldots, i-1]$也是以$A[i-1]$作结尾的最大子数组

反设$A[j', \ldots, i-1]$是以$A[i-1]$作结尾的最大子数组，$j' \ne j$，则$A[j', \ldots, i]$是更大的以$A[i]$作结尾的最大子数组

再看最大子数组

最优子结构性：设$A[j, \ldots, i-1]$是以$A[i]$作结尾的最大子数组，其中$j < i$，则$A[j, \ldots, i-1]$也是以$A[i-1]$作结尾的最大子数组

设$dp[i]$是以$A[i]$作结尾的最大子数组的和，两种情况：

1. 若$dp[i-1] \ge 0$，直接把$A[i]$接在$A[j, \ldots, i-1]$的后面即可
2. 若$dp[i - 1] < 0$，加上$A[j, \ldots, i-1]$反而更小了，不如从$A[i]$重新开始

$$ \begin{align*} \quad dp[i] = \begin{cases} dp[i-1] + A[i], & dp[i-1] \ge 0 \Longleftrightarrow j < i \\ A[i], & dp[i-1] < 0 \Longleftrightarrow j = i \end{cases} \end{align*} $$

再看最大子数组

def find_max_subarray_dp2():
    dp = [0] * n                    # dp[i]为以A[i]作结尾的最大子数组的和
    s = [0] * n                     # s[i]为以A[i]作结尾的最大子数组的起始索引
    dp[0], s[0] = A[0], 0
    for i in range(1, n):
        if dp[i-1] >= 0:
            dp[i] = dp[i-1] + A[i]  # A[i]接在A[j,...,i-1]后面
            s[i] = s[i-1]           # 继承其起始索引
        else:
            dp[i] = A[i]            # A[i]不接在A[j,...,i-1]后面 另起炉灶
            s[i] = i                # 起始索引就是当前位置
    max_index = np.argmax(dp)       # 遍历dp获取最大元的索引
    return s[max_index], max_index, dp[max_index]