一種快速和有效的方式從一系列產品

Question

AX，AY，AZ的創建矩陣：[N逐N]

B = AA（二進制）

這意味着：

B(i,j)= [Ax(i,j);Ay(i,j);Az(i,j)]*[Ax(i,j) Ay(i,j) Az(i,j)] 

B（i，j）：3×3矩陣。 一到構建體B的方法是：

N=2; 
Ax=rand(N); Ay=rand(N); Az=rand(N);  %# [N-by-N] 
t=1; 
F=zeros(3,3,N^2); 
for i=1:N 
    for j=1:N 
     F(:,:,t)= [Ax(i,j);Ay(i,j);Az(i,j)]*[Ax(i,j) Ay(i,j) Az(i,j)]; 
     t=t+1;       %# t is just a counter 
    end 
end 
%# then we can write 
B = mat2cell(F,3,3,ones(N^2,1)); 
B = reshape(B,N,N)'; 
B = cell2mat(B); 

是否有當N很大更快的方法。

編輯：

感謝您的回答。 （更快） 讓我們把： N = 2; Ax = [1 2; 3 4]; Ay = [5 6; 7 8]; Az = [9 10; 11 12];

B = 

1  5  9  4 12 20 
5 25 45 12 36 60 
9 45 81 20 60 100 
9 21 33 16 32 48 
21 49 77 32 64 96 
33 77 121 48 96 144 

運行：
???錯誤使用==> mtimes 內部矩陣維度必須一致。

如果我寫：P = Ai*Aj;然後

B = 

7 19 31 15 43 71 
23 67 111 31 91 151 
39 115 191 47 139 231 
10 22 34 22 50 78 
34 78 122 46 106 166 
58 134 210 70 162 254 

即從上方 甲不同模式（：，：，1）從[AX（1,1）Ay的（1,1）了Az deffer（1， 1）]

編輯：

N=100; 
Me  :Elapsed time is 1.614244 seconds. 
gnovice :Elapsed time is 0.056575 seconds. 
N=200; 
Me  :Elapsed time is 6.044628 seconds. 
gnovice :Elapsed time is 0.182455 seconds. 
N=400; 
Me  :Elapsed time is 23.775540 seconds. 
gnovice :Elapsed time is 0.756682 seconds. 
Fast! 

rwong: B was not the same. 

編輯：

經過一些修改我的applicati於： 通過gnovice碼

1st code : 19.303310 seconds 
2nd code: 23.128920 seconds 
3rd code: 13.363585 seconds 

似乎任何函數調用像小區，ind2sub ......讓THW循環緩慢，768,16儘可能避免。

symIndex很有趣！謝謝。

Answer 1

這是一個相當簡單和普通實現，它使用一個for循環執行linear indexing，避免了處理3維變量或重塑：

%# General solution: 
%# ---------------- 
B = cell(N); 
for index = 1:N^2 
    A = [Ax(index) Ay(index) Az(index)]; 
    B{index} = A(:)*A; 
end 
B = cell2mat(B); 

編輯＃1：

作爲對如何減少對稱矩陣B執行的計算次數的額外問題的迴應，您可以使用以下修改版本的th上述E的代碼：

%# Symmetric solution #1: 
%# --------------------- 
B = cell(N); 
for index = find(tril(ones(N))).' %'# Loop over the lower triangular part of B 
    A = [Ax(index) Ay(index) Az(index)]; 
    B{index} = A(:)*A; 
    symIndex = N*rem(index-1,N)+ceil(index/N); %# Find the linear index for the 
               %# symmetric element 
    if symIndex ~= index  %# If we're not on the main diagonal 
    B{symIndex} = B{index}; %# then copy the symmetric element 
    end 
end 
B = cell2mat(B); 

然而，在這樣的情況下，你可能會得到更好的性能（或者至少更簡單的查找代碼）由前述的線性索引和使用兩個用於與下標索引環路像這樣：

%# Symmetric solution #2: 
%# --------------------- 
B = cell(N); 
for c = 1:N %# Loop over the columns 
    for r = c:N %# Loop over a subset of the rows 
    A = [Ax(r,c) Ay(r,c) Az(r,c)]; 
    B{r,c} = A(:)*A; 
    if r ~= c   %# If we're not on the main diagonal 
     B{c,r} = B{r,c}; %# then copy the symmetric element 
    end 
    end 
end 
B = cell2mat(B); 

EDIT＃2：

第二對稱溶液可以進一步加快通過移動對角線計算的內循環之外（除去用於一個條件語句的需要）和overwrit因此我們可以使用A而不是B{r,c}來更新對稱單元格條目B{c,r}（即，使用A來代替B{r,c}）。更新一個細胞與另一出現的內容進行一些額外的開銷）：

%# Symmetric solution #3: 
%# --------------------- 
B = cell(N); 
for c = 1:N %# Loop over the columns 
    A = [Ax(c,c) Ay(c,c) Az(c,c)]; 
    B{c,c} = A(:)*A; 
    for r = c+1:N %# Loop over a subset of the rows 
    A = [Ax(r,c) Ay(r,c) Az(r,c)]; 
    A = A(:)*A; 
    B{r,c} = A; 
    B{c,r} = A; 
    end 
end 
B = cell2mat(B); 

這裏是使用下面的示例對稱矩陣Ax，Ay 3級對稱的解決方案的一些定時的結果，和Az：

N = 400; 
Ax = tril(rand(N));  %# Lower triangular matrix 
Ax = Ax+triu(Ax.',1); %'# Add transpose to fill upper triangle 
Ay = tril(rand(N));  %# Lower triangular matrix 
Ay = Ay+triu(Ay.',1); %'# Add transpose to fill upper triangle 
Az = tril(rand(N));  %# Lower triangular matrix 
Az = Az+triu(Az.',1); %'# Add transpose to fill upper triangle 

%# Timing results: 
%# -------------- 
%# Solution #1 = 0.779415 seconds 
%# Solution #2 = 0.704830 seconds 
%# Solution #3 = 0.325920 seconds 

解決方案3的大幅提速主要是通過用局部變量A更新B的單元格內容，而不是用另一個單元格的內容更新一個單元格。換句話說，使用像B{c,r} = B{r,c};這樣的操作複製單元格內容會帶來比我預期的更多的開銷。

Answer 2

A = cat(3, Ax, Ay, Az); % [N-by-N-by-3] 
F = zeros(3, 3, N^2); 
for i = 1:3, 
    for j = 1:3, 
    Ai = A(:,:,i); 
    Aj = A(:,:,j); 
    P = Ai(:) .* Aj(:); 
    F(i,j,:) = reshape(P, [1, 1, N^2]); 
    end 
end 

%# then we can write 
B = mat2cell(F,3,3,ones(N^2,1)); 
B = reshape(B,N,N)'; 
B = cell2mat(B);