Match Two Vectors
match_vectors.RdGiven two vectors x,y of length \(p \leq q\) respectively,
the routine match_vectors returns an integer vector
j, with unique elements, such that x matches y[j]
as best as possible. The procedure uses the "Munkres" algorithm for
solving this assignment problem. The procedure throws an error if the length of
x is larger than the length of y.
Usage
match_vectors(x, y = Conj(x))Examples
# Match the roots of two polynomials a1 and a2
p = 5
a1 = rnorm(p+1)
a2 = a1 + rnorm(p+1)*(1e-6) # a2 is a "noisy" copy of a1
z1 = polyroot(a1)
z2 = polyroot(a2)[order(stats::rnorm(p))] # reshuffle the roots of a2
j = match_vectors(z1, z2)
print(data.frame(z1 = z1, j = j, `z2[j]` = z2[j], d = z1-z2[j]))
#> z1 j z2.j. d
#> 1 -0.1895884+0.000000i 5 -0.1895887+0.000000i 2.975242e-07-0.000000e+00i
#> 2 -0.2888507+1.053496i 3 -0.2888506+1.053495i -9.844998e-08+3.192245e-07i
#> 3 -0.2888507-1.053496i 4 -0.2888506-1.053495i -9.844998e-08-3.192245e-07i
#> 4 1.2201385-0.998451i 2 1.2201376-0.998451i 8.943785e-07+1.775183e-08i
#> 5 1.2201385+0.998451i 1 1.2201376+0.998451i 8.943785e-07-1.775183e-08i
# A polynomial with real coefficients has pairs of complex conjugate roots.
# However, the roots returned by "polyroot" in general do not have this
# property!
# Match the roots and their complex conjugates
j = match_vectors(z1, Conj(z1))
print(data.frame(z = z1, j = j, `Conj(z[j])` = Conj(z1[j]),
d = z1-Conj(z1[j])))
#> z j Conj.z.j.. d
#> 1 -0.1895884+0.000000i 1 -0.1895884-0.000000i 0.000000e+00+6.955665e-16i
#> 2 -0.2888507+1.053496i 3 -0.2888507+1.053496i 2.775558e-16+0.000000e+00i
#> 3 -0.2888507-1.053496i 2 -0.2888507-1.053496i -2.775558e-16+0.000000e+00i
#> 4 1.2201385-0.998451i 5 1.2201385-0.998451i 0.000000e+00-3.330669e-16i
#> 5 1.2201385+0.998451i 4 1.2201385+0.998451i 0.000000e+00-3.330669e-16i