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.

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)

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.

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

@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.

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)

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.

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