MATLABでループなしの行列計算

Question

いくつかの配列演算を実行するコードに問題があります。ループを使用し、入力データが非常に大きいので、遅すぎます。それは私のための最も簡単な方法だったが、今は私はループよりも速いものを探している。私はコードを最適化または書き直そうとしていましたが、失敗しました。私はあなたの助けを本当に寄付します。私のコードで

私は3つのアレイx1、y1（グリッド内の点の座標）、g1（ポイント値）と、例えばそれらのサイズは300×300であり、Iは、9の組成物として各マトリックスを治療及びIが作るを有します中央の点の計算。たとえば、g1(101,101)で始まりますが、私はg1(1:201,1:201)=g2のデータを使用しています。 g1(1:201,1:201)からg1(101,101)（ll行列）の各点からの距離を計算する必要があります。次に、コード内のnnを計算し、次にg1(101,101)の値をnnから見つけ、N配列に入れます。その後、私はg1(101,102)に行き、g1(200,200)まで続きます。この最後の場合はg2=g1(99:300,99:300)です。

私が言ったように、このコードはあまり効率的ではありません。たとえ私がこの例で示したよりも大きな配列を使用しなければならなくても、あまり時間がかかります。私がコードから期待していることを十分に説明することを願っています。私はarrayfunを使うことを考えていましたが、この関数を使ったことは一度もありませんでしたので、どのように使うべきかわかりませんが、それは私には分かりません。たぶん、他の解決策があるかもしれませんが、私は適切なものを見つけることができませんでした。何が価値があるために

tic 
x1=randn(300,300); 
y1=randn(300,300); 
g1=randn(300,300); 
m=size(g1,1); 
n=size(g1,2); 
w=1/3*m; 
k=1/3*n; 
N=zeros(w,k); 
for i=w+1:2*w 
    for j=k+1:2*k 
     x=x1(i,j); 
     y=y1(i,j); 
     x2=y1(i-k:i+k,j-w:j+w); 
     y2=y1(i-k:i+k,j-w:j+w); 
     g2=g1(i-k:i+k,j-w:j+w); 
     ll=1./sqrt((x2-x).^2+(y2-y).^2); 
     ll(isinf(ll))=0; 
     nn=ifft2(fft2(g2).*fft2(ll)); 
     N(i-w,j-k)=nn(w+1,k+1); 
     end 
    end 
    czas=toc; 

Answer 1

、arrayfun()は、forループのための単なるラッパーなので、パフォーマンスの向上につながりません。また、おそらくx2の定義にタイプミスがあります。それはx1に依存していると仮定します。さもなければ余分な変数になります。また、i<->w/k、j<->k/wのペア設定が矛盾しているように見える場合は、それも確認する必要があります。 でも、ちょうどtic/tocでタイミングすることはめったに正確ではありません。あなたのコードをプロファイリングするときに、それを関数に入れて、そのタイミングを何度も実行し、そのタイミングから変数生成を除外します。さらに優れています：組み込みプロファイラを使用します。

免責事項：このソリューションは、膨大なメモリが必要となるため、実際の問題には役立たない可能性があります。 300x300行列を入力する場合は、サイズが300x300x100x100の配列で動作しますが、通常はノー・ゴーです。それでも、入力サイズが小さいほど参考になります。私はnlfilter（）をベースにしたソリューションを追加したいと思っていましたが、あなたの問題はそれを使用するには複雑すぎるようです。

いつもベクトル化では、メモリをスペアにすることができれば、より速く行うことができます。各[i,j]インデックスに対して、サイズ[2*k+1,2*w+1]の行列で作業しようとしています。これは、形状が[2*k+1,2*w+1,w,k]の4d配列を必要とします。各要素について[i,j]の場合は、x1とy1という対応する要素とともに扱うために、インデックスが[:,:,i,j]の行列があります。また、fft2が多次元配列を受け入れるのに役立ちます。ここで

は私の意味は次のとおりです。

tic 
x1 = randn(30,30); %// smaller input for tractability 
y1 = randn(30,30); 
g1 = randn(30,30); 
m = size(g1,1); 
n = size(g1,2); 
w = 1/3*m; 
k = 1/3*n; 

%// these will be indexed on the fly:  
%//x = x1(w+1:2*w,k+1:2*k);  %// size [w,k] 
%//y = x1(w+1:2*w,k+1:2*k);  %// size [w,k] 

x2 = zeros(2*k+1,2*w+1,w,k); %// size [2*k+1,2*w+1,w,k] 
y2 = zeros(2*k+1,2*w+1,w,k); %// size [2*k+1,2*w+1,w,k] 
g2 = zeros(2*k+1,2*w+1,w,k); %// size [2*k+1,2*w+1,w,k] 

%// manual definition for now, maybe could be done smarter: 
for ii=w+1:2*w  %// don't use i and j as variables 
    for jj=k+1:2*k %// don't use i and j as variables 
     x2(:,:,ii-w,jj-k) = x1(ii-k:ii+k,jj-w:jj+w); %// check w vs k here 
     y2(:,:,ii-w,jj-k) = y1(ii-k:ii+k,jj-w:jj+w); %// check w vs k here 
     g2(:,:,ii-w,jj-k) = g1(ii-k:ii+k,jj-w:jj+w); %// check w vs k here 
    end 
end 

%// use bsxfun to operate on [2*k+1,2*w+1,w,k] vs [w,k]-sized arrays 
%// need to introduce leading singletons with permute() in the latter 
%// in order to have shape [1,1,w,k] compatible with the first array 
ll = 1./sqrt(bsxfun(@minus,x2,permute(x1(w+1:2*w,k+1:2*k),[3,4,1,2])).^2 ... 
      + bsxfun(@minus,y2,permute(y1(w+1:2*w,k+1:2*k),[3,4,1,2])).^2); 
ll(isinf(ll)) = 0; 

%// compute fft2, operating on [2*k+1,2*w+1,w,k] 
%// will return fft2 for each index in the [w,k] subspace 
nn = ifft2(fft2(g2).*fft2(ll)); 

%// we need nn(w+1,k+1,:,:) which is exactly of size [w,k] as needed 
N = reshape(nn(w+1,k+1,:,:),[w,k]); %// quicker than squeeze() 
N = real(N); %// this solution leaves an imaginary part of around 1e-12 

czas=toc;