Theta(n**2) and Theta(n*lgn) algorithm perform unproperly

Question

I'm reading the Introduction to Algorithms and trying to finish the exercise in the book.

In Exercise 4.1-3

4.1-3 Implement both the brute-force and recursive algorithms for the maximum-subarray problem on your own computer. What problem size n0 gives the crossover point at which the recursive algorithm beats the brute-force algorithm? Then, change the base case of the recursive algorithm to use the brute-force algorithm whenever the problem size is less than n0. Does that change the crossover point?

I wrote the two algorithms according to the book's pseudo-code. However, there must be something wrong with my code because the second one, which is designed to be Theta(n*lgn) and supposed to run faster, always runs slower than the first Theta(n**2) one. My codes is shown below.

def find_maximum_subarray_bf(a):        #bf for brute force
    p1 = 0
    l = 0           # l for left
    r = 0           # r for right
    max_sum = 0
    for p1 in range(len(a)-1):
        sub_sum = 0
        for p2 in range(p1, len(a)):
            sub_sum += a[p2]
            if sub_sum > max_sum:
                max_sum  = sub_sum
                l = p1
                r = p2
    return l, r, max_sum

def find_maximum_subarray_dc(a):        #dc for divide and conquer

    # subfunction
    # given an arrary and three indics which can split the array into a[l:m]
    # and a[m+1:r], find out a subarray a[i:j] where l \leq i \less m \less j \leq r".
    # according to the definition above, the target subarray must
    # be combined by two subarray, a[i:m] and a[m+1:j]
    # Growing Rate: theta(n)

    def find_crossing_max(a, l, r, m):

        # left side
        # ls_r and ls_l indicate the right and left bound of the left subarray.
        # l_max_sum indicates the max sum of the left subarray
        # sub_sum indicates the sum of the current computing subarray      
        ls_l = 0
        ls_r = m-1
        l_max_sum = None
        sub_sum = 0
        for j in range(m+1)[::-1]:      # adding elements from right to left
            sub_sum += a[j]
            if sub_sum > l_max_sum:
                l_max_sum = sub_sum
                ls_l = j

        # right side
        # rs_r and rs_l indicate the right and left bound of the left subarray.
        # r_max_sum indicates the max sum of the left subarray
        # sub_sum indicates the sum of the current computing subarray                
        rs_l = m+1
        rs_r = 0
        r_max_sum = None
        sub_sum = 0
        for j in range(m+1,len(a)):
            sub_sum += a[j]
            if sub_sum > r_max_sum:
                r_max_sum = sub_sum
                rs_r = j

        #combine
        return (ls_l, rs_r, l_max_sum+r_max_sum)

    # subfunction
    # Growing Rate: should be theta(nlgn), but there is something wrong
    def recursion(a,l,r):           # T(n)
        if r == l:
            return (l,r,a[l])
        else:
            m = (l+r)//2                    # theta(1)
            left = recursion(a,l,m)         # T(n/2)
            right = recursion(a,m+1,r)      # T(n/2)
            crossing = find_crossing_max(a,l,r,m)   # theta(n)

            if left[2]>=right[2] and left[2]>=crossing[2]:
                return left
            elif right[2]>=left[2] and right[2]>=crossing[2]:
                return right
            else:
                return crossing

    #back to master function
    l = 0
    r = len(a)-1
    return recursion(a,l,r)

if __name__ == "__main__":

    from time import time

    a = [100,-10,1,2,-1,4,-6,2,5]
    a *= 2**10

    time0 = time()
    find_maximum_subarray_bf(a)
    time1 = time()
    find_maximum_subarray_dc(a)
    time2 = time()
    print "function 1:", time1-time0
    print "function 2:", time2-time1 
    print "ratio:", (time1-time0)/(time2-time1)

Answer 1

First, a mistake in the brute-force:

for p1 in range(len(a)-1):

that should be range(len(a)) [or xrange], as is, it wouldn't find the maximum subarray of [-12,10].

Now, the recursion:

def find_crossing_max(a, l, r, m):

    # left side
    # ls_r and ls_l indicate the right and left bound of the left subarray.
    # l_max_sum indicates the max sum of the left subarray
    # sub_sum indicates the sum of the current computing subarray      
    ls_l = 0
    ls_r = m-1
    l_max_sum = None
    sub_sum = 0
    for j in range(m+1)[::-1]:      # adding elements from right to left

You are checking all the indices to 0, but you should only check the indices to l. Instead of constructing the range list and reversing it, use xrange(m,l-1,-1)

        sub_sum += a[j]
        if sub_sum > l_max_sum:
            l_max_sum = sub_sum
            ls_l = j

For the sum to the right, the analogue holds, you should only check indices to r, so xrange(m+1,r+1).

Further, your intitial values for the sums resp. indices for the maximum subarray are dubious for the left part, and wrong for the right.

For the left part, we start with an empty sum but must include a[m]. That can be done by setting l_max_sum = None initially, or by setting l_max_sum = a[m] and letting j omit the index m. Either way, the initial value for ls_l should not be 0, and for ls_r it shouldn't be m-1. ls_r must be m, and ls_l should start as m+1 if the initial value for l_max_sum is None, and as m if l_max_sum starts as a[m].

For the right part, r_max_sum must start as 0, and rs_r should better start as m (though that isn't really important, it would only give you the wrong indices). If none of the sums on the right is ever non-negative, the right sum should be 0 and not the largest of the negative sums.

In recursion, we have a bit of duplication in

left = recursion(a,l,m)         # T(n/2)

the sums including a[m] have already been treated or majorised in find_crossing_max, so that could be

left = recursion(a,l,m-1)

But then one would have to also treat the possibility r < l in recursion, and the repetition is small, so I'll let that stand.

Since you always traverse the entire list in find_crossing_max, and that is called O(n) times, your divide-and-conquer implementation is actually O(n²) too.

If the range checked in find_crossing_max is restricted to [l,r], as it should be, you have (approximately) 2^k calls on ranges of length n/2^k, 0 <= k <= log_2 n, for a total cost of O(n*log n).

With these changes (and some random array generation),

def find_maximum_subarray_bf(a):        #bf for brute force
    p1 = 0
    l = 0           # l for left
    r = 0           # r for right
    max_sum = 0
    for p1 in xrange(len(a)):
        sub_sum = 0
        for p2 in xrange(p1, len(a)):
            sub_sum += a[p2]
            if sub_sum > max_sum:
                max_sum  = sub_sum
                l = p1
                r = p2
    return l, r, max_sum

def find_maximum_subarray_dc(a):        #dc for divide and conquer

    # subfunction
    # given an arrary and three indices which can split the array into a[l:m]
    # and a[m+1:r], find out a subarray a[i:j] where l \leq i \less m \less j \leq r".
    # according to the definition above, the target subarray must
    # be combined by two subarray, a[i:m] and a[m+1:j]
    # Growing Rate: theta(n)

    def find_crossing_max(a, l, r, m):

        # left side
        # ls_r and ls_l indicate the right and left bound of the left subarray.
        # l_max_sum indicates the max sum of the left subarray
        # sub_sum indicates the sum of the current computing subarray      
        ls_l = m+1
        ls_r = m
        l_max_sum = None
        sub_sum = 0
        for j in xrange(m,l-1,-1):      # adding elements from right to left
            sub_sum += a[j]
            if sub_sum > l_max_sum:
                l_max_sum = sub_sum
                ls_l = j

        # right side
        # rs_r and rs_l indicate the right and left bound of the left subarray.
        # r_max_sum indicates the max sum of the left subarray
        # sub_sum indicates the sum of the current computing subarray                
        rs_l = m+1
        rs_r = m
        r_max_sum = 0
        sub_sum = 0
        for j in range(m+1,r+1):
            sub_sum += a[j]
            if sub_sum > r_max_sum:
                r_max_sum = sub_sum
                rs_r = j

        #combine
        return (ls_l, rs_r, l_max_sum+r_max_sum)

    # subfunction
    # Growing Rate:  theta(nlgn)
    def recursion(a,l,r):           # T(n)
        if r == l:
            return (l,r,a[l])
        else:
            m = (l+r)//2                    # theta(1)
            left = recursion(a,l,m)         # T(n/2)
            right = recursion(a,m+1,r)      # T(n/2)
            crossing = find_crossing_max(a,l,r,m)   # theta(r-l+1)

            if left[2]>=right[2] and left[2]>=crossing[2]:
                return left
            elif right[2]>=left[2] and right[2]>=crossing[2]:
                return right
            else:
                return crossing

    #back to master function
    l = 0
    r = len(a)-1
    return recursion(a,l,r)

if __name__ == "__main__":

    from time import time
    from sys import argv
    from random import randint
    alen = 100
    if len(argv) > 1:
        alen = int(argv[1])
    a = [randint(-100,100) for i in xrange(alen)]

    time0 = time()
    print find_maximum_subarray_bf(a)
    time1 = time()
    print find_maximum_subarray_dc(a)
    time2 = time()
    print "function 1:", time1-time0
    print "function 2:", time2-time1 
    print "ratio:", (time1-time0)/(time2-time1)

We get something like we should expect:

$ python subarrays.py 50
(3, 48, 1131)
(3, 48, 1131)
function 1: 0.000184059143066
function 2: 0.00020382
ratio: 0.902923976608
$ python subarrays.py 100
(29, 61, 429)
(29, 61, 429)
function 1: 0.000745058059692
function 2: 0.000561952590942
ratio: 1.32583792957
$ python subarrays.py 500
(35, 350, 3049)
(35, 350, 3049)
function 1: 0.0115859508514
function 2: 0.00170588493347
ratio: 6.79175401817
$ python subarrays.py 1000
(313, 572, 3585)
(313, 572, 3585)
function 1: 0.0537149906158
function 2: 0.00334000587463
ratio: 16.082304233
$ python osubarrays.py 10000
(901, 2055, 4441)
(901, 2055, 4441)
function 1: 4.20316505432
function 2: 0.0381460189819
ratio: 110.186204655