BRICS

Basic Research in Computer Science

Faster Algorithms for Computing Longest Common Increasing Subsequences

Gerth Stølting Brodal Kanela Kaligosi Irit Katriel Martin Kutz

BRICS Report Series RS-05-37

ISSN 0909-0878 December 2005 BRICS RS-05-37 Brodal et al.: Faster Algorithms for Computing Longest Common Increasing Subsequences Copyright \begin{matrix}c \\ ◯ \end{matrix} 2005, Gerth Stølting Brodal & Kanela Kaligosi & Irit

Katriel & Martin Kutz. BRICS, Department of Computer Science University of Aarhus. All rights reserved.

Reproduction of all or part of this work is permitted for educational or research use on condition that this copyright notice is included in any copy.

See back inner page for a list of recent BRICS Report Series publications. Copies may be obtained by contacting:

BRICS Department of Computer Science University of Aarhus Ny Munkegade, building 540 DK–8000 Aarhus C Denmark

Telephone: +45 8942 3360
Telefax: +45 8942 3255
Internet: BRICS@brics.dk

BRICS publications are in general accessible through the World Wide Web and anonymous FTP through these URLs:

http://www.brics.dk
ftp://ftp.brics.dk

This document in subdirectory RS/05/37/

Faster Algorithms for Computing Longest Common Increasing Subsequences

Gerth Stølting Brodal^{∗}^{†} Kanela Kaligosi^{‡} Irit Katriel^{∗}^{§}

Martin Kutz^{‡}

December 14, 2005

Abstract

We present algorithms for finding a longest common increasing subsequence of two or more input sequences. For two sequences of lengths m and n, where m ≥ n, we present an algorithm with an output-dependent expected running time of O((m + nℓ) log log σ + Sort) and O(m) space, where ℓ is the length of a LCIS, σ is the size of the alphabet, and Sort is the time to sort each input sequence.

For k ≥ 3 length-n sequences we present an algorithm which improves the previous best bound by more than a factor k for many inputs. In both cases, our algorithms are conceptually quite simple but rely on existing sophisticated data structures.

Finally, we introduce the problem of longest common weakly-increasing (or non-decreasing) subsequences (LCWIS), for which we present an O(m+ n log n) time algorithm for the 3-letter alphabet case. For the extensively studied Longest Common Subsequence problem, comparable speedups have not been achieved for small alphabets.

1 Introduction

Algorithms that search for the longest common subsequence (LCS) of two input sequences or the longest increasing subsequence (LIS) of one input sequence date back several decades.

Formally, given two sequences A = (a_{1}, …, a_{n}) and B = (b_{1}, …, b_{m}) with elements from an alphabet Σ and with m ≥ n, a common subsequence of A and B is a subsequence (a_{j_{1}} = b_{κ_{1}}, a_{j_{2}} = b_{κ_{2}}, … a_{j_{ℓ}} = b_{κ_{ℓ}}), where j_{1} < j_{2} < ⋯ < j_{ℓ} and κ_{1} < κ_{2} < ⋯ < κ_{ℓ}. Given one sequence A = (a_{1}, …, a_{n}) where the a_{i}’s are drawn from a totally ordered set, an increasing subsequence of A is a subsequence (a_{j_{1}}, a_{j_{2}}, …, a_{j_{ℓ}}) such that j_{1} < j_{2} < ⋯ < j_{ℓ} and a_{j_{1}} < a_{j_{2}} < ⋯ < a_{j_{ℓ}}.

∗BRICS, Basic Research in Computer Science (www.brics.dk), funded by the Danish National Research Foundation, University of Aarhus, Århus, Denmark.

{gerth,irit}@daimi.au.dk.

†Partially supported by the Danish Research Agency.

‡Max-Plank-Institut für Informatik, Saarbrücken, Germany. {kaligosi,mkutz}@ mpi-sb.mpg.de.

§Supported by the Danish Research Agency (grant # 272-05-0081).

A classic algorithm by Wagner and Fischer [11] solves the LCS problem using dynamic programming in O(mn) time and space. Hirschberg [6] reduced the space complexity to O(n), using a divide-and-conquer approach. The fastest known algorithm by Masek and Paterson [8] runs in O(n^{2}/ log n) time. Faster algorithms are known for special cases, such as when the input consists of permutations or when the output is known to be very long or very short. Hunt and Szymanski [7] studied the LCS problem in terms of matching index pairs, i.e., they defined r to be the number of index-pairs (i, j) with a_{i} = b_{j} (such a pair is called a match) and designed an algorithm that finds the LCS of two sequences in O(r log n) time. For a survey on the LCS problem see [2].

Fredman [4] showed how to compute a LIS of a length-n sequence in optimal O(n log n) time. The expected length of a longest increasing subsequence of a random permutation has been shown (after successive improvments) to be 2\sqrt{n} − o(\sqrt{n}); for a survey see [1].

Note that after sorting both input sequences we can in linear time remove symbols that do not appear in both sequences and rename the remaining symbols to the alphabet {1, 2, …, σ}. We can therefore assume that this preprocessing stage was performed and hence the size of the alphabet, σ, is at most n. In the following we let Sort_{Σ}(m) denote the time required to sort a length-m input sequence drawn from the alphabet Σ.

Recently, Yang et al. [12] combined the two concepts, and defined a common increasing subsequence (CIS) of two sequences A and B, i.e., an increasing sequence which is a subsequence of both A and B. They designed a dynamic programming algorithm that finds a longest CIS (an LCIS, for short) of A and B using Θ(mn) time and space.

Subsequently, Chan et al. [3] obtained an upper bound of O(min{r log σ, mσ+ r} log log m + Sort_{Σ}(m)). The number of matches r is in the worst case Ω(mn), but in some important cases it is much smaller. For instance, when A and B are permutations of {1, …, n} then r = O(n).

Chan et al. proceeded to generalize their algorithm to find an LCIS of k ≥ 3 length-n sequences. They show that this can be done in

O(min{kr^{2}, kr log σ log^{k−1} r} + kSort_{Σ}(n))

time, where r is again the number of matches, i.e., k-coordinate vectors that contain an index from each input sequence, all with the same symbol.

1.1 Our results

In this paper we present three new upper bounds for the LCIS problem. The first is an output-dependent algorithm which runs in O((m+nℓ) log log σ+Sort_{Σ}(m)) expected time and O(m) worst-case space, where ℓ is the length of an LCIS. Whenever n = Ω(log log σ + Sort_{Σ}(m)/m) and either m = Ω(n log log σ) or ℓ = o(n/ log log n), it is faster than Yang et al.’s Θ(mn)-time algorithm.

Symbol Meaning
m, n Lengths of input sequences (we assume m ≥ n).
ℓ  Length of the LCIS/LCWIS.
k  Number of input sequences.
σ  Size of the alphabet (number of different symbols).
r  Number of matches in the input sequences.

Table 1: Parameters of the LCIS/LCWIS problems.

Previous Results  New
k = 2  O(mn) [12]  O((m + nℓ) log log σ + Sort_{Σ}(m))
O(min{r log σ, mσ + r} log log m  O(m) when σ = 2
+Sort_{Σ}(m)) [3]  O(m + n log n) when σ = 3

k ≥ 3 O(min{kr^{2}, kr log σ log^{k−1} r}  O(min{kr^{2}, r log^{k−1} r log log r}
+kSort_{Σ}(n)) [3]  +kSort_{Σ}(n))

Table 2: Previous and new results. The new upper bounds apply to both LCIS and LCWIS.

For a strictly-increasing subsequence we have ℓ ≤ σ. However, in the weaklyincreasing (i.e. non-decreasing) variant, the length of the output can be arbitrarily larger than the size of the alphabet. We show that a longest common weakly increasing subsequence (LCWIS) can be found in linear time for an alphabet of size two and in O(m + n log n) time for an alphabet of size three. These results are interesting because they pinpoint what seems to be a fundamental difference between the LCS and LWCIS problems. The approach we use cannot be applied to LCS, and to date, comparable speedups have not been achieved for LCS with small alphabets.

Finally, we consider the case of k ≥ 3 length-n sequences. The upper bound of Chan et al. is achieved by two algorithms; the first is a simple O(kr^{2} + kSort_{Σ}(n)) time algorithm and the second is a more complex implementation of the same approach, which runs in O(kr log σ log^{k−1} r + kSort_{Σ}(n)) time. We describe an algorithm which is significantly simpler than the latter and obtain a running time of O(min{kr^{2}, r log^{k−1} r log log r} + kSort_{Σ}(n)).

Table 1 provides a list of the symbols used in the paper and Table 2 summarizes the previous and new results.

The rest of the paper is organized as follows. In Section 2 we describe a dynamic programming algorithm that uses a data structure based on van Emde Boas trees and runs in expected O((m+nℓ) log log σ+Sort_{Σ}(m)) time and O(m) space. In Section 3 we present our results on LCWIS with small alphabets, which use different techniques. Finally, in Section 4 we describe how to use a data structure by Gabow et al. [5] to obtain an algorithm for finding an LCIS or LCWIS of k ≥ 3 sequences, which is simpler and faster than Chan et al.’s algorithm.

2 An Output-Dependent Upper Bound

2.1 Bounded heaps

In our output-dependent algorithm we use a data structure, in the following denoted a bounded heap (BH), that supports the following operations:

• Insert(H, k, p, d): Insert into the BH H the key k with priority p and associated data d.

• DecreasePriority(H, k, p, d): If the BH H does not contain the key k, perform Insert(H, k, p, d). Otherwise, set this key’s priority to min{p, p′}, where p′ is its previous priority.

• BoundedMin(H, k): Return the item that has minimum priority among all items in H with key smaller than k. If H does not contain any items with key smaller than k, return “invalid”.

The priority search tree (PST) of McCreight [9] supports each of these operations in O(log n) time. However, the PST also allows deletions, which the BH is not required to support. Using van Emde Boas trees, we obtain a faster BH for integer keys:

Lemma 1 There exists an implementation of bounded heaps that requires O(n) space and supports each operation in O(log log n) amortized time, where keys are drawn for the set {1, …, n}.

Proof. The data structure applies standard techniques, such as those described in Section 3 of [5].

We rely on the fact that a snapshot of the heap, at any point in time, can be represented as a decreasing step function. More precisely, let BM (s) be the value that would be returned by a BoundedMin(H, s) query. Then BM (s) ≤ BM(s′) whenever s > s′, i.e., the function BM is non-increasing in s (see Figure 1).

key k  1 2 3 4 5 6 7 8 9 10
priority 7 10 6 8 5 3 2 4 1 9
BM (k) ∞ 7 7 6 6 5 3 2 2 1

Figure 1: Example of BM values.

Assume that the keys are {s_{1}, s_{2}, …} with s_{i} ≤ s_{i+1} for all i. To answer BoundedMin queries, it suffices to maintain a search structure that contains the BM (s_{i}) value for every s_{i} at which the function BM changes, i.e., BM (s_{i}) < BM (s_{i−1}). Then, we answer a BoundedMin(s) by searching the data structure for the largest key which is at most s and returning its BM value. If the search structure is a van Emde Boas tree [10], this takes O(log log n) time.

It remains to show how to support Insert and DecreasePriority operations in O(log log n) amortized time. When the priority of a key s_{i} decreases to a new value of p, the following occurs:

1. s_{i} + 1 is inserted into the tree if p < BM (s_{i^{−}}), where s_{i^{−}} is the largest key in the tree which is smaller than s_{i}.

2. s_{j} is removed from the tree if j > i and BM (s_{j}) > p.

With van Emde Boas trees, the two steps are handled, respectively, as follows.

1. Searching for s_{i^{−}}, checking whether s_{i} should be inserted and inserting it if so, takes O(log log n) time.

2. Beginning at s_{i}, we repeatedly find the next item s_{j} in the tree (i.e., the smallest key larger than the current one) and remove it from the tree if BM (s_{j}) > p. The total time is O(k log log n), where k is the number of items that were removed. Since the total number of items deleted by DecreasePriority operations is upper bounded by the total number of Insert operations, we can charge the cost of each deletion to the insertion of the same item, and obtain that the amortized cost of each operation is O(log log n). □

2.2 An O((m + nℓ) log log σ + Sort_{Σ}(m)) time algorithm

Our output-dependent algorithm for the LCIS problem is shown in Figure 3 in the appendix. In a preprocessing step, it removes from each sequence all elements which do not appear in the other sequence; this is easy after the sequences are sorted. For every remaining element s, it generates a sorted list Occ_{s} that contains ∞ and the indices of all occurrences of s in B.

Then, the algorithm in n iterations identifies common increasing subsequences (CISs) of increasing lengths: In iteration i it identifies length-i CISs (using the results of iteration i − 1). More precisely, for every element a_{j} in A, it identifies the minimum index κ in B such that there is a length-i CIS which ends at a_{j} in A and at b_{κ} in B. The index κ is stored in L_{i}[j].

To compute the array L_{1}[1 … n], the algorithm traverses A and for each a_{j}, sets L_{1}[j] to be the minimum index in the list Occ_{a_{j}}, i.e., the earliest occurrence of a_{j} in B. Note that due to the preprocessing, there exists such an index in B.

For i > 1, the ith iteration proceeds as follows. The algorithm traverses A again, and for every a_{j}, it checks whether a_{j} (together with some b_{κ}) can extend a length-(i − 1) CIS to a length-i CIS, and if so, identifies the minimum such κ. For this purpose, the algorithm maintains a bounded heap H. When it begins processing a_{j}, H contains all elements a_{t} ∈ {a_{1}, …, a_{j−1}} for which L_{i−1}[t] ≠ ∞. The key of a_{t} in H is a_{t} itself and its priority is L_{i−1}[t], i.e., the minimum index of the endpoint in B of a length-(i − 1) CIS which ends, in A, at index t. The algorithm queries H to find the leftmost endpoint (in B) of a length-(i − 1) CIS which contains only elements smaller than a_{j}. Let κ′ be this endpoint. Then, L_{i}[j] is set to the first occurrence of a_{j} in B which lies behind κ′; we prove that this is the leftmost endpoint in B of a length-i CIS which ends, in A, at a_{j}.

We emphasize that H is built anew for every single pass. The only information saved between different scans of A and B is maintained in the arrays L_{i}.

The arrays Link_{1}, Link_{2}, … are used to save the information we need in order to construct the LCIS: Whenever we detect that the index pair (j, κ) can extend a length-(i − 1) CIS which ends at the index pair (j′, κ′), we set Link_{i}[j] = j′. Finally, if there is a length-(i − 1) CIS which ends at a_{j}, then a_{j} is inserted into H with priority L_{i−1}[a_{j}]; it may later be extended into a length-i CIS by some a_{j′} with j′ > j.

2.2.1 Correctness

The correctness of the algorithm relies on the following lemma, which states that if there is a solution then the algorithm finds it. It is straightforward to show that the algorithm will not produce an invalid sequence.

Lemma 2 Let A and B be two sequences that have a length-ℓ CIS which ends in A at index j and in B at index κ. Then at the end of the iteration in which i = ℓ, L_{ℓ}[j] ≤ κ.

Proof. By induction on ℓ. For ℓ = 1, the claim is obvious. Assume that it holds for any length-(ℓ − 1) CIS and that we are given A and B which have a length-ℓ CIS c_{1}, …, c_{ℓ}, which is located in A as a_{j_{1}}, …, a_{j_{ℓ}} and in B as b_{κ_{1}}, …, b_{κ_{ℓ}}.

By the induction hypothesis, at the end of the i = ℓ − 1 iteration, L_{i−1} contains entries which are not equal to ∞. Hence, the algorithm will proceed to perform iteration i = ℓ. Again by the induction hypothesis, L_{ℓ−1}[j_{ℓ−1}] ≤ κ_{ℓ−1}.

Since a_{j_{ℓ−1}} < a_{j_{ℓ}}, it is guaranteed that when j = j_{ℓ}, H contains an item with key a_{j_{ℓ−1}}, priority κ′ ≤ κ_{ℓ−1} and d = (j_{ℓ−1}, κ′). So the BoundedMin operation will return a valid value. If the value returned is (j_{ℓ−1}, κ_{ℓ−1}), then the smallest occurrence of a_{ℓ} in B after κ_{ℓ−1} is not beyond κ_{ℓ}. So the algorithm will set L_{ℓ}[j_{ℓ}] ≤ κ_{ℓ}. On the other hand, if the value returned is not (j_{ℓ−1}, κ_{ℓ−1}), then it is (j_{ℓ−1}, κ′) for some κ′ ≤ κ_{ℓ−1}. Since a_{j′} < a_{ℓ}, again we get that the smallest occurrence of a_{ℓ} in B after κ_{ℓ−1} is not beyond κ_{ℓ}. So the algorithm will set L_{ℓ}[j_{ℓ}] ≤ κ_{ℓ}. □

2.2.2 Time complexity

The preprocessing phase takes O(Sort_{Σ}(m)) time, to sort each of the sequences A and B. The construction of the Occ_{s}’s takes O(m) time.

The array A is traversed O(ℓ) times. During each traversal, O(n) operations are performed on the bounded heap, each of which takes O(log log σ) amortized time, and the Occ_{s} lists are queried at most n times. We now sketch a possible implementation of the Occ_{s} lists.

We partition the range {1, …, m} into m/σ blocks of σ consecutive locations and for every 1 ≤ i ≤ m/σ we denote by b_{i} the block containing locations (i − 1)σ + 1, …, iσ. For each i and each s ∈ Σ we create a data structure that represents occurances of s in the block b_{i} and is based on Willard’s y-fast tries. In addition, for each block we store the first occurance of s succeeding the block. To answer a query in Occ_{s}, we first identify the block containing the query point in constant time. We then search for the smallest index larger than the query point in the y-fast trie for this block in time O(log log σ). If we found one, we are done. Otherwise, we return the first s succeeding the block, using the stored information. Initializing the m y-fast tries with a total of m elements takes O(m log log σ) expected time. Note that this initialization step needs to be carried out only once.

In total, the main loop takes O(m + nℓ log log σ) time. Finally, Constructing the LCIS takes O(ℓ) time. We get that the total expected running time of the algorithm is O((m + nℓ) log log σ + Sort_{Σ}(m)).

2.2.3 Space complexity

As for space complexity, note that in the main loop we only use L_{i−1} and L_{i}. Therefore, we do not need to save the previous L’s. In order to construct the LCIS, the algorithm as described requires O(nℓ) space for the Link arrays.

However, we can reduce the space complexity to O(m) with the technique developed by Hirschberg [6] for LCS. First, we run the algorithm once to compute ℓ (without constructing the Link arrays). Then we run a recursive version of the algorithm that construct the LCIS. The top recursive level invokes the usual algorithm, except that this time we remember only some of the Link information: Each match in the second half of a CIS knows the location in A and B of the ⌊ℓ/2⌋-th match of the CIS that it was appended to. This information is found in the ⌊ℓ/2⌋-th iteration of the main loop and propagated by the later iterations while the L arrays are constructed. Then, we know for every LCIS the location (i, j) in A and B of the middle match. We select one LCIS and recursively run the same algorithm to find the length-⌊ℓ/2⌋−1 LCIS of (a_{1}, …, a_{i−1}) (b_{1}, …, b_{j−1}) and the length-⌈ℓ/2⌉ LCIS of (a_{i+1}, …, a_{n}) and (b_{j+1}, …, b_{m}). The base case is when we look for a constant-size LCIS. Then we run the original algorithm in linear space. To achive that the time complexity remains unchanged we need to limit the work done processing B during the recursion. For the preprocessing for the outermost recursion we need time Sort_{Σ}(m). For the remaining recursive calls we do not need to sort the arrays again and the preprocessing time is O(m). The computation of a middle match considers at most matches involving nℓ entries from B. These entries in B can be marked during the computation of the middle match, and only this subsequence of B is provided to the recursive calls. The thinning of B is done before each recursive call. Let T (m, n, ℓ) be the running time of the recursion on two sequences of lengths n and m with a length-ℓ LCIS and m ≤ nℓ. Assume that the middle match is (n_{1}, m_{1}). Then T (m, n, ℓ) ≤ nℓ log log σ + nℓ + T (m_{1}, n_{1}, ℓ/2) + T (m_{2}, n_{2}, ℓ/2), where n_{1} +n_{2} +1 = n and m_{1} +m_{2} +1 ≤ m. This recurrence solves to O(nℓ log log σ). The total running time becomes O((m + nℓ) log log σ + Sort_{Σ}(m)). It is easy to see that the amount of space we need is O(m). In conclusion, we have shown:

Theorem 1 An LCIS of two sequences of lengths m and n with m ≥ n can be found in O((m + nℓ) log log σ + Sort_{Σ}(m)) expected time and O(m) worst-case space where ℓ is the length of the output and Sort_{Σ}(m) is the time required to sort a length-m input sequence.

3 Weakly Increasing Subsequences

We now turn to longest common non-decreasing or weakly increasing subsequences (LCWIS) for small alphabets. By simply replacing < by ≤ in the BoundedMin operation in our algorithm for the LCIS problem, it is straightforward to verify that the algorithm solves the LCWIS problem in O((m + nℓ) log log σ + Sort(m)) time. But while the LCIS problem can be solved in linear time for alphabets of bounded size t, simply because the length of the solution is then also bounded by t, it is not clear how this fact should carry over to LCWIS, where the output size need not relate to t at all.

We show how to solve LCWIS for the 2– and the 3-letter alphabet in linear respectively O(m+n log n) time. This is in contrast to the classic LCS problem, where already the 2-letter case seems to be essentially as hard as the general problem. In fact, it seems that LCWIS behaves very different from both LCIS and LCS.

3.1 Preprocessing

Let us use as our alphabet the Greek letters Σ = {α, β, γ} in their standard order: α < β < γ. For both tasks, the 2-letter and 3-letter cases, we prepare arrays Num_{A,α},Num_{B,α},Num_{A,β}, …,Num_{B,γ} that count the number of αs, βs and γs, respectively, in prefixes of A and B. For example, Num_{A,γ}[9] contains the number of γs in A up to position 9 (inclusively). We also create arrays Pos_{A,α} through Pos_{B,γ}, which provide us with the position of the ith occurrence of α, β, or γ in A or B. E.g., Pos_{B,α}[5] contains the position of the 5th α within sequence B. These arrays can clearly be prepared in O(m) time. Note that they also provide constant-time lookup for the reversed mode, which counts elements from the end of the sequences A and B. For example, Pos_{A,β}[Num_{A,β}[n] − 23] tells us the position of the last but 23rd β in A.

3.2 The 2-letter case is simple

After the preprocessing, the 2-letter case becomes trivial. For each i, where 0 ≤ i ≤ min{Num_{A,α}[n], Num_{B,α}[m]}, we determine the position of the ith α in A and B and then the number of βs that come after those positions in the two sequences using the expression

min{ Num_{A,β}[n] − Num_{A,β}[Pos_{A,α}[i]] , Num_{B,β}[m] − Num_{B,β}[Pos_{B,α}[i]] } .

This gives us, for every i, the length of an LCWIS of type α^{i}β^{∗}. The longest of them over all i are the LCWIS’s of the two sequences. The total time is O(m).

3.3 Three-letter case – split diagrams

The naïve extension of the above approach to three letters would have to deal with a quadratic number of tentative exponent pairs (i, j) for subsequences of type α^{i}β^{j}γ^{∗}. We somehow need to avoid the testing of all such pairs. The basis of our almost-linear algorithm for a 3-letter alphabet are what we like to call “split-diagrams,” a data structure that stores information about parts of the given sequences in a compact way.

Assume we were only interested in subsequences of A that have all their αs up to some fixed position s and all their γs strictly after s. Likewise, we only consider subsequences in B with all αs up to some position t and all γs after that. We shall see that under these conditions, with a fixed split between αs and γs, it is possible to find an LCWIS in linear time.

Say, we try and see how long a sequence we can build if we started with exactly i many αs. We determine the ith pair of αs from the left and then count the number of βs in A and B up to the split (s, t). There are p = Num_{A,β}[s] − Num_{A,β}[Pos_{A,α}[i]] such βs in A and q = Num_{B,β}[t] − Num_{B,β}[Pos_{B,α}[i]] in B.

Assume p ≤ q for the moment. For the three values i, p, q, we define a piecewise-linear function f_{i}^{s,t} consisting of a slope-1 segment from (0, i + p) to (q − p, i + q) and a horizontal extension from that point to infinity as shown in the left diagram of Figure 2.

What is the purpose of this function? Assume we tried to find a long common subsequence by matching exactly j many γs in the two sequences. We would align these j pairs as far to the right as possible in order to gain as many βs as possible. So count the number of βs between position s and the leftmost matched γ in A and likewise in B. Say, there are x such βs in A and y in the respective part of B. We can now use our function f_{i}^{s,t} to obtain the length of an LCWIS of type α^{i}β^{∗}γ^{j}: Compute the surplus z = x − y of unmatchable βs in A on the right (assuming x ≥ y for the moment) and read off the function value of f_{i}^{s,t} for that argument. The value f_{i}^{s,t}(z) tells us exactly how long a subsequence we can build to the left of the split if we throw in a surplus of z βs into A.

For example, with no extra βs from the right, we only get min(p, q) = p many pairs of βs, which together with the i αs yield a sequence of length f_{i}^{s,t}(0) =

i + p i + q 0  q − p f_{i}^{s,t}

Figure 2: Split diagrams.

i + p. If we have q − p free βs on the right, we could get a sequence of length f(q − p) = i + q. More βs would not bring an advantage, which is expressed in the stagnation of the function f beyond q − p. The case q > p, which we had originally excluded for cleaner presentation, is simply covered by a function f̄_{i}^{s,t}, defined in the obvious way to handle free βs on the right of the split in sequence B.

Of course, we have not gained anything yet from the function f_{i}^{s,t}. The trick is now to draw the functions f_{i}^{s,t} for all values of i into one diagram. Their point-wise maximum f^{s,t}, the upper envelope of their plots, indicated in the right of Figure 2, gives us the best possible length to the left of the split for any surplus of βs from the right.

Lemma 3 Amongst all subsequences that have all their αs to the left and all γs to the right of a fixed split (s, t), we can find an LCWIS in linear time.

Proof. For 0 ≤ i ≤ min(Num_{A,α}[s], Num_{B,α}[t]), “draw” all functions f_{i}^{s,t} into one split diagram. One can build an array of function values of the upper envelope in O(n) time. After that, test for each right-aligned match of γs, how many βs match to the right of the split (s, t) and evaluate the envelope function for the respective surplus of βs. Actually constructing an LCWIS once its length is known is an easy task. □

In order to turn the split technique into a fast algorithm for the general case, where we do not have any pre-knowledge about good splits, we will have to refine it a little further. If we know that there is an LCWIS with many βs, we can apply Lemma 3 immediately.

Theorem 2 For two length-n sequences over three letters α < β < γ, we can find an LCWIS that contains at least rn many βs (r ∈ (0, 1)) in O(n/r^{2}) time.

Proof. Put a marker every rn positions in A and also in B. Test all ⌈1/r⌉^{2} candidate splits at marker pairs. Any α^{∗}β^{⌈rn⌉}β^{∗}γ^{∗} subsequence must cover at least one of those pairs with its β-section. Hence we will find it. □

3.4 A hierarchy of splits

In the general case, when we need to make sure that we identify subsequences with only a few βs, we need a few tricks to further reduce the number of splits. To this end, first note that we may restrict attention to splits (s, t) that are given by left-aligned α-matches: The collection S of all splits of the form (Pos_{A,α}[i], Pos_{B,α}[i]) suffices to find an LCWIS.

Note that S comes with a natural linear order since no two of its splits cross and hence, |S| = O(n). Yet, if we drew a complete split diagram for every split in S, we would still face a quadratic running-time. To reduce the work further, we avoid drawing complete diagrams for all splits but spread information over splits. Therefore, assign levels to the splits in S: let the level of the ith split (counting from left) be the index of the least significant bit equal to one in the binary representation of i. This scheme has the nice property that between any two splits on the same level there lies another split on a higher level.

Conceptually, our algorithm proceeds in two sweeps over the sequences. In the first sweep it constructs a split diagram for each of the splits in S. However, not all left-side configurations are entered into all diagrams. For each integer i, match the first i αs from A and B and enter the corresponding functions into the split diagram of the closest split (s, t) to the right on each level. This means that the effect of starting with exactly i αs is entered into O(log |S|) = O(log n) diagrams.

After all diagrams are prepared, the algorithm makes a second sweep of the sequences forming all right-aligned matches of γs. For each such partial subsequence we then query the split diagrams for the closest split to the left on each level to obtain the maximum length of an LCWIS with these many γs. A formal description of the algorithm is given in Figure 4 in the appendix.

Theorem 3 We can find an LCWIS of two three-letter sequences of lengths m and n, with m ≥ n, in O(m + n log n) time.

Proof. Consider an LCWIS embedded in A and B with its i αs left-aligned and j γs right-aligned so that the α/γ-free region is as large as possible. There is a unique highest S-split S within this region. Hence, the match of i αs is entered into the diagram of S and the match of j γs queries this diagram; so we are bound to detect the sequence. Again, it is easy to see how to construct an LCWIS once a suitable j is known.

The preprocessing phase takes O(m) time. The rest can be performed in O(m + n log n) time as follows. During the first sweep we simply create a list of O(n log n) quadruples (i, p, q, s) that represent the contents of the O(n) splitters: s is the identity of a splitter, and (i, p, q) are the parameters that define one of the functions illustrated in the left of Figure 2. Similarly, during the second sweep we construct a list of O(n log n) quadruples (i, p, q, s) where (i, p, q) is a query and s is the splitter on which it is to be performed. After bucket-sorting each list, all queries can be answered by a simultaneous linear scan of the lists.

□

4 Multiple Sequences

In this section we consider the problem of finding an LCIS of k length-n sequences, for k ≥ 3. We will denote the sequences by A^{1} = (a_{1}^{1}, …, a_{n}^{1}), A^{2} = (a_{1}^{2}, …, a_{n}^{2}), …, A^{k} = (a_{1}^{k}, …, a_{n}^{k}).

A match is a vector (i_{1}, i_{2}, …, i_{k}) of indices such that a_{i_{1}}^{1} = a_{i_{2}}^{2} = ⋯ = a_{i_{k}}^{k}. Let r be the number of matches. Chan et al. [3] showed that an LCIS can be found in O(min(kr^{2}, kr log σ log^{k−1} r) + kSort_{Σ}(n)) time (they show two algorithms, each corresponding to one of the terms in the min). We present a simpler solution which replaces the second term by O(r log^{k−1} r log log r).

We denote the ith coordinate of a vector v by v[i], and the alphabet symbol corresponding to the match described by a vector v will be denoted s(v). A vector v dominates a vector v′ if v[i] > v′[i] for all 1 ≤ i ≤ k, and we write v′ < v. Clearly, an LCIS corresponds to a sequence v_{1}, …, v_{ℓ} of matches such that v_{1} < v_{2} < ⋯ < v_{ℓ} and s(v_{1}) < s(v_{2}) < ⋯ < s(v_{ℓ}).

To find an LCIS, we use a data structure by Gabow et al. [5, Theorem 3.3], which stores a fixed set of n vectors from {1, …, n}^{k}. Initially all vectors are inactive. The data structure supports the following two operations:

1. Activate a vector with an associated priority.

2. A query of the form “what is the maximum priority of an active vector that is dominated by a vector p?”

A query takes O(log^{k−1} n log log n) time, and the total time for at most n activations is O(n log^{k−1} n log log n). The data structure requires O(n log^{k−1} n) preprocessing time and space.

Each of the r matches v = (v_{1}, …, v_{k}) corresponds to a vector. The priority of v will be the length of the longest LCIS that ends at the match v. We will consider the matches by non-decreasing order of their symbols. For each symbol s of the alphabet, we first compute the priority of every match v with s(v) = s. This is equal to 1 plus the maximum priority of a vector dominated by v. Then, we activate these vectors in the data structure with the priorities we have computed; they should be there when we compute the priorities for matches v with s(v) > s.

The algorithm applies to the case of a common weakly-increasing subsequence by the following modification: The matches will be considered by nondecreasing order of s(v) as before, but within each symbol also in non-decreasing lexicographic order of v. For each match, we compute its priority and immediately activate it in the data structure (so that it is active when considering other matches with the same symbol). The lexicographic order ensures that if v > v′ then v′ is in the data structure when v is considered.

Theorem 4 An LCIS or LCWIS of k length-n sequences can be computed in O(r log^{k−1} r log log r) time, where r counts the number of match vectors.

5 Outlook

The central question about the LCS problems is, whether it can be solved in O(n^{2−ϵ}) time in general. It seems that with LCIS we face the same frontier. Our new algorithms provide solutions fast in many situations, but in general, we do not obtain subquadratic running-time, either.

On the other hand, LCWIS seems to behave very different from the other two problems. Our result shows that it behaves somewhat like a mixture of LCS and LCIS. While already the 2-letter problem is unsolved for LCS, finite alphabets are trivial for LCIS. With LCWIS now, we present almost-linear solutions for alphabets with up to three letters, while it is unclear whether similar results can be obtained for all finite alphabets.

References

[1] David Aldous and Persi Diaconis. Longest increasing subsequences: From patience sorting to the Baik-Deift-Johansson theorem. Bulletin (New Series) of the American Mathematical Society, 36(4):413–432, 1999.

[2] L. Bergroth, H. Hakonen, and T. Raita. A survey of longest common subsequence algorithms. In SPIRE ’00: Proceedings of the Seventh International Symposium on String Processing Information Retrieval (SPIRE’00), pages 39–48. IEEE Computer Society, 2000.

[3] W.-T. Chan, Y. Zhang, S. P.Y. Fung, D. Ye, and H. Zhu. Efficient Algorithms for Finding A Longest Common Increasing Subsequence. In 16th Annual International Symposium on Algorithms and Computation (ISAAC 2005), 2005.

[4] M.L. Fredman. On computing the length of longest increasing subsequences. Discrete Mathematics, 11(1):29–35, 1975.

[5] Harold N. Gabow, Jon Louis Bentley, and Robert E. Tarjan. Scaling and related techniques for geometry problems. In STOC ’84: Proceedings of the sixteenth annual ACM symposium on Theory of computing, pages 135–143, New York, NY, USA, 1984. ACM Press.

[6] D. S. Hirschberg. A linear space algorithm for computing maximal common subsequences. Commun. ACM, 18(6):341–343, 1975.

[7] J. W. Hunt and T. G. Szymanski. A fast algorithm for computing longest common subsequences. Commun. ACM, 20(5):350–353, 1977.

[8] W.J. Masek and M.S. Paterson. A faster algorithm computing string edit distances. J. Comput. System Sci., 20:18–31, 1980.

[9] E. M. McCreight. Priority search trees. SIAM Journal on Computing, 14(2):257–276, 1985.

[10] P. van Emde Boas, R. Kaas, and E. Zijlstra. Design and implementation of an efficient priority queue. Mathematical Systems Theory, 10:99–127, 1977.

[11] Robert A. Wagner and Michael J. Fischer. The string-to-string correction problem. J. ACM, 21(1):168–173, 1974.

[12] I.-H. Yang, C.-P. Huang, and K.-M. Chao. A fast algorithm for computing a longest common increasing subsequence. Information Processing Letters, 93/5:249–253, 2005.

Appendix

Function LCIS(A = (a_{1}, … a_{n}), B = (b_{1}, … b_{m}))

Preprocess (* Clean A and B and build Occ_{s} for every s *) i ← 1

(* Compute L_{1}[1 … n] *) for j = 1 to n do L_{i}[j] ← MinimumKey(Occ_{a_{j}})

(* Main loop *) do

H ← [] (* Empty Bounded Heap *) i ← i + 1 for j = 1 to n do L_{i}[j] ← ∞ (j′, κ′) ← BoundedMin(H, a_{j}) if (j′, κ′) ≠ “invalid” then

L_{i}[j] ← min{κ : κ ∈ Occ_{a_{j}} ∧ κ > κ′}

Link_{i}[j] = j′

endif

if L_{i−1}[j] ≠ ∞ then (* Recall that DecreasePriority inserts a_{j} if it is not already there *) DecreasePriority(H, a_{j}, L_{i−1}[j], (j, L_{i−1}[j]))

endif endfor

while i < n and L_{i} ≠ ∞^{n}

(* Generate a LCIS in reverse order *) if L_{i} = ∞^{n} then i ← i − 1 j ← an index such that L_{i}[j] ≠ ∞ while i > 0 do

output a_{j}

j ← Link_{i}[j]
i ← i − 1

end while end

Figure 3: An O((m + nℓ) log log σ + Sort_{Σ}(m)) time LCIS algorithm for k = 2.

Function LCWIS_{3} (A = (a_{1}, …, a_{n}), B = (b_{1}, …, b_{m}))

Preprocess (* create arrays Num_{.,.}[] and Pos_{.,.}[] *)

t ← min(Num_{A,α}[n], Num_{B,α}[m] (* the size of S *) h := ⌊log_{2} t⌋ (* the highest level *) for i = 1 to t do create empty split diagrams D_{i} and D̄_{i}

for position (μ_{i}, ν_{i}) ← (Pos_{A,α}[i], Pos_{B,α}[i]) with level max{r : 2^{r}|i}

(* First sweep: filling the diagrams *) for i = 1 to t do

for r = 0 to h do

d ← index of closest level-r diagram to the right (≥)

of (Pos_{A,α}[i], Pos_{B,α}[i]) p ← Num_{A,β}[μ_{d}] − Num_{A,β}[Pos_{A,α}[i]] q ← Num_{B,β}[ν_{d}] − Num_{B,β}[Pos_{B,α}[i]] enter triple (i, p, q) into D_{d} and triple (i, q, p) into D̄_{d}

od od

Preprocess all D_{i} and D̄_{i} for quick look-up

(* Second sweep: reading the diagrams *) best ← 0 for i = 0 to min(Num_{A,γ}[n], Num_{B,γ}[m]) − 1 do

for r = 0 to h do

d ← index of closest level-r diagram to the left (<)

of (Pos_{A,γ}[n − i], Pos_{B,γ}[m − i]) x ← Num_{A,β}[Pos_{A,γ}[n − i]] − Num_{A,β}[μ_{d}] y ← Num_{B,β}[Pos_{B,γ}[m − i]] − Num_{B,β}[ν_{d}] length ← i + min(x, y) + max(D_{d}(x, y), D̄_{d}(y, x)) if length > best then best ← length

od od

return best

Figure 4: A 3-letter LCWIS algorithm.

Recent BRICS Report Series Publications

RS-05-37 Gerth Stølting Brodal, Kanela Kaligosi, Irit Katriel, and Martin Kutz. Faster Algorithms for Computing Longest Common Increasing Subsequences. December 2005. 16 pp.

RS-05-36 Dariusz Biernacki, Olivier Danvy, and Chung-chieh Shan. On the Static and Dynamic Extent of Delimited Continuations. December 2005. ii+33 pp. To appear in the journal Science of Computer Programming. Supersedes BRICS RS-05-13.

RS-05-35 Kristian Støvring. Extending the Extensional Lambda Calculus with Surjective Pairing is Conservative. November 2005. 19 pp.

RS-05-34 Henning Korsholm Rohde. Formal Aspects of Polyvariant Specialization. November 2005. 27 pp.

RS-05-33 Luca Aceto, Willem Jan Fokkink, Anna Ingólfsdóttir, and Sumit Nain. Bisimilarity is not Finitely Based over BPA with Interrupt. October 2005. 33 pp. This paper supersedes BRICS Report RS-04-24. An extended abstract of this paper appeared in Algebra and Coalgebra in Computer Science, 1st Conference, CALCO 2005, Swansea, Wales, 3–6 September 2005, Lecture Notes in Computer Science 3629, pp. 54–68, Springer-Verlag, 2005.

RS-05-32 Anders Møller, Mads Østerby Olesen, and Michael I. Schwartzbach. Static Validation of XSL Transformations. October 2005. 50 pp.

RS-05-31 Christian Kirkegaard and Anders Møller. Type Checking with XML Schema in Xact. September 2005. 20 pp.

RS-05-30 Karl Krukow. An Operational Semantics for Trust Policies. September 2005. 38 pp.

RS-05-29 Olivier Danvy and Henning Korsholm Rohde. On Obtaining the Boyer-Moore String-Matching Algorithm by Partial Evaluation. September 2005. ii+9 pp. To appear in Information Processing Letters. This version supersedes BRICS RS-05-14.

RS-05-28 Jiří Srba. On Counting the Number of Consistent Genotype Assignments for Pedigrees. September 2005. 15 pp. To appear in FSTTCS ’05.

RS-05-27 Pascal Zimmer. A Calculus for Context-Awareness. August 2005. 21 pp.