bugGNU Octave - Bugs: bug #52919, comparison of complex values

 
 

bug #52919: comparison of complex values

Submitted by:  Michael Leitner <mleitner>
Submitted on:  Wed 17 Jan 2018 11:06:17 AM UTC  
 
Category:  Interpreter Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  None Assigned to:  None
Originator Name:  mleitner Open/Closed:  Open
Release:  dev Operating System:  GNU/Linux

Add a New Comment (Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

( Jump to the original submission)

Tue 30 Jan 2018 04:02:39 PM UTC, comment #15: 

Just for information, R doc says:
Complex values are sorted first by the real part, then the imaginary part.

Looks like they have not given it much thought.

Michael Godfrey <godfrey>
Project Member
Mon 22 Jan 2018 05:41:49 PM UTC, comment #14: 

I think we all agree that sort and sortrows have performance issues that should be addressed.  That is really its own topic and should be segregated into a new bug report.  Could someone create that and port the comments over from this report?

As for Matlab's behavior, it really was quite bad.  As far as I know, the initial implementation was not a stable sort.  They added an extra option to commands like unique() so you could specify whether you wanted a "stable" or "sorted" output (http://www.mathworks.com/help/matlab/ref/unique.html).

Real numbers are not a special case of complex numbers in Matlab.  They are a different class.  The functions isreal, iscomplex, and isinteger all test the class that is used to represent the number, rather than testing the actual value of the number.  For example

x = complex (1,0)
x =  1 + 0i
octave:9> isreal (x)
ans = 0
octave:10> iscomplex (x)
ans = 1
octave:11> isinteger (x)
ans = 0

The number 1.0 is clearly an integer, but Octave says it isn't because it is being stored using 8 bytes and an IEEE double format.

The issue that you ran in to with indexing of a complex array producing a real result is known as narrowing.  Because complex numbers take twice as much memory to store, Octave and Matlab attempt to reduce complex arrays to real arrays whenever they can.  For example, you cannot creat a complex scalar with a zero magnitude imaginary component on the command line.  Try

octave:12> x = 1+0i
x =  1

In this case Octave's parser notices that you are trying to create a complex number, but the interpreter realizes that it can narrow this complex number to a real number and it does so.  The way around that is to use the complex() function which avoids automatic narrowing.

octave:13> x = complex (1,0)
x =  1 + 0i

Even so, Octave will narrow this to a real result whenever it has a chance.  An operation like indexing, or addition, will do this.  For example,

octave:13> x = complex (1,0)
x =  1 + 0i
octave:14> x(1)
ans =  1
octave:15> x
x =  1 + 0i
octave:16> x + 0
ans =  1

Matlab has some documentation to help out on this: https://www.mathworks.com/help/coder/ug/code-generation-for-complex-data.html.

Note:

For code generation, complex data that has all zero-valued imaginary parts remains complex. This data does not become real. This behavior has the following implications:

    In some cases, results from functions that sort complex data by absolute value can differ from the MATLAB® results. See Functions That Sort Complex Values by Absolute Value.

    For functions that require that complex inputs are sorted by absolute value, complex inputs with zero-valued imaginary parts must be sorted by absolute value. These functions include ismember, union, intersect, setdiff, and setxor.

Because Octave was second, we have to follow Matlab's rules which are that one sort algorithm is used for real numbers (value) and another algorithm for complex numbers (absolute value).

The extra complication around complex numbers on the negative real axis is not a buglet.  Matlab is (now) breaking ties in absolute value by using the phase angle which they get by calling the underlying atan2 library function.  This function returns an angle in the range [-pi, pi] (both endpoints are included).  But this causes problems because it can assign two values to the same underlying real number.  For example,

octave:21> x = [3-4i, complex(-5,0), 3+4i]
x =

   3 - 4i  -5 + 0i   3 + 4i

octave:22> abs (x)
ans =

   5   5   5

octave:23> angle (x)
ans =

  -0.92730   3.14159   0.92730

All numbers have the same magnitude, so the tie is broken by phase angle resulting in

[3-4i, 3+4i, -5]

But the IEEE standard for floating points includes both +0 and -0.  What happens when the imaginary part is -0

x(2) = complex (-5,-0)
x =

   3 - 4i  -5 - 0i   3 + 4i

octave:25> abs (x)
ans =

   5   5   5

octave:26> angle (x)
ans =

  -0.92730  -3.14159   0.92730

Now the sort order is

[-5, 3-4i, 3+4i]

That bothers most people, so Octave takes the additional step of remapping the atan2 function on to the range (-pi, pi] (first endpoint is not included).

The question to me is whether, given that we have to use abs() as the first test, it would be useful to prefer lexigraphic order during ties.  The current sort order puts negative real values after positive real values which is consistent, but doesn't generally accord with people's natural sorting algorithm.  For example,

sort ([complex(-5,0), complex(5,0), 3+4i, 3-4i])
ans =

   3 - 4i   5 + 0i   3 + 4i  -5 + 0i

We could have -5+0i appear first which might make more sense.  There are potentially other benefits as well.  It may be faster to simply compare magnitudes of real and imaginary parts rather than calculating the atan2 function.  Finally, there is an issue on Mac OSX platforms that makes the remapping of the range [-pi,pi] to (-pi, pi] cumbersome because they use a different value of pi than anyone else.  If the code were just comparing magnitudes it should be platform independent.

It appears that there could be some benefits to re-defining the sort algorithm, but first the Octave Maintainer's should decide whether it is worth deviating from Matlab compatibility in this regard.

Rik <rik5>
Project Administrator
Mon 22 Jan 2018 09:40:21 AM UTC, comment #13: 

To provisionally sum up the different issues that have crept up in this thread:

- the performance of sort() on complex values (compiled C++ is significantly slower than trivial m-code solution)
- should the order on complex values be changed (either for principal reasons, for least surprise, or for compatibility with Matlab, and if so, also for sort/min/max)?
- should perhaps a warning be emitted, as comparing complex values is arguably often unintentional?

And now a few thoughts on the comments over the weekend:

Re comment #12:
Yes, it is indeed slow. I think the difference is even larger than a factor 2.5. Consider this:

N=1e7;
r=floor(randn(N,1)+1i*randn(N,1))+.5+.5*1i;

tic;[~,j]=sortrows([abs(r) angle(r)]);rs1=r(j);toc
tic;rs2=sort(r);toc
isequal(rs1,rs2)

tic;[~,j]=unique([abs(r) angle(r)],"rows");u1=r(j);toc
tic;u2=unique(r);toc
isequal(u1,u2)

On my system, the difference is nearly a factor of four, and sortrows sorts two times, not just once according to abs(r). This is just the code that is also in sort_complex_test.m added in comment #6, using the fact that this is actually the algorithm used by sortrows. With the rounding, also some absolute values are equal, so that the decision on arg is relevant. Note that adding 0.5 to the numbers is necessary, as by conventional rounding you can end up with -1-0i and -1+0i, which seemingly are treated differently (which probably is also a small buglet).

In C++, one could definitely go below the time taken by the sortrows-solution by first sorting according to abs, and then make another pass and sort all stretches with same abs by arg. This approach should take N* log(N)+N* log(N* p), where p is the probability for two values having the same absolute value, and which is potentially very small. On the other hand, sortrows takes 2N* log(N) (disregarding the fact that also sort itself probably takes shortcuts when many arguments are equal). But I think, the most sensible solution would be just to change sort to do it in the sortrows-way for complex values.

I support comment #11, with the proviso that I am not quite convinced that it is a good idea to sort complex values (i.e., to rely on an order relation) in any case. But yes, if there is an order relation, it should be consistent.

Re comment #10:
As regards the Matlab behaviour: probably it was not that bad, because I would hope that the used sort algorithm is stable, that is, it keeps the input order for values that compare as equal. According to the documentation, Octave uses stable sort, which is also a prerequisite for sortrows.

As regards whether a real number is a special case of a complex number or not: I do not know what Matlab does, but for Octave I would argue that it is. Consider this:

a=[-1 1i];
iscomplex(a)
iscomplex(a(1))
iscomplex(a(2))

So without having told Octave to cast a(1) to real, it considers it to be real. And of course, it uses the real ordering relations on it, which gives the inconsistency that was my initial motivation for filing this bug report, because the following two are different:

(a<1)(1)
a(1)<1

So if you want to hold up this concept of "Octave has two categories of numbers", then you would need to have b=a(1);iscomplex(b) in the present example evaluate to true. Personally, I would not do this.

So, now for the question which ordering to use on the complex numbers: For the present (abs,arg)-sorting you can even add a real number to a complex number to reverse the order, compare the results of (a<1) and (a+1<1+1). When you do lexicographic (real,imag)-sorting, then this does not happen. Thus, with lexicographic sorting you could consider the real numbers as special case of complex numbers not only with respect to arithmetics, but also ordering, and these issues that real-valued entries taken out of a complex vector compare differently than the element-wise operation on the whole vector would not appear.

The pseudocode in the comment actually just implements lexicographic ordering, because if neither areal!=breal and aimag!=bimag, then of course a==b. I think there is no way to "have an ordering that would be consistent with the real numbers as the imaginary part of the number was shrunk to zero" but still sort primarily on the absolute value, because this would always make -2+1i* eps greater than 1+1i* eps.

Re comment #9:
Actually, when you want to efficiently sort complex numbers, you should rather do it not in C++ (at least not in the present C++), as I have shown. I have no idea at all what rlocus does and what typical input arguments are, but looking into its code I would be surprised if the performance of the internal sort() is actually deciding here. Are you sure that it is?

Michael Leitner <mleitner>
Fri 19 Jan 2018 10:53:11 PM UTC, comment #12: 

Irrespective of the sorting order, it does seem like sorting of complex numbers is rather slow.

I used the following test code which is also attached

N = 1e3;

x = [1e6:-1:1]';
y = x + 1i;
z = zeros (size (x));

bm1 = bm2 = bm3 = zeros (N,1);
for i = 1:N
  tic;
  z = sort (x);
  bm1(i) = toc;
  tic;
  z = sort (y);
  bm2(i) = toc;
  tic;
  z = sort (abs (y));
  bm3(i) = toc;
endfor

mean (bm1)
mean (bm2)
mean (bm3)

The results were

ans =  0.0028022
ans =  0.025930
ans =  0.011031

So sorting of complex numbers is about 10X slower than reals.  And more importantly, it is abot 2.5X slower than the equivalent sort (abs (y)) command.  The y values chosen all differ in abs() so it is not a question of extra processing of the phase angle.  It seems there are just a lot of hoops to go through before the correct code is engaged.

(file #42980)

Rik <rik5>
Project Administrator
Fri 19 Jan 2018 08:29:42 PM UTC, comment #11: 

Just adding, I think it's also important for max, min, sort, as well as eq, ge, gt, le, and lt to all be consistent. So that if someone writes their own simple bubble sort function using only lt, the results would be the same as the built in sort function. This is not true in Matlab.

Mike Miller <mtmiller>
Project Administrator
Fri 19 Jan 2018 08:23:49 PM UTC, comment #10: 

This is not an easy thing to address.  It should probably be discussed back on the Octave Maintainer's list if there is going to be a major change to ordering.

First, I agree with principle of least surprise.  To begin with that means operators and max, min, sort should all behave the same way.  Octave is not quite there yet.  The operators and sort are self-consistent.  And max/min are consistent when called with one argument.  When called with a two argument form such as max (x, 1+1i) the ordering is not consistent.  But we already have a bug report for that.  See bug #51307.

One thing that we are working against is Matlab's definition of complex.  iscomplex() does not mean the imaginary part is non-zero, it merely indicates that Matlab has reserved storage for a complex portion of the number.  I know that it may be tempting to think of real numbers as the limiting case of complex numbers as the imaginary part goes to zero.

xreal =   limit  (x + epsilon*i)
        epsilion -> 0

but this isn't what Matlab, and hence Octave, have implemented.  There is always a categorical distinction between a real number and a complex number.

So Octave has two categories of numbers and applies a different ordering system for each number type.  Incidentally, I believe Matlab used to sort real numbers by value and complex numbers by abs() without checking the phase angle for ties.  This is a horrible idea because the same program can produce different results based on the input data.  In the case of sort, equal magnitude numbers might appear in the output in the order they were in the input, or frankly any other permutation.  Not only does this create programs which fail based on differences in input data, but given that the output could depend on the internals of the sort algorithm, it also introduced a dependence on which version of Matlab you were using.  At least Matlab now uses all information (magnitude and phase angle) to unambiguously order things.  According to the documentation for sort

If A is complex, then by default, sort sorts the elements by magnitude. If more than one element has equal magnitude, then the elements are sorted by phase angle on the interval (−π, π].

This already seems to be what Octave does.

x = [-1+1i, 0, 1+1i]
x =

  -1 + 1i   0 + 0i   1 + 1i

angle (x)
ans =

   2.35619   0.00000   0.78540

sort (x)
ans =

   0 + 0i   1 + 1i  -1 + 1i

Hooray, Octave is Matlab-compatible as the number -1+1i with phase angle 2.356 is placed after 1+1i which only has a phase angle of 0.7854.

Octave could just stop there since we are compatible, but we do try to be better than Matlab where we can.

I agree that it is somewhat galling to see 1+1i ahead of -1+1i.  Intuitively, it feels like the number with a negative real part should be smaller.  In a perfect world, -1 + 1i is less than 1+1i.  In addition, I understand that it is counter-intuitive that adding a complex number to an existing complex number can reverse the sort order.

x = -1 + 1i;
y = 1+1i;
z = -2i;
x < y
ans = 0
(x+z) < (y+z)
ans = 1

One possibility would be to do a slightly more sophisticated ordering on the phase angle rather than just using the value.  If we first resolved whether the angles were in the left half of the complex plane (corresponding to negative real numbers) versus the right (corresponding to positive real numbers) then at least we would have an ordering that would be consistent with the real numbers as the imaginary part of the number was shrunk to zero.  For the test example above, the result of the '<' operator would also not flip based on the addition of an imaginary number.  That would be another nice gain.

Here is a sample of 4 relevant numbers with equal magnitude, but different phase angles.

x = [-1-1i, -1+1i, 1-1i, 1+1i]
x =

  -1 - 1i  -1 + 1i   1 - 1i   1 + 1i

octave:51> angle (x)
ans =

  -2.35619   2.35619  -0.78540   0.78540

That seems very distinguishable.  There may be an efficient way to calculate this without even calling the angle() function and relying on the signbits and magnitudes of the real and imaginary parts directly.

For reference, the code that implements operators is in liboctave/util/oct-cmplx.h

     template <typename T>                                                   \
     inline bool operator OP (const std::complex<T>& a,                      \
                              const std::complex<T>& b)                      \
     {                                                                       \
       OCTAVE_FLOAT_TRUNCATE const T ax = std::abs (a);                      \
       OCTAVE_FLOAT_TRUNCATE const T bx = std::abs (b);                      \
       if (ax == bx)                                                         \
         {                                                                   \
           OCTAVE_FLOAT_TRUNCATE const T ay = std::arg (a);                  \
           OCTAVE_FLOAT_TRUNCATE const T by = std::arg (b);                  \
           if (ay == static_cast<T> (-M_PI))                                 \
             {                                                               \
               if (by != static_cast<T> (-M_PI))                             \
                 return static_cast<T> (M_PI) OP by;                         \
             }                                                               \
           else if (by == static_cast<T> (-M_PI))                            \
             {                                                               \
               return ay OP static_cast<T> (M_PI);                           \
             }                                                               \
           return ay OP by;                                                  \
         }                                                                   \
       else                                                                  \
         return ax OPS bx;                                                   \
     }                                                                       \

I don't know how expensive the std::arg call is, but Octave is already doing a bit of work to re-map -PI to +PI that we might be able to do away with.

In pseduocode, after determining that the magnitudes are equal.

  areal = real (a);
  breal = real (b);
  if (areal != breal)
    return areal OP breal;   # could also just use signbit for comparison

  aimag = imag (a);
  bimag = imag (b);
  if (aimag != bimag)
    return aimag OP bimag;   # could also just use signbit for comparison

  ## In the same quadrant
  return (aimag / areal) OP (bimag / breal);

Rik <rik5>
Project Administrator
Thu 18 Jan 2018 04:56:44 PM UTC, comment #9: 

If it is going to change then it should be done in C++ for speed.
and it should not break rlocus.

I am concerned about the speed of rlocus already so any slowdown would be added slowness.

Doug Stewart <dastew>
Thu 18 Jan 2018 04:41:07 PM UTC, comment #8: 

It appears that Michael's comments #5 and #7 are a good
choice. Is there any objection?

Michael Godfrey <godfrey>
Project Member
Thu 18 Jan 2018 01:15:12 PM UTC, comment #7: 

Re comment #6:
I showed how the present behaviour of sort can be obtained when only comparisons of real values are used, and I admitted that this is less efficient than having the comparison function do it directly. And I doubted that this decreased efficiency would constitute a problem, which probably is also valid for your use case: most likely you sort the points that you plot afterwards, and even sorting two times will take much less time per point than plotting them.

Actually: in my system a sort of 1e7 real numbers takes 1.26 seconds, and a sort of 1e7 complex numbers takes 7.3 seconds. So it seems that the option to sort complex numbers is an afterthought that (nearly) nobody is using anyway, otherwise somebody would probably have come up with the idea to write a wrapper around sort for complex values to regain this factor of three by my trivial algorithm, which I now suggest. By the way, I think that for unique this is much more relevant.

For testing use the attached script. So it seems to me that my caveat of efficiency actually goes the other way. And it would not be any problem at all to implement Matlab's extended syntax for sort in this way.

(file #42955)

Michael Leitner <mleitner>
Thu 18 Jan 2018 11:44:09 AM UTC, comment #6: 

Michael Leitner said:
but when there is no use-case, who cares?

I believe there is a use-case, and I care.
If I remember correctly It was an issue in rlocus.
If it wasn't sorted this way then the plotting of the rlocus would switch paths.

Doug Stewart <dastew>
Thu 18 Jan 2018 08:23:40 AM UTC, comment #5: 

Yes, I definitely subscribe to the order relations being consistent with sort, min and max. As has been said, Matlab is (by default) not consistent in these three functions, and only for sort, the behaviour can be specified (they do not say anything about median and the like -- of course, Octave is again consistent).

However: I think, all of us are deeply convinced that real numbers are just a special case of complex numbers. All arithmetic operations when performed on complex numbers with imaginary part equal to zero reduce to those on the real numbers, and also all typical analytical functions like exp, sqrt, log and so on do. This is also the reason why Octave without being asked gives a real vector D

D=C(i)

for complex C, when all C(i) have imaginary part zero. However, only the ordering relations as Octave use it are inconsistent between C and D (that is, for very different x and y, when you add 1i*eps to x, they can order differently than without). Matlab's ordering relations would be consistent. That's why I like them better, and if there were no questions of backward compatibility, I would argue for changing it (whatever the compatibility with Matlab).

Further however: I think that in the vast majority of cases when complex numbers are compared (perhaps with the exception of the function unique), this is done erroneously. Apart from this, I see no use-case for ordering complex numbers. And when I would like to sort as Octave's sort does it now, I would do

A=complex(randn(10,1),randn(10,1))
B=sort(A)
[~,j]=sort(arg(A))
A=A(j)
[~,j]=sort(abs(A))
B1=A(j)
all(isequal(B,B1))

and analogously for Matlab's 'ComparisonMethod','real' (yes, vastly less efficient, but when there is no use-case, who cares?). In the case where I stumbled over this behaviour, I actually was aware that with small probability some entries in the two vectors I compared could be complex. These entries would have been sifted out afterwards, but in fact these few entries made the whole vector complex, thus the ordering of the whole vectors (for the negative values) was different.

Another example: do you think that x<y implies x+10<y+10? Yes, I do also, and I would not be surprised when this is assumed deep inside some internal Octave functions, but for Octave's complex comparison this is not the case. So as I think that comparing complex numbers for order (of course not for equality) is done in most cases unintendedly, I propose to emit a warning (the first time). I am someone who does read error messages, and this would really have been useful for me.

Michael Leitner <mleitner>
Wed 17 Jan 2018 08:21:40 PM UTC, comment #4: 

Oh, I see Matlab's sort and issorted functions now take a 'ComparisonMethod' option, so you can choose to sort by abs or by real. Default behavior stays the way it has been, inconsistent with eq, gt, lt. I'd rather Octave be consistent.

Mike Miller <mtmiller>
Project Administrator
Wed 17 Jan 2018 08:10:55 PM UTC, comment #3: 

Yeah, for operators like <, <=, etc., we use the same comparison for complex numbers as Matlab does for sort, min, and max.  That makes sense to me.  It's bizarre that Matlab would use a different sort ordering for those functions than it does for the operators.

But I guess consistency is not as important as copying the exact behavior of Matlab.  (Somewhat OT:  And now that Matlab has added double-quoted strings that are incompatible with Octave, people are thinking that we can keep Octave's behavior as it is now and do something different from Matlab?)

John W. Eaton <jwe>
Project Administrator
Wed 17 Jan 2018 07:54:08 PM UTC, comment #2: 

I believe the difference is Octave compares magnitude first, then phase, while Matlab compares only the real part (first? at all? not sure).

Here's a message from Jordi in 2013 confirming that Matlab only compares the real parts and we know Octave is different and we like our behavior better: https://lists.gnu.org/archive/html/octave-maintainers/2013-11/msg00281.html

Mike Miller <mtmiller>
Project Administrator
Wed 17 Jan 2018 07:00:11 PM UTC, comment #1: 

I thought we did try to match Matlab's behavior.  Is that not correct, or has the Matlab behavior changed?

John W. Eaton <jwe>
Project Administrator
Wed 17 Jan 2018 11:06:17 AM UTC, original submission:  

According to the documentation:
For complex numbers, the following ordering is defined: z1 < z2 if and only if

  abs (z1) < abs (z2)
  || (abs (z1) == abs (z2) && arg (z1) < arg (z2))

However, consider the following situation:

a=-rand(10,1);
b=-rand(10,1);
b(10)+=1i;

then

(a<b)(1)

gives a different result from

a(1)<b(1)

At least, this violates the principle of least surprise (and cost me a few hours tracking to another problem). Matlab's comparison of only the real part would of course solve this problem. I do not want to argue to make such a fundamental change.

Is there a viable solution? One could check in complex matrix comparisons whether the two numbers are actually real, but this would be less efficient and would also give inconsistencies within the matrix comparison.

Actually, I think that the Matlab way would be best, because then the ordering is transitive when real numbers are seen as special cases of complex numbers (I would rather redefine min, max and sort). But again, I see that this would be a very fundamental change.

Michael Leitner <mleitner>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #42980:  bm_sort.m added by rik5 (269B - text/x-matlab)
file #42955:  sort_complex_test.m added by mleitner (359B - text/x-matlab)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by amro_octave
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by dastew (Posted a comment)
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by mleitner (Submitted the item)
  • -email is unavailable- added by mleitner
  •  

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

     

     

    Follow 5 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-02-02 amro_octave Carbon-Copy- => Added amro_octave
    2018-01-22 rik5 Release4.0.3 => dev
    2018-01-19 rik5 Attached File- => Added bm_sort.m, #42980
    2018-01-18 mleitner Attached File- => Added sort_complex_test.m, #42955
    2018-01-17 mleitner Carbon-Copy- => Added mleitner

    Back to the top


    Powered by Savane 3.5