Obliczanie wyznacznika podczas rozwiązywania za pomocą CG

11

Rozwiązuję dla ogromnej rzadkiej dodatniej określonej macierzy za pomocą metody gradientu sprzężonego (CG). Czy można obliczyć wyznacznik podstawie informacji uzyskanych podczas rozwiązania? $Ax=b$ $A$ $A$

— Manuel Schmidt
źródło

Dlaczego chcesz obliczyć wyznacznik? Taki wynik z pewnością będzie albo niedomiarem, albo przepełnieniem ogromnej matrycy. Byłbym bardziej charytatywny, gdybyś poprosił o obliczenie numeru warunku, ale nie marnuj czasu na wyznacznik!

Prawdopodobnie już to wiesz, ale wartości Ritza podczas procesu gradientu sprzężonego są zbieżne z wartościami własnymi macierzy i na tej podstawie można uzyskać proste oszacowania dla wyznacznika.

— shuhalo,

10

Obliczanie wyznacznika rzadkiej macierzy jest zazwyczaj tak samo drogie jak bezpośrednie rozwiązanie, i jestem sceptyczny, że CG byłby bardzo pomocny w jej obliczeniu. Byłoby możliwe uruchomienie CG dla powtórzeń (gdzie jest ) w celu generowania informacji na całym spektrum , a następnie obliczyć determinantę jako iloczyn wartości własnych, ale byłoby to zarówno wolno i niestabilny numerycznie. $n$ $A$ $n \times n$ $A$

Byłoby lepszym pomysłem obliczyć rzadką bezpośrednią faktoryzację Choleskiego macierzy, powiedzmy , gdzie jest dolnym trójkątem. Następnie gdzie jest po prostu iloczynem przekątnych wejść dolnej trójkątnej matrycy ponieważ wartości własne macierzy trójkątnej leżą wzdłuż jej przekątnej. $A = L L^H$ $L$

det (ZA) = det (L.) det ({L.}^{H.}) = | det (L.) |^{2)},

$\text{det}(A) = \text{det}(L) \text{det}(L^H) = |\text{det}(L)|^2,$

det (L)

$\text{det}(L)$

L

$L$

W przypadku ogólnej macierzy niepodzielnej należy zastosować obrotowy układ LU, powiedzmy , gdzie jest macierzą permutacji, tak aby Ponieważ jest macierzą permutacji, , a ze względu na konstrukcję będzie zazwyczaj mieć przekątną wszystkich, co implikuje, że . W ten sposób można obliczyć jako $PA=LU$ $P$

det (ZA) = det ({P.}^{- 1}) \cdot det (L.) \cdot det (U) .

$\text{det}(A) = \text{det}(P^{-1}) \cdot \text{det}(L) \cdot \text{det}(U).$

P

$P$

det (P) = \pm 1

$\text{det}(P)=\pm 1$

L

$L$

det (L) = 1

$\text{det}(L)=1$

det (A)

$\text{det}(A)$

\pm det (U)

$\pm \text{det}(U)$ i ponownie uznają, że wyznacznik macierzy trójkątnej jest po prostu iloczynem jej przekątnych wejść. Tak więc koszt obliczenia wyznacznika jest zasadniczo tylko kosztem faktoryzacji.

— Jack Poulson
źródło

Byłaby to jedna z możliwości (choć użyłbym podziału na czynniki chłodnicze), gdyby matryca była niewielka, jednak

ma rozmiar ~

i dlatego nie można dokonać rozkładu

A

$A$

10^{6} x 10^{6}

$10^6 x 10^6$

— Manuel Schmidt

@ManuelSchmidt Rzadkie macierze o takim rozmiarze wynikające z dyskretyzacji typu elementów skończonych można zwykle łatwo rozłożyć na czynniki (na przykład) metodami wielopłaszczyznowymi. Zgadzam się, że należy zastosować faktoryzację Cholesky'ego, jeśli macierz ma HPD (a uogólnienie mojego powyższego argumentu jest oczywiste).

— Jack Poulson,

Dzięki za szybką odpowiedź i odpowiedź. Niestety matryca nie ma struktury spezialnej (co pozwoliłoby na łatwą faktoryzację).

— Manuel Schmidt,

2

Jestem ciekawy, dlaczego musisz obliczyć wyznacznik macierzy. Czy najwyższe i najniższe wartości własne są niewystarczające?

— Jack Poulson,

Jest to część złożonej funkcji rozkładu prawdopodobieństwa, a nie tylko stała normalizacyjna. Wiem, że rozkłady można uwzględnić (i to właśnie robimy w tej chwili), jednak mamy mnóstwo danych do modelowania i każdy z czynników staje się ogromny.

— Manuel Schmidt,

6

$A$ $B$ $\textrm{dim}A \ll \textrm{dim} B$ $\textrm{dim}B=\infty$

$B$ $A$ $B$ $A$ $B$ $\textrm{det}A$ $\textrm{det}B$ $A$ $B$

det ZA = \prod_{jot = 1 \dots ciemny ZA} λ_{ja} (ZA) \approx \prod_{jot = 1 \dots ciemny ZA} λ_{ja} (b) \prod_{jot = ciemny ZA + 1 \dots ciemny b} λ_{ja} (b)

$\textrm{det}A = \prod_{j=1\ldots \textrm{dim}A} \lambda_i(A) \approx \prod_{j=1\ldots \textrm{dim}A} \lambda_i(B) \prod_{j=\textrm{dim}A+1\ldots \textrm{dim}B} \lambda_i(B)$

B

$B$

A

$A$

dim B = \infty

$\textrm{dim}B=\infty$

det A

$\textrm{det}A$

det B

$\textrm{det}B$

— Wolfgang Bangerth
źródło

Okazuje się, że istnieją nasze naprawdę piękne i praktyczne algorytmy, które wymagają obliczenia sporej determinanty. Sprawdź www-m3.ma.tum.de/foswiki/pub/M3/Allgemeines/…

— Matt Knepley

2

Nie wchodząc (ponownie) w pytanie, dlaczego i jak determinanty są złe, załóżmy, że twój operator albo nie jest łatwo podzielny na czynniki, albo po prostu w ogóle nie jest dostępny jako matryca i że naprawdę musisz oszacować jego determinantę.

$A$ $A$

Prawdopodobnie możesz dokonać inżynierii wstecznej, w jaki sposób powstaje to oszacowanie wyznacznika w standardowej implementacji CG, postępując zgodnie z rozdziałem 6.7.3 książki.

— Dominique
źródło

2

re mi t (ZA) = \prod_{ja = 1}^{n} α_{k}^{- 1},

$det(A) = \prod_{i = 1}^n \alpha_k^{-1},$

α_{k} = \frac{r_{k}^{T} r_{k}}{p_{k}^{T} A p_{k}}

$\alpha_k = \cfrac{r_k^Tr_k}{p_k^TAp_k}$

r_{k} \neq 0 \forall k = 1, \dots, n

$r_k \ne 0 \;\,\forall k=1, \ldots, n$

R

$R$

r_{k}

$r_k$

P

$P$

p_{k}

$p_k$

p_{k} = - r_{k} + \sum_{ja = 1}^{k - 1} γ_{ja} r_{ja} .

$p_k = -r_k + \sum_{i=1}^{k-1}\gamma_i r_i.$

d e t (P) = (- 1)^{n} d e t (R)

$det(P) = (-1)^n det(R)$

r_{k}

$r_k$

p_{k}

$p_k$

A

$A$

\prod_{k = 1}^{n} α_{k} = \prod_{k = 1}^{n} \frac{r_{k}^{T.} r_{k}}{p_{k}^{T.} ZA p_{k}} = \frac{re mi t (R^{T.} R)}{re mi t ({P.}^{T.} ZA P.)} = \frac{re mi t (R^{T.} R)}{re mi t (ZA) re mi t ({P.}^{T.} P.)} = (re mi t (ZA))^{- 1} .

$\prod_{k = 1}^n \alpha_k = \prod_{k = 1}^n \frac{r_k^T r_k}{p_k^TAp_k} = \frac{det(R^TR)}{det(P^TAP)} = \frac{det(R^TR)}{det(A)det(P^TP)}= (det(A))^{-1}.$

— Jermek
źródło