动态规划之——最长公共子串和矩阵链乘法

1. 最长公共子串

最长公共子串与最长公共子序列有一些类似，只不过这里子串要求是连续的。

这里我们定义 lcs[i][j] 表示以 s[i] 与 t[j] 为末尾元素的最长公共子串长度，那么我们有：

$$lcs[i][j]=\{\begin{array}{ll}lcs[i-1][j-1]+1& \text{如果 }s[i]==t[j]\\ 0& \text{如果 }s[i]!=t[j]\end{array}$$

import numpy as np

s = 'mitcmu'
t = 'mtacmu'


def longest_commom_substr(s, t):

    m = len(s)
    n = len(t)
    lcs = np.zeros((m, n), 'uint8')
    max_lcs = 0

    for i in range(0, n):
        if s[0] == t[i]:
            lcs[0][i] = 1
        else:
            lcs[0][i] = 0

    for i in range(0, m):
        if s[i] == t[0]:
            lcs[i][0] = 1
        else:
            lcs[i][0] = 0

    for i in range(1, m):
        for j in range(1, n):
            if s[i] == t[j]:
                lcs[i][j] = lcs[i-1][j-1] + 1
                max_lcs = max(max_lcs, lcs[i][j])
            else:
                lcs[i][j] = 0

    print(lcs)
    return max_lcs


print(longest_commom_substr(s, t))

2. 矩阵链乘法

假设有一系列矩阵 $<A_1,A_2,\cdots ,A_n>$，计算 $A_1\times A_2\times \cdots \times A_n$ 的最小代价。若矩阵 $A$ 的大小为 $p\times q$，矩阵 $B$ 的大小为 $q\times r$，则 $A\times B$ 的代价是 $O(pqr)$。

由于矩阵乘法满足结合律，所以有：

$$A_1\times A_2\times A_3\times A_4=A_1\times (A_2\times A_3\times A_4)=(A_1\times A_2)\times (A_3\times A_4)=\cdots =A_1\times (A_2\times A_3)\times A_4$$

按照这几种顺序计算需要的代价是不一样的。比如：

我们定义状态 $state[i][j]$ 为计算矩阵 $A_i\times \cdots \times A_j$ 所需的最小代价，那么我们有：

$$state[i][j]=\{\begin{array}{ll}0& \text{如果 }i==j\\ mi{n}_{i⩽k<j}(state[i][k]+state[k+1][j]+p_{i-1}{p}_k{p}_j)& \text{如果 }i<j\end{array}$$

也就是我们把 $A_i\times \cdots \times A_j$ 在中间某个位置划分开，变成 $(A_i\times \cdots \times A_k)\times (A_{k+1}\times \cdots \times A_j)$ 两部分，这样总的计算代价也就等于两部分的代价之和加上将这两部分再乘起来需要的代价。

计算的时候，我们则需要自底向上地更新状态变量，最后逐渐得到问题的解。

算法伪代码如下所示：

import numpy as np

max_num = 2 ** 31 - 1
p = [(10, 100), (100, 5), (5, 50)]
p = [(30, 35), (35, 15), (15, 5), (5, 10), (10, 20), (20, 25)]


def matrix_chain_order(p):

    n = len(p)
    state = max_num * np.ones((n, n), 'int32')
    divide = np.zeros((n, n), 'uint8')

    for i in range(n):
        state[i][i] = 0

    for s in range(1, n):
        for i in range(n-s):
            j = i + s
            for k in range(i, j):
                temp = state[i, k] + state[k+1, j] + p[i][0] * p[k][1] * p[j][1]
                if temp < state[i][j]:
                    state[i][j] = temp
                    divide[i, j] = k

    print(state)
    print(divide)
    return state[0, n-1]


print(matrix_chain_order(p))

上述代码中有 6 个矩阵，运行结果如下，第一个输出为状态矩阵，第二个输出为划分矩阵，最终需要的计算次数为 $state[0][5]=15125$。

由 $divide[0][5]=2$ 可知计算所有 6 个矩阵时我们在矩阵 $A_2$ 处进行划分，由 $divide[0][2]=0$ 可知我们在计算前 3 个矩阵时我们在矩阵 $A_0$ 处进行划分，由 $divide[3][5]=0$ 可知我们在计算后 3 个矩阵时我们在矩阵 $A_4$ 处进行划分，所以最终的乘法顺序为：

$$A_0\times (A_1\times A_2)\times (A_3\times A_4)\times A_5$$

获取更多精彩，请关注「seniusen」!