Python implementation of the longest increasing subsequence

PEP8 is the python style guide. To improve your style a start is to run the pep8 script against your file. Other style "problems" you'll have to find them yourself (or ask here on core review :)

Comments and docstrings

Check the PEP8 session Comments.

This:

## Returns the longest subsequence (non-contiguous) of X that is strictly increasing.
def subsequence(X):

Should be a docstring:

def subsequence(X):
    """Returns the longest subsequence (non-contiguous) of X that is strictly
    increasing.

    """

Variable names

I know that you kept the wikipedia variable names, but you should try to make your script conistent with itself. So, squeeze your variable name imagination, and choose something short, sweet and meaningful. You can always leave a note at the beginning, linking to the wikipedia page and say which variable is which.

With this:

L = 1     ## length of longest subsequence (initially: just first element)

I'd do:

# length of longest subsequence (initially: just first element)
len_longest = 1

Since the comment was a bit too long to be an inline comment I moved it. A good name for X might be seq. For M maybe pos_smallest and for P pos_predecessor or something like that, maybe with a better understanding of the underling algorithm you can find better names :)

range() and xrange()

Check the link above and you'll see that one is a list and the other an itertor (genrator).

So, this:

for i in range(1,len(X)):

Should be:

for i in xrange(1, len(seq)):

Note that in Python 3 xrange() become range().

Reverse a list

This line:

return list(reversed(result))

Should be:

return result[::-1]

Other than being more clear, it will be Python 3 compatible.

Useless brackets are useless

Python is not C so don't put brackets everywhere.

Here:

while (trace_idx >= 0):

Could just be:

while trace_idx >= 0:

Do the same thing in another couple of places.

Binary search

Your use of bisect_left looks really strange:

j = bisect_left([X[M[idx]] for idx in range(L)], X[i])

And even more the previous version:

while True:
    if (start == end):
        if (X[M[start]] < X[i]):
            j = start
        break
    else:
        partition = 1 + ((end - start - 1) / 2)
        if (X[M[start + partition]] < X[i]):
            start += partition
            j = start
        else:
            end = start + partition - 1

Which could have easily been:

while start != end:
    partition = 1 + (end - start - 1) / 2
    if X[M[start + partition]] < X[i]:
        start += partition
        j = start
    else:
        end = start + partition - 1

if X[M[start]] < X[i]:
    j = start

But you shouldn't do it like that anyway. Take a look at Binary Search in Python, there are both solution, the accepted is in "pure" Python (the one that I like the most) and there's the one with bisect_left. If you care about performances time them both.

This is just a pointing and I would have really liked to give you the code of the binary search, but you didn't really explained how your algorithm works (and it's quite different from the one showed in the wikipedia article), so I wasn't able to improve yours without changing its core.

Usually what one should do is instantiate at the beginning the entire M and P list, something like:

M = [0] * len(seq)

So you won't have to deal with append and so on... that are quite hard to follow.
I couldn't also get around your -1 offestes, are they really necessary? That's the main different I've found with the "official" algorithm.

Note: some of the suggestions here concern the revised version of the code, that was posted as this other answer.

List comprehension in bisect

You could rewrite:

j = bisect_left([seq[head[idx]] for idx in xrange(len(head))], seq[i])

trivially as:

j = bisect_left([seq[k] for k in head], seq[i])

There's nothing to gain with the original approach.

Or, if you make seq a numpy.array, then you could use:

j = bisect_left(seq[head], seq[i])

If-statements in the for-loop

You don't really need to check if seq[i] < seq[head[j]].

From bisect_left line we know that all(seq[i] <= v for v in seq[head][j:]), so your check is equivalent to if seq[i] != seq[head[j]].

But I see no reason why you can't update head[j] = i in such case. Values in seq[head] stay the same, since seq[i] == seq[head[j]]. It also won't affect predecessor paths of longer subsequences.

On the other hand, j == len(head) implies that seq[i] == seq[head[j]], so maybe you just meant:

if j == len(head): head.append(i)
else: head[j] = i

to avoid unnecessary assignments. And that's ok. But in this particular case, I propose to ask forgiveness, not permission:

try: head[j] = i
except: head.append(j)

Parentheses

Also, you don't need the parentheses in while (trace_idx >= 0):, this should be just:

while trace_idx >= 0:

Generating random numbers

You don't need randomList if you use numpy:

from numpy.random import randint

l1 = list(randint(100, size=15))

Miscellaneous suggestions

• You may consider renaming the function to include the full name.

• You may consider handling the empty sequence correctly.

• I think it is possible to have a better name for the head variable.

• Using None gives better readability to mark that there are no predecessors.

• In some applications, it may be beneficial to return indices instead of actual values.

• It is possible to generate the indices in correct order directly, using recursion.

Here's the function with all remarks applied:

def longest_increasing_subsequence(seq):
  if not seq: return seq

  lastoflength = [0] # end position of subsequence with given length
  predecessor = [None] # penultimate element of l.i.s. ending at given position

  for i in range(1, len(seq)):
    # seq[i] can extend a subsequence that ends with a smaller element
    j = bisect_left([seq[k] for k in lastoflength], seq[i])
    # update existing subsequence of length j or extend the longest
    try: lastoflength[j] = i
    except: lastoflength.append(i)
    # remember element before seq[i] in the subsequence
    predecessor.append(lastoflength[j-1] if j > 0 else None)

  # return indices ..., p(p(p(i))), p(p(i)), p(i), i
  def trace(i):
    if i is not None:
      yield from trace(predecessor[i])
      yield i
  return list(trace(lastoflength[-1]))

Here's the revised code, based on the suggestions from code reviews, excluding @arekolek's answer, which appears to be a review of this revised code here.

from random import randrange
from bisect import bisect_left


def randomList(N, max):
    return [randrange(max) for x in xrange(N)]


def subsequence(seq):
    """Returns the longest subsequence (non-contiguous) of seq that is
    strictly increasing.

    """
    # head[j] = index in 'seq' of the final member of the best subsequence 
    # of length 'j + 1' yet found
    head = [0]
    # predecessor[j] = linked list of indices of best subsequence ending
    # at seq[j], in reverse order
    predecessor = [-1]
    for i in xrange(1, len(seq)):
        ## Find j such that:  seq[head[j - 1]] < seq[i] <= seq[head[j]]
        ## seq[head[j]] is increasing, so use binary search.
        j = bisect_left([seq[head[idx]] for idx in xrange(len(head))], seq[i])

        if j == len(head):
            head.append(i)
        if seq[i] < seq[head[j]]:
            head[j] = i

        predecessor.append(head[j - 1] if j > 0 else -1)

    ## trace subsequence back to output
    result = []
    trace_idx = head[-1]
    while (trace_idx >= 0):
        result.append(seq[trace_idx])
        trace_idx = predecessor[trace_idx]

    return result[::-1]


l1 = randomList(15, 100)