最大子数组查找算法
最近开始看《算法导论(第三版)》这本书,在第4章讲解分治策略的时候,举了一个例子:最大子数组问题(4.1节
),简单来说,就是在一个给定的数组中,找出一个非空连续子数组,使这个数组中元素的总和最大,我们称这样的子数组为最大子数组。
在书中讲解分析了两种解决方法,一种就是最容易想到的暴力求解方法,另一种自然就是主题相关的分治法;然后在结束的时候,提到了可以实现一种更快速的线性时间的算法,这里将给出我的思路(原书这里描述有误,不知道是作者/译者笔误还是印刷错误)。
下面先给出书中提到的两种方案,然后附上最后的线性时间的解决方法。
对于讨论的子数组问题,只有当被查找的数组中存在负数的时候才有意义,因为对于一个全部为正数的数组,最大子数组自然是其本身。
2016年2月10日更新:加入另一种线性时间的代码(来源于网络,地址见下文)和分析。
暴力求解方法
暴力求解,自然就是枚举所有的子数组,计算每个子数组的和,然后查找出和最大的那个子数组。对于一段长度为n的数组,我们需要检查n + (n-1) + ... + 1 = n(n+1)/2
个子数组,时间复杂度可以表示为O(n^2)
,另一方面,虽然计算一个子数组之和所需的时间是线性的,但是当计算所有子数组之和时,我们可以通过缓存之前计算出的子数组和来计算当前子数组的和,使得每个子数组和的计算时间为O(1)
,因此暴力求解方法所花费的时间为O(n^2)
。
暴力求解算法很简单,就不附加相关代码了。
分治策略求解方法
假定我们需要寻找子数组A[low..high]
的最大子数组。
下面的描述中,也将采用这种描述数组的方式,方括号中使用一对由
..
分开的整数(变量),分别指示数组下标的下限和上限(包含上下限)。
使用分治策略意味着需要将子数组划分为两个规模尽量相等的子数组,将子数组A[low..high]
的中央位置下标记为mid
,然后考虑求解两个子数组A[low..mid]
和A[mid+1..high]
,那么A[low..high]
的任何连续子数组A[i..j]
所处的位置必然是以下三种情况之一:
- 完全位于子数组
A[low..mid]
中,因此 low <= i <= j <= high。 - 完全位于子数组
A[mid+1..high]
中,因此 mid < i <= j <= high。 - 跨越了中点,因此 low <= i <= mid < j <= high。
因此,A[low..high]
的一个最大子数组所处位置必然是这三种情况之一。实际上,A[low..high]
的一个最大子数组必然是完全位于A[low..mid]
中、完全位于A[mid+1..high]
中或者跨越中点的所有子数组中和的最大者。我们可以递归的求解A[low..mid]
和A[mid+1..high]
的最大子数组,因为这两个子问题仍是最大子数组问题,只是规模更小。因此剩下的问题就是寻找跨越中点的最大子数组,然后在三种情况中选择最大者。
寻找跨越中点的最大子数组,并不是原问题的规模更小的实例,因为它加入了限制:求和的子数组必须跨越中点,因为任何跨越中点的子数组都是由两个子数组
A[i..mid]
和A[mid+1..j
组成,其中 low <= i <= mid 且 mid < j <= high。因此我们只需要找出形如A[i..mid]
和A[mid+1..j]
的最大子数组,然后将其合并即可。
FIND_MAX_CROSSING_SUBARRAY
使用过程FIND_MAX_CROSSING_SUBARRAY
描述这个方法,这个过程接收数组A和下标low、mid、high为输入,返回一个具有三个元素的元组,这三个元素分别表示最大子数组的上下边界以及值得和。
下面给出伪代码:
FIND_MAX_CROSSING_SUBARRAY(A, low, mid, high)
left_sum = -1
sum = 0
max_left = mid
for (i = mid; i >= low; --i)
sum += A[i]
if (sum > left_sum)
left_sum = sum
max_left = i
right_sum = -1
sum = 0
max_right = mid + 1
for (j = mid + 1; j <= high; ++j)
sum += A[j]
if (sum > right_sum)
right_sum = sum
max_right = j
return (max_left, max_right, left_sum + right_sum)
两个循环分别向左向右搜索使left_sum
和right_sum
取值最大的下标,然后返回结果:子数组A[max_left..max_right]
和相应的和left_sum + right_sum
。可以很容易看出上面两个循环迭代一共执行了(mid - lwo + 1) + (high - mid) = high - low + 1 = n
次,因此这个计算过程的时间复杂度可以记为O(n)
。
FIND_MAXIMUM_SUBARRAY
有了这个FIND_MAX_CORSSING_SUBARRAY
在手,就可以设计求解最大子数组的分治算法的伪代码了:
FIND_MAXIMUM_SUBARRAY(A, low, high)
if (high == low)
return (low, high, A[low]) // 递归底层,开始回归
else
mid = (low + high) / 2 // C语言整数计算方式,结果向下取整
(left_low, left_high, left_sum) = FIND_MAXIMUM_SUBARRAY(A, low, mid)
(right_low, right_high, right_sum) = FIND_MAXIMUM_SUBARRAY(A, mid+1, high)
(cross_low, cross_high, cross_sum) = FIND_MAX_CROSSING_SUBARRAY(A, low, mid, high)
if (left_sum >= right_sum && left_sum >= cross_sum)
return (left_low, left_high, left_sum)
if (right_sum >= left_sum && right_sum >= cross_sum)
return (right_low, right_high, right_sum)
else
return (cross_low, cross_high, cross_sum)
初始调用FIND_MAXIMUM_SUBARRAY(A, 0, A.length-1)
。
时间复杂度不再细述,由递归深度lg(n)
可以得出过程FIND_MAXIMUM_SUBARRAY
的时间复杂度为O(n*lg(n))
。
一个线性时间的算法
在书中分治算法分析结束的时候,提到了对于这个求解最大子数组的问题,存在一个线性时间的算法,并在附加的练习题目中给出了提示:
使用如下思想为最大子数组问题设计一个非递归的、线性时间的算法。从数组的左边界开始,由左至右处理,记录到目前为止已经处理过的最大子数组。若已知
A[0..j]
的最大子数组,那么基于如下性质将解扩展为A[0..j+1]
的最大子数组:A[0..j+1]
的最大子数组要么是A[0..j]
的最大子数组,要么是某个子数组A[i..j+1]
(0 <= i <= j+1)。在已知A[0..j]
的最大子数组的情况下,可以在线性时间内找出形如A[i..j+1]
的最大子数组。
上面引用中最后一句话中被标记为加粗斜体的文字,是我认为书中记录失误的地方,很显然,如果找出最大子数组的时间仍然是线性的,那么总的时间必然是O(n^2)
,可以想象为一个双重循环样式。而实际上,通过缓存一些中间结果的方式,我们可以做到在常量时间内从前一个最大子数组推出下一个最大子数组。
分析
根据题目的提示,我们很容易可以推断出这个迭代过程中,第一个最大子数组就是A[0]
,然后开始下一步迭代。在我们将数组扩展到下一个元素的时候(从A[0..j]
扩展到A[0..j+1]
),那么下一个最大子数组将只有以下两种形式:
- 包含原来的最大子数组。
- 不包含原来的最大子数组。
而不会出现跨越原来的最大子数组的情况。我们将原来的最大子数组记为A[m..n]
,假设可以跨越原来的最大子数组的话,我们将新的最大子数组记为A[x..y]
,其中m < x <= n且n <= y,那么很显然我们可以构造一个新的子数组A[m..y]
,这个子数组完全包含了原来的最大子数组A[m..n]
,A[x..y]
和A[m..y]
可表示如下:
A[x..y] = A[x..n] + A[n..y] // 去掉重合的A[n]
A[m..y] = A[m..n] + A[n..y] // 去掉重合的A[n]
由于A[x..n]
被最大子数组A[m..n]
所包含,那么很容易得知A[m..n]
的和大于A[x..n]
的和,所以可推得A[m..y]
的和大于A[x..y]
的和,那么A[x..y]
就不是最大子数组了,与假设矛盾,因此不存在部分包含原来的最大子数组的情况。
现在我们考虑推导下一个最大子数组的问题。首先要想通下面几条规则:
- 如果包含原来的最大子数组的话,那么下一个最大子数组起始下标必然等于原来的最大子数组起始下标。
- 如果不包含原来的最大子数组的话,那么下一个最大子数组起始下标必然是原来的最大子数组结束下标之后的某一个正数的位置。不可能是负数,那是在做负功。
求解
为了方便说明,从现在开始,我们将原来的最大子数组记为A[p..q]
,而下一个最大子数组记为A[r..s]
。请丢掉上面反证法中定义的各种东西。
那么现在第一个问题是如何确定下一个最大子数组是否要包含原来的最大子数组呢?先考虑一种最简单的情况,根据最大子数组的定义,如果包含的话(p == r && s > q
),很显然必须有A[q+1..s]
的和大于0,否则岂不是做了负功,让A[r..s]
的和比A[p..q]
的和还要小。
因此可以顺便提一下,当我们遍历数组元素的时候,就可以顺便计算出上面需要的
A[q+1..s]
的和并缓存了,而不用等到需要的时候又再去重复遍历计算。
而在这里,我们也可以确定触发计算下一个最大子数组的几乎必要的条件:当前元素必须为正数,如果为负数,我们只需要简单的将其用来计算A[q+1..s]
的和;否则如果计算确定需要包含的话,那么可以很简单的将原最大子数组“扩容”即可,也就是将结束下标更改为当前遍历时刻的元素下标,并重置各种中间缓存变量的值,比如我们缓存的从当前最大子数组结束下标到下一个可能的最大子数组开始坐标之间的元素之和,就需要重置为0了。
那么第二个问题也来了,当我们确定不包含原来的最大子数组之和,那么下一个最大子数组的起始下标该如何获取?按照上面说的,这个位置是个正数元素,但未必是原来的最大子数组结束之后的第一个正数位置,考虑下面的序列:
8, -10, 3, -4, 9
如果当前最大子数组为第一个元素,那么当遍历到第三个元素(值为3)的时候,由于-10 + 3 = -7 < 0
,因此新的最大子数组不可能为8, -10, 3
,而3 < 8
,所以新的最大子数组也不可能为3
,我们需要继续遍历,但是要记录当前这个正数位置,万一它后面就是个非常大的正数呢。当然在这里,为了说明情况,我给它后面弄了个比-3还小的元素,因此我们就可以丢掉刚刚记录的正数位置了,因为没有用,如果新的最大子数组不包含原来的话,那么也不应该从这个记录的正数位置开始,因为它的值被后面的负数“中和”了,如果带上它,我们将会被迫同时带上它后面的那个负数,结果就很容易想到了。
所以继续遍历,第五个元素是9,是个正数,我们可以更新记录位置的那个变量了。很显然,-10 + 3 - 4 + 9 = -2 < 0
,新的可能的最大子数组开始位置就应该是刚刚记录的正数位置了,而9 > 8
,因此,新的最大子数组起始位置应该为第五个元素了(当然,目前结束位置也应该是)。
现在考虑一点例外。
把上面的结论推广,也就是从当前最大子数组开始考虑,如果为了方便思考的话,可以将当前最大子数组等价替换为当前位置的一个元素,其值等于最大子数组的和。那么有时候即使上面的A[q+1..s]
计算为正,我们也不能直接包含当前最大子数组。考虑下面这个序列:
6, -8, 10
最开始的最大子数组为6
,然后我们按照前面的计算方式,可以得到-8 + 10 = 2 > 0
,但是这个时候,很明显可以看出正确的最大子数组应该为10
,而不是6, -8, 10
。至于原因,也很容易根据上面计算需要的正数位置的方法来解释,因为6 - 8 = -2
,所以我们的下一个可能的最大子数组不可以从6
开始,这个正数指针(此处不一定代表C语言中的那种指针,只是意义类似)应该被更新到下一个正数,然后继续进行相似的验证。
到这里,恰好统一了,初始情况下我们的正数指针应该指向第一个最大子数组的起始位置。而任何时候我们更新了当前最大子数组,也需要同步设置这个正数指针的值。而只要这个指针的值不等于当前最大子数组的起始地址,说明如果需要更新最大子数组的信息的话,需要从当前这个正数地址开始。
最后,还有个例外也需要考虑,那就是数组中全负数的情况,这时候就完全不需要上面这么麻烦了,只需要找出整个序列的最大值即可。
总结
- 当正数指针值等于当前最大子数组的起始地址时,如果从当前最大子数组结束后被遍历元素之和大于0,则需要“扩容”当前最大子数组。
- 当正数指针值不等于当前最大子数组的起始地址时,如果从这个正数指针开始的一个序列元素的和大于当前子数组的和,则需要将这个序列设置为当前最大子数组。
- 当开头一段序列全部是负数的时候,可以忽略掉这些元素;作为一个特例的情况,序列如果全部是负数,那么只需要返回其中的最大值。
好了,现在基本分析可以说已经完毕了,剩下的还是通过代码来理顺吧。
伪C++
代码
tuple<unsigned int, unsigned int, int> find_maximum_subarray(vector<int> &array)
{
/* 不考虑异常等情况,比如array为空之类的 */
unsigned int start = 0;
unsigned int end = 0;
int cur_sum = array[0];
unsigned int i = 1;
if (cur_sum <= 0)
{
for (; i < array.size() && array[i] <= 0; ++i) // 略过开始的非正数序列
{
if (cur_sum < array[i])
{
cur_sum = array[i];
start = i;
}
}
if (i == array.size()) // 全负数特例
{
return tuple<unsigned int, unsigned int, int>(start, start, cur_sum);
}
cur_sum = array[i];
start = end = i;
}
int temp_sum = 0; // 缓存当前最大子数组后面被遍历元素之和
unsigned int next_start = start; // 下一个候选最大子数组,以正数开始
int next_sum = cur_sum; // 下一个候选最大子数组的累积和
for (; i < array.size(); ++i)
{
temp_sum += array[i];
if (array[i] > 0 && next_sum <= 0)
{
/* 出现正数被负数抵消情况 */
next_start = i;
next_sum = array[i];
}
else
{
next_sum += array[i];
}
if (next_start == start && temp_sum > 0)
{
end = i;
cur_sum += temp_sum;
temp_sum = 0;
next_sum = cur_sum;
}
else if (next_start != start && next_sum > cur_sum)
{
start = next_start;
end = i;
cur_sum = next_sum;
temp_sum = 0;
}
}
return tuple<unsigned int, unsigned int, int>(start, end, cur_sum);
}
大概代码最终如上,还没有进行测试,也许考虑仍有不周。但是就上述代码来说,运行时间为线性的,第一个循环是针对开始的负数序列,同时也包括了全负数的形式。而当进入第二个循环的时候,两个保存元素和的变量cur_sum
、next_sum
均为正值,因此这个循环中的代码可以只针对正数/正负数混合优化,而不用考虑有全负数序列,也不用测试cur_sum
等是否为负。而通过变量i
的值可知,迭代次数一共为n(==array.size()
)次,因此,整个过程的时间复杂度为O(n)
。
另一种解法
在clrs.skanev.com上又发现了另外一种解决方法,使用C
语言实现;相比较上面的方法,思路更加通俗易懂,可以很容易的从简单的思考中想到并优化获得最终的优化过的代码。
下面贴出页面的代码,然后简单分析一下。
typedef struct
{
unsigned left;
unsigned right;
int sum;
} max_subarray;
max_subarray find_maximum_subarray(int A[], unsigned low, unsigned high)
{
max_subarray suffixes[high - low];
suffixes[0].left = low;
suffixes[0].right = low + 1;
suffixes[0].sum = A[low];
for (int i = low + 1; i < high; i++)
{
if (suffixes[i - 1].sum < 0)
{
suffixes[i].left = i;
suffixes[i].right = i + 1;
suffixes[i].sum = A[i];
}
else
{
max_subarray *previous = &suffixes[i - 1];
suffixes[i].left = previous->left;
suffixes[i].right = i + 1;
suffixes[i].sum = previous->sum + A[i];
}
}
max_subarray *max = &suffixes[0];
for (int i = low + 1; i < high; i++)
{
if (max->sum < suffixes[i].sum)
{
max = &suffixes[i];
}
}
return *max;
}
这段代码中一共需要对数组遍历两次,但是每次迭代过程中运行的代码也相对更少一些,与C++
版本相比,在时间效率上相差不大,对于全正数序列,这段C
代码相对要快一点,而如果是一个全负数序列,上面的C++
版本运行应该更快,而其它情况下,则需要视实际序列情况而定,但总的来说,差别不大。而在空间复杂度上,C
版本的代码很明显通常就需要更多的空间了,因为在第十行声明了一个长度等于n
(==high-low
)、大小为sizeof(max_subarray) * n
的数组,相比而言,C++
版本中总共只使用了大约7 * sizeof(int)
字节的空间。
除此之外,这段C代码由于使用了不定长数组语法,因此需要能支持
C99
标准的编译器来进行编译,当然,这在现在基本不是问题,只不过目前很多编译器即使支持部分或完整的C99
标准,但默认编译选项依然是使用C89
(ISO C
/ANSI C
)标准,需要注意查看。
回到分析代码上,首先定义了一个结构体用于存储子数组的信息,left
指示了子数组开始位置,right
代表了子数组结束位置的下一位,注意这个位置是个哨兵位置,并不应该参与子数组和sum
的计算。在方法体中声明的suffixes
数组中,suffixes[i]
存储在原数组A
中以A[i]
为结束位置的最大子数组的信息(所以suffixes[i].right = i + 1
)。
因此,从最简单粗暴的方法入手,要获取以A[i]
作为结束位置的最大子数组信息,只需要向前遍历计算就可以,但是这样一来,时间复杂度就达到了O(n^2)
了。所以我们就要想办法利用缓存,利用上次计算的数据,在O(1)
时间内计算出下一个最大子数组信息,分析一下上面的过程,每次都需要前向遍历,而实际上,每次遍历只比上一次同样的过程多分析一个元素:第一次是A[0]
,第二次则是A[0..1]
,第三次则是A[0..2]
……一直重复,那么我们应该想办法利用上一次的结果以及当前的元素,计算出本次结果。借用上面的分析,如果一个最大子数组开始位置要不同于前一次的结果,那么只有当前一次的结果小于0,否则应该直接在前一次的结果上继续累加,只有前一次结果小于0的时候,才需要以当前位置作为新的开始位置。
需要始终记住的是:当前元素必须作为结束位置,因此其实我们只有两种选择,要不加入前一次的结果(因为前一次的结果其实就相当于不包括当前元素的遍历结果),要不单独起始。而很容易想到,如果前一次结果已经被加到为负数了,当前元素无论为正为负,都肯定比加一个负数所得结果要大。
因此出了第一次我们直接把起始元素的值作为sum
的值,之后每次只需要利用前一次的结果经过常量时间的运算即可得到本次的结果。当所有结果都计算出来了之后,就很简单了,只需要一次寻找最大值得操作即可。
这段C
代码很容易就可以证明其正确性,而有这段分析,类比之下,也许上面的C++
算法的正确性也可以更容易验证了。
尾注
心血来潮,偶然为之。在最后的完成线性时间的算法时,也是经历了不断的修正,也许结果仍然不完善,待测试。最后所为算法并没有完全按照题目的提示,但是第二个循环中最后的if
和else if
语句中,由end = i
可知,这两个判断都是为了从当前最大子数组获取下一个具有A[i..j+1]
形式的最大子数组,而在条件判断之外则没有更新过最大子数组信息,也就是依然等于A[0..j]
的最大子数组。