最大子序列和求解过程

问题来历

该问题最早出现在布朗大学Grenander面对的一个模式匹配的问题，问题的最初形式是给定nxn的实数数组，求出矩形子数组的最大总和。最大综合子数组是数字图像中某种特定模式的最大似然估计量，因为二维问题求解需要太多时间，所以将它简化为一维问题从而深入了解其结构。这个问题的解决经历了四种复杂度的算法最终得到解决，一开始的算法的时间复杂度为 O(n3),第一次改良之后的时间复杂度为O(n2)，再改良为 O(nlogn),最终改良为 O(n)。从一开始解决规模为100 000的问题需要15天到最后只需要5毫秒，提升是显而易见的。
在这篇博客中将会从 O(n3) 开始逐步得到 O(n)，在这里先定义一组测试数据x[0..9]如下

31?415926?535897?96?93?2384

最终的结果是x[2..6]，和为187。知道数据跟结果我们就开始我们的解题吧

O(n3)解法

O(n3)解法很简单，就是我们看到这道题的时候潜意识中想象到的方法，也叫做暴力破解法。用语言描述就是对所有满足0≤i≤j的(i,j)整数对进行迭代，对每一个整数对计算x[i..j]的和，然后检验总和是否大于迄今为止的最大综合，用伪代码表示就是

tmax=0
for i=[0,n]
    for j=[i,n]
        sum=0
        for k=[i,j]
            sum+=x[k]
        tmax=max(tmax,sum)

tmax=0
x=[31,-41,59,26,-53,58,97,-93,-23,84]
for i in range(0,10):
    for j in range(i,10):
        sum=0
        for k in range(i,j+1):
            sum+=x[k]
        tmax=max(tmax,sum)
print tmax

O(n2)的两种解法

一般来说都不会想到O(n3)的解法，都是从O(n2)开始，从O(n3)到O(n2)有两种改良的方法，都是为了去除第三套循环。第一个改良的关注点是，x[i..j]的总和与前面已计算出的总和x[i..j-1]的总和密切相关，改良后的伪代码

tmax=0
for i=[0,n)
    sum=0
    for j=[i,n)
        sum+=x[j]
        tmax=max(tmax,sum)

可执行的python代码

tmax=0
x=[31,-41,59,26,-53,58,97,-93,-23,84]
for i in range(0,10):
    sum=0
    for j in range(i,10):
        sum+=x[j]
        tmax=max(tmax,sum)
print tmax

第二种改良方法是通过访问在外循环执行之前就已构建的数据结构的方式在内循环中计算总和。假设一个数组sumarr的第i个元素包含了x[0..i]中各个数的累加和，所以x[i..j]中各个数的总和可以通过计算sumarr[j]-sumarr[i-1]得到，该方法改良后的伪代码为

sumarr[-1]=0
for i=[0,n)
    sumarr[i]=sumarr[i-1]+x[i]
tmax=0
for i=[0,n)
    for j=[i,n)
        sum+=x[j]
        tmax=max(tmax,sum)

这里转化为实际代码的时候存在取数组第-1个数的问题，解决的方法就是把原来的数组加长一个数。
可运行的python代码

tmax=0
x=[31,-41,59,26,-53,58,97,-93,-23,84]
sumarr=[0]*11
for i in range(0,10):
    sumarr[i+1]=sumarr[i]+x[i]
for i in range(0,10):
    for j in range(i,10):
        sum=sumarr[j+1]-sumarr[i]
        tmax=max(tmax,sum)
print tmax

O(nlogn)解法

O(nlogn)算法利用的是分治策略，也就是，要解决规模为n的问题，可递归地解决两个规模近似为n/2的子问题，然后对它们的答案进行合并以得到整个问题的答案。
使用分治策略解决的思路是，把数组分为ab两部分，然后求出a部分以及b部分的最大子序列和为ma ，mb，还有一个最重要的是求出横跨ab部分的最大子序列的和为mc，然后比较三个之中最大的值。根据这个指导思想可以得出以下代码

float mss(l,u)
    if (l>u)
        return 0
    if (l==u)
        return max(0,x[1])
    m=(l+u)/2
    lmax=sum=0
    for (i=m;i>=1;i--)
        sum+=x[i]
        lmax=max(lmax,sum)

    rmax=sum=0
    for i=(m,u]
        sum+=x[i]
        rmax=max(rmax,sum)
    return max(lmax+rmax,mss(l,m),mss(m+1,u))

转换为可执行的python代码为

x=[31,-41,59,26,-53,58,97,-93,-23,84]
def mss(l,u):
    if (l>u):
    return 0
    if (l==u):
        return max(0,x[1])
    m=(l+u)/2

    lmax=sum=0
    for i in (range(0,m))[::-1]:
        sum+=x[i]
        lmax=max(lmax,sum)

    rmax=sum=0

    for i in range(m,u+1):
        sum+=x[i]
        rmax=max(rmax,sum)
    return max(lmax+rmax,mss(l,m),mss(m+1,u))

if __name__=='__main__':
    print mss(0,9)

O(n)解法

这种解法应该是最广为流传的最佳解法，也就是一般人说的扫描算法，解题的思路是，从数组最左端开始扫描，一直到最右端，并几下所遇到的总和最大的子序列。最重要的原理就是，前i个元素中，最大总和子序列要么在前i-1个元素中（设为tmax），要么其结束位置为i（设为fmax）。
计算fmax有一个巧妙的方法，不需要从头开始计算结束位置为i的最大子序列，而是利用结束位置为i-1的最大子序列进行计算，具体实现如下

tmax=0
fmax=0
for i=[0,n)
    fmax=max(fmax+x[i],0)
    tmax=max(tmax,fmax)

该代码需要细细体会。关键点在于fmax，fmax肯定是包含最后一个数的最大的子序列，而tmax是目前位置总和最大的子序列。
很容易转换为python代码

tmax=0
fmax=0
x=[31,-41,59,26,-53,58,97,-93,-23,84]


for i in range(0,10):
    fmax=max(fmax+x[i],0)
    tmax=max(tmax,fmax)
print tmax

总结

这是我看《编程珠玑》的时候第一次感受到一个好的算法居然可以让程序提速得如此夸张，还有一个原因是这也是面试或者机试的常考题，也许以前在一些竞赛书上能够看到最快速的解法，但能够看到一个算法从O(n3)到O(n)的进化简直是一种享受。而且在这其中还有一个小故事是，当得到O(nlogn)算法的时候，《编程珠玑》的作者和他身边的人都以为这是最佳解法，在一个会议上公布这个成果的时候，一个与者这只花了几分钟就想到了O(n)的解法，可以看出算法设计也是要看灵感的。好的算法往往更加简洁而不容易出错。