Vectorization of multiple embedded for loops

Question

I have the following code that includes 3 iterated for loops in order to create an upper diagonal matrix, I plan on performing on large data set many times and want to make as computationally efficient as possible.

data = magic(3);
n = size(data,1);
W = zeros(n,n);
for i = 1:n
    for j = i:n
        if i==j
            W(i,j)=0;
        else
            for k = 1:n
                temp(1,k) = (data(i,k)-data(j,k))^2;
                sumTemp = sumTemp + temp(1,k);
            end
          W(i,j)=sqrt(sumTemp);
        end
        temp = 0;
        sumTemp = 0;        
    end
end

Answer should look like:

[0 6.4807 9.7980
0 0 6.4807
0 0 0]

I am working it hard right now, but figure I would throw it out there in case anyone has any suggestions that would save me hours of fiddling around.

What is `m`? You better explain what is the inner operation, it's gon na be easier to propose a solution. — Oleg, Oct 24 '13 at 18:32
Good pickup, that was typo - should also be n - I fixed it in original post. — Keith D, Oct 24 '13 at 18:49
you can try it yourself, read this one for example: http://stackoverflow.com/questions/17540681/improving-matlab-matrix-construction-code-or-code-vectorization-for-begginers/17585725#17585725 — bla, Oct 24 '13 at 19:13
Daniel - thanks. Let me examine this. The www was old variable name - I have since update original post. I like the explanation link above!! I'd vote you up but I don't have enough reputation points yet. — Keith D, Oct 24 '13 at 19:45
@user2913719 - I think the loop `k = 1:n` should be `k = 1:size(data,2)` to work for a non-square `data` matrix.` My solution will handle square or rectangular with dimensions of size(data,2). — chappjc, Oct 24 '13 at 20:29

Daniel · Accepted Answer · 2013-10-24T20:22:09.960

1

This is hat I have at the moment:

data = magic(3);
n = size(data,1);
W = zeros(n,n);
for i = 1:n
  for j = i+1:n
    W(i,j)= norm(data(i,:)-data(j,:))
    %W(i,j)= sqrt(sum((data(i,:)-data(j,:)).^2));      
  end
end

What I did:

vecorized the inner loop
removed www, which is unused
changed 2nd loop, start at i+1 because nothing is done for i=j
Replaced sqrt((a-b).^2) with norm(a-b)

And now the "full" vectorization:

data = magic(3);
n = size(data,1);
W = zeros(n,n);
tri=triu(ones(n,n),1)>0;
[i,j]=find(tri);
W(tri)=arrayfun(@(i,j)norm(data(i,:)-data(j,:)),i,j)

edited Oct 24 '13 at 20:22

answered Oct 24 '13 at 18:58

Daniel

36,610
3
36
69

2

That sqrt-of-sum-of-squares looks very much like a 2-norm. – Ben Voigt Oct 24 '13 at 20:05

score 1 · Answer 2 · answered Oct 24 '13 at 20:23

Here is a straightforward solution with bsxfun:

Wfull = sqrt(squeeze(sum(bsxfun(@minus,data,permute(data,[3 2 1])).^2,2)))
W = triu(Wfull)

Use this where data is N-by-D, where N is the number of points and D is dimensions. For example,

>> data = magic(3);
>> triu(sqrt(squeeze(sum(bsxfun(@minus,data,permute(data,[3 2 1])).^2,2))))
ans =
         0    6.4807    9.7980
         0         0    6.4807
         0         0         0
>> data = magic(5); data(:,end-1:end)=[]
data =
    17    24     1
    23     5     7
     4     6    13
    10    12    19
    11    18    25
>> triu(sqrt(squeeze(sum(bsxfun(@minus,data,permute(data,[3 2 1])).^2,2))))
ans =
         0   20.8087   25.2389   22.7376   25.4558
         0         0   19.9499   19.0263   25.2389
         0         0         0   10.3923   18.3576
         0         0         0         0    8.5440
         0         0         0         0         0
>>

Thanks - great stuff, I am new to vectorization but loving it. — Keith D, Oct 25 '13 at 18:32
@user2913719 - Great! It can be lots of fun. Note that [`arrayfun` is usually slower than for loops](http://stackoverflow.com/questions/12522888/arrayfun-can-be-significantly-slower-than-an-explicit-loop-in-matlab-why) and not really vectorization IMO. Try to use `bsxfun` or just linear algebra to vectorize. — chappjc, Oct 25 '13 at 19:17
I guess one line was too simple. That's just how it goes at SO sometimes, I suppose. — chappjc, Oct 26 '13 at 18:45

Vectorization of multiple embedded for loops

2 Answers2