Fill matrix efficiently without loop

9 views (last 30 days)
Markus Mikschi
Markus Mikschi on 21 Feb 2017
Edited: Jan on 21 Feb 2017
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 (=
  1 Comment
Jan
Jan on 21 Feb 2017
This is a nice example for investigations in the efficiency of loops and vectorization.

Sign in to comment.

Answers (1)

Jan
Jan on 21 Feb 2017
Edited: Jan on 21 Feb 2017
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.
  1 Comment
Jan
Jan on 21 Feb 2017
Edited: Jan on 21 Feb 2017
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.

Sign in to comment.

Categories

Find more on Matrices and Arrays in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!