complex comparison ops vs. Matlab

[Jaroslav Hajek-2]

hi all,

by this changeset:
http://hg.savannah.gnu.org/hgweb/octave/rev/7dafdb8b062f

I finished my project of refactoring dense binary operator
implementations in liboctave. The approach is similar to what I've
done before for reduction operations:
The innermost loops on pointer arrays are encapsulated as overloaded
template functions, providing good performance; appliers are provided
that take array arguments and return array, applying a given loop.
Using this machinery, the (MS|SM|MM|NDS|SND|NDND)_(BIN|CMP|BOOL)_OP(S)
macros in mx-op-defs were changed into trivial one-line wrappers, and
the *OPS1, *OPS2 alternatives are gone.

The total of mx-op-defs.h, mx-inlines.cc and MArray-defs.h is now 200
lines shorter than 3.2.x, and that also includes the new stuff for
in-place binary ops. It will be even more if also the VS_,SV_,and
VV_OPS are replaced (which was not immediately possible because
*Vector classes are not constructible from dim_vector).

The new implementation also probably creates less object code; as the
loop bodies are passed by pointer and thus can be shared amongst the
matrix and ndarray wrappers. I'm not sure whether the GNU linker
actually merges the instances, but it certainly has the possibility
(the last change cut down liboctave size by 330KB. It's still 400KB
larger than 3.2.x, but that's not fair comparison because there's a
lot of new stuff since 3.2.x).

The logical binary ops were sped up significantly by this change,
although the code is logically equal - probably GCC was not powerful
enough to optimize the previous code:

n = 4e3;
a = rand (n) < 0.5;
b = rand (n) < 0.5;
tic; a & b; toc
tic; a | b; toc
tic; a & !b; toc
tic; a | !b; toc
tic; !a | b; toc
tic; !a & b; toc

GCC 4.3.1, -O3 -march=native, Core 2 Duo @ 2.83 GHz:

3.2.x:

Elapsed time is 0.0800021 seconds.
Elapsed time is 0.0852289 seconds.
Elapsed time is 0.080343 seconds.
Elapsed time is 0.0791581 seconds.
Elapsed time is 0.0736232 seconds.
Elapsed time is 0.074126 seconds.

tip:

Elapsed time is 0.018801 seconds.
Elapsed time is 0.0197818 seconds.
Elapsed time is 0.026242 seconds.
Elapsed time is 0.020916 seconds.
Elapsed time is 0.0185721 seconds.
Elapsed time is 0.0178182 seconds.

One important point is that I changed the way comparing complex number works.
I think that Matlab compares only real parts of complex arrays when
involved in comparison operators; at least the previous Octave
implementation did so.
OTOH, Matlab defines a rigorous ordering on complex numbers as the
lexicographical ordering of [abs(z), arg(z)], which is used in sort,
min and max, and is shared with Octave.
I couldn't find the Matlab behavior documented anywhere, so I presume
it's a relic of the past and probably stemming from the fact that
Matlab is widely known to store the real and imaginary parts of arrays
separately - so it's the lazy do-nothing approach. However, it
simplifies the implementation significantly to consider one universal
ordering of complex numbers (defined in oct-cmplx.h) that gets
automatically used everywhere. As a result, Octave now uses the strict
abs/arg ordering of complex numbers for comparing complex arrays:

octave:1> 1+5i > 2+2i
ans =  1
octave:2> 1+i > 1-i
ans =  1

This means that various "obvious" invariants, such as

octave:3> a = sort (rand(100, 1) + i*rand(100,1)); all (a(1:99) <= a(2:100))
ans =  1
octave:4> a = rand(100, 1) + i*rand(100,1); all (a <= max (a(:)))
ans =  1

now hold for complex arrays as well (previously, this was not true).
Personally, I think this new behavior is more logical and consistent
and I vote for a change. I suppose few code depends on this and broken
scripts can be trivially fixed by wrapping operands in real () (which
I'd suggest anyway).

If I'm outvoted on this, I'd still like to keep that Matlab
compatibility mess out of liboctave, so I'll probably add a hack into
the interpreter to apply "real" to the operands first. It could also
be made conditional (for instance, keep the old behavior in the
braindamage mode), but I think John already said he doesn't like
global variables fundamentally influencing code flow.

Opinions? Comments?

regards

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz

[John W. Eaton-3]

On 27-Aug-2009, Jaroslav Hajek wrote:

| by this changeset:
| http://hg.savannah.gnu.org/hgweb/octave/rev/7dafdb8b062f
|
| I finished my project of refactoring dense binary operator
| implementations in liboctave.

Thanks for taking on this project.

| One important point is that I changed the way comparing complex number works.
| I think that Matlab compares only real parts of complex arrays when
| involved in comparison operators; at least the previous Octave
| implementation did so.

Yes.  Not that it makes much sense (at least to me) but Matlab is
documented to work this way.

| OTOH, Matlab defines a rigorous ordering on complex numbers as the
| lexicographical ordering of [abs(z), arg(z)], which is used in sort,
| min and max, and is shared with Octave.
| I couldn't find the Matlab behavior documented anywhere,

Look at the alphabetical list of the documentation for Matlab
functions on the web.  Near the top of the page for relational
operators, it says

  The operators <, >, <=, and >= use only the real part of their
  operands for the comparison. The operators == and ~= test real and
  imaginary parts.

| However, it
| simplifies the implementation significantly to consider one universal
| ordering of complex numbers (defined in oct-cmplx.h) that gets
| automatically used everywhere.

Yes, I think it makes much more sense to have one definition for the
ordering.  Can anyone think of a case where it is really useful to
compare complex numbers using only the real part?good reason, then

| now hold for complex arrays as well (previously, this was not true).
| Personally, I think this new behavior is more logical and consistent
| and I vote for a change. I suppose few code depends on this and broken
| scripts can be trivially fixed by wrapping operands in real () (which
| I'd suggest anyway).

Yes, I agree that if you want to compare complex numbers by real part
only, then you should say so explicitly.

| If I'm outvoted on this, I'd still like to keep that Matlab
| compatibility mess out of liboctave, so I'll probably add a hack into
| the interpreter to apply "real" to the operands first. It could also
| be made conditional (for instance, keep the old behavior in the
| braindamage mode), but I think John already said he doesn't like
| global variables fundamentally influencing code flow.

What about an optional "matlab-incompatibility" warning?  Especially
since we are changing the behavior of Octave too, and this is the kind
of change that will otherwise silently change the results obtained
from of previously "working" code.

jwe

[Jaroslav Hajek-2]

On Thu, Aug 27, 2009 at 10:02 PM, John W. Eaton<jwe@...> wrote:

> On 27-Aug-2009, Jaroslav Hajek wrote:
>
> | by this changeset:
> | http://hg.savannah.gnu.org/hgweb/octave/rev/7dafdb8b062f
> |
> | I finished my project of refactoring dense binary operator
> | implementations in liboctave.
>
> Thanks for taking on this project.
>
> | One important point is that I changed the way comparing complex number works.
> | I think that Matlab compares only real parts of complex arrays when
> | involved in comparison operators; at least the previous Octave
> | implementation did so.
>
> Yes.  Not that it makes much sense (at least to me) but Matlab is
> documented to work this way.
>
> | OTOH, Matlab defines a rigorous ordering on complex numbers as the
> | lexicographical ordering of [abs(z), arg(z)], which is used in sort,
> | min and max, and is shared with Octave.
> | I couldn't find the Matlab behavior documented anywhere,
>
> Look at the alphabetical list of the documentation for Matlab
> functions on the web.  Near the top of the page for relational
> operators, it says
>
>  The operators <, >, <=, and >= use only the real part of their
>  operands for the comparison. The operators == and ~= test real and
>  imaginary parts.
>
> | However, it
> | simplifies the implementation significantly to consider one universal
> | ordering of complex numbers (defined in oct-cmplx.h) that gets
> | automatically used everywhere.
>
> Yes, I think it makes much more sense to have one definition for the
> ordering.  Can anyone think of a case where it is really useful to
> compare complex numbers using only the real part?good reason, then
>

Right now I can't think of any (which doesn't mean there isn't any),
but even if it is, it can be done with real (), which will also be
more readable. On the contrary, getting the complex ordering is not
trivial. For instance, how would you check that a complex array is
strictly upwards sorted in Matlab? Compare
za = abs (z); zf = arg (z); all (za(1:n-1) < za(2:n) || (za(1:n-1) ==
za(2:n) && zf(1:n-1) < zf(2:n))
to the simple
z(1:n-1) < z(2:n)

But of course, it's primarily the inconsistency with max/min/sort that
makes the Matlab behavior nutty. Also, the complex ordering defined by
sort is consistent with complex equality (one of <, >, == is always
true, except for NaNs) but comparison by real parts isn't. It' also
consistent with the normal ordering of reals.

>
> | now hold for complex arrays as well (previously, this was not true).
> | Personally, I think this new behavior is more logical and consistent
> | and I vote for a change. I suppose few code depends on this and broken
> | scripts can be trivially fixed by wrapping operands in real () (which
> | I'd suggest anyway).
>
> Yes, I agree that if you want to compare complex numbers by real part
> only, then you should say so explicitly.
>
> | If I'm outvoted on this, I'd still like to keep that Matlab
> | compatibility mess out of liboctave, so I'll probably add a hack into
> | the interpreter to apply "real" to the operands first. It could also
> | be made conditional (for instance, keep the old behavior in the
> | braindamage mode), but I think John already said he doesn't like
> | global variables fundamentally influencing code flow.
>
> What about an optional "matlab-incompatibility" warning?  Especially
> since we are changing the behavior of Octave too, and this is the kind
> of change that will otherwise silently change the results obtained
> from of previously "working" code.
>
> jwe
>

OK, that would probably work best. And the warning on in traditional
mode, I suppose. I'll try making a patch.

regards

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz