14

私は初めての大規模な MATLAB プログラムを構築しています。平射図法でベクトル密度を表す画像を作成しようとするまで、すべての元のベクトル化されたコードを書くことができました。数回失敗した後、私は Mathworks ファイル交換サイトにアクセスし、Malcolm Mcleanの厚意により、私のニーズに合ったオープン ソース プログラムを見つけました。テストマトリックスを使用すると、彼の関数は次のようなものを生成します。

平射図法におけるベクトル密度

これは私が望んでいたこととほとんど同じですが、彼のコードは 3 重にネストされた for ループに依存しています。私のワークステーションでは、サイズ 25000x2 のテスト データ マトリックスは、コードのこのセクションで 65 秒かかりました。プロジェクトでサイズ 500000x2 のデータ マトリックスにスケールアップするため、これは受け入れられません。

これまでのところ、最も内側のループ (最長/最悪のループ) をベクトル化できましたが、続行して、可能であればループを完全に取り除きたいと考えています。ベクトル化する必要がある Malcolm の元のコードは次のとおりです。

dmap = zeros(height, width); % height, width: scalar with default value = 32
for ii = 0: height - 1          % 32 iterations of this loop
    yi = limits(3) + ii * deltay + deltay/2; % limits(3) & deltay: scalars
    for jj = 0 : width - 1      % 32 iterations of this loop
        xi = limits(1) + jj * deltax + deltax/2; % limits(1) & deltax: scalars
        dd = 0;
        for kk = 1: length(x)   % up to 500,000 iterations in this loop
            dist2 = (x(kk) - xi)^2 + (y(kk) - yi)^2;
            dd = dd + 1 / ( dist2 + fudge); % fudge is a scalar
        end
        dmap(ii+1,jj+1) = dd;
    end
end

ここでは、最も内側のループに加えた変更を示します (これが効率の最大の浪費でした)。これにより、私のマシンでは同じテスト マトリックスの時間が 65 秒から 12 秒に短縮されました。 

     dmap = zeros(height, width);
    for ii = 0: height - 1
        yi = limits(3) + ii * deltay + deltay/2;
        for jj = 0 : width - 1
            xi = limits(1) + jj * deltax + deltax/2;
            dist2 = (x - xi) .^ 2 + (y - yi) .^ 2;
            dmap(ii + 1, jj + 1) = sum(1 ./ (dist2 + fudge));
        end
    end
私の主な質問は、このコードを最適化するためにできる変更はありますか? または、問題にアプローチする別の方法でさえありますか? プログラムのこのセクションでは、MATLAB の代わりに C++ または F# を使用することを検討しました。MATLAB コードで妥当な効率レベルに達しない場合は、そうするかもしれません。

また、この時点で追加のツールボックスを持っていないことにも注意してください。もしそうなら、これは些細なことであることがわかります（たとえば、統計ツールボックスの hist3 を使用します）。

4

3 に答える 3

1

あるいは、この問題は、カーネル密度推定技術を使用して解決できます。これは Statistics Toolbox の一部です。または、Zdravko Botev による KDE 実装があります (ツールボックスは必要ありません)。

以下のコード例では、N = 500000 の場合は 0.3 秒、N = 1000000 の場合は 0.7 秒になります。

N = 500000;
data = [randn(N,2); rand(N,1)+3.5, randn(N,1);];  % 2 overlaid distrib
tic; [bandwidth,density,X,Y] = kde2d(data); toc;
imagesc(density);

分布密度の出力例

于 2013-05-08T23:17:15.417 に答える