0

I have two vectors, say A of size nx1, and B of size 1xm. I want to create a result matrix C (nxm) from a non-linear formula so that

C(i,j) = A(i)/(A(i)+B(j)).

I know I can do this with a loop such as:

for i=1:n,
    C(i,1:m)=A(i)./(A(i)+B(1:m));
end

but is it possible in some way to do this without using a loop?

EDIT: Thanks for your answers! As a small addition, before I saw them a friend suggested the following solution:

C = A*ones(1,m)./(ones(n,1)*B+A*ones(1,m))
yuyu2809
  • 27
  • 1
  • 7

2 Answers2

6

If you are on MATLAB R2016a or earlier you'll want to use bsxfun to accomplish this

result = bsxfun(@rdivide, A, bsxfun(@plus, A, B));

If you are on R2016b or newer, then there is implicit expansion which allows you to remove bsxfun and just apply the element-operators directly

result = A ./ (A + B);

Benchmark

Here is a real benchmark using timeit to compare the execution speed of using bsxfun, repmat, implicit broadcasting and a for loop. As you can see from the results, the bsxfun and implicit broadcasting methods yield the fastest execution time.

function comparision()

    sz = round(linspace(10, 5000, 30));

    times1 = nan(size(sz));
    times2 = nan(size(sz));
    times3 = nan(size(sz));
    times4 = nan(size(sz));

    for k = 1:numel(sz)
        A = rand(sz(k), 1);
        B = rand(1, sz(k));
        times1(k) = timeit(@()option1(A,B));
        A = rand(sz(k), 1);
        B = rand(1, sz(k));
        times2(k) = timeit(@()option2(A,B));
        A = rand(sz(k), 1);
        B = rand(1, sz(k));
        times3(k) = timeit(@()option3(A,B));
        A = rand(sz(k), 1);
        B = rand(1, sz(k));
        times4(k) = timeit(@()option4(A,B));
    end

    figure
    p(1) = plot(sz, times1 * 1000, 'DisplayName', 'bsxfun');
    hold on
    p(2) = plot(sz, times2 * 1000, 'DisplayName', 'repmat');
    p(3) = plot(sz, times3 * 1000, 'DisplayName', 'implicit');
    p(4) = plot(sz, times4 * 1000, 'DisplayName', 'for loop');
    ylabel('Execution time (ms)')
    xlabel('Elements in A')
    legend(p)
end

function C = option1(A,B)
    C = bsxfun(@rdivide, A, bsxfun(@plus, A, B));
end

function C = option2(A,B)
    m = numel(B);
    n = numel(A);
    C = repmat(A,1,m) ./ (repmat(A,1,m) + repmat(B,n,1));
end

function C = option3(A, B)
    C = A ./ (A + B);
end

function C = option4(A, B)
    m = numel(B);
    n = numel(A);
    C = zeros(n, m);
    for i=1:n,
        C(i,1:m)=A(i)./(A(i)+B(1:m));
    end
end

enter image description here

See this answer for more information comparing implicit expansion and bsxfun.

Graham
  • 7,431
  • 18
  • 59
  • 84
Suever
  • 64,497
  • 14
  • 82
  • 101
  • 1
    I'm not sure `B` should be transposed? I get a "Non-singleton dimensions not equal" error unless I run it as `bsxfun(@rdivide, A, bsxfun(@plus, A, B));` – Wolfie Mar 24 '17 at 14:32
  • @Wolfie You're right, I didn't realize that `B` was a row vector. – Suever Mar 24 '17 at 14:33
1

Implicit expansion is the way to go if you have 2016b or newer, as mentioned by Suever. Another approach without that is to do element-wise operations after making A and B the same size as C, using repmat...

A1 = repmat(A,1,m);
B1 = repmat(B,n,1);
C = A1 ./ (A1 + B1);

Or in 1 line...

C = repmat(A,1,m) ./ (repmat(A,1,m) + repmat(B,n,1));

As a benchmark, I ran your loop method, the above repmat method, and Suever's bsxfun method for m = 1000, n = 100, averaging over 1,000 runs each:

Using for loop 0.00290520 sec
Using repmat   0.00014693 sec
Using bsxfun   0.00016402 sec

So for large matrices, repmat and bsxfun are comparable with repmat edging it. For small matrices though, just looping can be quicker than both, especially with a small n as that's your loop variable!

It's worth testing the given methods for your specific use case, as the timing results seem fairly variable depending on the inputs.

Wolfie
  • 27,562
  • 7
  • 28
  • 55
  • 1
    When benchmarking it's best to use `timeit` as that guarantees a relatively stable estimate of execution time. See my updated post for a comparison of all methods. – Suever Mar 24 '17 at 14:48
  • @Suever nice one, I knew `timeit` existed but have never really used it, nice to see your example. I assume that the `repmat` method also uses more memory than `bsxfun`, so I was expecting it to be slower! – Wolfie Mar 24 '17 at 16:53