动态规划系列(1): 线性DP

从斐波那契数列说起

能你从来没有接触过动态规划,只是隐约听别人提起过DP是种特别厉害的算法;又或许你已经看过动态规划的相关定义和教程,能够做一些像LCS、数字三角形之类的简单DP题目。

但无论你对DP了解如何,阅读本文都会使你对动态规划有更为清晰的认识。不同于你可能看过的很多资料,本文将会以一种更为崭新的视角,从更为本质而数学化的角度出发,尽力发掘出隐藏在所谓的“状态转移方程”、“最优子结构”等重重迷雾下动态规划的真实一面。

首先要搞明白的一点是:类似递归、二分甚至深度优先搜索等,动态规划是一种策略而不是一种特定的算法;也就是说,动态规划提供一种分析问题的方式。 它通过将问题分解成一个个子问题代表的状态,通过逐渐求解这些子问题,在不同状态之间不断转移,最终计算得出问题的解。也就是说,动态规划的理论性和实践性都比较强,既要求你能够正确地对问题进行建模,又要求你在获得问题的理论分析之后,能够正确构造出解决问题的算法。可以这么说:对动态规划的掌握程度,就是你对问题分析和建模能力的直接体现。

我们先从一个常见的例子说起:斐波那契数列的计算。在学习分治法时,你可能已经见过这个问题,并且还知道这个问题是使用递归的一个不良例子。我们再来回顾一下这个例子:

	// 使用递归计算斐波那契数列的第n项
int Fib(int n) {
	if(n == 0)
		return 0;
	else if(n == 1)
		return 1;
	else
		return Fib(n - 1) + Fib(n - 2);
}

大概计算到n = 40时就要开始等几秒才会出现结果了,随着n不断增大,这个程序的计算速度会迅速降低,最后甚至引发崩溃。原因也很简单:如果我们将每次求解Fib(n)视作一个问题,那么当我们求解这个问题时,我们先要求解子问题Fib(n - 1),在这个过程中我们已经递归求解了子问题Fib(n - 2);而之后我们又把子问题Fib(n - 2)重复求解了一遍!因为无论是Fib(n - 1)还是Fib(n - 2),其求解都涉及更多下属的子子问题,当n很大时这些子子问题重复求解的次数加起来就会十分可观,最终这个算法的时间复杂度将会达到惊人的$O(2^{n})$。

因此,我们说:斐波那契数列的计算问题具有重复子问题性质。分治法面对这类问题是无效的,因为它的目的是将问题划分为互不相交的子问题,之后分别求解。因此,分治法并不考虑是不是做了重复工作;它只是简单地求解给定的子问题。因此,如果我们试图使用分治法来求解具有重复子问题性质的问题,分治法将陷入无休止地对重复子问题的反复求解中。

但是将问题划分成子问题是一个非常有效的简化问题的途径,我们不想也不能抛弃这一方法;那么有什么方法可以改善我们的处境呢?回想上文关于动态规划的描述:我们将每个子问题视作一个状态。而我们可以将不同的状态记录下来。当我们需要知道某个状态的值时,我们先查看之前是否已经记录这一状态:如果有的话,就直接提取;否则计算这一状态的值,之后保存下来。通过这种做法我们便可避免对子问题的重复计算。下面的代码描述了这一思想的实现:

#define MAX 1000
int f[MAX];

memset(f, -1, sizeof(f));
f[0] = 0; f[1] = 1;

int Fib(int n) {
	if(f[n] != -1)
		return f[n];
	else
		return (f[n] = Fib(n - 1) + Fib(n - 2));
}

我们称这一方式为记忆化搜索。顾名思义,它将问题中已经“搜索”过的部分记忆下来,之后每次需要时调用。

但是我们仍觉得,这个问题可以有更符合直觉的解法:我们观察上述解法可以发现,我们在求解一个子问题前,总要先求解比其更小的子问题。在上面的解法中,我们是从最大的总问题出发去求解比其小的问题,这样我们仍避免不了对子问题的计算。如果我们将顺序反过来:从最小的子问题开始递推,逐步由小到大自底向上地解决问题,这样当我们需要求解一个新的子问题时,所用到的子子问题就总是已经解决好的。通常,这种思路更符合我们的直觉,因而也更常用;同时也可以避免记忆化搜索中仍然存在的函数调用的开销。

#define MAX 1000
int f[MAX];

memset(f, -1, sizeof(f));
f[0] = 0; f[1] = 1;

int Fib(int n) {
	for(int i = 2; i <= n; i++)
		f[i] = f[i - 1] + f[i - 2];
	
	return f[i];
}

从这个问题中我们可以看出:动态规划利用问题的可划分性以及子问题之间的相似性来进行归纳,从而降低求解的复杂度。动态规划将一个问题视作若干个重复子问题的逐层递进,并且将每个子问题的求解过程视作一个阶段;只有完成前一个阶段的计算后,才会开始后一个阶段的计算。

但是这里还没有展现出动态规划的重要用途:给出问题的一个最优解(因为问题的最优解通常不止一个)。通常情况下,问题的某个状态可能并不来源于固定的子问题,而是由很多子问题中的一个或一些导出,而这些子问题的解包含足够的信息来决定何种选择可使这个状态达到最优。并且,我们还没有系统地阐明如何使用动态规划方法来对一个问题进行分析。因此,我们来看下一个例子:

初探状态设计:数塔 (HDU.2084)

Problem Description

在讲述DP算法的时候,一个经典的例子就是数塔问题,它是这样描述的:

有如下所示的数塔,要求从顶层走到底层,若每一步只能走到相邻的结点,则经过的结点的数字之和最大是多少?

已经告诉你了,这是个DP的题目,你能AC吗?

Input

输入数据首先包括一个整数C,表示测试实例的个数,每个测试实例的第一行是一个整数N(1 <= N <= 100),表示数塔的高度,接下来用N行数字表示数塔,其中第i行有个i个整数,且所有的整数均在区间[0,99]内。

Output

对于每个测试实例,输出可能得到的最大和,每个实例的输出占一行。

Sample Input

1
5
7
3 8
8 1 0 
2 7 4 4
4 5 2 6 5

Sample Output

30 

分析

初看这道题时,dfs是很容易想到的方法,但是使用dfs遍历整个数塔,其时间复杂度则又会是指数级的;我们需要进行动态决策来缩小此题的状态空间。因此,我们仍然尝试使用动态规划来解决此题。

回想上文关于动态规划的描述:我们首先需要将原问题分解为很多逐层递进的子问题,每个子问题代表了问题的一个状态,而每个子问题的求解过程则构成求解过程的一个阶段。很显然,一个问题的状态表示和阶段划分直接相关。

那么,我们又如何定义数塔问题的状态呢?我们首先要从精确地定义问题,也即明确目标入手。当问题被精确地定义时,子问题的划分也就变得显然,从而问题中状态的定义也就自然显现。

我们还需要从之前的dfs的想法中获取思路。我们目前需要获知的信息有以下三点:

  1. 问题的开始条件(也即边界)在何处?
  2. 问题的结束条件是什么?也即,何时我们停止一次搜寻,获得问题的一个解?
  3. 我们怎样从问题的解空间中获取一个最优解?

在数塔问题中,这三个问题的答案十分显然。我们要从顶层走到底层,因此搜索肯定从顶层开始——这意味着问题的边界是顶层。同样,搜索在底层结束,这意味着一次搜寻也在到达底层时结束,因此问题的结束条件是到达底层。而每次问题的求解都是在一条特定的路径上进行的,这意味着问题求解的每个阶段是一条开始于顶层、结束于底层的路径。最后,我们将遍历这些路径获得的解收集起来,从中取出最优的一个作为最优解。

但是这样的想法仍然并没有良好地定义这个问题。如果将不同的路径做为问题求解的阶段,则这样定义的子问题显然无助于问题的最终解决——我们无法从已经求解完的子问题中获得任何信息来指导下一阶段的计算,也必须逐个遍历完全部的子问题后才能获知问题最终的最优解。但是我们仍然获得了重要信息——问题的边界和结束条件。我们给出这个问题的良好定义时需要用到它们。

现在,我们试图从更为抽象的角度来定义这个问题:我们是否可以摒弃对遍历数塔的具体过程的关注,转而直接记录经过数塔底层某个点的路径集合在这个点上的最优解?**也即,我们模仿前面处理斐波那契数列时的方式,采用一个一维数组F[n]来记录从数塔顶部走到数塔底层的第i个点所经过路径的最大和。**现在,我们暂且将此作为关于数塔问题状态的表示。

相应地,根据重复子问题的原则,我们可以很容易地想到与原问题相似并且大量重复的子问题形式——从数塔顶部走到数塔内第(i, j)点时所经过路径的最大和。显然此时我们需要拓展状态的维度,将存储状态的数组由一维扩展到二维。这样我们就可以给出数塔问题的最终状态表示: $f(i, j)$表示从数塔顶部走到第$i$行第$j$列,所经过路径的最大和

现在我们需要考虑问题的求解具体如何进行:我们已经从考虑遍历数塔的具体细节转为考虑数塔中不同的点代表的(子)问题最优解,但我们仍然需要逐个求解这些子问题。我们怎样才能从求解某些子问题,转移到求解其下一阶段的子问题呢?在获得了数塔问题的状态表示后,这个问题不难考虑:$(i, j)$ 点只能从 $(i-1, j)$ 或 $(i-1, j-1)$ 两个点到达,因此 $f(i,j)$也只能由 $f(i-1, j)$ 和 $f(i-1, j-1)$ 两个状态转移而来。更直接地,假设存放数塔三角形的矩阵为 $A$,我们可以列出状态转移方程: $$ f(i,j) = A(i,j)+max[f(i-1, j-1), f(i-1, j)] $$ 凭借这个方程,结合递推我们便可以在问题的不同状态中逐次推进,最终得出我们需要的最优解。

但是我们此时又面临一个问题:为什么我们在之前定义数塔问题时,通过定义原问题是求“最大和”,直接就能定义其子问题也是求“最大和”?也即:为什么我们能够保证,子问题取到最优,即可以保证其下一阶段的子问题、乃至原问题都取得最优解

我们可以通过反证法证明这一点:在上述状态转移方程中,如果我们不取$max[f(i-1,j-1), f(i-1,j)]$ 而是取两者中最小值,所得的 $f(i,j)$ 一定小于第一种情况所得,这与 $f(i,j)$ 应取得最优解不符,因为此时存在一个更优的解。因此原命题得证。我们将这一性质称为最优子结构性质。它保证了动态规划能够对状态的抽象和子问题的重叠递进起到优化作用。

至此,我们终于分析完了数塔问题。现在我们对上面的分析结果做一个总结:

问题数塔
状态表示$f(i,j)$ 表示从数塔顶部走到第$i$行第$j$列,所经过路径的最大和
阶段划分某条路径的结尾位置$(i,j)$
转移方程$f(i,j) = A(i,j)+max[f(i-1, j-1), f(i-1, j)]$
边界$ f(1,1)=A(1,1) $
目标$\underset{1 \leq i \leq n}{max}[f(n,i)]$

题解

#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;

const int INF = 0x3f3f3f3f;
const int MAXN = 100;

int map[MAXN + 10][MAXN + 10];
int dp[MAXN + 10][MAXN + 10];

int main() {
    ios::sync_with_stdio(false);
    
    int c = 0;
    cin >> c;
    
    while(c--) {
        int n = 0;
        cin >> n;
        
        for(int i = 1; i <= n; i++)
            for(int j = 1; j <= i; j++)
                cin >> map[i][j];
        
        memset(dp, 0, sizeof(dp));
        dp[1][1] = map[1][1];
        
        for(int i = 2; i <= n; i++)
            for(int j = 1; j <= i; j++)
                dp[i][j] = map[i][j] + max(dp[i - 1][j], dp[i - 1][j - 1]);
            
        int ans = 0;
        for(int i = 1; i <= n; i++)
            ans = max(ans, dp[n][i]);
        
        cout << ans << endl;
    }
}

状态设计进阶:最长公共子序列(HDU.1159)

Problem Description

A subsequence of a given sequence is the given sequence with some elements (possible none) left out. Given a sequence X = <x1, x2, …, xm> another sequence Z = <z1, z2, …, zk> is a subsequence of X if there exists a strictly increasing sequence <i1, i2, …, ik> of indices of X such that for all j = 1,2,…,k, xij = zj. For example, Z = <a, b, f, c> is a subsequence of X = <a, b, c, f, b, c> with index sequence <1, 2, 4, 6>. Given two sequences X and Y the problem is to find the length of the maximum-length common subsequence of X and Y.  The program input is from a text file. Each data set in the file contains two strings representing the given sequences. The sequences are separated by any number of white spaces. The input data are correct. For each set of data the program prints on the standard output the length of the maximum-length common subsequence from the beginning of a separate line.

Sample Input

abcfbc abfcab
programming contest 
abcd mnp

Sample Output

4
2
0

分析

不同于上一道数塔问题尚且具有搜索解法,如果我们仍执着于依据题意处理这道题的细节,将会变得非常困难:暴力枚举此题显然不可行,假设字符串X长度为 $n$ ,字符串Z长度为 $m$,则字符串X有子序列 $2^{n}$ 个,字符串Z有子序列 $2^{m}$ 个,此时需要处理的字串对多达 $2^{m+n}$个,显然远远超出任何我们能处理的数据范围。而若不采用暴力枚举的方式,因为最长公共子序列(LCS)逐个分布不均地排列在两个字符串中,使得我们难以找到一个确定的问题开始条件,求解仿佛无从下手。因此,我们仍需借助抽象度更高的动态规划方法来解决这道题。

通过设计数塔问题状态的经验,我们现在试着设计这样一个状态:$f(i,j)$ 表示字符串X的 前缀子串$[1..i]$ (这里,我们令字符串下标从1开始,为了方便后期直接写出边界条件,实际编程中处理也很简单)和字符串Z的前缀子串 $[1..j]$ 内取得的LCS的长度。那么,这样设计的状态是否具有最优子结构性质呢?我们又如何在不同状态之间进行转移呢?

我们再次考虑如何将问题缩小一阶,即考虑 $f(i,j)$ 与 $f(i-1, j-1)$ 等相邻阶段状态的关系。很显然 $X[i]$ 和 $Z[j]$ 的关系在这里成为了重点。

如果 $X[i] = Z[j]$ ,那么状态转移方程很显然会是 $$ f(i, j) = f(i-1, j-1) + 1 $$ 因为 $X[i]$ 或 $Z[j]$ 可以直接接在原有的LCS后方形成新的LCS。

但是,如果 $X[i] \neq Z[j]$ 呢?这时要考虑 $f(i-2, j-2)$ 了吗?先别急。我们令 $LCS(i, j)$ 表示字符串X的 $[1..i]$ 区间和字符串Z的 $[1..j]$ 区间内取得的LCS,则 $LCS(i, j)$ 的最后一项还有可能是 $X[i]$ 或 $Z[j]$ !这是显然的,虽然 $X[i] \neq Z[j]$ ,但是并没有排除这两个字符仍作为 $LCS(i, j)$ 的最后一项的可能性,因此我们不能直接在状态的两个维度上都直接向后缩小一阶,而是应分别考虑两个维度。因此我们可以写出这种情况下的状态转移方程 $$ f(i,j)=max[f(i-1,j),f(i,j-1)] $$ 写出这个方程时我们实际上已经默认这个问题具有最优子结构性质,根据上面的分析过程证明这一性质也十分容易,在这里我们不再详述。

现在,我们对LCS问题的分析结果做一个总结:

问题描述LCS问题
状态表示$f(i, j)$ 表示A的前缀子串 $[1..i]$ 和B的前缀子串 $[1..j]$ 的LCS长度
阶段划分已经处理的前缀在两个字符串中的位置(一个二维坐标)
转移方程$\begin{cases} & f(i,j)=f(i-1,j-1)+1,X[i]=Z[j]\ & f(i,j)=max[f(i-1,j),f(i,j-1)],X[i] \neq Z[j] \end{cases}$
边界$f(i,0)=f(0,j)=0,1 \leq i \leq n, 1 \leq j \leq m$
目标$f(n,m)$

题解

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

const int MAXN = 1000;
char X[MAXN + 10], Z[MAXN + 10];
int f[MAXN + 10][MAXN + 10];

int main() {
    while(scanf("%s %s", X + 1, Z + 1) != EOF) {
        int n = (int)strlen(X + 1);
        int m = (int)strlen(Z + 1);

        memset(f, 0, sizeof(f));
        for(int i = 1; i <= n; i++) {
            for(int j = 1; j <= m; j++) {
                if(X[i] == Z[j]) {
                    f[i][j] = f[i - 1][j - 1] + 1;
                } else {
                    f[i][j] = max(f[i - 1][j], f[i][j - 1]);
                }
            }
        }

        printf("%d\n", f[n][m]);
    }

    return 0;
}

我们还能做些什么

在前面的两道例题中,我们的目标都是输出问题的一个最优解的值而非这个最优解本身。但是如果现在要求我们输出一个LCS的具体内容,我们又该怎么做?

输出一个LCS的内容

我们再次将目光投向LCS问题的状态转移方程: $$ \begin{cases} & f(i,j)=f(i-1,j-1)+1,X[i]=Z[j]\ & f(i,j)=max[f(i-1,j),f(i,j-1)],X[i] \neq Z[j] \end{cases} $$ 可以看到,这里的状态转移过程显得很有“方向性”:$f(i,j)$ 总是从 $f(i-1,j-1)$ , $f(i-1,j)$ , $f(i, j-1)$ 三个状态转移过来,并且当且仅当 $f(i,j)$ 从 $f(i-1,j-1)$ 转移而来时,$X[i]$ 或 $Z[j]$ 才是所求LCS中的一个元素。那么,我们如果在计算的中途将 $f(i,j)$ 具体由哪个状态转移而来记录下来,当LCS的长度计算完成时,我们便可以从 $f(n,m)$ 出发一路追踪下去,最终打印出整个LCS序列。

如果我们将三种状态视作三种“方向”:$f(i-1,j-1)$代表从“左上方”转移而来,$f(i-1,j)$ 代表从“左侧”转移而来,$f(i, j-1)$代表从“上方”转移而来,我们便可以在已有的二维数组 $f$ 内继续维护有关方向的信息,之后输出时根据方向向后追踪并逆序输出结果。

在原有的代码上稍加改动,我们便可以得到下面的程序:

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

#define UL 1
#define L 2
#define U 3

const int MAXN = 1000;

struct Node {
	int status, arrow;
};

char X[MAXN + 10], Z[MAXN + 10];
Node f[MAXN + 10][MAXN + 10];

int LCSlength(char* X, char* Z, Node* f, int n, int m);
void PrintLCS(char* X, char* Z, Node* f, int p, int q);
int main() {
    while(scanf("%s %s", X + 1, Z + 1) != EOF) {
        int n = (int)strlen(X + 1);
        int m = (int)strlen(Z + 1);

        printf("%d\n", LCSlength(X, Z, f, n, m));
	    PrintLCS(X, Z, f, n, m);
	 	putchar();
    }

    return 0;
}

int LCSlength(char* X, char* Z, Node* f, int n, int m) {
	memset(f, 0, sizeof(f));

    for(int i = 1; i <= n; i++) {
        for(int j = 1; j <= m; j++) {
            if(X[i] == Z[j]) {
                f[i][j].status = f[i - 1][j - 1] + 1;
				f[i][j].arrow = UL;
            } else {
                if(f[i - 1][j].status > f[i][j - 1].status) {
					f[i][j].status = f[i - 1][j].status;
					f[i][j].arrow = U;
				} else {
					f[i][j].status = f[i][j - 1].status;
					f[i][j].arrow = L;
				}
            }
        }
    }

	return f[n][m].status;
}

void PrintLCS(char* X, char* Z, Node* f, int p, int q) {
	if(p == 0 || q == 0)
		return;
	
	if(f[p][q].arrow == UL) {
		PrintLCS(X, Z, f, p - 1, q - 1);
		printf("%c ", X[p]);
	} else if (f[p][q].arrow == U) {
		PrintLCS(X, Z, f, p - 1, q);
	} else {
		PrintLCS(X, Z, f, p, q - 1);
	}
}
updatedupdated2023-03-112023-03-11