# Longest increasing subsequence

### Question

Given an input sequence, what is the best way to find the longest (not necessarily continuous) non-decreasing subsequence.

``````0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15 # sequence

1, 9, 13, 15 # non-decreasing subsequence

0, 2, 6, 9, 13, 15 # longest non-deceasing subsequence (not unique)
``````

I'm looking for the best algorithm. If there is code, Python would be nice, but anything is alright.

1
31
10/21/2010 11:04:17 PM

I just stumbled in this problem, and came up with this Python 3 implementation:

``````def subsequence(seq):
if not seq:
return seq

M = [None] * len(seq)    # offset by 1 (j -> j-1)
P = [None] * len(seq)

# Since we have at least one element in our list, we can start by
# knowing that the there's at least an increasing subsequence of length one:
# the first element.
L = 1
M = 0

# Looping over the sequence starting from the second element
for i in range(1, len(seq)):
# Binary search: we want the largest j <= L
#  such that seq[M[j]] < seq[i] (default j = 0),
#  hence we want the lower bound at the end of the search process.
lower = 0
upper = L

# Since the binary search will not look at the upper bound value,
# we'll have to check that manually
if seq[M[upper-1]] < seq[i]:
j = upper

else:
# actual binary search loop
while upper - lower > 1:
mid = (upper + lower) // 2
if seq[M[mid-1]] < seq[i]:
lower = mid
else:
upper = mid

j = lower    # this will also set the default value to 0

P[i] = M[j-1]

if j == L or seq[i] < seq[M[j]]:
M[j] = i
L = max(L, j+1)

# Building the result: [seq[M[L-1]], seq[P[M[L-1]]], seq[P[P[M[L-1]]]], ...]
result = []
pos = M[L-1]
for _ in range(L):
result.append(seq[pos])
pos = P[pos]

return result[::-1]    # reversing
``````

Since it took me some time to understand how the algorithm works I was a little verbose with comments, and I'll also add a quick explanation:

• `seq` is the input sequence.
• `L` is a number: it gets updated while looping over the sequence and it marks the length of longest incresing subsequence found up to that moment.
• `M` is a list. `M[j-1]` will point to an index of `seq` that holds the smallest value that could be used (at the end) to build an increasing subsequence of length `j`.
• `P` is a list. `P[i]` will point to `M[j]`, where `i` is the index of `seq`. In a few words, it tells which is the previous element of the subsequence. `P` is used to build the result at the end.

How the algorithm works:

1. Handle the special case of an empty sequence.
3. Loop over the input sequence with index `i`.
4. With a binary search find the `j` that let `seq[M[j]` be `<` than `seq[i]`.
5. Update `P`, `M` and `L`.
6. Traceback the result and return it reversed.

Note: The only differences with the wikipedia algorithm are the offset of 1 in the `M` list, and that `X` is here called `seq`. I also test it with a slightly improved unit test version of the one showed in Eric Gustavson answer and it passed all tests.

Example:

``````seq = [30, 10, 20, 50, 40, 80, 60]

0    1   2   3   4   5   6   <-- indexes
``````

At the end we'll have:

``````M = [1, 2, 4, 6, None, None, None]
P = [None, None, 1, 2, 2, 4, 4]
result = [10, 20, 40, 60]
``````

As you'll see `P` is pretty straightforward. We have to look at it from the end, so it tells that before `60` there's `40,`before `80` there's `40`, before `40` there's `20`, before `50` there's `20` and before `20` there's `10`, stop.

The complicated part is on `M`. At the beginning `M` was `[0, None, None, ...]` since the last element of the subsequence of length 1 (hence position 0 in `M`) was at the index 0: `30`.

At this point we'll start looping on `seq` and look at `10`, since `10` is `<` than `30`, `M` will be updated:

``````if j == L or seq[i] < seq[M[j]]:
M[j] = i
``````

So now `M` looks like: `[1, None, None, ...]`. This is a good thing, because `10` have more chanches to create a longer increasing subsequence. (The new 1 is the index of 10)

Now it's the turn of `20`. With `10` and `20` we have subsequence of length 2 (index 1 in `M`), so `M` will be: `[1, 2, None, ...]`. (The new 2 is the index of 20)

Now it's the turn of `50`. `50` will not be part of any subsequence so nothing changes.

Now it's the turn of `40`. With `10`, `20` and `40` we have a sub of length 3 (index 2 in `M`, so `M` will be: `[1, 2, 4, None, ...]` . (The new 4 is the index of 40)

And so on...

For a complete walk through the code you can copy and paste it here :)

29
5/23/2017 12:10:37 PM

Here is how to simply find longest increasing/decreasing subsequence in Mathematica:

`````` LIS[list_] := LongestCommonSequence[Sort[list], list];
input={0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15};
LIS[input]
-1*LIS[-1*input]
``````

Output:

``````{0, 2, 6, 9, 11, 15}
{12, 10, 9, 5, 3}
``````

Mathematica has also LongestIncreasingSubsequence function in the Combinatorica` libary. If you do not have Mathematica you can query the WolframAlpha.

## C++ O(nlogn) solution

There's also an O(nlogn) solution based on some observations. Let Ai,j be the smallest possible tail out of all increasing subsequences of length j using elements a1, a2, ... , ai. Observe that, for any particular i, Ai,1, Ai,2, ... , Ai,j. This suggests that if we want the longest subsequence that ends with ai + 1, we only need to look for a j such that Ai,j < ai + 1 < = Ai,j + 1 and the length will be j + 1. Notice that in this case, Ai + 1,j + 1 will be equal to ai + 1, and all Ai + 1,k will be equal to Ai,k for k!=j+1. Furthermore, there is at most one difference between the set Ai and the set Ai + 1, which is caused by this search. Since A is always ordered in increasing order, and the operation does not change this ordering, we can do a binary search for every single a1, a2, ... , an.

## Implementation C++ (O(nlogn) algorithm)

``````#include <vector>
using namespace std;

/* Finds longest strictly increasing subsequence. O(n log k) algorithm. */
void find_lis(vector<int> &a, vector<int> &b)
{
vector<int> p(a.size());
int u, v;

if (a.empty()) return;

b.push_back(0);

for (size_t i = 1; i < a.size(); i++) {
if (a[b.back()] < a[i]) {
p[i] = b.back();
b.push_back(i);
continue;
}

for (u = 0, v = b.size()-1; u < v;) {
int c = (u + v) / 2;
if (a[b[c]] < a[i]) u=c+1; else v=c;
}

if (a[i] < a[b[u]]) {
if (u > 0) p[i] = b[u-1];
b[u] = i;
}
}

for (u = b.size(), v = b.back(); u--; v = p[v]) b[u] = v;
}

/* Example of usage: */
#include <cstdio>
int main()
{
int a[] = { 1, 9, 3, 8, 11, 4, 5, 6, 4, 19, 7, 1, 7 };
vector<int> seq(a, a+sizeof(a)/sizeof(a));
vector<int> lis;
find_lis(seq, lis);

for (size_t i = 0; i < lis.size(); i++)
printf("%d ", seq[lis[i]]);
printf("\n");

return 0;
}
``````

Example data is: `{ 1, 9, 3, 8, 11, 4, 5, 6, 4, 19, 7, 1, 7 }` and answer: `1 3 4 5 6 7`.