对齐多个有序序列的最有效算法

Most Efficient Algorithm to Align an Multiple Ordered Sequences

我有一种奇怪的感觉,这是一个非常容易解决的问题,但我找不到不使用蛮力或动态编程的好方法。开始了:

Given N arrays of ordered and monotonic values, find the set of positions for each array i1, i2 ... in that minimises pair-wise difference of values at those indexes between all arrays. In other words, find the positions for all arrays whose values are closest to each other. Multiple solutions may exist and arrays may or may not be equally sized.

If A denotes the list of all arrays, the pair-wise difference is given by the sum of absolute differences between all values at the given indexes between all different arrays, as so:

一个例子,3个数组abc:

a = [20 29 30 32 33]
b = [28 29 30 32 33]
c = [10 12 28 31 32 33]

此数组的最佳对齐方式是 a[3] b[3] c[4]a[4] b[4] c[5],因为 (32,32,32) 和 (33,33,33) 都是相等的值,因此具有最小值两两之间的差异。 (假设数组索引从 0 开始)

这是生物信息学中的一个常见问题,通常用 Dynamic Programming 解决,但由于这是一个有序序列,我认为有某种方法可以利用这种顺序概念。我首先考虑成对地这样做,但这并不能保证全局最优,因为最好的局部答案可能不是最好的全局答案。

这意味着与语言无关,但我真的不介意针对特定语言的答案,只要不失一般性即可。我知道 Dynamic Programming 是一个选项,但我觉得有更简单的方法吗?

棘手的事情是解析数组,以便在某些时候保证您考虑实现成对最小值的索引集。对值使用最小堆不起作用。具有 4 个数组的反例:[0,5]、[1,2]、[2]、[2]。我们从 d(0,1,2,2) = 7 开始,最优是 d(0,2,2,2) = 6,但是最小堆将我们从 7 移动到 d(5,1,2,2) = 12,则 d(5,2,2,2) = 9.

我相信(但尚未证明)如果我们总是增加最大程度地改善成对距离(或最小程度地降低成对距离)的索引,我们可以保证访问每个局部最小值和全局最小值。

假设 k 个数组中共有 n 个元素:

简单的方法:我们重复获得成对距离增量(delta wrt。增加每个索引),增加最好的一个,并且任何时候这样做都会将我们从改进切换到退化(即局部最小值)我们计算成对的距离。所有这些都是每个增量的 O(k^2),总共 运行 时间 O((n-k) * (k^2)).

使用 O(k^2) 存储,我们可以保留一个数组,其中 (i,j) 存储通过增加数组 i wrt 的索引实现的成对距离增量。数组 j。我们还存储列总和。然后在递增索引时,我们可以在 O(k) 中更新适当的行&列&列总和。这给了我们 O((n-k)*k)

的 运行 时间

为了完成 Dave 的回答,这里是 delta 算法的伪代码:

initialise index_table to 0's where each row i denotes the index for the ith array
initialise delta_table with the corresponding cost of incrementing index of ith array and keeping the other indexes at their current values
cur_cost <- cost of current index table
best_cost <- cur_cost
best_solutions <- list with the current index table
while (can_at_least_one_index_increase)
    i <- index whose delta is lowest
    increment i-th entry of the index_table
    if cost(index_table) < cur_cost
        cur_cost = cost(index_table)
        best_solutions = {} U {index_table}
    if cost(index_table) = cur_cost
         best_solutions = best_solutions U {index_table}
    update delta_table

重要说明:在迭代期间,某些 index_table 条目可能已经达到该数组的最大值。每当更新 delta_table 时,永远不要选择这些值,否则这将导致 Array Out of BoundsSegmentation Fault 或未定义的行为。一个巧妙的技巧是简单地检查哪些索引已经达到最大值并设置一个 足够大的 值,这样它们就不会被选中。如果没有索引可以再增加,循环将结束。


这是 Python 中的一个实现:

def align_ordered_sequences(arrays: list):
    def get_cost(index_table):
        n = len(arrays)
        if n == 1:
            return 0
        sum = 0
        for i in range(0, n-1):
            for j in range(i+1, n):
                v1 = arrays[i][index_table[i]]
                v2 = arrays[j][index_table[j]]
                sum += math.sqrt((v1 - v2) ** 2)
        return sum

    def compute_delta_table(index_table):
        # Initialise the delta table: we switch each index element to 1, call
        # the cost method and then revert the change, this avoids having to
        # create copies, which decreases performance unnecessarily
        delta_table = []
        for i in range(n):
            if index_table[i] + 1 >= len(arrays[i]):
                # Implementation detail: if the index is outside the bounds of
                # array i, choose a "large enough" number
                delta_table.append(999999999999999)
            else:
                index_table[i] = index_table[i] + 1
                delta_table.append(get_cost(index_table))
                index_table[i] = index_table[i] - 1
        return delta_table

    def can_at_least_one_index_increase(index_table):
        answer = False
        for i in range(len(arrays)):
            if index_table[i] < len(arrays[i]) - 1:
                answer = True
        return answer

    n = len(arrays)

    index_table = [0] * n
    delta_table = compute_delta_table(index_table)
    best_solutions = [index_table.copy()]
    cur_cost = get_cost(index_table)
    best_cost = cur_cost

    while can_at_least_one_index_increase(index_table):
        i = delta_table.index(min(delta_table))
        index_table[i] = index_table[i] + 1

        new_cost = get_cost(index_table)
        # A new best solution was found
        if new_cost < cur_cost:
            cur_cost = new_cost
            best_solutions = [index_table.copy()]
        # A new solution with the same cost was found
        elif new_cost == cur_cost:
            best_solutions.append(index_table.copy())
        # Update the delta table
        delta_table = compute_delta_table(index_table)

    return best_solutions

这里有一些例子:

>>> print(align_ordered_sequences([[0,5], [1,2], [2], [2]]))
[[0, 1, 0, 0]]

>> print(align_ordered_sequences([[3, 5, 8, 29, 40, 50], [1, 4, 14, 17, 29, 50]]))
[[3, 4], [5, 5]]

注2:这里输出的索引不是每个数组的实际值。