Finding the maximum overall subsequence
In this article, I would like to review popular algorithms for solving the problem of finding the maximum common subsequence or LCS (longest common sequense). Since the emphasis is on learning, and not on real use, Python is chosen as the language for implementation, which will reduce the amount of code and concentrate on the main ideas.
A sequence is an ordered set of elements. A string is a special case of a sequence, further examples will be considered just for simplicity strings, but without changes they can be used for arbitrary text or anything else consistent.
Let there be a sequence x consisting of elements x 1 x 2 ... x m and a sequence y consisting of elements y 1 y 2 ... y n . z is a subsequence of x if there is a strictly increasing set of indices of elements xfrom which z is obtained .
A common subsequence for x and y is a sequence z that is both a subsequence of x and a subsequence of y .
The maximum total subsequence is the total subsequence with the maximum length. Further in the text we will use the abbreviation LCS .
As an example, let x = HA B R AHA BR , y = HARB OU R , in this case LCS (x, y) = HARBR. You can already go directly to the LCS calculation algorithm , but it would be nice to understand why we might need this.
The most common use is the use in file comparison programs, such as GNU diff. Having found the LCS for two texts, compiling a list of elementary changes to turn x into y or vice versa is a trivial task. As a bonus, based on the length of the overall subsequence, you can specify a metric to determine the similarity of the two sequences. That's it, now you can definitely get down to business.
To start a couple of observations:
These observations are already enough to implement recursion:
Now you might think that this is not the case with this implementation. In the worst case, when there are no identical elements between x and y, LCS_RECURSIVE is called 2 len (x) + len (y) times, at the recursion level len (x) + len (y). Even if you close your eyes to performance, the code:
without an additional call, sys.setrecursionlimit will fail. Not the most successful implementation.
The algorithm in question is also known as the Needleman-Wunsch algorithm.
The whole approach boils down to phasing in the matrix, where the rows are x elements and the columns are y elements. When filling in, two rules follow from the observations already made:
1. If the element x i is equal to y j, then in the cell (i, j) the value of the cell (i-1, j-1) is written with the addition of one
2. If the element x i is not equal to y j, then the maximum of the values (i-1, j) and (i, j-1) is recorded in the cell (i, j).
Filling takes place in a double cycle in i and j with increasing values by one, so at each iteration the cell values needed at this step are already calculated:
Illustration of what is happening:
Cells in which an increase in the value directly occurred are highlighted. After filling in the matrix, connecting these cells, we get the desired LCS. In this case, you need to connect while moving from maximum to minimum indices, if the cell is highlighted, then add the corresponding element to the LCS, if not, then move up or to the left, depending on where the larger value is in the matrix:
The complexity of the algorithm is O (len (x) * len (y)), the same memory score. Thus, if I want to compare line by line two files of 100,000 lines, then I will need to store a matrix of 10 10 elements in memory . Those. actual use threatens to receive a MemoryError almost out of the blue. Before moving on to the next algorithm, we note that when filling in the matrix L at each iteration over the elements x, we need only the row obtained in the previous move. Those. if the task were limited only to finding the length of the LCS without the need to calculate the LCS itself, then it would be possible to reduce the memory usage to O (len (y)), saving only two rows of the matrix L.
The idea behind this algorithm is simple: if you divide the input sequence x = x 1 x 2 ... x m into two arbitrary parts for any boundary index i by xb = x 1 x 2 ... x i and xe = x i + 1 x i + 2 ... x m , then there is a way to similarly partition y (there is such an index j dividing y into yb = y 1 y 2 ... y j and ye = y j + 1 y j + 2 ... y n ) such that LCS (x, y) = LCS (xb, yb) + LCS (xe, ye).
In order to find this partition y is proposed:
This can be represented as filling the matrix L from two opposite ends:
Note that since there is a need only in the last row of the matrix L, then when calculating the whole matrix is not needed. Having slightly changed the implementation of fill_dyn_matrix we get:
Now directly, about the algorithm itself. It is a classic divide and conquer: while there are elements in the sequence x, we divide x in half, find the appropriate partition for y and return the sum of the recursive calls for pairs of sequences (xb, yb) and (xe, ye). Note that in the trivial case, if x consists of one element and occurs in y, we simply return a sequence from this single element x.
Now we have O (len (y)) memory requirements, the time required for calculation has doubled, but O (len (x) len (y)) is still asymptotically.
The first thing we need is to create a hash of the matchlist table, which will contain the indices of the elements of the second sequence. The time required for this is O (len (y)). The following python code does this:
For the HARBOR sequence, the hash will be as follows {'A': [1], 'B': [3], 'H': [0], 'O': [4], 'R': [2, 6] , 'U': [5]}.
Further, sorting out the elements of the sequence x, fill the THRESH array with the corresponding indices from the prepared matchlist, so that the value of the k-th element of THRESH should be the index y_index, provided that THRESH [k-1] <y_index and y_index <THRESH [k]. Thus, at any moment in time, the THRESH array is sorted and we can use binary search to find a suitable k. When updating the THRESH element, we also add the sequence element corresponding to the y_index index to our LCS. The following code can clarify:
After that, we can only collect the sequence itself from LINK_x.
The running time of this algorithm is O ((r + n) log n), where n is the length of the larger sequence and r is the number of matches, in the worst case, with r = n * n, we get a worse working time than in the previous solution. Memory requirements remain of the order of O (r + n).
There are a lot of algorithms that solve this problem, asymptotically, the most efficient algorithm (Masek and Paterson) has an O (n * n / log n) time estimate. Given the overall small efficiency in LCS calculations, in practice, the simplest preparations are performed before the algorithm runs, such as discarding identical elements at the beginning and at the end of sequences and searching for trivial differences between sequences. Also, there are some optimizations using bit operations that do not affect the asymptotic behavior of the operating time.
download all code with examples
Definitions
A sequence is an ordered set of elements. A string is a special case of a sequence, further examples will be considered just for simplicity strings, but without changes they can be used for arbitrary text or anything else consistent.
Let there be a sequence x consisting of elements x 1 x 2 ... x m and a sequence y consisting of elements y 1 y 2 ... y n . z is a subsequence of x if there is a strictly increasing set of indices of elements xfrom which z is obtained .
A common subsequence for x and y is a sequence z that is both a subsequence of x and a subsequence of y .
The maximum total subsequence is the total subsequence with the maximum length. Further in the text we will use the abbreviation LCS .
As an example, let x = HA B R AHA BR , y = HARB OU R , in this case LCS (x, y) = HARBR. You can already go directly to the LCS calculation algorithm , but it would be nice to understand why we might need this.
Practical application
The most common use is the use in file comparison programs, such as GNU diff. Having found the LCS for two texts, compiling a list of elementary changes to turn x into y or vice versa is a trivial task. As a bonus, based on the length of the overall subsequence, you can specify a metric to determine the similarity of the two sequences. That's it, now you can definitely get down to business.
First approach or folk art
To start a couple of observations:
- If for sequences x and y we already calculated LCS (x, y) = z, then LCS for sequences obtained from x and y by adding the same element will consist of z and this added element
- If we add one different element to the sequences x and y, then for the resulting xa and yb, LCS should be the larger of the two: LCS (x, yb) or LCS (xa, y)
These observations are already enough to implement recursion:
defLCS_RECURSIVE(x, y):if len(x) == 0or len(y) == 0:
return []
if x[-1] == y[-1]:
return LCS_RECURSIVE(x[:-1], y[:-1]) + [x[-1]]
else:
left = LCS_RECURSIVE(x[:-1], y)
right = LCS_RECURSIVE(x, y[:-1])
return left if len(left) > len(right) else right
Now you might think that this is not the case with this implementation. In the worst case, when there are no identical elements between x and y, LCS_RECURSIVE is called 2 len (x) + len (y) times, at the recursion level len (x) + len (y). Even if you close your eyes to performance, the code:
from uuid import uuid4
x = [uuid4().hex for _ in xrange(500)]
y = [uuid4().hex for _ in xrange(500)]
print LCS_RECURSIVE(x,y)
without an additional call, sys.setrecursionlimit will fail. Not the most successful implementation.
Dynamic programming or 64 kb is enough for everyone
The algorithm in question is also known as the Needleman-Wunsch algorithm.
The whole approach boils down to phasing in the matrix, where the rows are x elements and the columns are y elements. When filling in, two rules follow from the observations already made:
1. If the element x i is equal to y j, then in the cell (i, j) the value of the cell (i-1, j-1) is written with the addition of one
2. If the element x i is not equal to y j, then the maximum of the values (i-1, j) and (i, j-1) is recorded in the cell (i, j).
Filling takes place in a double cycle in i and j with increasing values by one, so at each iteration the cell values needed at this step are already calculated:
deffill_dyn_matrix(x, y):
L = [[0]*(len(y)+1) for _ in xrange(len(x)+1)]
for x_i,x_elem in enumerate(x):
for y_i,y_elem in enumerate(y):
if x_elem == y_elem:
L[x_i][y_i] = L[x_i-1][y_i-1] + 1else:
L[x_i][y_i] = max((L[x_i][y_i-1],L[x_i-1][y_i]))
return L
Illustration of what is happening:
Cells in which an increase in the value directly occurred are highlighted. After filling in the matrix, connecting these cells, we get the desired LCS. In this case, you need to connect while moving from maximum to minimum indices, if the cell is highlighted, then add the corresponding element to the LCS, if not, then move up or to the left, depending on where the larger value is in the matrix:
defLCS_DYN(x, y):
L = fill_dyn_matrix(x, y)
LCS = []
x_i,y_i = len(x)-1,len(y)-1while x_i >= 0and y_i >= 0:
if x[x_i] == y[y_i]:
LCS.append(x[x_i])
x_i, y_i = x_i-1, y_i-1elif L[x_i-1][y_i] > L[x_i][y_i-1]:
x_i -= 1else:
y_i -= 1
LCS.reverse()
return LCS
The complexity of the algorithm is O (len (x) * len (y)), the same memory score. Thus, if I want to compare line by line two files of 100,000 lines, then I will need to store a matrix of 10 10 elements in memory . Those. actual use threatens to receive a MemoryError almost out of the blue. Before moving on to the next algorithm, we note that when filling in the matrix L at each iteration over the elements x, we need only the row obtained in the previous move. Those. if the task were limited only to finding the length of the LCS without the need to calculate the LCS itself, then it would be possible to reduce the memory usage to O (len (y)), saving only two rows of the matrix L.
Fighting for a place or Hirschberg Algorithm
The idea behind this algorithm is simple: if you divide the input sequence x = x 1 x 2 ... x m into two arbitrary parts for any boundary index i by xb = x 1 x 2 ... x i and xe = x i + 1 x i + 2 ... x m , then there is a way to similarly partition y (there is such an index j dividing y into yb = y 1 y 2 ... y j and ye = y j + 1 y j + 2 ... y n ) such that LCS (x, y) = LCS (xb, yb) + LCS (xe, ye).
In order to find this partition y is proposed:
- Fill dynamic matrix L with fill_dyn_matrix for xs and y. L1 we equate the last row of the calculated matrix L
- Fill in the matrix L for * xe and * y (inverse sequences for xe and y, in terms of Python: list (reversed (x)), list (reversed (y))). We equate * L2 the last row of the calculated matrix L, L2 is the river from * L2
- Take j equal to the index for which the sum L1 [j] + L2 [j] is maximum
This can be represented as filling the matrix L from two opposite ends:
Note that since there is a need only in the last row of the matrix L, then when calculating the whole matrix is not needed. Having slightly changed the implementation of fill_dyn_matrix we get:
deflcs_length(x, y):
curr = [0]*(1 + len(y))
for x_elem in x:
prev = curr[:]
for y_i, y_elem in enumerate(y):
if x_elem == y_elem:
curr[y_i + 1] = prev[y_i] + 1else:
curr[y_i + 1] = max(curr[y_i], prev[y_i + 1])
return curr
Now directly, about the algorithm itself. It is a classic divide and conquer: while there are elements in the sequence x, we divide x in half, find the appropriate partition for y and return the sum of the recursive calls for pairs of sequences (xb, yb) and (xe, ye). Note that in the trivial case, if x consists of one element and occurs in y, we simply return a sequence from this single element x.
defLCS_HIRSHBERG(x, y):
x_len = len(x)
if x_len == 0:
return []
elif x_len == 1:
if x[0] in y:
return [x[0]]
else:
return []
else:
i = x_len // 2
xb, xe = x[:i], x[i:]
L1 = lcs_length(xb, y)
L2 = reversed(lcs_length(xe[::-1], y[::-1]))
SUM = (l1 + l2 for l1,l2 in itertools.izip(L1, L2))
_, j = max((sum_val, sum_i) for sum_i, sum_val in enumerate(SUM))
yb, ye = y[:j], y[j:]
return LCS_HIRSHBERG(xb, yb) + LCS_HIRSHBERG(xe, ye)
Now we have O (len (y)) memory requirements, the time required for calculation has doubled, but O (len (x) len (y)) is still asymptotically.
Hunt-Szymanski Algorithm
The first thing we need is to create a hash of the matchlist table, which will contain the indices of the elements of the second sequence. The time required for this is O (len (y)). The following python code does this:
y_matchlist = {}
for index, elem in enumerate(y):
indexes = y_matchlist.setdefault(elem, [])
indexes.append(index)
y_matchlist[elem] = indexes
For the HARBOR sequence, the hash will be as follows {'A': [1], 'B': [3], 'H': [0], 'O': [4], 'R': [2, 6] , 'U': [5]}.
Further, sorting out the elements of the sequence x, fill the THRESH array with the corresponding indices from the prepared matchlist, so that the value of the k-th element of THRESH should be the index y_index, provided that THRESH [k-1] <y_index and y_index <THRESH [k]. Thus, at any moment in time, the THRESH array is sorted and we can use binary search to find a suitable k. When updating the THRESH element, we also add the sequence element corresponding to the y_index index to our LCS. The following code can clarify:
x_length, y_length = len(x), len(y)
min_length = min(x_length, y_length)
THRESH = list(itertools.repeat(y_length, min_length+1))
LINK_s1 = list(itertools.repeat(None, min_length+1))
LINK_s2 = list(itertools.repeat(None, min_length+1))
THRESH[0], t = -1, 0for x_index,x_elem in enumerate(x):
y_indexes = y_matchlist.get(x_elem, [])
for y_index in reversed(y_indexes):
k_start = bisect.bisect_left(THRESH, y_index, 1, t)
k_end = bisect.bisect_right(THRESH, y_index, k_start, t)
for k in xrange(k_start, k_end+2):
if THRESH[k-1] < y_index and y_index < THRESH[k]:
THRESH[k] = y_index
LINK_x[k] = (x_index, LINK_x[k-1])
if k > t:
t = k
After that, we can only collect the sequence itself from LINK_x.
The running time of this algorithm is O ((r + n) log n), where n is the length of the larger sequence and r is the number of matches, in the worst case, with r = n * n, we get a worse working time than in the previous solution. Memory requirements remain of the order of O (r + n).
Total
There are a lot of algorithms that solve this problem, asymptotically, the most efficient algorithm (Masek and Paterson) has an O (n * n / log n) time estimate. Given the overall small efficiency in LCS calculations, in practice, the simplest preparations are performed before the algorithm runs, such as discarding identical elements at the beginning and at the end of sequences and searching for trivial differences between sequences. Also, there are some optimizations using bit operations that do not affect the asymptotic behavior of the operating time.
download all code with examples
