DP的五类优化(2) - 快速幂，四边形不等式

在上一章中，我们介绍了基于单调队列和二进制DP的优化。
今天我们来看另外3类，斜率优化，四边形不等式，快速幂优化。

斐波那契数列

一般大学的DP课，都会从这个有名的数列讲起。通常会给你们演示的递归写法，发现在算接近40的菲波那切项的时候就长时间返回不出值了。这种做法被证明是指数级的复杂度。随后便开始讲解递归过程中，比如F(K) 在很多递归树的分支里都被展开进行了重复计算。如果我们可以保存一个已经算好的结果，之后相同K的计算其实是可以复用之前这个算好的结果的。这就是记忆化搜索，也称自顶向下的DP。这种做法可以把原来的指数级别的时间复杂度给优化到线性的。

其实这个数列还可以更快，这里就要用到矩阵乘法的思想了。在讲这个之前，我们先来介绍一下快速幂是什么？

我们来看一道LEETCODE的题目

https://leetcode.com/problems/powx-n/

我们在算X 的 N次方时，最基本的做法是X 乘以N次。其实我们也可以用二进制的思想来用LOG N 的时间把它算出来。
比如我们算4次方，我们可以用X^2 的结果直接再平方。同理算8次方的话，我们可以先把x^2 给算好，接下来只要算x^2 的四次方了。然后我们把x^4算好，只要算它的平方了。
那么如果是奇数怎么办，我们可以把当前的值预先乘进答案里来解决。比如三次方我们发现是奇数，我们可以先把X 乘一次到答案，然后再算X^2即可。
所以会有如下代码

public double myPow(double x, int n) {
    if (x == 0) return 0;
    if (n == Integer.MIN_VALUE) return (1 / x) * myPow(1 / x, Integer.MAX_VALUE);
    if (n < 0) return myPow(1 / x, -n);
    double res = 1, p = x;
    while (n > 0) {
        if (n % 2 == 1) res *= p; // 是奇数，把答案先乘进结果
        p = p * p; // 把x 的基础 变成 x^2
        n >>= 1; // 之后只需要求原来的一半次幂
    }
    return res;
}

上述是快速幂的基本思想。这里假设小伙伴们已经知道了矩阵乘法是如何做的。以及1个 M * N 的矩阵乘以一个 N * K的矩阵结果是 M * K的矩阵。如果不知道的，可以看我的这篇博客或上网查阅资料。
博客里也介绍了只有方阵才有矩阵的幂。

那么矩阵乘法快速幂的DP优化的核心思想如下：
一组DP状态，其实等价于一个向量。而DP状态的转移方程，可以是对一个向量做变形的矩阵。那么本质上从1个向量到另一个状态的向量，是可以通过一个矩阵来做到。矩阵具有结合律，我们可以先对右半部分矩阵用快速幂得到一个终极的变形矩阵，再乘以向量，就可以把O(N)的计算优化到 O (LOG (N))
第一次接触这个思想的小伙伴一定会觉得非常陌生，不过我们就拿斐波那契数列来下手。
我们可以知道斐波那契的递推公式为 dp[i] = dp[i-1] + dp[i - 2]; 那么每一个新的数的计算依赖于前2个，所以我们结果可以构建这么一个向量为【dp[n], dp[n-1]】
那么怎么转移呢，其实就是找用什么样的矩阵和这个向量做乘法后，可以让N ++
dp[n + 1] = dp[n] * 1 + dp[n - 1] *1; dp[n] = dp[n] * 1 + dp[n - 1] * 0;
我们可以发现，只需要用【[1,1],[1,0]】这个矩阵对向量【dp[n], dp[n-1]】做乘法即可得到【dp[n+1], dp[n]】。

image.png

那么有了上述公式，如果要从 dp[1], dp[2] 求到 dp[n], dp[n-1] 中间需要有N-2个同样的变形矩阵的乘法。

image.png

综上我们可以实现如下代码
https://leetcode.com/problems/fibonacci-number

public int fib(int N) {
    if (N == 0) return 0;
    if (N <= 2) return 1;
    int[][] dp = {{1, 1}}; // dp[2], dp[1]
    int[][] ma = {{1,1},{1,0}};
    N -= 2;
    while (N > 1) {
        if ((N & 1) == 1) dp = mul(dp, ma);
        ma = mul(ma, ma);
        N >>= 1;
    }
    dp = mul(init, ma);
    return dp[0][0]; // dp[n]
}
int[][] mul(int[][] a, int[][] b) { // 矩阵乘法 (m * n)  X (n * k) = m * k 
    int m = a.length, n = a[0].length, k = b[0].length;
    int[][] c = new int[m][k];
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < k; j++) {
            for (int p = 0; p < n; p++) {
                c[i][j] += a[i][p] * b[p][j];
            }
        }
    }
    return c;
}

另一道LEETCODE

https://leetcode.com/problems/knight-dialer/
这道题主要是讲马子跳跃法，然后再一个键盘上可能打出多少不同的数字，在跳N步之后。比如跳1步，那么就是10. 因为第一个键可以按在任何一个位置。跳2步，是20.
这道题可以直接TOP DOWN去做，比如我的目标是最后在1这个位置，跳完K步。那么能跳到1 的，只有

image.png

我们这里定K=2,也就是求跳完2步的数量
所以就有 dfs(1, K=2） = dfs (6, K=1) +dfs(8, K=1)
K = 1 是递归出口，返回1就好了。
有了这个思路，我们其实只要把每个点可以由哪些点跳过来的信息保存好。
就可以搜索了。

NEIGHBORS_MAP = {
  0: (4, 6),
  1: (6, 8),
  2: (7, 9),
  3: (4, 8),
  4: (3, 9, 0),
  5: tuple(), # 5 has no neighbors
  6: (1, 7, 0),
  7: (2, 6),
  8: (1, 3),
  9: (2, 4),
}

因为到达一个点之后K步，和之前怎么跳过来的是不相关的。所以我们可以人为无论前面怎么跳，你现在跳到了M的数字，并且还有K步，余下到1的可能性都是不变的。那么就可以引入记忆化搜索来避免重复的递归展开计算。这样时间复杂度就是状态数量 * 每次递归函数要做的操作。状态数量是 10 * K（K为跳的步数）
其实这道题就解决了。
其实我们可以发现这道题也是当前的状态是通过上一步的状态，根据固定的公式去转移的，最后是去求个数，我们就可以使用矩阵来为向量做变换的思想把它优化到LOG (N)。
上面这个邻居信息表的含义其实就是dp[i][1] = dp[i-1][4] + dp[i-1][6];
那么我们把需要的位置给设置成系数1，不需要的位置设置成系数0，上面的MAP等价于下面的矩阵

NEIGHBORS_MAP = {
  0: (0, 0, 0, 0, 1, 0, 1, 0, 0, 0),
  1: (0, 0, 0, 0, 0, 0, 1, 0, 1, 0),
  2: (0, 0, 0, 0, 0, 0, 0, 1, 0, 1),
  3: (0, 0, 0, 1, 0, 0, 0, 0, 1, 0),
  4: (1, 0, 0, 1, 0, 0, 0, 0, 0, 1),
  5: (0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
  6: (1, 1, 0, 0, 0, 0, 0, 1, 0, 0),
  7: (0, 0, 1, 0, 0, 0, 1, 0, 0, 0),
  8: (0, 1, 0, 1, 0, 0, 0, 0, 0 ,0),
  9: (0, 0, 1, 0, 1, 0, 0, 0, 0, 0),
}

接下来的事情似乎就是用快速幂，来求终极变换方案。然后把初始向量和终极矩阵方案直接相乘。然后把1~10的值求和即可。
如果理解的斐波那契那里的代码，下面其实是一样一样的。
注意因为向量是行向量，所以做乘法的时候，是和矩阵的列去乘，所以上面的MAP，再转矩阵的时候，应该从行映射到矩阵的列。比如上面的MAP的第一行其实等价下面矩阵的第一列。当然你定义初始向量为列向量，就可以不用做这个变换。

int M = 1000000007;
public int knightDialer(int n) {
    long[][] m = {{0,0,0,0,1,0,1,0,0,0},
                  {0,0,0,0,0,0,1,0,1,0},
                  {0,0,0,0,0,0,0,1,0,1},
                  {0,0,0,0,1,0,0,0,1,0},
                  {1,0,0,1,0,0,0,0,0,1},
                  {0,0,0,0,0,0,0,0,0,0},
                  {1,1,0,0,0,0,0,1,0,0},
                  {0,0,1,0,0,0,1,0,0,0},
                  {0,1,0,1,0,0,0,0,0,0},
                  {0,0,1,0,1,0,0,0,0,0}};
    long[][] res = {{1,1,1,1,1,1,1,1,1,1}};
    n--;
    while (n > 0) {
        if (n % 2 == 1) res = mul(res, m);
        m = mul(m , m);
        n /= 2;
    }
    long sum = 0;
    for (long i : res[0]) {
        sum += i;
    }
    return (int) (sum % M);
}
// A[m][p] * B[p][n] = C[m][n]
long[][] mul(long[][] a, long[][] b) {
    int m = a.length, n = b[0].length, p = a[0].length;
    long[][] res = new long[m][n];
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            for (int k = 0; k < p; k++) {
                res[i][j] += a[i][k] * b[k][j];
                res[i][j] %= M;
            }
        }
    }
    return res;
}

什么样的问题可以用矩阵快速幂优化

我们看到这里发现这里有个状态机转化的思想，我们把它写成矩阵和向量的乘法形式，这类DP都可以使用快速幂; 当然这种题目会要求去求个数，而不是MIN/ MAX
但并不是只要是状态机变化的计数DP都可以用矩阵乘法快速幂。比如经典的DECODE ways，https://leetcode.com/problems/decode-ways/

虽然他是从DP[N-1] 和 DP[N-2] 过来，但是里面涉及到了条件分支，这种题目无法写成矩阵变换的形式。

下面我们再看一道可以用快速幂的思想去解的题，然后看我们怎么定义不同的状态使得可以用不同的矩阵转换来表示的情况。

https://leetcode.com/problems/student-attendance-record-ii/

题目中要求一个学生最多只能出现1个A，和2个连续的L（也就是说不要求L的总数，只要没有3个连续的L即可）

同时我们发现这道题也是求个数，我们可以之后思考是否可以用快速幂来优化。

我们看怎么来定义DP的状态。首先A和L一定会在状态机里的。不可以定dp[n][i][j] 表示第N天时，这个学生已经连续i天是L了，且历史上发生了j次A

那么根据这个状态，我们可以知道如果i >0 那么其实只能从dp[n-1][i-1][j]转移过来，因为你希望连续2天是L，必然要从连续1天是L转过来。

如果i = 0的时候，可以从其他所有状态转过来，因为只要再结尾加个P或加个A，就可以破坏掉连续i天是L的连续。根据这个思路，我们也可以列出关系MAP。

NEIGHBORS_MAP = {
  k,0,0: (k-1, 1, 0),  (k-1, 2, 0),  (k-1, 0, 0) [最后加P]
  k,1,0: (k-1, 0, 0)[最后加L]
  k,2,0: (k-1, 1, 0)[最后加L]
  k,0,1: (k-1, 1, 0),  (k-1, 2, 0),  (k-1, 0, 0) [最后加A]  (k-1, 1, 1),  (k-1, 2, 1),  (k-1, 0, 1) [最后加P]
  k,1,1: (k-1, 0, 1)[最后加L]
  k,2,1: (k-1, 1, 1)[最后加L]
}

有了上述的MAP，我们就可以很方便的转换成矩阵。这里我们要对I,J 做编码。因为J最多2种取值。所以编码之后的数为 i * 2 + j
那么矩阵就是

0: 1,1,1,0,0,0
1: 1,0,0,0,0,0
2: 0,1,0,0,0,0
3: 1,1,1,1,1,1
4: 0,0,0,1,0,0
5: 0,0,0,0,1,0

有了这些,下面我们来思考初始向量是什么，根据定义，在最开始只有可能是没有A,没有L，所以只有dp[0][0][0] = 1,其他都为0.
最后记得把这6种状态结尾的个数做一个求和即是题目的答案

int M = 1000000007;
public int checkRecord(int n) {
    long[][] m = new long[][]{
        {1,1,0,1,0,0},
        {1,0,1,1,0,0},
        {1,0,0,1,0,0},
        {0,0,0,1,1,0},
        {0,0,0,1,0,1},
        {0,0,0,1,0,0}
    };
    long[][] res = new long[][]{{1,0,0,0,0,0}};
    while (n > 1) {
        if (n % 2 == 1) res = mul(res, m);
        m = mul(m, m);
        n >>= 1;
    }
    res = mul(res, m);
    long sum = 0;
    for (int i = 0; i < 6; i++) sum += res[0][i];
    return (int) (sum % M);
}

我们再来看一个状态的表示方法，之前我们定义的是结尾有多少个L。这里我们可以定义结尾最多可能有多少个L。这样定义的好处是最后不用作那个求和。因为我们的状态是最多有多2个L，所以也包含了1个L和0个L的情况了。为了把A也给不求和，所以我们把状态也转成历史上最多有了多少个A
这样我们最后只要返回dp[n][2][1] 就是所有结果。
那么因为都改为最多，所以第一个变化的就是初始向量，原来除了0,0 其他都不合法。现在因为是最多，也就是L和A可有可无。所以求变得全合法了。

long[][] res = new long[][]{{1,1,1,1,1,1}};

然后状态转移是如何呢，我们知道0,0 现在只能从前一个2,0过来了，不然就破坏了最多的定义。因为只要加一个P，就可以使得结尾最多又恢复到0个L。
1,0 也可以从2,0转过来（通过最后加P）
也可以从（0,0）转过来，通过最后加L。但是不能最后加A，因为当前定义是历史上最多0个A。（注意历史上最多和只看结尾上最多还是有区别的）

所以我们发现任何状态都可以从（0,2）转过来，因为最后都可以加P。

只有当I >0时，可以从（i-1， x)转过来，通过加L, 注意这里的X和上一个状态要一致，因为这里是历史上最多。

当J >0，可以从（2， j -1)转过来, 通过加A。这里前面可以直接取最大值2，因为加了一个A我们就不会让最多2个L的性质不合法，所以可以取2.

那么下面就是写转移矩阵了。因为所有通过加P 都可以从(2,0)转移过来，我们可以看到第二列全是1

0: 0,0,1,0,0,0
1: 1,0,1,0,0,0
2: 0,1,1,0,0,0
3: 0,0,1,0,0,1
4: 0,0,1,1,0,1
5: 0,0,1,0,1,1

下面只要改一下初始矩阵，和变换矩阵，其余代码不用动，最后直接返回即可。

int M = 1000000007;
public int checkRecord(int n) {
    long[][] m = new long[][]{
        {0,1,0,0,0,0},
        {0,0,1,0,0,0},
        {1,1,1,1,1,1},
        {0,0,0,0,1,0},
        {0,0,0,0,0,1},
        {0,0,0,1,1,1}
    };
    long[][] res = new long[][]{{1,1,1,1,1,1}};
    while (n > 1) {
        if (n % 2 == 1) res = mul(res, m);
        m = mul(m, m);
        n >>= 1;
    }
    res = mul(res, m);
    return (int) res[0][5];
}

讲了这么多，我们最后再来总结一下快速幂优化的思想，就是把计数的状态机转换的DP，通过把初始状态表示为初始向量，转移方程表示为变换矩阵。通过矩阵快速幂的方式优化时间复杂度从O(n) 到 O(log(n))的一类技巧。

讲到这里矩阵快速幂优化就要告一段落，再开始新的篇章前，我给你们留一道思考题。

上题中L和A 都是定值，如果L和A，是可变的（假设2者之和不超过10）,你该如何实现LOG N的算法呢？

四边形不等式优化

四边形不等式DP理论非常复杂，编码还是比较简单。
先说下他的由来。我们都知道一个东西叫最优树。还记得我们在学编码时的哈夫曼数吗，因为每个字母的出现频率不一样，所以我们希望频率高的编码尽可能短，就有了哈夫曼树的思想。他就是贪心的去合并权值最小的2个树，最后合到一颗为止。该树即为所求的哈夫曼树。
随后计算机鼻祖高纳德在解决最优二叉搜索树时发明的一个算法，随后姚期智的夫人，做了深入研究，扩展为一般性的DP优化方法。可以把一些时间复杂度O(n^{3)的DP问题优化到O(n}2), 所以这个方法又被成为 (Knuth-Yao speedup theorem)

最优二叉搜索树问题：

现有 n 个排序过的数，记为数组 a。为叙述方便使 a 的下标从 1 开始。已知给定两组概率 P1...PN 和 Q0...QN，Pi 为“每一次搜索的目标正好为 a i的概率， Qi
为“每一次搜索的目标正好在a (i) 和 a (i+1) 之间”的概率，其中设边界值 a (0)
为负无穷，边界值 a (n+1)为正无穷。求根据这些概率组成一个高效的二叉搜索树，使得每次搜索的平均时间最小化。只需要返回该树的最优平均时间，不需要表达或者返回树的结构。

我们来思考下因为二叉搜索树需要保持节点本身有序的特性，所以我们不能像哈夫曼树那样贪心的取2个概率最小的子树去合并，因为会破坏搜索树的特性。其实这里等价于只有相邻的子树可以合并。这样我们可以把最小的子问题给求好，然后依据最小求次小。次小的时候我们需要枚举决策点，然后再所有决策点里找最小。这样的做法是O(n^3)的。
因为要遍历N层，每一层我们要遍历N个窗口，每个窗口我们又要枚举最优决策点。
我们来看下N^3的代码如实写

double calculateOptimalBST(double[] recordProbability, double[] gapProbability) {
    int n = gapProbability.length;
    double[][] dp = new double[n+1][n+1];
    double[][] subTreeProbabilitySum = new double[n+1][n+1];
    for (int i = 1;i <= n; i++) {
        dp[i][i - 1] = gapProbability[i-1];
        subTreeProbabilitySum[i][i - 1] = gapProbability[i-1];
    }
    for (int len = 1; len < n; len++) { // 枚举节点数为LEN的子树的最优解
        for (int i = 1; i < n + 1 - len; i++) { // 滑动每一个窗口i~j
            int j = i + len - 1;
            subTreeProbabilitySum[i][j] =
                    subTreeProbabilitySum[i][j - 1] + recordProbability[j] + gapProbability[j];
            dp[i][j] = Double.MAX_VALUE;
            for (int k = i; k <= j; k++) { // 枚举决策点,K是根节点
                if (dp[i][j] > dp[i][k-1] + dp[k+1][j]) {
                    dp[i][j] = dp[i][k-1] + dp[k+1][j]; // 左右子树的搜索代价各加一层
                }
            }
            dp[i][j] += subTreeProbabilitySum[i][j]; // 有这个树的搜索代价加一层
        }
    }
    return dp[1][n - 1];
}

上面的代码我们可以跑2个简单的例子论证下正确性
比如，只含一个节点的树，他的最优解就是本身，但是左右的GAP 因为在遍历的时候是需要搜到NULL节点才能确定这个值不存在，所以搜的层数都为2.

Assert.assertTrue(0.2 + (0.3 + 0.5) * 2 ==
        calculateOptimalBST(new double[]{-1, 0.2}, new double[]{0.3, 0.5}));

我们再来验证2个节点的情况，加入第一个GAP区间概率很高，我们应该拿左边的节点为根节点更优。

Assert.assertTrue(0.25 + (0.4 + 0.15) * 2 + (0.08 + 0.12) * 3
        == calculateOptimalBST(new double[]{-1, 0.25, 0.15}, new double[]{0.4, 0.08, 0.12}));

我们再来验证2个节点的情况，加入最后一个GAP区间概率很高，我们应该拿右边的节点为根节点更优。

Assert.assertTrue(0.15 + (0.4 + 0.25) * 2 + (0.08 + 0.12) * 3
        == calculateOptimalBST(new double[]{-1, 0.25, 0.15}, new double[]{0.12, 0.08, 0.4}));

下面我们来讲四边形不等式优化。
这个优化的证明过程非常繁琐，我这里只讲技巧，具体证明我给大家一些不错的资料，有兴趣的朋友大家可以根据资料去学习。比如B站的这个视频

这类的优化过程通常是这样的，比如原来的O（N^3）的写法是这样

image.png

优化之后的代码会长这样

image.png

上面我们引入一个关键的S表，他代表我们之前求过的最优决策。
这个决策表会有一些初始值我们需要赋值，之后就是用最优决策来锁定第三层循环的范围。被证明可以让第三层在第二层下的总次数是O （N）的。也就是表面上看是2层循环，但是实际只有O（N）的遍历次数，具体可以通过打CNT，每遍历一次
CNT++来直观感受。

下面我们来分析什么时候可以用四边形不等式。

1、如果上述的w函数同时满足区间包含单调性和四边形不等式性质，那么函数dp也满足四边形不等式性质
我们再定义s(i,j)表示 dp(i,j) 取得最优值时对应的下标（即 i≤k≤j 时，k 处的 dp 值最大，则 s(i,j)=k此时有如下定理

2、假如dp(i,j)满足四边形不等式，那么s(i,j)单调，即 s(i,j)≤s(i,j+1)≤s(i+1,j+1)

所以也就是说只要W函数，有2个性质，我们可以知道S[I,J]单调，那么就可以套模板来优化。也就是第三层循环由原来的I~J，变成S[I][J-1] ~S[I+1][J]

我们来看下什么是包含单调性和四边形不等式性。

区间包含单调性：如果对于任意 a<=b<=c<=d ，均有w(b,c) <= w(a,d) 成立，则称函数 w 对于区间包含关系具有单调性。
四边形不等式：如果对于任意 a <= b <= c <= d ，均有 w(a,c) + w(b,d) <= w(a,d) + w(b,c) 成立，则称函数满足四边形不等式（简记为“交叉小于包含”）。若等号永远成立，则称函数 w 满足四边形恒等式。

我们回过头来看上面最优二叉搜索树的W函数。
本质上是subTreeProbabilitySum，因为加的都是>=0很容易得出大区间一定>=小区间，所以满足区间包含单调性
第二个四边形不等式，有时可以直接证明出来是满足的，有些时候不太好想，我们可以直接对W函数打表，然后验证所有的ABCD，如果是满足的。那么大概率可以用四边形不等式优化。
我们来写个打表函数。

for (int i = 1; i < n; i++) {
    for (int j = i; j < n; j++) {
        for (int k = j; k < n; k++) {
            for (int m = k; m < n; m++) {
                double contain = subTreeProbabilitySum[i][m] + subTreeProbabilitySum[j][k];
                double cross = subTreeProbabilitySum[i][k] + subTreeProbabilitySum[j][m];
                assert contain >= cross;
            }
        }
    }
}

如果没有报错，我们可以尝试用一下四边形不等式优化。

double calculateOptimalBST(double[] recordProbability, double[] gapProbability) {
    int n = gapProbability.length;
    double[][] dp = new double[n+1][n+1];

    double[][] subTreeProbabilitySum = new double[n+1][n+1];
    for (int i = 1;i <= n; i++) {
        dp[i][i - 1] = gapProbability[i-1];
        subTreeProbabilitySum[i][i - 1] = gapProbability[i-1];
    }
    int[][] s = new int[n+1][n+1]; // step 1. 引入决策表
    for (int i = 1; i <= n; i++) // step 4. 给s 赋初始值
        s[i][i - 1] = i;
    for (int len = 1; len < n; len++) {
        for (int i = 1; i < n + 1 - len; i++) {
            int j = i + len - 1;
            subTreeProbabilitySum[i][j] =
                    subTreeProbabilitySum[i][j - 1] + recordProbability[j] + gapProbability[j];
            dp[i][j] = Double.MAX_VALUE;
            int st = s[i][j-1], ed = Math.min(j, s[i+1][j]); // step 3. 用决策表更新搜索范围
            for (int k = st; k <= ed; k++) {
                if (dp[i][j] > dp[i][k-1] + dp[k+1][j]) {
                    dp[i][j] = dp[i][k-1] + dp[k+1][j];
                    s[i][j] = k; // step2. 记录最优决策
                }
            }
            dp[i][j] += subTreeProbabilitySum[i][j];
        }
    }
    return dp[1][n - 1];
}

大致分为4步。

第一步，引入决策表。
第二步，在更新时更新决策表。
第三步，在第三层循环时，用决策表的数据来循环。
第四步，是需要思考的，如何赋初始值。这道题因为最开始算的是DP[I][I]，也就是1的单位，那么对S[I][j]来说，他的前驱和后继s[i][j-1] 和 s[i+1][j] 都是0 的单位，所以要在赋初始值时用s[i][i-1]，其次就是在第一次循环时ed 这个位置因为只有一个长度好遍历，所以这里要加个MIN

我们再来研究下这个优化到底在遍历什么？

首先我们会发现，他是根据长度从小到大在遍历每个窗口。在每层长度中，每个窗口要遍历的范围则是S[i][j-1] 到 S[i+1][j]. 具体一下再LEN=1, I = 1时，其实就是 s[2][1] - s[1][0] 个决策点. 随后I到2了，变成s[3][2] - s[2][1] 个决策点。那么发现前项和后项有正负号可以抵消。所以最终一个LEN的遍历次数就是 s[n-1][n-len+1] - s[len-1][0]。根据S数组的定义，就是决策点，所有决策定都不会超过N，所以对于一个LEN来说，内部2个循环和为s[n-1][n-len+1] - s[len-1][0] <= n.
就证明是O(n ^ 2) 了

我们可以看一道类似的题目
https://www.lintcode.com/problem/stone-game
这道题其实和最优二叉树还是比较像的，区别就是不用考虑GAP区间。其实也就更加简单。
因为不考虑GAP区间，我们甚至可以直接证明COST[I,J] 是满足四边形恒等式的。
cost[a,d] = cost[a,b] + cost[b,c] + cost[c,d]

cost[a,c] = cost[a,b] + cost[b,c]
cost[b,d] = cost[b,c] + cost[c,d]
这样变形之后带入原式，就会发现左右两边相等。
另外一个最优二叉搜索树的DP[I,J] 其实是包含了gap[i-1] ~gap[j] 的。
所以枚举最优决策不让左右最优子树的GAP有重叠需要从dp[i][k-1] 和 dp[k+1][j]来转移。因为他们代表了gap[i-1]~gap[k-1] 和 gap[k]~gap[j]的2颗最优子树。刚好排除了决策节点的全覆盖。
而这道题因为不存在GAP，DP[I,J]的定义也发生了变化，指的是合并stone[i] ~ stone[j], 所以枚举K的时候，是dp[i][k] 和 dp[k+1][j]转移过来，意思是最后是由[i,k]这堆石头和[k+1, j]这堆石头合并。
我们来看下直接用四边形，因为石子这个问题，循环是从长度为2的情况开始（1是终态不用算的），所以在初始化的时候是初始化1，而不是像上一题初始化0.这些都是要注意的细节。下面上代码

public int stoneGame(int[] A) {
    int l = A.length;
    if (l == 0) return 0;
    int[][] dp = new int[l][l];
    int[][] s = new int[l][l];
    int[] presum = new int[l + 1];
    for (int i = 0; i < l; i++) {
        presum[i + 1] = presum[i] + A[i];
        s[i][i] = i;
    }
    for (int len = 2; len <= l; len++) {
        for (int i = 0; i < l - len + 1; i++) {
            int j = i + len - 1;
            dp[i][j] = Integer.MAX_VALUE / 2;
            int st = s[i][j-1], ed = Math.min(j - 1, s[i+1][j]);
            for (int k = st; k <= ed; k++) {
                if (dp[i][k] + dp[k + 1][j] < dp[i][j]) {
                    dp[i][j] = dp[i][k] + dp[k + 1][j];
                    s[i][j] = k;
                }
            }
            dp[i][j] += presum[j + 1] - presum[i];
        }
    }
    return dp[0][l - 1];
}

通过四边形不等式我们又把一道N^3的区间DP问题给优化到了N ^ 2
这道题还有一种n log n的解法，叫GarsiaWachs算法的最优解法, 学有余力的小伙伴可以自行搜索学习

我们再来看一道这周周赛的一个问题

https://leetcode.com/problems/allocate-mailboxes/
这道题暴力的思路是，我们只分析前K个房子，假设此时要建M个邮局。我们已经知道M-1个邮局，前0~K-1个房子下的最优解。我们只要枚举最后一个邮局建造的位置管辖的屋子，然后不管辖的屋子靠前面算过的最优子问题直接获得答案即可得到当前问题的答案。
因为这里最优子问题是由2个维度组成N和K。那么枚举决策的时候还是要根据N来枚举最后一步能管的房子数量。所以这道题DP 是O(N^3）的
DP方程为
dp[i][j]=min(dp[i][j],dp[k][j-1]+w[k+1][i]);
其中dp[i][j]表示前i个村庄建j个邮局的最小距离和，k枚举1到i的所有村庄；w[k+1][i]表示第k+1个村庄到第i个村庄建一个邮局的最小距离和，有一个显然的性质：在某一段区间上建一个邮局，最小距离和为在其中点村庄上建

那么我先来思考怎么快速的把这个W数组给求出来
我们知道如果只有一个房子w[i][i] = 0.
如果有2个房子，我们是取中点建是最优的。所以我们新加了一个房子就是加上他到中点的距离。

int mid=(i+j)/2;
w[i][j] = w[i][j-1]+abs(x[j]-x[mid]);

有读者可能会问，除了最后一个点不算，为什么w[i][j] = w[i][j-1]，中点不是会随着加了一个屋子，而往后移吗？
这个我们可以分类讨论看，如果是从奇数房子增加到偶数房子。中点是不会移动的，所以可以等价。
如果是偶数房子增加到奇数房子，中点是需要往后移动一格。但是原来小于等于中点的房子数量是偶数的一半，我把中点往后移，这一半的房子的距离都要加上中点移动的距离。同时剩下一半的房子的距离都会减去中点移动的距离。因为2半的房子数量相同，所以值还是不变的。综上这个等式是成立的。
那么我们就可以在O(N ^ 2）的时间求完W数组

接下来主要的时间瓶颈就是DP这个O(n * n * K)的复杂度。
这个式子dp[i][j] = min{ dp[i-1][k] + w(k+1,j) | i-1 <= k < j }乍一看和我们之前说到可以用四边形不等式的式子dp[i][j] = min{ dp[i][k] + dp[k+1][j] + w(i,j) | i <= k < j }似乎不太一样，但有一个相似点是他也要枚举最优的决策，这个时候有一个技巧就是，我们把决策表打印出来，看他是不是每行单调递增(允许>=），同时每列也单调递增(允许>=）。如果满足这个性质，大概率这个式子也是满足决策单调性，就可以用四边形不等式的套路进行优化
所以基础代码如下

int inf = Integer.MAX_VALUE / 2;
private int[][] dis(int[] a) {
    int l = a.length;
    int[][] dis = new int[l][l];
    for (int i = 0; i < l; i++) {
        for (int j = i + 1; j < l; j++) {
            dis[i][j] = dis[i][j - 1] + a[j] - a[(j + i)/2];
        }
    }
    return dis;
}
public int minDistance(int[] houses, int k) {
    Arrays.sort(houses);
    int[][] dis = dis(houses);
    int n = houses.length;
// DP I J 表示前J个屋子用了I个邮局的最小距离和
    int[][] dp = new int[k + 1][n];
    int[][] s = new int[n+1][n+1];
    for (int[] i : s) Arrays.fill(i, -1);
    for (int i = 0; i < n; i++) {
        dp[1][i] = dis[0][i];
    }
    for (int l = 2; l <= k; l++) {
        for (int i = l; i < n; i++) {
            dp[l][i] = inf;
            for (int j = l-1; j <= i; j++) { // 枚举最后一个邮局COVER 多少房子
                if (dp[l - 1][j - 1] + dis[j][i] < dp[l][i]) {
                    dp[l][i] = dp[l - 1][j - 1] + dis[j][i];
                    s[l][i] = j; 
                }
            }
        }
    }
    // 验证行单调
    for (int[] i : s) {
        boolean seeingMinusOne = true;
        for (int j = 1; j < i.length; j++) {
            if (seeingMinusOne && i[j-1] != -1) seeingMinusOne = false;
            if (i[j] == -1 && !seeingMinusOne ) i[j] = inf;
            assert i[j] >= i[j-1];

        }
    }
    // 验证列单调
    for (int i = 0; i < n; i++) {
        boolean seeingMinusOne = true;
        for (int j = 1; j < n; j++) {
            if (seeingMinusOne && s[j-1][i] != -1) seeingMinusOne = false;
            if (s[j][i] == -1 && !seeingMinusOne ) s[j][i] = inf;
            assert s[j][i] >= s[j-1][i];
        }
    }
    return dp[k][n - 1];
}

用这个代码去跑一下评测系统，如果出现ASSERTION ERROR，那么就代表决策不具备单调性，如果没出问题，那么我们可以去尝试用四边形不等式优化。

下面我就是要思考这里我们要算的是DP[L][I]，那么更新的决策点就是S[L][I]，这个时候S[L-1][I] 是已经求好了，另外一侧需要S[L][I+1]，那么就需要第二层循环从大到小，再检查一下DP只需要上一层的元素，所以从大到小是没问题的。

因为I都是从I-1开始，那么初始的第三层循环的右侧<=值最大应该是N-1，所以S[L][N] = N - 1.同时因为L是从2开始的，我们需要构建好左边的初始值，S[1][i] = 1. (因为L =2 , 然后J从L-1开始，所以=1)

用四边形不等式优化的代码如下：

int inf = Integer.MAX_VALUE / 2;
private int[][] dis(int[] a) {
    int l = a.length;
    int[][] dis = new int[l][l];
    for (int i = 0; i < l; i++) {
        for (int j = i + 1; j < l; j++) {
            dis[i][j] = dis[i][j - 1] + a[j] - a[(j + i)/2];
        }
    }
    return dis;
}
public int minDistance(int[] houses, int k) {
    Arrays.sort(houses);
    int[][] dis = dis(houses);
    int n = houses.length;
    int[][] dp = new int[k + 1][n];
    int[][] s = new int[n+1][n+1];  // step.1 
    for (int[] i : s) Arrays.fill(i, -1);
    for (int i = 0; i < n; i++) {
        dp[1][i] = dis[0][i];
        s[1][i] = 1; // step 4.
    }
    
    for (int l = 2; l <= k; l++) {
        s[l][n] = n - 1; // step 4.
        for (int i = n - 1; i >= l - 1; i--) {
            dp[l][i] = inf;
            int st = s[l-1][i], ed = s[l][i+1]; // step 3.
            for (int j = st; j <= ed; j++) {
                if (dp[l - 1][j - 1] + dis[j][i] < dp[l][i]) {
                    dp[l][i] = dp[l - 1][j - 1] + dis[j][i];
                    s[l][i] = j; // step 2.
                }
            }
        }
    }
    
    return dp[k][n - 1];
}

到这里我已经把四边形不等式如何优化的思想已经介绍完了，当然四边形不等式的优化还可以运用再一维DP里，不过本文已经相当长了。而且我还有斜率优化也还没写，所以我们1维DP的四边形不等式优化和斜率优化放在下一章讲。因为这2个算法，目前LC还没有对应的题目，所以算是超纲讲授。不过既然都开始写了，就写写完完整整。所以小伙伴们，我们下章见。

总结

这次我们主要学习了矩阵快速幂来优化基于状态机转移的DP计数类问题。原理就是把初始状态设计成向量，转移方程设计为矩阵。转移过程就是向量和矩阵的乘法，然后因为矩阵乘法具有结合律，所以我们可以先算矩阵乘法通过快速幂的方式达到优化效果。

随后是四边形不等式的优化，原理就是在区间类DP中需要枚举最优决策点时，我们可以通过判断代价函数是否满足区间包含单调性，和四边形不等式来得知决策点是否具备单调性，如果具备单调性就可以用4步法来把O(N^3)的复杂度优化至 O (N^2）。

老惯例，给大家留2道思考题。

在矩阵快速幂中，那道同学上课缺席迟到的题目，如果L和A是动态传入，应该如何通过代码来构建转移矩阵呢？

石子合并那道题，如果是求最大COST，应该怎么做呢？

请研究https://leetcode.com/problems/minimum-cost-to-merge-stones/ 怎么做，能否用四边形不等式。

上期思考题时间

在一棵有根多叉树中，如何使用二进制优化，来找最近公共祖先呢？

这道题分为3步。

第一步，同样我们对这颗树的每个节点构建工具包，使得每个节点向上的一步，二步，四步。。。的节点编号都直接被存下来。同时把每个节点的深度给存下来。

随后深度深的那个节点，开始从大工具（步数大的开始跳）开始使用，使用的前提是用完之后，依然没有比深度浅的节点更浅。那么就用，继续换小工具。这样做的目标是使得2个节点在同一深度。

随后如果2个节点是一个点了，那么就直接返回。如果不是，就来到第三步。

第三步就是2个节点使用工具一起往上跳，用该工具的前提是，2个节点用完不会使得他们来到了同一个节点（因为可能跳过头了）；我们的目标要找到最浅的一个不同的父节点，那么他们上面一个就是最近公共祖先。

总共有 n 道题目要抄，编号 0,1,…,n，抄第 i 题要花 ai 分钟。老师要求最多可以连续空的题为K，求消耗时间最少满足老师要求的方案。

首先我们可以在最后加1个时间为0 的题。然后这样就可以用DP[N+1]来得到答案。

DP[N+1]表示的是N+1这道题写了的话花的最小时间。

随后我们就可以知道转移方程就是，因为最多K道空，那么前K道里面必然需要一个题是写的，那么就从这K个DP转移过来，求最小。那么这里也是维护区间最小值，我们可以用单调队列来解决。