bsxfun doesn't work as I expect on a constant function

Question

In Matlab R2016a, I have a large set of small X-vectors and Y-vectors which are paired (e.g. 10,000 1x3 X-vectors paired with 10,000 1x3 Y vectors). For each {X,Y} pair, I want to calculate a 2-scalar-argument function for every pairwise combination of the elements in X and Y, (so in my example I would get 10,000 3x3 matrices).

I thought I could use bsxfun to perform these calculations, but it doesn't work when I try to do some simple tests. bsxfun(@(x,y) x*y,[1 2],[1 2]') returns:

ans =

     1     2
     2     4

Which is what I would expect. However, bsxfun(@(x,y) 1,[1 2],[1 2]') returns:

Error using bsxfun
Specified function handle produces invalid output dimensions. The function handle
must be a binary elementwise function.

Which makes no sense. The function handle is a binary elementwise function that always returns the scalar 1, so bsxfun should give the same result as ones(2,2), unless I'm not understanding how bsxfun works.

Answer 1

The inputs to the function handle that are passed to bsxfun are not scalars. In versions prior to R2016b, the inputs are either scalar or they are the same size.

FUNC can also be a handle to any binary element-wise function not listed above. A binary element-wise function in the form of C = FUNC(A,B) accepts arrays A and B of arbitrary but equal size and returns output of the same size. Each element in the output array C is the result of an operation on the corresponding elements of A and B only. FUNC must also support scalar expansion, such that if A or B is a scalar, C is the result of applying the scalar to every element in the other input array.

In releases since R2016b, they do not have to be equal sizes, but should be compatible sizes

In the example you have shown, the first input to the function handle is a scalar and the second is a vector (y) and the function is evaluated for every element of x and the output is expected to be the size of y

In the case you've posted, the call to bsxfun is essentially the equivalent of:

x = [1 2];
y = [1 2].';

yourfunc = @(x,y)x * y;

for k = 1:numel(x)
    output(:,k) = yourfunc(x(k), y)
end

If you want to return a 1 for every entry, you need to replace your function with something that yields the appropriately sized output.

bsxfun(@(x,y)ones(max(size(x), size(y))), [1 2], [1 2]')

How you formulate the function handle really depends upon your specific problem