Longest Increasing Subsequence in O(n log n)

Given an array, find the length of the longest strictly increasing subsequence, the set of values that appear in order, though not necessarily next to each other. In [10, 9, 2, 5, 3, 7, 101, 18] the answer is 4, from [2, 3, 7, 18].

It's easy to come up with a quadratic solution using two for loops, but there exists an O(n log n) solution.

Why the obvious dynamic program is too slow

Let best[i] be the length of the longest increasing subsequence ending exactly at index i. To fill it you look back at every earlier index j with a smaller value and take the best of them, plus one.

def length_of_lis(nums):           # O(n^2) - correct, too slow
    best = [1] * len(nums)
    for i in range(len(nums)):
        for j in range(i):
            if nums[j] < nums[i]:
                best[i] = max(best[i], best[j] + 1)
    return max(best) if best else 0

This is easy to justify and easy to get right. It is also two nested loops, and at 100,000 elements that is ten billion comparisons.

The idea: track the best ending value per length

Turn the question around; instead of asking how long a run ending at each position can be, ask: for each possible length, what is the smallest value that a run of that length could end on?

Call that array tails, so tails[k] is the smallest value any increasing subsequence of length k + 1 ends with.

Observe:

1. tails is always sorted ascending. A subsequence one element longer must end on a larger value, because its final element extends a run that was already shorter.
2. Its length is always the answer so far. Position k only ever gets filled once a run of that length has genuinely been seen.

Why sorted matters: a sorted array is searchable in logarithmic time. That single fact is what removes the inner loop.

For each value x, find the first entry in tails that is not smaller than x.

• If there is no such entry, every run can be extended by x, so append it and the longest length grows by one.
• If there is one, overwrite it with x. A run of that length can now end on something smaller, which can only make future extensions easier.

Example

Look at [10, 9, 2, 5, 3, 7, 101, 18]. The right-hand column is the state after each step.

The code

Python has the binary search in the standard library; the other two spell it out. All three are the same algorithm.

import bisect


def length_of_lis(nums):
    tails = []
    for n in nums:
        # First index whose value is >= n. bisect_left, not bisect_right:
        # the subsequence must be strictly increasing, so an equal value
        # has to replace rather than extend.
        i = bisect.bisect_left(tails, n)
        if i == len(tails):
            tails.append(n)
        else:
            tails[i] = n
    return len(tails)

Why it is O(n log n)

Each of the n values is handled once. The work per value is a binary search over tails, which never holds more than n entries, so each search is O(log n). Appending and overwriting are both O(1). Overall, O(n log n) time and O(n) space.

The quadratic version does O(n) work per element instead of O(log n). At n = 100,000 that is the difference between roughly 1.7 million operations and 10 billion, which is ~5800x more efficient!

The mistakes that actually happen

• Using upper bound instead of lower bound. That accepts equal values as extensions and silently solves the non-decreasing problem. [7, 7, 7, 7] should answer 1, not 4 — it is the fastest way to catch this.
• Returning the contents of tails as the subsequence. The length is always right; the array usually is not a real subsequence.
• Appending when the search lands on the last index. Appending is only correct when the index is past the end, meaning no existing entry was large enough.
• Forgetting the empty input. An empty array has length 0, and tails is naturally empty, so the base case needs no special handling — but the quadratic version's max() call does.