6

私は3つの次元Y（i、j、w）を持つ行列を持っています。各数値が行列 Y(:,:,w) の行列式となる行列式ベクトル d(w) を取得したいと考えています。

それにはエレガントな構文がありますか、それともループを使用する必要がありますか?

ありがとう

4

2 に答える 2

8

まず第一に、行列式を計算したいと思うことはほとんどありません。実際、決定要因のスケーリングが不十分なため、これが良いことになることはほとんどありません。行列の特異点状態を推測するために使用されることが多すぎますが、これは数値解析の観点からはひどいことです。

一般的な決定要因に対する私の小さな暴言を述べました...

オプション1：

3 次元配列を、配列の各平面を 1 つのセルとして、正方行列のセル配列に変換します。mat2cell は、このトリックを簡単かつ効率的に実行します。

次に、cell 配列に対して cellfun を使用します。cellfun は関数 (@det) をすべてのセルに適用でき、行列式のベクトルを返します。これは信じられないほど効率的ですか？ループを作成するときにベクトルを事前に割り当てている限り、ループで det を適用することよりもおそらく大きなメリットはありません。

オプション 2:

行列が小さい場合、たとえば 2x2 または 3x3 行列の場合、行列式の乗算を明示的なベクトルの乗算として展開します。私が書いているので、これは明確ではないと思うので、Y が 2x2xn である 2x2 の場合:

d = Y(1,1,:).*Y(2,2,:) - Y(1,2,:).*Y(2,1,:);

確かに、これは行列 Y のすべての平面に対して 2x2 の行列式のベクトルを形成することがわかります。3x3 の場合は、項の 6 つの 3 元積としても簡単に記述できます。以下の 3x3 のケースは注意深くチェックしていませんが、近いはずです。

d = Y(1,1,:).*Y(2,2,:).*Y(3,3,:) + ...
    Y(2,1,:).*Y(3,2,:).*Y(1,3,:) + ...
    Y(3,1,:).*Y(1,2,:).*Y(2,3,:) - ...
    Y(3,1,:).*Y(2,2,:).*Y(1,3,:) - ...
    Y(2,1,:).*Y(1,2,:).*Y(3,3,:) - ...
    Y(1,1,:).*Y(3,2,:).*Y(2,3,:);

ご覧のとおり、OPTION 2 は非常に高速で、ベクトル化されています。

編集: Chris への応答として、所要時間に大きな違いがあります。1e5 行列のセットに必要な時間を考えてみましょう。

p = 2;
n = 1e5;
Y = rand(p,p,n);

tic,
d0 = squeeze(Y(1,1,:).*Y(2,2,:) - Y(2,1,:).*Y(1,2,:));
toc

Elapsed time is 0.002141 seconds.

tic,
X = squeeze(mat2cell(Y,p,p,ones(1,n)));
d1= cellfun(@det,X);
toc

Elapsed time is 12.041883 seconds.

2 つの呼び出しは、同じ値を浮動小数点トラッシュ内に返します。

std(d0-d1)
ans =
   3.8312e-17

ループは良くはありません。実際には、確実に悪くなります。したがって、配列内の多数のそのような行列の行列式を生成するタスクに直面するコードを書くとしたら、2x2 および 3x3 行列のコードを特別なケースにします。4x4 行列用に書き出すことさえあるかもしれません。はい、書き出すのは面倒ですが、所要時間に大きな違いがあります。

理由の 1 つは、MATLAB の det が LU の呼び出しを使用して、行列を因数分解することです。これは理論的には中規模の行列の乗算よりも優れていますが、2x2 または 3x3 の場合、余分なオーバーヘッドが致命的です。(損益分岐点がどこに落ちるかはわかりませんが、簡単にテストできます。)

于 2012-06-02T11:43:10.917 に答える
4

私はarrayfunを使用します：

d = arrayfun(@(w) det(Y(:, :, w)), 1 : size(Y, 3));

編集：スピードテスト：

p = 10;
n = 1e4;
Y = rand(p,p,n);

テスト 1:

>> tic, d1 = arrayfun(@(w) det(Y(:, :, w)), 1 : size(Y, 3)); toc
Elapsed time is 0.139030 seconds.

テスト 2 (木材チップによる):

>> tic, X = squeeze(mat2cell(Y,p,p,ones(1,n))); d2= cellfun(@det,X); toc
Elapsed time is 1.318396 seconds.

テスト 3 (素朴なアプローチ):

>> p = 10;
>> n = 1e4;
>> Y = rand(p,p,n);
>> tic; d = nan(n, 1); for w = 1 : length(d), d(w) = det(Y(:, :, w)); end; toc
Elapsed time is 0.069279 seconds.

結論: 素朴なアプローチが最速です。

于 2012-06-03T15:34:37.263 に答える