Matlab: optimal-distance matching of the elements of two sets
I have two floating-point number vectors which contain the same values up to a small error, but not necessarily sor开发者_开发技巧ted in the same way; for instance, A=rand(10);a=eig(A);b=eig(A+1e-10);
(remember that eig
outputs eigenvalues in no specified order).
I need to find a permutation p
that matches the corresponding elements, i.e. p=mysterious_function(a,b)
such that norm(a-b(p))
is small.
Is there an existing function that does this in a sane and safe way, or do I really need to roll out my own slow and poorly-error-checked implementation?
I need this only for test purposes for now, it need not be excessively optimized. Notice that the solution which involves sorting both vectors with sort
fails in case of vectors containing complex equal-modulus arguments, such as the typical output of eig()
.
You seem to want to solve the linear assignment problem. I haven't tested it myself, but this piece of code should help you.
I believe that the sort()
solution you discarded might actually work for you; The criteria you have defined minimize norm(a-b)
is, by definition, considering the modulus (absolute value) of the complex number: norm(a-b) == sqrt(sum(abs(a-b).^2))
And as you know, SORT
orders complex numbers based on their absolute value: sort(a)
is equivalent to sort(abs(a))
for complex input.
%# sort by complex-magnitude
[sort(a) sort(b)]
As long as the same order is applied to both, you might as well try lexicographic ordering (sort by real part, if equal, then sort by imaginary part):
%# lexicographic sort order
[~,ordA] = sortrows([real(a) imag(a)],[1 2]);
[~,ordB] = sortrows([real(b) imag(b)],[1 2]);
[b(ordB) a(ordA)]
If you are too lazy to implement the Hungarian algorithm that @AnthonyLabarre suggested, go for brute-forcing:
A = rand(5);
a = eig(A);
b = eig(A+1e-10);
bb = perms(b); %# all permutations of b
nrm = sqrt( sum(abs(bsxfun(@minus, a,bb')).^2) ); %'
[~,idx] = min(nrm); %# argmin norm(a-bb(i,:))
[bb(idx,:)' a]
Beside the fact that eigenvalues returned by EIG are not guaranteed to be sorted, there is another difficulty you have to deal with if you to match eigenvectors as well: they are not unique in the sense that if v
is an eigenvector, then k*v
is also one, especially for k=-1
. Usually you would enforce a sign convention like: multiply by -/+1 so that the largest element in each vector have a positive sign.
精彩评论