非常特殊矩阵的最大子矩形

Question

在处理图像处理任务时，我遇到了以下问题：单位正方形中有n个点，分别具有坐标$x_i$和$y_i$，每个点都赋予一个正或负权重$w_i$。找到一个矩形，使得该矩形内的所有点的权重之和为正且最大。

通过定义一个适当的网格，可以将问题重新表述为在一个n×n矩阵A中找到一个子矩阵，其元素之和最大。这也被称为“最大子矩形问题”，并已经在SO上讨论过。虽然暴力方法的运行时间为O(n^5)，但有一种巧妙的解决方法，其运行时间为O(n^3)。它利用了对应一维问题的解决方案，称为"最大子数组问题", 其运行时间为O(n)。

我已经在R中实现了这两种算法，并且可以在几秒钟内解决数百个点。但是对于成千上万个点，它将太慢了，即使将循环外包给一些Fortran或C代码。

现在看一下矩阵A。假设（没有丧失一般性），所有点的x或y坐标都不同，A具有特殊形式：在A的每一行和列中，恰好有一个非零元素。对于具有这种特殊属性的矩阵，我认为应该有一种算法可以在O(n^2)时间内执行任务，甚至更好。

这里是添加最佳矩形的示例：

set.seed(723)
N <- 50; w <- rnorm(N)
x <- runif(N); y <- runif(N)
clr <- ifelse (w >= 0, "blue", "red")
plot(x, y, pch = 20, col = clr, xlim = c(0, 1), ylim = c(0, 1))
rect(0.075, 0.45, 0.31, 0.95, border="gray")

你会发现在最优矩形中可能出现红色，即负数点。这也表明仅解决 x 和 y 坐标的一维情况是不够的。

我将把标准解法翻译成Fortran，但我肯定希望能掌握更有效的算法。

- Hans W.

Answer 1

这些人（从维基页面中发现）声称已经找到了二维情况下更简单的亚立方解决方案。这可能是你已经知道的那个。

- phs

---

谢谢。请记住，我只对矩阵每行和每列仅有一个非零元素的这种特殊情况感兴趣。这个属性应该可以让我们实现更高效的解决方案。 - Hans W.

啊，我错过了。每行或每列只有一个非零元素的方阵是对角矩阵和置换矩阵的乘积。很可能可以找到一维情况的变体，应用于对角矩阵，并将其扩展到可由某个特定置换类达到的矩阵。 - phs

是的，有点儿。不幸的是，仅将一维情况应用于对角线并不起作用。在阅读文章和思考一些实际几何应用之后，这看起来更像一个研究问题。一个合理的解决方案可能会被发表在某个地方。 - Hans W.

Answer 2

请参考“稀疏矩阵中的最大总和子矩形”中的被接受答案。对于具有m个非零元素的nxn矩阵，该解决方案需要O(nm logn)的时间。因此，对于您而言，由于恰好有n个非零元素，这将需要O(n^2 logn)的时间。可能您将能够处理n比标准O(n^3)解决方案大50倍或更多的情况。

- user2566092

Answer 3

我能做到的最好时间复杂度是O(n^2 log n)。

如果我们看一下Kadane的2D算法对于你这种类型的输入所做的n+1个选择2次调用Kadane的1D算法，除了O(n)个连续的对之外，其余都在仅相差一个元素的1D数组上。我将介绍Kadane 1D的分治变体；通过缓存每个递归调用的结果，只需要重新计算涉及更改的数组元素的O(log n)次，将内部循环的（平摊）运行时间从Theta(n)降低到Theta(log n)。

def maxsubarray(arr, a, b):
    # this function returns a 4-tuple
    # element 0 is the max over intervals of the form [i, j)
    # element 1 is the max over intervals of the form [i, b)
    # element 2 is the max over intervals of the form [a, j)
    # element 3 is the max over intervals of the form [a, b), i.e., sum(arr[a:b])
    n = b - a
    if n == 0:
        return (0, 0, 0, 0)
    elif n == 1:
        x = arr[a]
        y = max(x, 0)
        return (y, y, y, x)
    else:
        m = a + n // 2
        l = maxsubarray(arr, a, m)
        r = maxsubarray(arr, m, b)
        return (max(l[0], r[0], l[1] + r[2]),
                max(r[1], l[1] + r[3]),
                max(l[2], l[3] + r[2]),
                l[3] + r[3])

- Stacy

---

注意：如果您知道arr不为空，则可以删除n == 0的情况。 - Stacy

抱歉，这是什么样的回答呢：（1）你的代码没有正确运行；你有试过吗？例如使用 a=[1,-1.5,2,-1,3,-2,4,-5,6] 进行测试吗？（2）声称可以在 O(log n) 的时间内解决一维情况是错误的。算法必须至少触及每个元素，因此最小值为 O(n)。（3）递归会减慢速度，也可能导致长数组的“无限递归”。（4）我的问题不是一维情况（我可以在不到0.1秒的时间内解决包含数百万元素的数组），而是如何利用矩阵的特殊形式。 - Hans W.

那只是一个愚蠢的例子。我在一些更长的例子（100个元素）上尝试了您的代码，但它没有给出正确的结果。顺便说一句：要应用这样的函数，它需要返回索引，即最大子数组从哪里到哪里运行。 - Hans W.

(1) 我之前没有测试过，但它运行良好：使用maxsubarray（arr，0，len（arr））调用。 (2) 第一个需要线性时间；更新为O（log n）。 (3) 然后将其变成迭代形式。 (4) 是的，我意识到了，这就是为什么我告诉你将这个一维算法插入到二维Kadane中。 (5) 修改最大值以记住索引并不难。 (6) 那是一系列粗鲁的评论，所以我将不会帮助你完成任何事情。 - Stacy

如果我说话听起来很严厉，我向大家道歉。对于一维情况，我有一个非常快速的编译例程，并且我将其应用于矩阵情况。使我的解决方案变慢的是通过矩阵行进行双重循环（O（n ^ 2）部分）。因此，我仍然对如何利用矩阵的非常特殊形式的想法感兴趣。谢谢。 - Hans W.