2012-07-27 75 views
4

我希望我的問題有一個非常簡單的解決方案。我只是無法找到它:第三維向量的Matlab產品

假設你有兩個向量(一個列向量,一個行向量)A,B:

A = [1,2,3] 
B = [4;5;6] 

如果我們乘他們如下,我們得到了一個矩陣:

>> B*A 
ans = 
4  8 12 
5 10 15 
6 12 18 

現在我的問題是:我有大小m × n × p的兩個3D矩陣和m × n × q

沿尺寸想象n我們有像素,每個像素有一個向量(長度爲p或q)。 現在我想要的是,爲每個對應的像素乘以這兩個圖像的向量,使得對於每個像素我得到一個矩陣,因此最終總共是一個4D矩陣。

我該如何有效地做到這一點?

回答

4

循環在Matlab不再是令人擔心的事情,或避免本身

當然,在使用它們時應該非常小心,但是JIT可以處理多種循環,甚至可以超越內建函數來提高性能。

請看下面的測試案例:

clc 

m = 512; n = 384; 
p = 5;  q = 3; 

A = rand(m,n,p); % some sample data 
B = rand(m,n,q); % some sample data 

%% non-loop approach 

tic 
A2 = reshape(A,[],p); 
B2 = reshape(B,[],q); 
C2 = arrayfun(@(ii) A2(ii,:)'*B2(ii,:),1:m*n,'uni',false); 
C0 = permute(reshape(cell2mat(C2),p,q,m,n),[3 4 1 2]); 
toc 

%% looped approach, simplest 

tic 
C = zeros(m,n,p,q); 
for mm = 1:m 
    for nn = 1:n   
     C(mm,nn,:,:) = ... 
      squeeze(A(mm,nn,:))*squeeze(B(mm,nn,:)).'; 
    end 
end 
toc 

% check for equality 
all(C0(:) == C(:)) 

%% looped approach, slightly optimized 

tic 
C = zeros(m,n,p,q); 
pp = zeros(p,1); 
qq = zeros(1,q); 
for mm = 1:m 
    for nn = 1:n 
     pp(:) = A(mm,nn,:); 
     qq(:) = B(mm,nn,:); 
     C(mm,nn,:,:) = pp*qq; 
    end 
end 
toc 

% check for equality 
all(C0(:) == C(:)) 

%% looped approach, optimized 

tic 
C = zeros(p,q,m*n); 
A2 = reshape(A,[],p); 
B2 = reshape(B,[],q); 
for mn = 1:m*n 
    C(:,:,mn) = A2(mn,:).'*B2(mn,:); 
end 
C = permute(reshape(C, p,q,m,n), [3,4,1,2]); 
toc 

% check for equality 
all(C0(:) == C(:)) 

結果:

Elapsed time is 3.955728 seconds. 
Elapsed time is 21.013715 seconds. 
ans = 
    1 
Elapsed time is 1.334897 seconds. 
ans = 
    1 
Elapsed time is 0.573624 seconds. 
ans = 
    1 

無論表現,我也覺得最後一種情況很多比非循環情況下更爲直觀性和可讀性。

+0

哇,這是一個很酷的事情。我不知道MATLAB能做這種事情。我會嘗試一下,如果它實際上以這種方式工作,速度更快,我會用它;) – 2012-07-29 11:31:46

+0

@rody_o:+1進行全面比較,再加上我無法就JIT加速達成一致 – Amro 2012-07-29 22:05:06

3

對於某些reshapingarrayfunpermute

m=5; 
n=4; 
p=3; 
q=2; 
A=randi(10,m,n,p); %some sample data 
B=randi(10,m,n,q); %some sample data 

A2=reshape(A,[],p); 
B2=reshape(B,[],q); 
C2=arrayfun(@(ii) A2(ii,:)'*B2(ii,:),1:m*n,'uni',false); 

C=permute(reshape(cell2mat(C2),p,q,m,n),[3 4 1 2]); 

擊穿:

  • 前兩個重塑改變A和B mxnx(p or q)矩陣轉換成(m*n)x(p or q)格式
  • 使得arrayfun可以很容易地循環通過他們計算行的向量乘積
  • 然後cell2mat,重塑和置換結果改回mxnxpxq格式
0

我用rody_o的解決方案,並修改它擺脫了重塑和置換的:

C = zeros(m*n, p, q); 
A2 = reshape(A,[],p); 
B2 = reshape(B,[],q); 
for mn = 1:m*n 
    C(mn,:,:) = A2(mn,:).' * B2(mn,:); 
end 
+0

哦,我仍然會需要重塑,但我沒有意識到,導致「所有(C0(:) == C(:))」不檢查此 – 2012-07-30 07:57:34