最大子矩阵和

摘要: 最大子矩阵和问题，有两种处理二维情况的方式。

【对算法，数学，计算机感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：算法题刷刷
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

---

在文章 最大子数组和的三种解法 中，我们详细拆解了最大子数组和这个问题，并且了解到这个问题有三种解法，都非常主流，并且各个方法都有它适用的变种问题。

在文章 环形数组上的最大子数组和 中，我们看了最直接的一个变种：把数组变成环形数组，那道题的重点是怎样把环形数组上的问题转换为普通数组上的问题。只要完成了正确的转换，就可以套用普通数组上我们已经知道的解法。

本文我们看另一个比较直接的变种：一维数组变成矩阵。重点也是如何做问题的转换：将矩阵上的问题转换为数组上的问题，进而使用数组上的已知解法去解决。本题的两个算法对应了二维问题的两种思考方式：直接考虑二维情况看能不能找到解法，以及先固定其中一维，看另一维能不能直接用一维情况的解法来做。

---

题目

面试题 17.24. 最大子矩阵

给定一个正整数、负整数和 0 组成的 N × M 矩阵，编写代码找出元素总和最大的子矩阵。

返回一个数组 [r1, c1, r2, c2]，其中 r1, c1 分别代表子矩阵左上角的行号和列号，r2, c2 分别代表右下角的行号和列号。若有多个满足条件的子矩阵，返回任意一个均可。

说明：

1 <= matrix.length, matrix[0].length <= 200

样例:

输入：
[
[-1,0],
[0,-1]
]
输出：[0,1,0,1]
解释：输入中标粗的元素即为输出所表示的矩阵

题解

算法1: 直接在二维做

这是一个二维的矩阵上的问题。因此我们第一感觉是棋盘 DP。关于棋盘 DP 可以参考文章 棋盘DP。

如果定义 dp[i][j] 表示在矩形 matrix[0..i][0..j] 上，以 (i, j) 为右下角的子矩形中的最大的子矩形和，那么我们有没有办法找到状态转移的路径呢。

如下图，假如现在状态推导到了 dp[i][j]，下一个状态要计算 dp[i][j+1]。以 (i, j) 为右下角的最大子矩阵为图中蓝色阴影部分，它的高度 l 为 4，那么 dp[i][j+1] 怎么算呢。如果说以 (i, j+1) 为右下角的子矩形和最大时，它的高度也是 4 的话，那么我们就可以立刻知道 dp[i][j+1] 为 dp[i][j] 加上图中红色阴影部分的和。

但是以 (i, j) 为右下角的最大子矩阵高度为 4 的时候，以 (i, j+1) 为右下角的最大子矩阵高度可不一定是 4，而是所有高度都有可能。因此我们想到高度 l 是不是可以作为一个状态，我们把以 (i, j) 为右下角，且高度为 l 的最大子矩形都求一下，记为 dp[i][j][l]，这样对于每个 (i, j)，我们把所有可能的高度 l=1,2,...,j+1 下的最大子矩形都转移一次。

这样我们就找到了增加一个高度维度 l 的棋盘 DP 算法，如下：

状态定义：
dp[i][j][l] := 以 (i, j) 为右下角，高度为 l 的最大矩形和, 其中 1 <= l <= j + 1

答案
min(dp[i][j][l])

初始化
dp[0][0][1] = matrix[0][0]

转移
dp[i][j][l] = rec_sum(i, j, i, j - l + 1) + max(0, dp[i - 1][j][l])

仔细一看，其中状态转移部分 dp[i][j][l] = rec_sum(i, j, i, j - l + 1) + max(0, dp[i - 1][j][l]) 就是一维的最大子数组和的动态规划算法的转移方程 dp[i] = nums[i] + max(dp[i - 1], 0) 的二维版本。

此外，rec_sum(i, j, i, j - l + 1) 表示子矩阵 matrix[i][j..j-l+1] 的和。这部分可以用二维前缀和进行预处理，参考 【模板】前缀和与差分。

由于返回的是具体的 [r1, c1, r2, c2]，因此需要另一个 dp2 用于输出决策：

dp2[i][j][l] := 以 (i, j) 为右下角，高度为 l 的最大矩形和对应的上界下标

代码(C++)

class Solution {
public:
    vector<int> getMaxMatrix(vector<vector<int>>& matrix) {
        if(matrix.empty()) return {};
        int n = matrix.size(), m = matrix[0].size();

        vector<vector<int>> sums(n + 1, vector<int>(m + 1, 0));
        for(int i = 1; i <= n; ++i)
            for(int j = 1; j <= m; ++j)
                sums[i][j] = sums[i - 1][j] + sums[i][j - 1] - sums[i - 1][j - 1] + matrix[i - 1][j - 1];

        vector<vector<vector<int>>> dp(n, vector<vector<int>>(m, vector<int>(m + 1, -1)));
        // dp[i][j][l] := 以 (i, j) 为右下角，高度为 l 的最大矩形和
        vector<vector<vector<int>>> dp2(n, vector<vector<int>>(m, vector<int>(m + 1, 0)));
        // dp2[i][j][l] := 以 (i, j) 为右下角，高度为 l 的最大矩形和对应的上界下标

        dp[0][0][1] = matrix[0][0];
        for(int j = 1; j < m; ++j)
            for(int l = 1; l <= j + 1; ++l)
                dp[0][j][l] = rec_sum(0, j, 0, j - l + 1, sums);
        for(int i = 1; i < n; ++i)
        {
            dp[i][0][1] = matrix[i][0] + max(dp[i - 1][0][1], 0);
            if(dp[i - 1][0][1] > 0)
                dp2[i][0][1] = dp2[i - 1][0][1];
            else
                dp2[i][0][1] = i;
        }

        for(int i = 1; i < n; ++i)
            for(int j = 1; j < m; ++j)
            {
                for(int l = 1; l <= j + 1; ++l)
                {
                    dp[i][j][l] = rec_sum(i, j, i, j - l + 1, sums) + max(dp[i - 1][j][l], 0);
                    if(dp[i - 1][j][l] > 0)
                        dp2[i][j][l] = dp2[i - 1][j][l];
                    else
                        dp2[i][j][l] = i;
                }
            }

        int max_sum = INT_MIN;
        int mi = -1, mj = -1, ml = -1;
        for(int i = 0; i < n; ++i)
            for(int j = 0; j < m; ++j)
                for(int l = 1; l <= j + 1; ++l)
                    if(dp[i][j][l] > max_sum)
                    {
                        max_sum = dp[i][j][l];
                        mi = i;
                        mj = j;
                        ml = l;
                    }
        return vector<int>{dp2[mi][mj][ml], mj - ml + 1, mi, mj};
    }

private:
    int rec_sum(int i1, int j1, int i2, int j2, const vector<vector<int>>& sums)
    {
        return sums[i1 + 1][j1 + 1] - sums[i1 + 1][j2] - sums[i2][j1 + 1] + sums[i2][j2];
    }
};

算法2: 枚举一个维度，在另一个维度上以一维来做

前面我们的思考方式是直接从二维的角度考虑状态如何转移，发现直接转移是转移不了的，如果固定了高度 l 则可以转移，于是我们增加额外维度 l，进行状态的定义和转移，写出了相应的棋盘 DP 算法，然后我们发现转移的过程其实可以对应到一维情况下的最大子数组和。

这里我们尝试先从一维的算法出发，看能不能直接得到二维情况下的算法。

二维的情况一共有两个轴向，枚举其中一个轴向的所有区间，例如 axis=1 方向的区间 [j1, j2]。当固定 axis=1 方向的区间 [j1, j2] 后，在 axis=0 方向就可以按照类似于一维的的方式做了。

回顾一维数组的最大子数组和的动态规划算法(参考 最大子数组和的三种解法)，如下：

状态定义：
dp[i] := [0..i] 中，以 nums[i] 结尾的最大子串和

答案:
max(dp[i])

初始化：
dp[0] = nums[0]

状态转移：
dp[i] = nums[i] + max(dp[i - 1], 0)

因为求 dp[i] 时只用到了 dp[i - 1]，所以dp数组可以压缩掉，只用一个变量 sum 维护状态，状态转移过程的代码可以写成:

for(int i = 0; i < n; ++i)
{
    sum = nums[i] + max(sum, 0);
    ans = max(ans, sum);
}

回到二维的情况，如果我们固定 axis=1 轴向的 [j1, j2] 这一段，就可以仿照一维情况的算法来写了，如下:

状态定义：
dp[i] := axis=1 轴向选定区间 [j1, j2] 后，axis=0 轴向考虑 [0..i] 时，以 nums[i][j1..j2] 结尾的最大子串和

答案:
max(dp[i])

初始化：
dp[0] = sum(nums[0][j1..j2])

状态转移：
dp[i] = sum(nums[i][j1..j2]) + max(dp[i - 1], 0)

sum(nums[i][j1..j2]) 同样用前缀和预处理。

因为求 dp[i] 时只用到了 dp[i - 1]，所以dp数组可以压缩掉，只用一个变量 sum 维护状态，状态转移过程的代码可以写成:

for(int i = 0; i < n; ++i)
{
    sum = sum(nums[i][j1..j2]) + max(sum, 0);
    ans = max(ans, sum);
}

因为本题要返回具体的坐标，因此两个 max 都需要打开。一个 sum 的信息就不够了，还需要一个 pre 表示求和的起点，即当前在 axis=0 轴向上考虑求和的区间是 [pre, i]。

for(int i = 0, pre = 0; i < n; ++i)
{
    int val = sums[i + 1][j2 + 1] - sums[i + 1][j1] - sums[pre][j2 + 1] + sums[pre][j1];
    if(val > max_sum)
    {
        max_sum = val;
        result = {pre, j1, i, j2};
    }
    if(val <= 0)
        pre = i + 1;
}

代码 (c++)

class Solution {
public:
    vector<int> getMaxMatrix(vector<vector<int>>& matrix) {
        if(matrix.empty()) return {};
        int n = matrix.size(), m = matrix[0].size();

        vector<vector<int>> sums(n + 1, vector<int>(m + 1,0));
        for(int i = 1; i <= n; ++i)
            for(int j = 1; j <= m; ++j)
                sums[i][j] = sums[i - 1][j] + sums[i][j - 1] - sums[i - 1][j - 1] + matrix[i - 1][j - 1];

        int max_sum = INT_MIN;
        vector<int> result;
        for(int j1 = 0; j1 < m; ++j1)
        {
            for(int j2 = j1; j2 < m; ++j2)
            {
                // 枚举 axis=1 方向的区间 [j1, j2]
                // 固定 axis=1 方向的区间后，在 axis=0 方向做一维的最大子数组和
                for(int i = 0, pre = 0; i < n; ++i)
                {
                    int val = sums[i + 1][j2 + 1] - sums[i + 1][j1] - sums[pre][j2 + 1] + sums[pre][j1];
                    if(val > max_sum)
                    {
                        max_sum = val;
                        result = {pre, j1, i, j2};
                    }
                    if(val <= 0)
                        pre = i + 1;
                }
            }
        }
        return result;
    }
};

---