2012-06-02 23 views
5

我有一個3維矩陣Y(i,j,w)。 我想得到一個行列式向量d(w),其中每個數字都是矩陣行列式的行列式(:,:,w)。返回行列式向量 - Matlab

它是否有一個優雅的語法,或者我只需要使用循環?

感謝

回答

7

嗯,首先,你幾乎從來沒有真正要計算一個決定因素,你只是覺得你這樣做。事實上,它幾乎從來都不是一件好事,因爲行列式的縮放比例很小。通常他們被用來推斷矩陣的奇點狀態,這在數值分析方面是一件可怕的事情。

在闡述對一般決定我的小咆哮......

選項1:

轉換您的3-d陣列成方陣的單元陣列,與陣列爲一體的每個平面細胞。 mat2cell將輕鬆高效地完成這項技巧。

接下來,在cell陣列上使用cellfun。 cellfun可以爲每個單元格應用一個函數(@det),然後返回一個行列式向量。這是非常有效的嗎?在循環中應用det可能不是一個巨大的收益,只要您在編寫循環時提前預先分配vector。

選項2:

如果基質是小的,因此說2×2 3×3或基質,然後展開出來的行列式作爲明確的矢量相乘的乘法。我覺得這是我寫的,所以對於一個2×2的情況下,其中Y是2x2xn不明確:

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

當然,你看,這形成2x2的決定因素的向量矩陣Y的每一個面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,:); 

正如你所看到的,選項2會很快,而且是矢量化的。

編輯:作爲對克里斯的迴應,所需時間存在顯着差異。考慮一組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. 

這兩個調用將相同的值返回到浮點垃圾內。

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

循環不會更好,事實上肯定會更糟。所以我要寫一段代碼來面對在數組中產生許多這樣的矩陣的決定因素的任務,我會特別說明2x2和3x3矩陣的代碼。我甚至可以寫出4x4矩陣。是的,寫出來很麻煩,但所需時間有很大差異。

一個原因是MATLAB的det使用了一個調用LU的因子分解矩陣。理論上這比中等大矩陣的乘法更好,但對於2x2或3x3,額外的開銷是一個殺手。 (我不會猜測盈虧平衡點落在哪裏,但可以輕鬆測試。)

+0

我更喜歡選項1或者甚至選項2的循環。如果循環上解釋器的速度損失是重要的,我寧願使用更靈活的代碼。想象一下,用這種方式做一個5x5的行列式 –

+2

@ChrisA - 是12.041883/0.002141 = 5624.4一個令人震驚的區別? – 2012-06-02 17:08:16

+0

+1是!好的迴應。也許有一個更靈活的方法來做到這一點...因爲它使得術語的數量變得更大,所以'p!' –

3

我會使用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. 

結論:幼稚的做法是最快的。

+0

您可以使用@woodchips答案中的測試用例來比較此語法和cellfun和特殊情況選項之間的速度嗎? – tmpearce

+0

@tmpearce,沒問題,我添加了測試結果 – Serg

+0

結果將取決於'p'的大小:在p = 2的情況下,@ woodchips方式可能相當快。有趣的是,儘管p = 10,但循環比arrayfun快。 – tmpearce