bugGNU Octave - Bugs: bug #61912, Slow performance of complex()

 
 

bug #61912: Slow performance of complex()

Submitter:  Rik <rik5>
Submitted:  Mon 24 Jan 2022 10:00:31 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Performance
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sun 30 Jan 2022 11:43:56 AM UTC, comment #13: 

Should we open a new report for the apparently large overhead for function calls? (At least, so we don't forget about it...)

Markus Mützel <mmuetzel>
Group administrator
Fri 28 Jan 2022 11:41:51 PM UTC, comment #12: 

I checked in some small improvements here: http://hg.savannah.gnu.org/hgweb/octave/rev/1182fcd6cdaf.

Re-benchmarking, I find that the complex() code takes 0.38 seconds.  But just calling a function body (with no code) took 0.32 seconds.  So the actual operations inside Fcomplex are fast but the function call and loop overhead seems to be quite high.

Marking as Fixed and closing report.

Rik <rik5>
Group administrator
Wed 26 Jan 2022 05:41:28 PM UTC, comment #11: 

I tried one more syntax hoping that avoiding the Complex constructor would save some time since the Complex NDArray already has zero-initialized Complex() values.  But, it made maybe a 2% difference so not really worth it.


              for (octave_idx_type i = 0; i < re_val.numel (); i++)
                {
                result.xelem (i).real (re_val.xelem (i));
                result.xelem (i).imag (im_val.xelem (i));
                //result.xelem (i) = Complex (re_val(i), im_val(i));
                }


Rik <rik5>
Group administrator
Wed 26 Jan 2022 04:13:52 PM UTC, comment #10: 

Responses to comments #7 and #8

@jwe: I agree, no need to use reinterpret_cast.  I was looking for any technique that might show an improvement.  I guess I'll prepare a changeset that adds in a comment about why the "new" operator is being used and make the other small changes to not initialize the array and to use xelem.

@hg200: I noticed the same thing, specifically that function calling seems to be quite slow.  As a test, I replaced the entire Fcomplex function in data.cc with


return octave_value (1.0);


I then used this test code


N = 1e5;

if (! exist ("val"))
  val = rand (1e2, 1);
endif

tic;
for i = 1 : N
#  z = complex (val, val);
  z = val + 1i*val;
endfor
bm = toc


Timings were

complex(): 0.32
multiplication: 0.23

This is remarkable given that in the multiplication case actual work was being done and it still beat the no-op version of complex().  That observation probably deserves its own bug report.

Rik <rik5>
Group administrator
Wed 26 Jan 2022 04:11:50 PM UTC, comment #9: 

@hg200: Interesting observation. That probably means that the function body already executes pretty fast. But the overhead from calling any(?) function is surprisingly high.

Markus Mützel <mmuetzel>
Group administrator
Wed 26 Jan 2022 06:11:04 AM UTC, comment #8: 

Maybe this is a too simple argument but: If i take the loop-code from the original posting its runtime is 0.190 over here. The fastest way to exit from data.cc is when only one argument is passed, that is complex:


  if (nargin == 1)
    {
      octave_value arg = args(0);
      if (arg.iscomplex ())
        retval = arg;
      else



N = 1e5;
val = 1i*pi;
tic;
for i = 1 : N
  z = complex (val);
endfor
bm = toc


With this change the runtime is 0.183. Not much faster.

Additionally you can exchange complex() with other functions. For example "isempty". Then the runtime is 0.168.

Hg200 <hg200>
Wed 26 Jan 2022 06:07:39 AM UTC, comment #7: 

The constructors like


octave_value (const ComplexNDArray&)


will transform the value from complex to real if the imaginary part is zero and we don't want that to happen in the complex function.  The constructor that takes a pointer to the rep object doesn't do that.  There should be a an explanatory comment for that because it wasn't obvious to me either at this point.

My guess is that at the time that the complex function was added, there was no way to create an uninitialized array.  If there is now, then yes, we should be using it if the next step will be to redefine every element.

Using xelem or even just pointers to the data seems fine to me, but I'd rather avoid things like using "reinterpret_cast<double*>".

John W. Eaton <jwe>
Group administrator
Tue 25 Jan 2022 10:27:56 PM UTC, comment #6: 

Adding jwe to the CC list as I'd like his opinion.

The code for Fcomplex is in data.cc.  It is unusual in that it uses the new operator from C++ whereas all but one other function in data.cc do not.  Why is that necessary?  It doesn't seem to be the cause of the performance issues (although one could imagine that a lot of new/delete operations would be slow).  For reference, I quote an m-file case of


val = rand (1e4, 1);
z = complex (val, val);


The code in data.cc begins at line 3683.


          const NDArray re_val = re.array_value ();

          if (im.numel () == 1)
            {
              double im_val = im.double_value ();

              ComplexNDArray result (re_val.dims (), Complex ());

              for (octave_idx_type i = 0; i < re_val.numel (); i++)
                result.xelem (i) = Complex (re_val(i), im_val);

              retval = octave_value (new octave_complex_matrix (result));
            }
          else
            {
              const NDArray im_val = im.array_value ();

              if (re_val.dims () != im_val.dims ())
                error ("complex: dimension mismatch");

              ComplexNDArray result (re_val.dims (), Complex ());

              for (octave_idx_type i = 0; i < re_val.numel (); i++)
                result.xelem (i) = Complex (re_val(i), im_val(i));

              retval = octave_value (new octave_complex_matrix (result));
            }


It doesn't speed things up, but why not write the simpler code below?


diff -r 20b05c5c41a1 libinterp/corefcn/data.cc
--- a/libinterp/corefcn/data.cc Mon Jan 24 16:50:27 2022 +0100
+++ b/libinterp/corefcn/data.cc Tue Jan 25 13:59:12 2022 -0800
@@ -3705,7 +3705,7 @@ complex ([1, 2], [3, 4])
               for (octave_idx_type i = 0; i < re_val.numel (); i++)
                 result.xelem (i) = Complex (re_val(i), im_val(i));

-              retval = octave_value (new octave_complex_matrix (result));
+              retval = octave_value (result);
             }
         }
     }


Similarly, it doesn't make a timing difference but why fill the array with a known value when the next for loop is guaranteed to overwrite each value?  It seems better to just initialize the array with a dimension vector and then fill it.


-              ComplexNDArray result (re_val.dims (), Complex ());
+              ComplexNDArray result (re_val.dims ());


Lastly, since the for loop is guaranteed to run over only the existing elements it theoretically should be faster us use xelem() rather than operator () which calls elem.  In practice, it seems to have slowed things down by ~2% but there is noise there.


               for (octave_idx_type i = 0; i < re_val.numel (); i++)
-                result.xelem (i) = Complex (re_val(i), im_val(i));
+                result.xelem (i) = Complex (re_val.xelem (i), im_val.xelem (i));


In the past avoiding elem() and xelem() and using direct array [] access has been faster by about 10%, but it wasn't the case here.  I tried


              const double *re_data = re_val.data ();
              const double *im_data = im_val.data ();
              Complex *result_data = result.fortran_vec ();
              for (octave_idx_type i = 0; i < re_val.numel (); i++)
                result_data[i] = Complex (re_data[i], im_data[i]);


Even going to great lengths didn't seem to produce a measurable improvement


              const double *re_data = re_val.data ();
              const double *im_data = im_val.data ();
              double *result_data = reinterpret_cast<double *> (result.fortran_vec ());
              octave_idx_type n = re_val.numel ();
              for (octave_idx_type i = 0; i < n; i++)
                {
                result_data[2*i] = re_data[i];
                result_data[2*i + 1] = im_data[i];
                }


So, I don't really have a good way to speed this up but I wonder whether we can at least get rid of the new operator since it is so different from all the other functions in data.cc

Rik <rik5>
Group administrator
Tue 25 Jan 2022 05:48:46 PM UTC, comment #5: 

So, the original test case was unfair to complex().  The code was


z = complex (0, val);
OR
z = i*val;


In the multiplication situation I wasn't creating a real part while complex had to deal with that even though the value was 0.  A proper comparison is


z = 0 + i*val;


Nevertheless, on my machine the vector has to be pretty big, ~5,500, before complex is faster.

Rik <rik5>
Group administrator
Tue 25 Jan 2022 12:41:11 AM UTC, comment #4: 

Finally found a situation, vectors for both RE and IM parts, where complex is faster.


N = 1e5;

val = rand (1e4, 1);

tic;
for i = 1 : N
  #z = complex (val, val);
  z = val + 1i*val;
endfor
bm = toc


Timings:
complex() : 2.5 seconds
multiply  : 3.3 seconds


(file #52722, file #52723)

Rik <rik5>
Group administrator
Tue 25 Jan 2022 12:33:30 AM UTC, comment #3: 

I was specifically testing the scalar case where the scalar is a real that needs to be made complex.

I modified to try a vector using the code below


N = 1e5;

val = rand (1e4, 1);

tic;
for i = 1 : N
  z = complex (0, val);
  #z = 1i * val;
endfor
bm = toc


Results still favor using multiplication
complex() : 3.1 seconds
multiply  : 1.7 seconds

@Arun: I think you need to run the test multiple times in a for loop in order to get adequate statistics.  If you are executing only one time it is possible that an interrupt happens and the CPU goes off for a while working on something else and then returns to the Octave thread.  If that happens only for one of the test conditions then it is unfairly penalized.

Rik <rik5>
Group administrator
Mon 24 Jan 2022 10:48:42 PM UTC, comment #2: 

I am getting the opposite result: calling complex (re, im) is faster than doing the calculation re + i*im.

Scalar, N = 1:

>> N = 1; re = rand(N,1); im = rand(N,1); tic; z1 = re + i*im; toc, tic; z2 = complex (re, im); toc, assert(all(z1==z2))
Elapsed time is 1.21593e-05 seconds.     ## time to do re + i*im
Elapsed time is 5.96046e-06 seconds.     ## time to do complex (re, im)


Vector, N = 1e7:

>> N = 1e7; re = rand(N,1); im = rand(N,1); tic; z1 = re + i*im; toc, tic; z2 = complex (re, im); toc, assert(all(z1==z2))

Elapsed time is 0.0973542 seconds.     ## time to do re + i*im
Elapsed time is 0.059479 seconds.      ## time to do complex (re, im)


Arun Giridhar <arungiridhar>
Group Member
Mon 24 Jan 2022 10:16:11 PM UTC, comment #1: 

I guess the difference is when val is already complex:



octave:1> val = 1+pi*i
val =  1.0000 + 3.1416i
octave:2> val*i
ans = -3.1416 + 1.0000i
octave:3> complex(val)
ans =  1.0000 + 3.1416i


Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Mon 24 Jan 2022 10:00:31 PM UTC, original submission:  

I recently benchmarked the complex() function for a scalar and found it very slow.  I don't know why that is, but either it should be improved or we should add something to the documentation about alternative constructs to create a complex number.

Benchmark code was


N = 1e5;

val = pi;

tic;
for i = 1 : N
  z = complex (0, val);
  #z = 1i * val;
endfor
bm = toc


Timing on my machine were

complex() : 0.33 seconds
multiply  : 0.08 seconds

It is 4x faster to use multiply to create a complex scalar than the alternative.

Rik <rik5>
Group administrator

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #52722:  bm_complex2.m added by rik5 (110B - text/x-matlab)
file #52723:  bm_complex3.m added by rik5 (116B - text/x-matlab)
file #52720:  bm_complex.m added by rik5 (99B - 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 mmuetzel (Posted a comment)
  • -email is unavailable- added by hg200 (Posted a comment)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by rik5
  • -email is unavailable- added by arungiridhar (Posted a comment)
  • -email is unavailable- added by dasergatskov (Posted a comment)
  • -email is unavailable- added by rik5 (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-01-28 rik5 StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2022-01-25 rik5 Carbon-Copy- Added jwe
    2022-01-25 rik5 Attached File- Added bm_complex2.m, #52722
        Attached File- Added bm_complex3.m, #52723
    2022-01-24 rik5 Attached File- Added bm_complex.m, #52720

    Back to the top

    Powered by Savane 3.13-758e.
    Corresponding source code