Norm of the difference of each row of two matrices of different dimensions

Hello, I have two matrices A and B of dimensions m-by-3 and n-by-3 respectively where n < m (they are basically RGB values of an image). For sake of simplicity assume m = 8 and n = 4. For each row in m, I need to calculate the the norm of the difference of that row with the each row of n. The way I am doing it now is using a nested for loop, but my m and n are way higher. So this is taking forever to run. I was wondering if there is any other efficient way to do it.

Just for the sake of completion here is my code for performing this operation:

for i = 1:size(m, 1)
 for j = 1:size(n, 1)
  d(j, 1) = norm(m(i, 1) - n(j, 1))
 end
 % I do some stuff with d here, but its not important to the 
   question.
end

Thanks.

Answer 1

[Matt Fig]

Did you properly pre-allocate d before the loop? This would be your largest time killer if not...

Also, your example loops don't look like what you described. You are only storing the value of the norm for the final value of i. Also, I thought you were doing this by rows, yet you are only calculating the norm of the difference between two scalars. Your verbal description sounds more like this:

m = 8;
n = 4;
A = round(rand(m,3)*10);
B = round(rand(n,3)*10);
cnt = 1;
d = zeros(32,1);  % Pre-allocate d

for ii = 1:m 
    for jj = 1:n 
        d(cnt) = norm(A(ii,:) - B(jj,:)) ;
        cnt = cnt + 1;
    end 
end

If not, please modify this example and/or your description. Perhaps you are not meaning to store all of the norms? But then why only use the first column when you said you were using rows?

.

.

EDIT

O.k., So if your original code looks like this:

for ii = 1:m
    for jj = 1:n 
        d(jj, 1) = norm(A(ii, :) - B(jj, :)) ;
    end 
end

Then this is much faster:

t = ones(1,size(B,1),'single');  % Use as an index into A
for ii = 1:m
    Ar = A(ii,:);
    d = sqrt(sum((Ar(t,:) - B).^2,2));
end

Answer 2

[Sudarshan]

Hi,

1)Yes. d is preallocated with d = zeros(size(n, 1), 1).

2) If you look at the code, I have a comment indicating that I do something with d, so for each iteration of the outer loop I use the calculated d vector. I only need it for each iteration of the outer loop.

3) You are right. That was a typo on my part. It is a difference between 2 three element vectors.

Here's is my actual code for that loop. Please ignore the names of the variables, they pertain to my project. But the general idea is the same:

for pix = 1:npixels
      for m = 1:size(means, 1)
          dist(m, 1) = norm(yClustered(pix, 1:end-1) - means(m, :));
      end
      [mDist, mIdx] = min(dist);
      yClustered(pix, 4) = mIdx;
      means(mIdx, :) = mean(yClustered(yClustered(:, end) == mIdx, 1:end-1));
end

Like I mentioned, here dist is preallocated with dist = zeros(size(means, 1), 1). Thanks.

[Matt Fig]

See my response above. Also note that indexing into an array with 'end' in a loop is slower than passing the size if known...

Answer 3

[Teja Muppirala]

What's wrong with BSXFUN here?

for pix = 1:npixels
      dist = sum(bsxfun(@minus,means,yClustered(pix, 1:end-1)).^2,2);
      [mDist, mIdx] = min(dist);
      yClustered(pix, 4) = mIdx;
      means(mIdx, :) = mean(yClustered(yClustered(:, end) == mIdx, 1:end-1));
end

[Sudarshan]

Speed-wise I wonder whether this would be slower than Matt's response, since you need call a function every iteration. But I could be wrong.

Answer 4

[Sudarshan]

I used Matt's answer in my code. It is definitely faster but not as fast as I hoped as 110600 and the max number of rows for the other matrix is 256. To slow things further, I need to do about 100 iterations of the pixel loop (or terminate on an error condition). For 1 iteration of the complete pixel loop takes about 2 mins with the other matrix length of 2.

I was thinking maybe there would be a way where instead of having dist be a vector, I could have it be a matrix where each element of each row would correspond to the norm of the current pixel with all the means and the rows themselves will be the pixels. In other words, I have a matrix A (110600 x 3) and B (256 x 3) and get d (110600 x 256). d(1, 1) would be norm(A(:, 1) - B(:, 1), d(1, 2) would be norm(A(:, 1 - B(:, 2)) and so on. I have been cracking my head to do this without loops but have come up with no idea. Any help is appreciated.

Thanks.