动态编程

算法的力量真强大。

这几天在耍一耍 rosalind 上面的题目，遇到了几道所谓“斐波那契兔子”的问题，还有一道求最长公共子串的问题。根据题目给的提示，这类“斐波那契”问题都可以归属于动态编程的范畴。那么动态编程具体指的是什么呢？

动态规划(dynamic programming)是运筹学的一个分支，是求解决策过程(decision process)最优化的数学方法。20世纪50年代初美国数学家R.E.Bellman等人在研究多阶段决策过程(multistep decision process)的优化问题时，提出了著名的最优化原理(principle of optimality)，把多阶段过程转化为一系列单阶段问题，利用各阶段之间的关系，逐个求解，创立了解决这类过程优化问题的新方法——动态规划。

以上是百度百科上的定义，看不懂啊看不懂。不过没关系，经过下面的问题就应该会好懂一点。

以下问题的相关解思路参考于 IBM - developerworks，含示例图片。有关原始实现，可看原文。

典型的斐波那契数列

众所周知，典型的斐波那契数列问题的通解是 F(n) = F(n-1) + F(n-2)，这样一个式子，明显就是可以用递归来进行计算的。比如下面的例子：

1
2
3
4
5
6
7
8
def fibonacci_ord(n):
	def calc(n):
		if n <= 2:
			return 1
		else:
			return calc(n-1) + calc(n-2)
	res = calc(n)
	return res

就上面的代码而言，fibonacci_ord(n) 中的某些数总是会被重复计算。这时候有两个办法，一个是使用额外的列表或者字典存储已经计算后的 f(n)。
另外一个是把 fibonacci_ord(n) 的递归计算改成迭代计算，从最小的开始算起，这里还是需要额外的列表。下面使用第二种方式进行改写一下：

1
2
3
4
5
6
7
8
9
10
def fibonacci_iter(n):
	values = dict()
	values[0] = 0
	values[1] = 1
	for i in range(n+1):
		if i <= 1:
			count = values[i]
		else:
			values[i] = values[i-1] + values[i-2]
	return values[n]

可以看到，就计算 F(n) 的形式而言，仍旧是使用原来的公式，但是把递归改造成使用 for 循环进行。而这就涉及动态编程了。

动态编程

在斐波那契数列这类问题中，问题的解来源于子问题，子问题的解来源于子问题。这类问题的解都可以使用递归，但是里面却有重复计算子问题解的缺陷存在。
在这时，在这类可以使用递归的问题中，因为递归重复解决相同的子问题造成效率低下，则可以采用动态编程进行改写。像上面的使用迭代的方式计算斐波那契数列的思路就是动态编程的工作原理，把递归分解成迭代。动态编程能够解决这类的问题，它主要有下面三个特征：

  1. 每个问题的解都能用递归关系表示。
  2. 用递归方法对这些递归关系的直接实现会造成解决方案效率低下，因为其中包含了对子问题的多次计算。
  3. 一个问题的最优解能够用原始问题的子问题的最优解构造得到。

最长子序列问题

下面我们来看常见的一个问题，最长子序列问题（longest common subsequence, LCS）。对于两个序列，记为 C 和 S。我们把这两个序列的左侧各切除 1 个字符，得到 c1 和 s1. 那么现在，LCS 的解依赖于 3 个方面：

  • (C-c1) 与 S 的解 LCS1
  • (S-s1) 与 C 的解 LCS2
  • (C-c1) 和 (S-s1) 的解 LCS3，如果 c1 与 s1 相同，这里的解就是 LCS3 + 1

LCS 是这三个解之间的最大值。基于这样一个形式，就十分适合于使用动态规划的打分矩阵来进行。当前 cell 的打分依赖于左侧 cell、上侧 cell 和左上 cell 的打分，分别对应那三个方面得到的三个解。对于左上 cell，考虑 +1 的问题。

打分

在得到打分矩阵之后，从得分最大的点开始回溯。如果当前 cell 的打分比左上 cell 打分大 1，表明在该位置处 c1 和 s1 的值相同，添加到最大子序列里面去。
由于这里要回溯到前一个，所以构建打分矩阵的时候要添加连接信息。下面是示例代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
def LCS(C, S):
	# 构建打分矩阵, x=C,y=S
	score = []
	max_score = 0
	for i in range(len(C) + 1): # 含 0 ，即 C 与 S 都是 “” “”
		temp = []
		for j in range(len(S) + 1):
			if j == 0 or i == 0:
				temp.append([0, None])	# None is prev
			else:
				max_score = score[i-1][j-1][0] # 暂定左上
				prev = "lefttop"
				if C[i-1] == S[j-1]:	# c1 和 s1 相同
					max_score = max_score + 1
				if temp[j-1][0] > max_score:	# 此时的 左边的 cell 值还在 temp 中
					max_score = temp[j-1][0]
					prev = "left"
				if score[i-1][j][0] > max_score:
					max_score = score[i-1][j][0]
					prev = "top"
				temp.append([max_score, prev])
		score.append(temp)

	# 回溯
	i = len(C)
	j = len(S)
	shared = []
	while i >= 1 and j >= 1:
		cell = score[i][j]
		# 由于前面默认是左上作为初始max分数
		# 所以有可能出现不是左上但分数多了1
		if cell[0] == score[i-1][j-1][0] + 1 and cell[1] == "lefttop":
			shared.append(S[j-1])
			i -= 1
			j -= 1
		elif cell[1] == "left":
			j -= 1
		else:
			i -= 1
	return "".join(shared)

print(LCS("GCCCTAGCG", "GCGCAATG"))

最长子串问题

最长子串和最长子序列问题是不同的，对于子序列来说，只要某个字符在两个序列中出现，那么就视为是共同子序列，最长子序列就是这些共同子序列的并集。而最长子串则是需要连续的共同子序列组成，其重点在于连续和共同。我们前面已经实现了寻找最长子序列。在此基础上，我们只要稍作更改即可。

试想一下，基于前面的代码，如果某个 cell 回溯的前一个 cell 位于其 “lefttop” 则表明，这两个 cell 所对应的两个序列的字符是相同且连续的。那现在要找到连续最长的，我们只需在打分时，只取两个碱基相等时的打分，其余为 0，这时候得到的分数则是连续共同子串的长度。这样做以后，我们就可以得到得分最高的那个 cell，由那个 cell 进行回溯就可以得到最长子串。

在 rosalind 里面的最长子串是 100 个序列的最长子串，采用这里的代码去跑的话，需要 20s 的时间。这说明，对于 n 个序列的最长公共子串这一问题，有更好的算法。事实也在 rosalind 的 solution 里面有展示，不到 30 行代码，跑完是时间不到 0.1s。不得不说，算法的厉害就在于此。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
def LCSS(C, S):# longest common substring
	# 构建打分矩阵, x=C,y=S
	score = []
	max_score = 0
	max_score_pos = [0, 0]
	for i in range(len(C) + 1): # 注意这里的 i 和 j 的取值范围[0, len(C/S)] 是闭区间
		temp = []
		for j in range(len(S) + 1):
			if j == 0 or i == 0:
				temp.append(0)	# None is prev
			else:
				if C[i-1] == S[j-1]:
					current_score = score[i-1][j-1] + 1
					if current_score > max_score:
						max_score = current_score
						max_score_pos = [i, j]
				else:
					current_score = 0
				temp.append(current_score)
		score.append(temp)

	# 回溯
	i, j = max_score_pos
	shared = []
	while max_score > 0:
		shared.append(C[i-1])	# i-1是因为构建打分矩阵多引入了 0 的初始行
		max_score = score[i-1][j-1]
		i = i - 1
		j = j - 1
	shared.reverse() # 因为是回溯，要翻转
	return "".join(shared)

print(LCSS("GCCCTAGCG", "GCGCAATG"))

序列比对

有了最长子串和最长子序列的问题，自然而然就联系到序列比对中去了。序列比对不就是在 2 条 DNA 序列中寻找子串或子序列吗。有关序列比对这里展示 2 种类似的算法，NeedlemanWunsch 和 SmithWaterman 算法。

NeedlemanWunsch 算法

对于 NeedlemanWunsch 算法，在计算打分矩阵的时候，会根据两个序列字符是否匹配和 gap 添加对应的加分和罚分，以此来寻找最优解。由于这里我们会考虑 gap 和 错配的情况，所以 NeedlemanWunsch 算法也是一种全局比对算法。其基本思路如下：

  • 比对时加入空格 -2 分，错配 -1 分，比对上加 2 分。这里可以自由指定。
  • 对于来自上面的单元格，代表将左侧的字符与空格比对。
  • 对于来自左侧的单元格，代表将上面的字符与空格比对。
  • 对于来自左上侧的单元格，代表与左侧和上面的字符比对（可能匹配也可能不匹配）

所以，根据上面这个思路，打分矩阵的第一行和第一列的分数是逐渐递减的。下面是我的实现。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
def NeedlemanWunsch(C, S, space=-2, match=2, mismatch=-1):
	# 同样是构建矩阵
	score = []
	# 使用 0,1,-1 分别代表lefttop，left，top
	for i in range(len(C) + 1):
		temp = []
		for j in range(len(S) + 1):
			if i == 0 and j == 0:
				temp.append([0, None])			# 起点得分 0
			elif i == 0:
				temp.append([j*space, None])	# 第一行得分为 space * column_index
			elif j == 0:
				temp.append([i*space, None])	# 第一列得分为 space * row_index
			else: # 其余行列
				max_score = score[i-1][j-1][0]
				prev = "lefttop"
				if C[i-1] == S[j-1]:			# 能比对上
					max_score += match
				else:
					max_score += mismatch		# 错配
				if temp[j-1][0] + space > max_score:
					max_score = temp[j-1][0] + space
					prev = "left"
				if score[i-1][j][0] + space > max_score:
					max_score = score[i-1][j][0] + space
					prev = "top"
				temp.append([max_score, prev])
		score.append(temp)

	# 回溯
	i = len(C)
	j = len(S)
	C_string = []
	S_string = []
	while i >= 1 and j >= 1:
		cell = score[i][j]
		prev = cell[1]
		if prev == "lefttop":
			C_string.append(C[i-1])
			S_string.append(S[j-1])
			i -= 1
			j -= 1
		elif prev == "top": # top
			S_string.append("-")	# 向下延伸是在 水平 上添加 gap，即 S
			C_string.append(C[i-1])
			i -= 1
		else:
			C_string.append("-")
			C_string.append(S[j-1])
			j -= 1
	C_string.reverse()
	S_string.reverse()
	return ["".join(C_string), "".join(S_string)]

print(NeedlemanWunsch("GCCCTAGCG", "GCGCAATG"))

Smith-Waterman 算法

在 Smith-Waterman 算法中，不必比对整个序列。两个零长字符串即为得分为 0 的局部比对，这一事实表明在构建局部比对时，不需要使用负分。这样会造成进一步比对所得到的分数低于通过 “重设” 两个零长字符串所能得到的分数。而且，局部比对不需要到达任何一个序列的末端，所以也不需要从右下角开始回溯：可以从得分最高的单元格开始回溯。这一事实表明在构建局部比对时，不需要使用负分。这样会造成进一步比对所得到的分数低于通过 “重设” 两个零长字符串所能得到的分数。而且，局部比对不需要到达任何一个序列的末端，所以也不需要从右下角开始回溯：可以从得分最高的单元格开始回溯。这导致 Smith-Waterman 算法与 Needleman-Wunsch 算法存在着三个区别。

  1. 首先，在初始化阶段，第一行和第一列全填充为 0（而且第一行和第一列的指针均为空）。
  2. 第二，在填充表格时，如果某个得分为负，那么就用 0 代替，只对得分为正的单元格添加返回指针。
  3. 最后，在回溯的时候，从得分最高的单元格开始，回溯到得分为 0 的单元格为止。
    除此之外，回溯的方式与 Needleman-Wunsch 算法完全相同。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def SmithWaterman(C, S, space=-2, match=1, mismatch=-1):
	# 同样是构建矩阵
	score = []
	max_score = 0
	max_score_pos = [0, 0]
	# 使用 0,1,-1 分别代表lefttop，left，top
	for i in range(len(C) + 1):
		temp = []
		for j in range(len(S) + 1):
			if i == 0 or j == 0:
				temp.append([0, None])			# 起点得分 0
			else: # 其余行列
				current_score = score[i-1][j-1][0]
				prev = "lefttop"
				if C[i-1] == S[j-1]:			# 能比对上
					current_score += match
				else:
					current_score += mismatch		# 错配
				# 找到最大的分数
				if temp[j-1][0] + space > current_score:
					current_score = temp[j-1][0] + space
					prev = "left"
				if score[i-1][j][0] + space > current_score:
					current_score = score[i-1][j][0] + space
					prev = "top"
				# 得分小于 0 就重置为0
				if current_score < 0:
					current_score = 0
					prev = None
				# 不为 0，就看是否是当前最大得分
				elif current_score > max_score:
					max_score = current_score
					max_score_pos = [i, j]
				temp.append([current_score, prev])
		score.append(temp)
        
	# 回溯
	i, j = max_score_pos
	C_string = []
	S_string = []
	while i >= 1 and j >= 1:
		cell = score[i][j]
		prev = cell[1]
		if prev == "lefttop":	# lefttop
			C_string.append(C[i-1])
			S_string.append(S[j-1])
			i -= 1
			j -= 1
		elif prev == "top": # top
			S_string.append("-")	# 向下延伸是在 水平 上添加 gap，即 S
			C_string.append(C[i-1])
			i -= 1
		elif:
			C_string.append("-")
			C_string.append(S[j-1])
			j -= 1
	C_string.reverse()
	S_string.reverse()
	return ["".join(C_string), "".join(S_string)]

print(SmithWaterman("GCCCTAGCG", "GCGCAATG"))

参考

Comments

▲