4

このコードは、2次元行列値関数の1次元畳み込みを実行するために作成しました（kは私の時間インデックス、kendは10e3のオーダーです）。おそらく組み込み関数を使用して、これを行うためのより高速またはよりクリーンな方法はありますか？

for k=1:kend
  C(:,:,k)=zeros(3);
  for l=0:k-1
      C(:,:,k)=C(:,:,k)+A(:,:,k-l)*B(:,:,l+1);
  end
end
4

3 に答える 3

4

新しいソリューション：

これは、以前に与えられた式を解決した古いソリューションに基づいて構築された新しいソリューションです。問題のコードは、実際にはその式を修正したものであり、3次元の2つの行列間のオーバーラップが繰り返しシフトされます（データの3次元に沿った畳み込みに似ています）。私が与えた前の解決策は、質問のコードの最後の反復の結果のみを計算しました（つまりk = kend）。したがって、これが1000のオーダーの質問のコードよりもはるかに効率的であるはずの完全なソリューションです。kend

kend = size(A,3);                         %# Get the value for kend
C = zeros(3,3,kend);                      %# Preallocate the output
Anew = reshape(flipdim(A,3),3,[]);        %# Reshape A into a 3-by-3*kend matrix
Bnew = reshape(permute(B,[1 3 2]),[],3);  %# Reshape B into a 3*kend-by-3 matrix
for k = 1:kend
  C(:,:,k) = Anew(:,3*(kend-k)+1:end)*Bnew(1:3*k,:);  %# Index Anew and Bnew so
end                                                   %#   they overlap in steps
                                                      %#   of three

を使用した場合でもkend = 100、このソリューションは、質問のソリューションよりも約30倍高速であり、純粋なforループベースのソリューション（5つのループが含まれる！）よりも約4倍高速であることがわかりました。以下の浮動小数点精度の説明が引き続き適用されることに注意してください。したがって、相対的な浮動小数点精度のオーダーでソリューション間にわずかな違いが見られるのは正常であり、予想されます。

古いソリューション：

コメントでリンクしたこの式に基づいて：

ここに画像の説明を入力してください

質問で提供したコードとは異なることを実際に実行したいようです。が3x3 x kの行列であるとすると、結果Aは3 x 3の行列になり、ネストされたforループのセットとして書き出されたリンクの数式は次のようになります。BC

%# Solution #1: for loops
k = size(A,3);
C = zeros(3);
for i = 1:3
  for j = 1:3
    for r = 1:3
      for l = 0:k-1
        C(i,j) = C(i,j) + A(i,r,k-l)*B(r,j,l+1);
      end
    end
  end
end

これで、再形成と再編成を適切に行うことで、forループなしでこの操作を実行できます。AB

%# Solution #2: matrix multiply
Anew = reshape(flipdim(A,3),3,[]);        %# Create a 3-by-3*k matrix
Bnew = reshape(permute(B,[1 3 2]),[],3);  %# Create a 3*k-by-3 matrix
C = Anew*Bnew;                            %# Perform a single matrix multiply

質問にあるコードを作り直して、3行3列の行列の行列乗算を実行する単一のループを持つソリューションを作成することもできます。

%# Solution #3: mixed (loop and matrix multiplication)
k = size(A,3);
C = zeros(3);
for l = 0:k-1
  C = C + A(:,:,k-l)*B(:,:,l+1);
end

さて、質問：これらのアプローチのどれがより速く/よりクリーンですか？

ええと、「クリーナー」は非常に主観的であり、私は正直に言って、上記のコードのどれが操作が何をしているのかを理解しやすくするのかをあなたに伝えることができませんでした。最初のソリューションのすべてのループと変数により、何が起こっているのかを追跡するのが少し難しくなりますが、それは明らかに式を反映しています。2番目のソリューションは、すべてを単純な行列演算に分解しますが、元の数式とどのように関連しているかを確認するのは困難です。3番目の解決策は、2つの中間点のように見えます。

それでは、スピードをタイブレーカーにしましょう。のいくつかの値に対して上記のソリューションの時間をk計ると、次の結果が得られます（与えられたソリューション、MATLAB R2010bの10,000回の反復を実行するのに必要な秒数）。

 k   |  loop  | matrix multiply |  mixed
-----+--------+-----------------+--------
 5   | 0.0915 |      0.3242     | 0.1657
 10  | 0.1094 |      0.3093     | 0.2981
 20  | 0.1674 |      0.3301     | 0.5838
 50  | 0.3181 |      0.3737     | 1.3585
 100 | 0.5800 |      0.4131     | 2.7311   * The matrix multiply is now fastest
 200 | 1.2859 |      0.5538     | 5.9280

さて、小さい値k（約50以下）の場合、forループソリューションが実際に勝ちます。これは、forループが以前のバージョンのMATLABで考慮されていたほど「悪」ではないことをもう一度示しています。特定の状況下では、巧妙なベクトル化よりも効率的です。ただし、の値がk約100より大きい場合、ベクトル化された行列乗算ソリューションが優先され始めk、forループソリューションよりも増加に伴ってはるかに適切にスケーリングされます。for-loop / matrix-multiplyの混合ソリューションは、私が正確に確信していない理由で、ひどくスケーリングします。

したがって、k大きいと予想される場合は、ベクトル化された行列乗算ソリューションを使用します。覚えておくべきことの1つは、各ソリューションで実行される加算と乗算の順序が異なるため、各ソリューション（行列）から得られる結果が（浮動小数点の精度Cのレベルで）わずかに異なることです。丸め誤差の累積に違いが生じます。つまり、これらのソリューションの結果の違いはごくわずかであるはずですが、注意する必要があります。

于 2011-04-25T15:15:09.443 に答える
2

Matlabのconvメソッドを調べましたか？

k=1Aの0番目の要素にアクセスしようとすると問題が発生するため、提供されたコードと比較することはできませんk-1=0。

于 2011-04-23T16:38:13.117 に答える
1

FFTを使用して畳み込みを行うことを検討しましたか？畳み込み演算は、周波数領域でのポイントごとの乗算です。注意しないと巡回畳み込みが発生するため、有限シーケンスでは注意が必要です（ただし、これは簡単です）。

これは、1Dケースの簡単な例です。

>> a=rand(4,1);
>> b=rand(3,1);
>> c=conv(a,b)

c =

    0.1167
    0.3133
    0.4024
    0.5023
    0.6454
    0.3511

FFTを使用して同じ

>> A=fft(a,6);
>> B=fft(b,6);
>> C=real(ifft(A.*B))

C =

    0.1167
    0.3133
    0.4024
    0.5023
    0.6454
    0.3511

ポイントベクトルとポイントベクトルの畳み込みは、MポイントベクトルになりNますM+N-1。したがって、FFTを取得する前に、各ベクトルaとゼロをパディングしました（これは、ポイントFFTbを取得するときに自動的に処理されます）。4+3-1=6

編集

あなたが示した方程式は巡回畳み込みに似ていますが、それは正確ではありません。conv*したがって、FFTアプローチと組み込み関数を捨てることができます。あなたの質問に答えるために、これは明示的なループなしで行われる同じ操作です：

dim1=3;dim2=dim1;
dim3=10;
a=rand(dim1,dim2,dim3);
b=rand(dim1,dim2,dim3);


mIndx=cellfun(@(x)(1:x),num2cell(1:dim3),'UniformOutput',false);
fun=@(x)sum(reshape(cell2mat(cellfun(@(y,z)a(:,:,y)*b(:,:,z),num2cell(x),num2cell(fliplr(x)),'UniformOutput',false)),[dim1,dim2,max(x)]),3);
c=reshape(cell2mat(cellfun(@(x)fun(x),mIndx,'UniformOutput',false)),[dim1,dim2,dim3]);
  • mIndxこれがセルで、ithセルにはベクトルが含まれてい1:iます。これはあなたのインデックスです（他の人が指摘しているように、変数名としてl使用しないでください）。l
  • 次の行は、畳み込み演算を実行する無名関数であり、インデックスが反転されkたインデックスであるという事実を利用しています。l操作は個々のセルで実行され、次に組み立てられます。
  • 最後の行は、実際に行列の演算を実行します。

答えは、ループで得られたものと同じです。ただし、ループされたソリューションは実際には1桁高速であることがわかります（コードの平均は0.007秒、ループの平均は0.0006秒です）。これは、ループが非常に単純であるのに対し、この種のネストされた構造では、関数呼び出しのオーバーヘッドが多く、再形成が繰り返されるため、ループが遅くなるためです。

MATLABのループは、ループが恐ろしい初期の頃から長い道のりを歩んできました。確かに、ベクトル化された操作は非常に高速です。ただし、すべてをベクトル化できるわけではなく、ループはそのような複雑な無名関数よりも効率的である場合があります。構造を最適化することで（または別のアプローチをとることで）、あちこちで10分の数を削ることができるかもしれませんが、そうするつもりはありません。

優れたコードは読みやすくなければならず、読みやすさを犠牲にして効率的でマイナーな最適化は誰にも役立たないことを忘れないでください。私は上記のコードを書きましたが、1か月後に再訪した場合、それが何をするのかを確実に解読することはできません。ループされたコードは明確で、読みやすく、高速でした。それを堅持することをお勧めします。

于 2011-04-23T19:42:48.940 に答える