Most Efficient Algorithm to Align an Multiple Ordered Sequences

The tricky thing is parsing the arrays so that at some point you're guaranteed to be considering the set of indices that realize the pairwise min. Using a min heap on the values doesn't work. Counterexample with 4 arrays: [0,5], [1,2], [2], [2]. We start with a d(0,1,2,2) = 7, optimal is d(0,2,2,2) = 6, but the min heap moves us from 7 to d(5,1,2,2) = 12, then d(5,2,2,2) = 9.

I believe (but haven't proved) that if we alway increment the index that improves pairwise distance the most (or degrades it the least), we're guaranteed to visit every local min and the global min.

Assuming n total elements across k arrays:

Simple approach: we repeatedly get the pairwise distance deltas (delta wrt. incrementing each index), increment the best one, and any time doing so switch us from improvement to degradation (i.e. a local minimum) we calculate the pairwise distance. All this is O(k^2) per increment for a total running time of O((n-k) * (k^2)).

With O(k^2) storage, we could keep an array where (i,j) stores the pairwise distance delta achieve by increment the index of array i wrt. array j. We also store the column sums. Then on incrementing an index we can update the appropriate row & column & column sums in O(k). This gives us a running time of O((n-k)*k)

To just complete Dave's answer, here is the pseudocode of the delta algorithm:

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

Important Note: During an iteration, some index_table entries might have already reached the maximum value for that array. Whenever updating the delta_table, it is necessary to never pick those values, otherwise this will result in a Array Out of Bounds,Segmentation Fault or undefined behaviour. A neat trick is to simply check which indexes are already at max and set a sufficiently large value, so they are never picked. If no index can increase anymore, the loop will end.

Here's an implementation in 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

And here are some examples:

>>> 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]]

Note 2: this outputs indexes not the actual values of each array.