bugGNU Octave - Bugs: bug #38628, bsxfun slow for complex

 
 

bug #38628: bsxfun slow for complex

Submitter:  John Hunt <huntj>
Submitted:  Fri 29 Mar 2013 07:17:37 PM UTC
   
 
Category:  Interpreter Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Performance
Status:  Confirmed Assigned to:  None
Originator Name:  Open/Closed:  * Open
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

Mon 19 Sep 2022 09:43:36 PM UTC, comment #14: 

Well, it might get messy in the same way that cellfun gets messy.  In cellfun certain functions are detected and accelerated, and the rest are processed normally.

I think it might be similarly useful if functions that are accelerated via broadcasting, such as @plus, were simply forwarded to the broadcasting operator.

But, we also need to be able to handle arbitrary binary functions.  For example, this code


fun = @(a,b) a - exp(b);
a = 1:7;
b = pi*[0 1/4 1/3 1/2 2/3 3/4 1].';
C = bsxfun(fun,a,b)


For those, we should do everything possible to speed things up.

It is not, in general, possible to know the real/complex nature of the result without calculating the function.  For specific operations, say @plus, it is possible to work out the return type in advance given the inputs.  But something like


f = @(a,b) = sqrt (a - b);


Can produce real or complex output even though the inputs a,b are both real.

Rik <rik5>
Group administrator
Sat 17 Sep 2022 05:50:05 PM UTC, comment #13: 

Question. Why are we doing this pattern for each element?

if (! tmp(0).isfloat ())
  {
    have_FloatNDArray = false;
    C = result_FloatNDArray;
    C = cat_op (C, tmp(0), ra_idx);
  }
else if (tmp(0).isreal ())
  result_FloatNDArray.insert
    (tmp(0).float_array_value (), ra_idx);
else
  {
    result_FloatComplexNDArray
      = FloatComplexNDArray (result_FloatNDArray);
    result_FloatComplexNDArray.insert
      (tmp(0).float_complex_array_value (), ra_idx);
    have_FloatNDArray = false;
    have_FloatComplexNDArray = true;
  }

This is casting the partial result generated so far into a new type, then appending the current element's result, and declaring that to be the new augmented result.

From what I've seen so far, some things like sparse/full or single/double or float/integer can be decided up-front, so the above checks could be done once before starting the loop. Does that extend to real/complex as well? If the first few elements of a complex array have no imaginary part, are they being treated as real, and is that why the above code pattern exists to allow the whole result to become complex?

Arun Giridhar <arungiridhar>
Group Member
Sat 17 Sep 2022 11:03:30 AM UTC, comment #12: 

I started work on this today.

Digging through bsxfun.cc, I found that the slow code path happens only if the following conditions are all met:
1. A and B are different types.
2. A and B are different sizes.
3. Neither A nor B is empty in any dimension.

If any of the conditions is not met, there's a faster code path available.

In the slow code path, there's a 150-line for-loop with a large if-else tree inside it, which is done for each element of the result.

The first big question is whether that whole loop can be eliminated and replaced with a call to broadcasting. (I don't know that yet and would like to hear your thoughts on the viability of that approach.)

If we do need that big loop, the next question is whether we can do all the type-checking up front and then end up with a bunch of different tight loops with no decisions inside, one of which will be chosen based on the type.

I instrumented the code with some cout statements to trace the code paths. I also added comments for myself as I work on this.

The hg diff is attached.

Test based on comment #10, packed in one line to ease command line recall:

>> Nc = 8; F = complex (randn(7,11189), randn(7,11189)); tic; G1 = bsxfun(@times, eye(Nc, Nc), shiftdim(F, -2)); toc, tic; Gr = bsxfun(@times, eye(Nc, Nc), shiftdim(real(F), -2)); Gi = bsxfun(@times, eye(Nc, Nc), shiftdim(imag(F), -2)); G2 = complex(Gr, Gi); toc, assert (all (G1(:) == G2(:)))


Output that was used to deduce the above rules about which code path is used where:

maybe_optimized_builtin with name = times
        retval.numel == 0
Entering special:
        dva.numel == 64
        dvb.numel == 78323
Code path 3
Elapsed time is 2.46602 seconds.
maybe_optimized_builtin with name = times
        retval.numel == 5012672
maybe_optimized_builtin with name = times
        retval.numel == 5012672
Elapsed time is 0.068238 seconds.



(file #53708)

Arun Giridhar <arungiridhar>
Group Member
Thu 19 Nov 2020 02:57:22 PM UTC, comment #11: 

I completely agree.  A slight difference in coding style should not cause a 10X difference in performance.  Seem my comment #9 (admittedly from 2018) where I tracked the problem down somewhat.

Rik <rik5>
Group administrator
Thu 19 Nov 2020 10:49:37 AM UTC, comment #10: 

Hi, here is a piece of code that I found to be very slow in octave:

G = bsxfun(@times, eye(Nc, Nc), shiftdim(F, -2));

(dimensions for F (complex array) were 7x11189 and Nc was equal to 8, computation time was 8.96s)

Then I found this discussion. I effectively managed to recover a normal execution time using the following trick:


Gr = bsxfun(@times, eye(Nc, Nc), shiftdim(real(F), -2));
Gi = bsxfun(@times, eye(Nc, Nc), shiftdim(imag(F), -2));
G = complex(Gr, Gi);

(computation time down to 0.75s !)

That would be great however if the issue could be fixed in a later revision of octave.





Vincent Gras <gvinc>
Thu 22 Mar 2018 06:28:50 PM UTC, comment #9: 

I just reviewed and updated bsxfun for a failure to operate on float complex values (http://hg.savannah.gnu.org/hgweb/octave/rev/bff4a7c7bc39).

I think we should, in the future, try to get rid of the duplicate code paths and simply forward calls to bsxfun on to the code used for broadcasting.

From what I could determine, the issue seems to be that Octave calls complex_array_value() for every single element of the array.  This is slower than calling complex_array_value() once on the array, and then operating on each of the elements.

Rik <rik5>
Group administrator
Sun 20 Nov 2016 09:58:41 PM UTC, comment #8: 

This (performance) issue is still present in Octave 4.2.0.

Hartmut <hardy>
Tue 13 Oct 2015 05:07:58 PM UTC, comment #7: 

As a follow-up to comment #2, I think the reason it is slower is that bsxfun doesn't have any special optimizations in place for mixed-type operations.

John W. Eaton <jwe>
Group administrator
Tue 13 Oct 2015 03:07:36 PM UTC, comment #6: 

@Jordi: Did you write the automatic broadcasting code?  I haven't followed this much, but having two different code paths does seem problematic.  Certainly we need to keep bsxfun to be Matlab-compatible, but maybe it should call to the automatic broadcasting code since it seems to be faster.

Rik <rik5>
Group administrator
Tue 13 Oct 2015 11:17:01 AM UTC, comment #5: 

The broadcast warning is gone, but it was just replaced by a different warning tag (that change is the subject of another bug report).

Also, since @times is handled as a special case in bsxfun, why is it slow for complex?

I was looking at bsxfun last night to see why it is slow for @bitxor, and I got a bit lost.  It seems quite complicated.  Also bad that we now have two ways of doing these operations.  I thought that broadcasting operators were just syntactic sugar for bsxfun, but it seems there are actually two different code paths for the same thing?

John W. Eaton <jwe>
Group administrator
Tue 13 Oct 2015 04:30:11 AM UTC, comment #4: 

I replaced bsxfun in all cases in Octave core where there might be complex inputs.  See this cset (http://hg.savannah.gnu.org/hgweb/octave/rev/a3b9ee5c040a)

Rik <rik5>
Group administrator
Fri 09 Oct 2015 04:50:09 PM UTC, comment #3: 

The bug is still there in 4.0.0.

Anyway, now that the "Octave:broadcast" warning is gone, trapz could use broadcasting ? (patch attached for trapz).


(file #35126)

ederag <ederag>
Thu 04 Apr 2013 02:38:11 AM UTC, comment #2: 

Exceedingly strange.  This is, perhaps, a clue.  If dx is promoted to a complex matrix (with 0 for the imaginary part) then the bsxfun call goes back to being extremely fast.


dxc = complex (dx);
tic; r = bsxfun (@times, dxc, ym); toc
Elapsed time is 0.005 seconds.


Rik <rik5>
Group administrator
Wed 03 Apr 2013 07:01:30 PM UTC, comment #1: 

Hm, this is quite bizarre. I'll try to remember to look into this later.

Jordi GutiƩrrez Hermoso <jordigh>
Group Member
Fri 29 Mar 2013 07:17:37 PM UTC, original submission:  

bsxfun is slow for complex array. Much slower than broadcasting.
trapz which uses bsxfun is impacted for instance.


octave:1> dx = ones(1, 1, 1000);
octave:2> ym = ones(1, 1, 1000, 1, 201)+i;
octave:3> tic; r = bsxfun(@times, dx, ym); toc
Elapsed time is 1.45173 seconds.
octave:4> tic; r = bsxfun(@times, dx, ym); toc
Elapsed time is 1.46083 seconds.
octave:5> tic; r1 = dx .* ym; toc
warning: product: automatic broadcasting operation applied
Elapsed time is 0.000883102 seconds.
octave:6> assert(r1, r)
octave:7>


these outputs are from dev
hg id -in
a695ee2dc17e+ 16392+

but it holds also for 3.6.3 and 3.6.4

John Hunt <huntj>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #35126:  octave-4.0.0_trapz.m.patch added by ederag (456B - text/x-patch)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by arungiridhar (Updated the item)
  • -email is unavailable- added by gvinc (Posted a comment)
  • -email is unavailable- added by hardy (Posted a comment)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by ederag (Updated the item)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by jordigh (Posted a comment)
  • -email is unavailable- added by huntj (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 3 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-09-17 arungiridhar Attached File- Added bsxfun.commented.patch, #53708
    2015-10-09 ederag Attached File- Added octave-4.0.0_trapz.m.patch, #35126
    2013-04-03 jordigh StatusNone Confirmed

    Back to the top

    Powered by Savane 3.13-758e.
    Corresponding source code