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

2 次查看（过去 30 天）
tlawren
tlawren 2011-2-18
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
  2 个评论
Jan
Jan 2011-2-19
Is "B(jj:jj+size(A,1)-1,ii:ii+size(A,2)-1, ..." a typo? Do you mean "ii"<->"jj" here?
tlawren
tlawren 2011-2-21
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Bruno Luong
Bruno Luong 2011-2-20
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)

更多回答（6 个）

Jan
Jan 2011-2-20
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
Doug Hull
Doug Hull 2011-2-18
Is there a way that convn can be used to do this? I don't know the details, but I think it might help.
tlawren
tlawren 2011-2-18
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.
tlawren
tlawren 2011-2-21
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!
  1 个评论
Bruno Luong
Bruno Luong 2011-2-21
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.

请先登录，再进行评论。

tlawren
tlawren 2011-2-21
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.
  2 个评论
Bruno Luong
Bruno Luong 2011-2-21
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
tlawren 2011-2-21
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.

请先登录，再进行评论。

tlawren
tlawren 2011-2-21
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.

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by