np.median on a complex array return incorrect value for the imaginary part

[flomertens]

Calling np.median on a complex array return a complex number with correct real part, but incorrect imaginary part. Calling np.median on the real and imaginary part of the complex array separately return the correct values.

Reproducing code example:

d = (1 + 1j) + np.random.randn(10000) + 1j * np.random.randn(10000)
print 'median(d) : %.3f + %.3f i' % (np.median(d).real, np.median(d).imag)
print 'median(d.real): %.3f' % np.median(d.real)
print 'median(d.imag): %.3f' % np.median(d.imag)

median(d) : 1.002 + 0.124 i
median(d.real): 1.002
median(d.imag): 0.998

Numpy version: 1.15.4

[charris]

The result for complex types depends on the sort order of complex, which will not give you what you want. I am not convinced that returning the medians of the real and imaginary parts as you desire makes much more sense as a default. For multidimensional points, the geometric median might be a better choice.

[seberg]

Hmm, frankly, it seems to me that median on complex is probably just broken (at least unless the complex numbers have some additional properties)? Should we deprecate it, or put a user warning?

[charris]

@seberg Deprecation seems reasonable. Should probably raise the issue on the list first.

[WarrenWeckesser]

I agree that, outside the context of NumPy, taking the median of a set of complex numbers is not well-defined. However, my understanding is that NumPy has imposed an ordering on the complex numbers by using lexicographic ordering of the real and imaginary parts. This affects all functions that are based on the ordering of the inputs, including the functions and methods min, argmin, minimum, fmin, max, argmax, maximum, fmax, sort, argsort, partition, argpartition, percentile, quantile, searchsorted, clip, all the comparison functions & ufuncs, and probably some others that I can't think of. So it seems consistent that median should also use the same ordering. In that context, the results given above and in an example such as:

In [80]: z
Out[80]: array([10.+0.j, 1.+2.j, -0.-3.j, -3.+3.j, 4.-5.j])

In [81]: np.sort(z)
Out[81]: array([-3.+3.j, -0.-3.j, 1.+2.j, 4.-5.j, 10.+0.j])

In [82]: np.median(z)
Out[82]: (1+2j)

are correct.

[WarrenWeckesser]

... NumPy has imposed an ordering on the complex numbers by using lexicographic ordering of the real and imaginary parts.

There is some interest in deprecating this: #15981