Fill matrix efficiently without loop

Hey guys!

I need to construct a matrix and want to do it the "right" way. I have heard that the use of loops in Matlab is very inefficient and even frowned upon. However, since I have a programming background it is hard for me to not jump straight to nested loops when it comes to dealing with matrices.

In this particular case I need to construce a nxn Matrix C witch is defined by: c_ij = sqrt((x_i-x_j)^2 + (y_i-y_j)^2 + G). The matrix is obviously symmetrical. I have saved the x and y values in seperat vectors aMaybe some of you recognize this as the Hardy interpolation.

Is there a way to do this elegantly without loops? Any help is highly appreciated (=

[Jan]

This is a nice example for investigations in the efficiency of loops and vectorization.

Answer 1

[Jan]

It is a rumor, that loops are very inefficient in general, because Matlab's JIT acceleration works successfully since version 6.5, R13 (2002). Try this:

function main
x = rand(1, 1000);
y = rand(1, 1000);
G = rand;

tic;
for k = 1:100
  c1 = testA(x, y, G);
end
toc

tic;
for k = 1:100
  c2 = testB(x, y, G);
end
toc
isequal(c1, c2)
end

function c = testA(x, y, G)
% c = sqrt((x - x.') .^ 2 + (y - y.') .^ 2 + G);  % Matlab >= 2016b
c = sqrt(bsxfun(@minus, x, x.') .^ 2 + ...
         bsxfun(@minus, y, y.') .^ 2 + G);  % Matlab < 2016b
end

function c = testB(x, y, G)
nx = numel(x);
ny = numel(y);
c  = zeros(nx, ny);
sG = sqrt(G);
for i2 = 1:ny
  c(i2, i2) = sG;
  for i1 = i2+1:nx
    t = sqrt((x(i1) - x(i2)) ^ 2 + (y(i1) - y(i2)) ^ 2 + G);
    c(i1, i2) = t;
    c(i2, i1) = t;
  end
end
end

Under my old Matlab R2009a on a i5 (I'm going to run this under R2016b in the evening):

 Elapsed time is 2.967072 seconds.  % Vectorized
 Elapsed time is 2.357717 seconds.  % Loop
 1                                  % Results are equal

Here loops can avoid almost the half of the expensive SQRT evaluations and the creation of the large temporary arrays x-x.' and y-y.'. This saves more time than the loop overhead costs in this case.

[Jan]

Well, the timings on my other machine are slightly different:

 2009a/64 on a Core2Duo:
 Elapsed time is 7.994782 seconds.
 Elapsed time is 4.876842 seconds.

 2016b/64:
 Elapsed time is 4.267039 seconds.
 Elapsed time is 4.472175 seconds.

And when I use the automatic expanding of R2016b:

c = sqrt((x - x.') .^ 2 + (y - y.') .^ 2 + G);

this needs 5.41 seconds. This means, that bsxfun rules.