Given two arrays, find the permutations that give closest distance between two arrays

Question

Let's say I have two arrays of the same length n, named A and B.

These two arrays contain real values. We define the distance between two arrays as the mean squared distance.

dist(A,B) = sqrt( sum((A - B)2) )

I want to find the permutation of A that gives the min distance to B. The naive method is to try every permutation of A and record the min distance. However, this method is of complexity O(n!).

Is there an algorithm of complexity less than O(n!)?

Neptunian

For what purpose is this. Do you need the best solution or just a good enough solution? – Thirion

Can you eleborate a bit with an example problem and the answer you expect? – Mowery

I am not certain if I fully get the issue you are describing, sounds like you want to try all the possible permutations of a vector of coordinates to use Pythagoras and find which combination provides the shortest distance. So you are basically building a tree of all possible combinations, for this you would need to consider what is it you are after, degree of precision or if it needs to be exact, how many results you want (can be many), speed [of the algorithm]? Having this in mind what you could do to speed up such process is to trim your tree; again how exact does it need to be? – Saimon

because the other answer has better references i guess. I always appreciate good references. – Neptunian

Answer 1

The problem you describe is equivalent to the Minimum Cost Perfect Matching Problem which can be solved efficiently (and exactly) using The Hungarian Algorithm. In the Minimum Cost Perfect Matching Problem you have an input weighted bipartite graph where the two sets have the same size n, and each edge has a non-negative cost. The goal is to find a perfect matching of minimum cost.

In your case, the bipartite graph is a biclique. That is, every vertex in one set is connected to every vertex in the other set, and the cost of edge (i,j) is (A[i] - B[i])^2 (where i corresponds to index i in array A and j corresponds to index j in array B).

EDIT: This is not the best solution for the problem. Ivo Merchiers came up with a better solution both in terms of efficiency and simplicity. The reason I'm not deleting my answer is because my suggested solution is valuable for distance measures to which Ivo's solution does not apply (as his approach works by exploiting a property of the Euclidean distance).

Pennywise

Correct but inefficient for the problem. See Ivo's solution. – Mouseear

@Dave, thanks for the comment. I've added a reference to Ivo's better solution. – Pennywise

Answer 2

You can just sort both A and B. In that case, the Euclidean distance is minimal.

If B has to remain fixed, then you just need to invert the permutation needed to sort B and apply that to the sorted version of A.

This solution does assume that you want to find just a permutation and not the most simple permutation (since sorting and unsorting through permutations will not be incredibly efficient).

---

Proof: Let S,T be our pair of arrays. We can assume S to be sorted without loss of generality, since all that matters is the mapping between the two sets of elements.

Let T be the permutation that minimizes the distance between the two arrays, and let d be that distance.

Suppose that T is not sorted. Then there exist elements i,j s.t. T_i > T_j

S_i + k1 = S_j
T_i = T_j + k2
where k1,k2 > 0

Let x be the total distance of all elements except i and j.

d = x + (S_i - T_i)^2 + ((S_i + k1) - (T_i - k2))^2

If we swap the order of T_i and T_j, then our new distance is:

d' = x + (S_i - (T_i - k2))^2 + ((S_i + k1) - T_i)^2

Thus: d - d' = 2 * k1 * k2, which contradicts our assumption that T is the permutation that minimizes the distance, so the permutation that does so must be sorted.

Sorting the two arrays can be done in O(n log n) using a variety of methods.

Embezzle

If this approach is wrong, I'd be glad to hear a counter-example, since I'm not sure what the issue with it is and I'm glad to learn something new. – Embezzle

It is mathematically sound. I don't have time to write up a proof, but it's along the lines of: If the best solution is not sorted, find two indices x, y such that A[x] < A[y] and B[x] > B[y]. Swap the two in A, and you have a better solution. – Thirion

Seems correct in my use case. would be great if you have any reference to a proof though. – Neptunian

@WangDuo: It's easy to prove with simple algebra. But I'm not going to try again without MathJax. Ask on Mathematics (or Computer Science) where MathJax is available. – Ninetta

If this can be proven correct, it's superior to the Hungarian Algorithm in the accepted answer; that one's O(n^3). – Starling

Yair's idea works. I've submitted an edit with a proof along those lines. – Starling

Please note, that this is true for the Euclidean distance, but it's not necessarily true for an other distance measure. Just take sqrt(A[i] - B[i]) as a distance measure and take the elements [1 4] , [4 5]. Sorting yields error sqrt(3) + sqrt(1) which is bigger than sqrt(4) + sqrt(0) (which is the error of the ideal permutation). – Conferral

Answer 3

You can just sort A and B and match up the corresponding elements.

Imagine that there are two elements of A, Ai and Aj, corresponding to Bi and Bj.

The error contribution from these matches is:

(Ai-Bi)^2 + (Aj-Bj)^2

= Ai^2 + Bi^2 + Aj^2 + Bj^2 - 2(AiBi + AjBj)

Is it better to swap the matches, or to keep them the way they are?

Well, the difference in the error if we swap them is:

2(AiBi + AjBj) - 2(AiBj + AjBi)

~ AiBi - AiBj + AjBj - AjBi

= Ai(Bi - Bj) - Aj(Bi - Bj)

= (Ai - Aj)(Bi - Bj)

So, if the As and Bs are in the same order, then this product is positive and the error will go up if you swap them. If they are not in the same order, then this product is negative an the error will go down if you swap them.

If you repeatedly swap any pairs that are out of order until there are no such pairs, your error will keep going down, and you'll end up with with the nth largest A matched up with the nth largest B all the way through the array.

Just sorting them and matching them up is therefor optimal, and of course it's way faster than the Hungarian algorithm.

Electric

Answer 4

Construct a bipartite graph from the vectors. Find minimum weight perfect matching in this graph.

How to construct the graph.

1. Let A, B be two parts of the graph. Each with n nodes.
2. Connect i in A to j in B with an edge of weight abs(A[i] - B[j]).

I believe this can be done in O(n^2).

See http://www.cse.iitd.ernet.in/~naveen/courses/CSL851/lec4.pdf

---

If each number in A has only one closest number in B then you can do this in O(n \log n). This probably might be the case since you have the real numbers.

How?

1. Sort A O(n \log n)
2. Binary search for the closest number for each number in B. O(n \log n).

If the numbers are coming from the real world and have even a little bit of randomness to them, then the differences between each pair of the numbers is probably unique. You can verify if this is the case by running experiments on the input vectors. Then the problem becomes easy to solve yay!

Appropriate

Thanks, yours answer is correct! But @Pennywise replied earlier with good references so I decided to accept his. – Neptunian

Answer 5

I needed this in python so I'll just share my solution based on Ivo Merchiers answer here:

target = [12, 14, 4512, 123, 4412]
source = [12, 14, 120, 4413, 5512]

permutationToSortTarget = [i[0] for i in sorted(enumerate(target), key=lambda x: x[1])] # get permutation that would sort target
permutationNeeded = [i[0] for i in sorted(enumerate(permutationToSortTarget), key=lambda x: x[1])] # get needed permutation

source.sort()
source = [source[i] for i in permutationNeeded] # apply permutation to sorted source

Sagunto

Answer 6

I was interested in some variations of this problem, in a different context, and was thinking of solutions along similar lines as Ivo's answer.

So here I give a solution (in JavaScript) for the extended problem which seeks the alignment that minimizes total given distance metric between sequences and which can also handle sequences of unequal length.

Examples
note: alignment is like the permutation (and/or warping) of (parts of) b sequence that minimizes total given distance metric with a sequence

a:0,1,2 b:0,1,2 alignment:0,1,2
a:0,1,2 b:2,1,0 alignment:2,1,0
a:2,1,0 b:0,1,2 alignment:2,1,0
a:2,1,0 b:2,1,0 alignment:0,1,2
a:0,1,2 b:0,1 alignment:0,1,1
a:0,1,2 b:1,0 alignment:1,0,0
a:2,1,0 b:0,1 alignment:1,1,0
a:2,1,0 b:1,0 alignment:0,0,1
a:2,1,0,3,4 b:1,0,2 alignment:2,0,1,2,2
a:2,3,1,4,0 b:1,3,2 alignment:2,1,0,1,0
a:-1,2,3,1,4,0,5 b:1,3,2 alignment:0,2,1,0,1,0,1
a:3,4,2,1,0 b:1,2 alignment:1,1,1,0,0
a:2,1,0,3,4 b:0,2 alignment:1,1,0,1,1
a:0,1,2 b:-2,-1,0,1,2 alignment:2,3,4
a:0,1,2 b:2,1,0,4,3 alignment:2,1,0
a:2,1,0 b:-2,-1,0,1,2 alignment:4,3,2
a:2,1,0 b:2,1,0,4,3 alignment:0,1,2
a:2,1,0 b:0,2,4,1,3 alignment:1,3,0

JS Code:

function cmp_with_indices(cmp)
{
    return function(a, b) {
        var res = cmp(a[0], b[0]);
        return 0 > res ? -1 : (0 < res ? 1 : (a[1] - b[1]));
    };
}
function with_indices(arr)
{
    var n = arr.length, vi = new Array(n), i;
    for (i=0; i<n; ++i) vi[i] = [arr[i], i];
    return vi;
}
function get_indices(vi, inv)
{
    var n = vi.length, i = new Array(n), j;
    if (inv)
    {
        for (j=0; j<n; ++j) i[vi[j][1]] = j;
    }
    else
    {
        for (j=0; j<n; ++j) i[j] = vi[j][1];
    }
    return i;
}
function cmp(a, b)
{
    return a < b ? -1 : (a > b ? 1 : 0);
}
function dist(a, b)
{
    return Math.abs(a - b);
}
function align(A, B, dist_AB, cmp_AA, cmp_BB)
{
    var n = A.length, m = B.length, i, j, k, s, sm, sM, km, jm, perm_A, perm_B, iperm_A, alignment;
    if (n && m)
    {
        // O(NlogN), N = max(n,m)
        // assume that cmp_AA, cmp_BB, dist_AB are "compatible" in that:
        // (0 > cmp_AA(A, A') and 0 > cmp_BB(B, B')) ==> (dist_AB(A, B) + dist_AB(A', B')) <= (dist_AB(A', B) + dist_AB(A, B'))
        // in other words, minimum overall distance is when alignment implies similarly sorted sequences
        dist_AB = dist_AB || dist;
        perm_A = get_indices(with_indices(A).sort(cmp_with_indices(cmp_AA || cmp)));
        perm_B = get_indices(with_indices(B).sort(cmp_with_indices(cmp_BB || cmp)));
        alignment = new Array(n);
        if (n > m)
        {
            iperm_A = new Array(n);
            for (i=0; i<n; ++i)
            {
                iperm_A[perm_A[i]] = i;
            }
            sm = Infinity; km = 0;
            for (k=0; k+m<=n; ++k)
            {
                s = 0;
                for (i=0; i<m; ++i) s += dist_AB(A[perm_A[i+k]], B[perm_B[i]]);
                if (s < sm) {sm = s; km = k;}
            }
            for (i=0; i<m; ++i)
            {
                alignment[perm_A[i+km]] = perm_B[i];
            }
            for (i=0; i<n; ++i)
            {
                if (null == alignment[i])
                {
                    // pad/interpolate
                    k = iperm_A[i]-km;
                    j = k<=0 ? 0 : (k>=m ? (m-1) : k);
                    s = dist_AB(A[i], B[perm_B[j]]);
                    sm = j-1>=0 ? dist_AB(A[i], B[perm_B[j-1]]) : Infinity;
                    sM = j+1<m ? dist_AB(A[i], B[perm_B[j+1]]) : Infinity;
                    jm = j;
                    if (sm < s)
                    {
                        s = sm;
                        jm = j-1;
                    }
                    if (sM < s)
                    {
                        s = sM;
                        jm = j+1;
                    }
                    alignment[i] = perm_B[jm];
                }
            }
        }
        else if (n < m)
        {
            sm = Infinity; km = 0;
            for (k=0; k+n<=m; ++k)
            {
                s = 0;
                for (i=0; i<n; ++i) s += dist_AB(A[perm_A[i]], B[perm_B[i+k]]);
                if (s < sm) {sm = s; km = k;}
            }
            for (i=0; i<n; ++i)
            {
                alignment[perm_A[i]] = perm_B[i+km];
            }
        }
        else// if (n === m)
        {
            for (i=0; i<n; ++i)
            {
                alignment[perm_A[i]] = perm_B[i];
            }
        }
        return alignment;
    }
    return [];
}

function test(a, b, dist)
{
    console.log('a:'+a.join(',')+' b:'+b.join(',')+' alignment:'+align(a, b, dist).join(','));
}

test([0, 1, 2], [0, 1, 2], dist);
test([0, 1, 2], [2, 1, 0], dist);
test([2, 1, 0], [0, 1, 2], dist);
test([2, 1, 0], [2, 1, 0], dist);
test([0, 1, 2], [0, 1], dist);
test([0, 1, 2], [1, 0], dist);
test([2, 1, 0], [0, 1], dist);
test([2, 1, 0], [1, 0], dist);
test([2, 1, 0, 3, 4], [1, 0, 2], dist);
test([2, 3, 1, 4, 0], [1, 3, 2], dist);
test([-1, 2, 3, 1, 4, 0, 5], [1, 3, 2], dist);
test([3, 4, 2, 1, 0], [1, 2], dist);
test([2, 1, 0, 3, 4], [0, 2], dist);
test([0, 1, 2], [-2, -1, 0, 1, 2], dist);
test([0, 1, 2], [2, 1, 0, 4, 3], dist);
test([2, 1, 0], [-2, -1, 0, 1, 2], dist);
test([2, 1, 0], [2, 1, 0, 4, 3], dist);
test([2, 1, 0], [0, 2, 4, 1, 3], dist);

Monck