Can I speed this code up, looking for similarity between two 3-d matrices.

Is it possible to vectorize the following code in order to speed it up?

This little routine gets called tens of thousands of times in some code I'm working on, and it is listed as the major bottle neck in profile viewer.

What it does is shift a small 3D matrix A through a large 3D matrix B looking for the best match in B to A. My current thoughts are that it can be vectorized by reshaping and replicating A and B, but I can't figure out how to do it.

temp_position = zeros(3,1); 
new_position = zeros(3,1);

ad_current = Inf;
for ii = 1:size(B,1)-size(A,1)+1
  for jj = 1:size(B,2)-size(A,2)+1
    for kk = 1:size(B,3)-size(A,3)+1
      ad_new = sum(reshape(abs(B(ii:ii+size(A,1)-1,jj:jj+size(A,2)-1,kk:kk+size(A,3)-1) - A),[],1));
      if ad_new < ad_current
        ad_current = ad_new;                 
        temp_position = [ii,jj,kk];
      end
    end
  end
end
new_position = ... something + temp_position

[Jan]

Is "B(jj:jj+size(A,1)-1,ii:ii+size(A,2)-1, ..." a typo? Do you mean "ii"<->"jj" here?

[tlawren]

Yes, this is "kind of" a typo too. I say kind of, because it actually isn't. The data I'm working on is stored in a weird way, so my code is written to where it makes sense within the context of the data. It is hard to explain, but it works. I tried to change it to make it more normal and readable for this post. Given that I've had to edit my post twice, I apparently failed at doing this.

Answer 1

[Bruno Luong]

The trick here is loop on A (small size) and vectorized on B

% Test data
B=rand(50,60,100);
A=rand(2,3,4);

%%Engine
szA = size(A);
szB = size(B);
szC = szB-szA+1;
C = zeros(szC);

szBs = ones(1,6);
szBs([1 3 5]) = szA;
AA = reshape(A, szBs);
for n=1:numel(A)
    first = cell(1,3);
    [first{:}] = ind2sub(szA,n);
    first = cat(2,first{:});
    nb = floor((szB-first+1)./szA);
    lgt = nb.*szA;
    last = first + lgt - 1;
    iB = arrayfun(@(i,j) i:j, first, last, 'Unif', false);
    Bs = B(iB{:});
    szBs([2 4 6]) = nb;
    Bs = reshape(Bs, szBs);
    BmA = bsxfun(@minus, Bs, AA);
    d = sum(sum(sum(abs(BmA),1),3),5);
    iC = arrayfun(@(i,s,j) i:s:j, first, szA, last, 'Unif', false);
    C(iC{:}) = reshape(d, nb);
end
[mindiff loc] = min(C(:));
[ii jj kk] = ind2sub(szC, loc);
loc = [ii jj kk];

% Check
disp(mindiff)
disp(loc)

Answer 2

[Jan]

Just some marginal changes:

sizeA = size(A);
sA3m1 = sizeA(3) - 1;
sizeB = size(B);
ad_current = Inf;
for ii = 1:sizeB(1)-sizeA(1)+1
  Q = B(ii:ii+sizeA(1)-1, :, :);
  for jj = 1:sizeB(2)-sizeA(2)+1
    P = Q(:, jj:jj+sizeA(2)-1, :);
    for kk = 1:sizeB(3)-sA3m1
      ad_new = sum(reshape(abs( ...
               P(:, :, kk:kk+sA3m1) - A),[],1));
      if ad_new < ad_current
        ad_current = ad_new;                 
        temp_position = [ii,jj,kk];
      end
    end
  end
end

EDITED: Q(:,:, kk:...) -> P(:, :, kk:...) Thanks Bruno

Answer 3

[Doug Hull]

Is there a way that convn can be used to do this? I don't know the details, but I think it might help.

Answer 4

[tlawren]

Sorry! Jan, you're absolutely right. The loop is currently commented out in my code (to restrict my searches), and the kk calculations are removed. In general, they will be there. I edited my original post to correct this.

[Jan]

I've deleted my intermediate answer, such that this posting is not useful anymore. Consider to delete it.

Answer 5

[tlawren]

Thanks for all the replies!!!

@Jan - Your suggestions are easy to implement and the changes are giving me about a 25% decrease in my computation times.

@Bruno - Your little routine blows my current routine out of the water in terms of computational time. For example, my routine takes 70 s to do what your routine does in 3 s. However, when I run it on multiple As on the same B, your routine takes a lot longer. The problem lies in my larger code, so I need to figure out how to change it (vectorize it).

If you have any suggestion on how to adapt your above routine to do multiple As on the same B, I would love to see them.

Thanks again!

[Bruno Luong]

Are the As same size?

I don't understand how you run multiple A in your code. If you run multiple As with a for loop, I can't see how you function can be suddenly faster.

When asking a question, attach a little code is better than 1000 words.

Answer 6

[tlawren]

Yes, the A's are all the same size and they are being looped over. The logic behind my code goes as follows: get an A, get a B (the same for a set of As), do the mindiff routine, and move on to the next A. Once all the As in a given set are done, move on to a new set. The new set will have a new B. I have about 700 As per set and about 300 sets, so 300 Bs too.

I'm sure there is a way to do a whole set of A's at once, but I don't know how. I'm new to vectorization, so thinking about problems in a vectorized way is hard for me right now.

[Bruno Luong]

Something is fishy. If you loop over A, the compute time is just proportional with the number of A. It did not make sense to me why your code is more efficient when running with many As.

[tlawren]

You are right. I created a test scenario and ran the two routines to see what the times would be. Your routine finished in 27 seconds and mine never finished. I killed it at 575 seconds.

However, when I plug your routine into my project code, it runs slower. I haven't figured out why, so I need to study it more.

I can say that...

iC = arrayfun(@(i,s,j) i:s:j, ...

iB = arrayfun(@(i,j) i:j, firs...

[first{:}] = ind2sub(szA,n);

are all taking the most time.

Answer 7

[tlawren]

A test scenario would be running the code on the following...

B = rand(10,200,150,100);
A = rand(10,100,5,5,5);

where I've got 10 sets of As and 100 As per set. There maybe a better way to store the data, but this is the first thing that came to mind. Your code is significantly faster than mine on this test.