How can I use bsxfun to replace for-loop in this code?

Dear all,

I want to use the function 'bsxfun' instead of inner for-loop. My code is as following. However, the result, matrix B, calculated by 'bsxfun' is error. I feel very confused with how 'bsxfun' works. Could anyone help me to understand the result and use 'bsxfun' correctly instead of inner for-loop?

Thanks very much!

sequence = [1,2,3,4,5];
avg = @(i, j) mean(sequence(i:j));

for i = 1:4

      for j = i+1:5
          A(i,j) = avg(i,j);
      end

      B(i, i+1:5) = bsxfun(avg, i, i+1:5);

end

[Adam]

What is the end result you are looking for? You should state that clearly ahead of being focused entirely on using a specific function.

[Shi Shi]

Thanks for your answer. The result calculated by inner for-loop is just what I want. For each pair (i,j), where i<j, I want to calculate the average of sequence(i:j).

Answer 1

[Guillaume]

"I want to use the function 'bsxfun' instead of inner for-loop" Rather than fixating on a function that may not be a solution to your problem, tell us what you want to do. See XY problem.

Your avg function cannot be used with bsxfun. The documentation of bsxfun is very explicit: fun must support scalar expansion. There is no way to support scalar expansion with an expression that includes a colon. So bsxfun is not an option at all.

There might be some obscure way to generate your A array without a loop (some weird combination of cumsum, toeplitz, and some other functions maybe) but honestly, you're better off with your loop approach. It will be a lot clearer and most likely, just as fast.

[Shi Shi]

Thanks a lot for your reply. What I want is the average of sequence(i:j) for each pair (i,j), where i<j. I used for-for-loop before, anyway, it runs too slowly when the size of sequence was too large, in size of millions. Therefore, I try 'bsxfun' and failed as you see :-( Would you know any way to fast my code? Thank you advanced.

[Guillaume]

An obscure non-loopy version to generate that in R2016b or later, which doesn't need bsxfun would be:

triu((cumsum(sequence) - [0; cumsum(sequence(1:end-1)).']) ./ toeplitz(1:numel(sequence)), 1)

in versions < R2016b, bsxfun is required, so this becomes:

triu(bsxfun(@minus, cumsum(sequence), [0; cumsum(sequence(1:end-1)).']) ./ toeplitz(1:numel(sequence)), 1)

You'd better write some comment in your code if you use that, explaining exactly what the expression does. The first cumsum sort of slides along the end of the sequence, the second cumsum sort of slides along the beginning. The difference of both is the sum of sequence(i:j). toeplitz calculates the number of elements in the sequence, so the ratio is the mean. triu(..., 1) to keep the values above the main diagonal.

[Shi Shi]

Thanks again. I will try the code later. And I still wonder the process detail of 'bsxfun' in the code. For example, when i = 3, bsxfun(avg, 3, 4:5) gets the result 3.5000, which is avg(3,4), rather than [avg(3,4), avg(3,5)]. Why doesn't bsxfun replicate 3 along [4:5] here? Maybe I have used bsxfun uncorrectly before?

[Walter Roberson]

bsxfun(avg, 3, 4:5) fits case #2 that I gave: that if either of the two values are scalars then the function is called on the arguments as-is, with the expectation that the function will do scalar expansion. bxsfun requires that the given function do scalar expansion. scalar expansion means that given f(scalar,vector) then f must act like f(repmat(scalar, size(vector)), vector) . So in your case, your avg would be required to treat avg(3, 4:5) the same as avg([3 3], [4 5])

Replacement code:

B(i, i+1:5) = arrayfun(avg, repmat(i, 1, 5-i), i+1:5);

[Shi Shi]

Thanks very much for your explain. Why does avg([3,3], [4,5]) get 3.5, which means [3,3]:[4,5] = [3,4]? I feel so confused about this result.

[Guillaume]

The reason you can't use your avg function with bsxfun and the reason it returns 3.5 are all down to the : (colon) operator.

As per the documentation of colon, "If you specify nonscalar arrays, then MATLAB interprets j:i:k as j(1):i(1):k(1)." So it interprets [3, 3] : [4, 5] simply as 3 : 4, whose mean is indeed 3.5.

In effect your function ignores all but the first element of any vector that it is passed. Hence, why it does not work with bsxfun.

[Shi Shi]

well, I understand now. Thanks very very much.

Answer 2

[Walter Roberson]

https://www.mathworks.com/help/matlab/ref/bsxfun.html#inputarg_fun

"Binary function to apply, specified as a function handle. fun must be a binary (two-input) element-wise function of the form C = fun(A,B) that accepts arrays A and B with compatible sizes. For more information, see Compatible Array Sizes for Basic Operations. fun must support scalar expansion, such that if A or B is a scalar, then C is the result of applying the scalar to every element in the other input array."

In practice what this means for bxfun(@fun, A, B) is:

1. if A and B are the same size(), fun(A,B) is called directly and size(A) = size(B) outputs are expected
2. if either A or B are scalars and the other is not, then fun(A,B) is called directly and size() of the non-scalar is the expected output size
3. otherwise fun(A(:,K), B(K)) is called once for each column K in A with size(A,1)x1 expected output size

Your case matches the first of those, so avg(i, i+1:5) is going to be called -- but your avg code expects the second input to be a scalar rather than a vector.

[Walter Roberson]

For R2016b or later:

S = cumsum([0; sequence(:)]);
result = tril((S - S.') ./ ((1:length(S)).' - (1:length(S))),-1);

For earlier:

S = cumsum([0; sequence(:)]);
result = tril( bsxfun(@minus, S, S.') ./ bsxfun(@minus, (1:length(S)).', (1:length(S))), -1);

[Shi Shi]

Thanks very much for your reply. When avg(i, i+1:5) is called, I hope the result would be the same as arrayfun(avg, repmat(i, length(i+1:5)), i+1:5), because the document says "Whenever a dimension of A or B is singleton (equal to one), bsxfun virtually replicates the array along that dimension to match the other array". But the result is not as I thought. Actually,I have no idea what happens here of my bsxfun.

[Shi Shi]

Thanks very much. It works on small sample very well. I will try it on a large sample later, in size of million.