2012-06-02 22 views
5

Mam macierz z 3 wymiarami Y (i, j, w). Chcę uzyskać wektor determinujący d (w), w którym każda liczba byłaby wyznacznikiem macierzy Y (:,:, w).powracający wektor determinant - Matlab

Czy jest to elegancka składnia, czy po prostu muszę użyć pętli?

dzięki

Odpowiedz

7

Cóż, przede wszystkim, praktycznie nigdy tak naprawdę nie chcemy obliczyć wyznacznik, wystarczy pomyśleć ty. W rzeczywistości prawie nigdy nie jest to dobre, ponieważ determinanty są tak słabo skalowane. Zbyt często są one używane do wnioskowania o statusie osobliwości macierzy, co jest straszne, jeśli chodzi o analizę numeryczną.

Stwierdziwszy mój mini-rant przeciw uwarunkowań w ogóle ...

OPCJA 1:

Konwersja 3-d tablicy do tablicy komórek macierzy kwadratowych, z każdej płaszczyzny tablicy jako jeden komórka. Mat2cell z łatwością i sprawnie załatwi sprawę.

Następnie użyj cellfun w macierzy komórek. cellfun może zastosować funkcję (@det) do każdej komórki, a następnie zwróci wektor determinant. Czy to jest niesamowicie wydajne? Prawdopodobnie nie jest to ogromna korzyść z zastosowania detu w pętli, o ile wcześniej wstępnie przydzielisz wektor podczas pisania pętli.

Opcja 2:

przypadku matryce są małe, zatem powiedzieć, 2x2 lub 3x3 macierzy, wówczas rozwinąć mnożenia dla wyznacznika jako wyraźne przemnożenia wektora. Myślę, że to nie jest jasne, jak piszę, więc dla przypadku 2x2, gdzie Y jest 2x2xn:

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

pewnością widać, że to tworzy wektor determinant 2x2 dla każdej płaszczyzny matrycy Y. Sprawa 3x3 jest na tyle prosta, że ​​można ją napisać również jako sześć trójdrożnych produktów terminów. Nie sprawdziłem dokładnie poniższego przykładu 3x3, ale powinno być blisko.

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,:); 

Jak widać, OPCJA 2 będzie dość szybka i wektoryzowana.

Edytuj: jako odpowiedź dla Chrisa istnieje ZNACZA różnica w wymaganym czasie. Rozważ czas potrzebny na zestaw macierzy 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. 

Te dwie rozmowy zwracają te same wartości w obrębie kosza zmiennoprzecinkowego.

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

Pętla nie będzie lepsza, w rzeczywistości na pewno gorsza. Miałem więc napisać fragment kodu, który miałby za zadanie wytworzyć determinanty dla WIELU takich macierzy w tablicy, w szczególnym przypadku byłby to kod dla macierzy 2x2 i 3x3. Mogę nawet napisać to dla macierzy 4x4. Tak, pisanie jest kłopotliwe, ale jest duża różnica w wymaganym czasie.

Jednym z powodów jest to, że detektor MATLAB używa wywołania do LU, rozkładając macierz na czynniki pierwsze. Jest to lepsze w teorii niż mnożenie dla nawet średnich macierzy, ale dla 2x2 lub 3x3 dodatkowe obciążenie jest mordercze. (Nie będę zgadywał, gdzie mieści się punkt równości szans, ale można to łatwo sprawdzić.)

+0

Preferuję opcję 1 lub nawet pętlę nad opcją 2. Byłbym zszokowany, gdyby utrata prędkości z interpretera w pętli była znacząca i wolałbym bardziej elastyczny kod. Wyobraźmy sobie, że wyznaczamy w ten sposób determinantę 5x5. –

+2

@ChrisA - Czy 12.041883/0.002141 = 5624.4 jest szokującą różnicą? –

+0

+1 Tak! Dobra odpowiedź. Może istnieje bardziej elastyczny sposób na zrobienie tego ... ponieważ robi się naprawdę nieprzyjemnie, ponieważ liczba terminów zwiększa się dzięki 'p!' –

3

I użyłby arrayfun:

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

Edycja: Test prędkości:

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

testowy 1:

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

Test 2 (przy zrębków)

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

testowy 3 (naiwne podejście)

>> 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. 

Wniosek: naiwne podejście jest najszybszy.

+0

Czy możesz użyć testu w odpowiedzi @woodchips, aby dokonać porównania prędkości między tą składnią a opcjami cellfun i special-case? – tmpearce

+0

@tmpearce, bez problemu, dodałem wynik testu – Serg

+0

Wyniki będą zależały od wielkości 'p': Przy p = 2, sposób @woodchips jest prawdopodobnie dość szybki. Interesujące jest to, że przy p = 10 pętla jest szybsza niż arrayfun. – tmpearce

Powiązane problemy