1

Axe、Ay、Az：[N-by-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）：3x3行列。Bを作成する1つの方法は次のとおりです。

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

これは、上記のA（：、：、1）とは異なります[Ax（1,1）Ay（1,1）Az（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.

編集：

私のアプリケーションのためにいくつかの変更を加えた後：gnoviceコードによる

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

ceil、ind2sub ...のような関数呼び出しは、thwループを遅くし、可能であれば回避する必要があるようです。

symIndex面白かったです！ありがとうございました。

4

2 に答える 2

2

これは、単一のforループを使用して線形インデックスを実行し、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、上記のコードの次の修正バージョンを使用できます。

%# 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);

ただし、このような場合は、線形インデックスを作成し、次のように添え字付きのインデックスを使用して2つのforループを使用することで、パフォーマンスが向上する（または少なくともコードが単純になる）場合があります。

%# 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);

編集＃2：

2番目の対称ソリューションは、対角計算を内部ループの外側に移動し（条件ステートメントの必要性を排除）A、結果を上書きして、代わりに（つまり、1つのセルを更新する）を使用しA(:)*Aて対称セルエントリを更新できるようにすることで、さらに高速化できます。別のコンテンツを使用すると、余分なオーバーヘッドが発生するようです）：B{c,r}AB{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);

そして、以下のサンプル対称行列、、、を使用した3つの対称解のタイミング結果を次にAx示しAyます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の大幅な高速化は、あるセルを別のセルの内容で更新するのではなく、主Bにローカル変数でのセルの内容を更新することによってもたらされAます。つまり、セルの内容をなどの操作でコピーすると、B{c,r} = B{r,c};予想以上にオーバーヘッドが発生します。

于 2011-06-06T13:49:29.010 に答える
1 
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);
于 2011-06-05T04:45:09.827 に答える