Próbkowanie z wielowymiarowego Gaussa z grafem kowariancji Laplaciana (odwrotna)

Wiemy np. Z Koutis-Miller-Peng (na podstawie pracy Spielmana i Tenga), że możemy bardzo szybko rozwiązać układy liniowe dla macierzy które są wykresem macierzy Laplaciana dla niektórych rzadkich wykresów z nieujemnymi wagami krawędzi . $A x = b$ $A$

Teraz (pierwsze pytanie) rozważ użycie jednej z tych grafów macierzy Laplaciana jako kowariancji lub (drugie pytanie) odwrotnej macierzy kowariancji o zerowym średnim wielowymiarowym rozkładzie normalnym lub . W każdym z tych przypadków mam dwa pytania: $A$ $\mathcal{N}(\boldsymbol{0}, A)$ $\mathcal{N}(\boldsymbol{0}, A^{-1})$

A. Jak skutecznie możemy pobrać próbkę z tej dystrybucji? (Zazwyczaj w celu narysowania próbki obliczamy rozkład Cholesky'ego $A = LL^T$ , rysujemy standardowy normalny $y \sim \mathcal{N}(\boldsymbol{0}, I)$ , a następnie obliczamy próbkę jako $x = L^{-1} y$ ).

B. Jak skutecznie możemy obliczyć wyznacznik $A$ ?

Zauważ, że oba z nich można łatwo rozwiązać, biorąc pod uwagę rozkład Cholesky'ego, ale nie widzę od razu, jak wyodrębnić $L$ bardziej wydajnie niż po prostu za pomocą standardowego rzadkiego algorytmu Cholesky'ego, który nie użyłby technik przedstawionych w wyżej wspomnianym działa i który miałby złożoność sześcienną dla wykresów rzadkich, ale wysokich.

— dan_x
źródło

Myślę, że może być nieco bardziej szczegółowe określenie tego, co w obu przypadkach można uznać za „wydajne”. Czy „efektywny” jest tym samym, co „nie zależy od rozkładu Choleskiego”?

— Suresh Venkat

Dzieki za sugestie. Możliwe, że odpowiedź na wszystkie pytania brzmi: „musisz obliczyć rozkład Choleskiego i nie ma struktury, którą można by wykorzystać poza rzadkością matrycy”. Chciałbym wiedzieć, czy to prawda (ale mam nadzieję, że tak nie jest). W odniesieniu do „wydajnie” w ostatnim akapicie, tak, mam na myśli przede wszystkim bardziej efektywnie niż standardowe rzadkie algorytmy Choleskiego. Chociaż gdyby istniał sposób wykorzystania technik powyższej pracy, aby obliczyć Cholesky równie szybko, jak można to zrobić za pomocą innych środków, byłoby to również interesujące.

— dan_x

Jeśli chcesz próbkować z , możesz użyć tego , gdzie jest macierzą padania wykresu. W ten sposób można próbki ze standardowej krzywej Gaussa na ( są krawędziami) i stosuje się do transformacji liniowej . Nie wiem, jak to porównać z poniższymi sugestiami, ale nie trzeba obliczać rozkładu Choleskiego.

N (0, A)

$N(0,A)$

A = B^{T} B

$A = B^T B$

B

$B$

R^{E}

$\mathbb{R}^E$

E

$E$

B

$B$

— Lorenzo Najt,

Istnieją tutaj dwa osobne problemy.

Jak stosować wydajne solwery dla , aby zastosować . $Ax=b$ $A^{1/2}b$
Jak obliczyć wyznacznik.

Krótkie odpowiedzi to 1) użyj przybliżeń funkcji racjonalnej macierzy i 2) nie musisz, ale i tak nie musisz. Obie te kwestie omawiam poniżej.

Przybliżenia pierwiastka kwadratowego macierzy

Chodzi tutaj o konwersję aproksymacji funkcji wymiernej dla funkcji skalarnych na aproksymację funkcji wymiernej dla funkcji macierzowych.

Wiemy, że istnieją funkcje wymierne, które mogą bardzo dobrze przybliżać funkcję pierwiastka kwadratowego, dla pozytywnego . Rzeczywiście, aby uzyskać wysoką dokładność w przedziale , potrzebujesz wyrażeń w serii. Aby uzyskać odpowiednie wagi ( ) i bieguny ( ), po prostu wyszukaj racjonalne przybliżenie funkcji online lub w książce.

\sqrt{x} \approx r (x) := \frac{a_{1}}{x + b_{1}} + \frac{a_{2}}{x + b_{2}} + \dots + \frac{a_{N}}{x + b_{N}},

$\sqrt{x} \approx r(x) := \frac{a_1}{x+b_1} + \frac{a_2}{x+b_2} + \dots + \frac{a_N}{x+b_N},$

b_{i}

$b_i$

[m, M]

$[m,M]$

O (\log \frac{M}{m})

$O(\log \frac{M}{m})$

a_{i}

$a_i$

- b_{i}

$-b_i$

Rozważmy teraz zastosowanie tej funkcji wymiernej do macierzy:

r (A) = a_{1} (A + b_{1} I)^{- 1} + a_{2} (A + b_{2} I)^{- 1} + \dots + a_{N} (A + b_{N} I)^{- 1} .

$r(A) = a_1(A + b_1 I)^{-1} + a_2(A + b_2 I)^{-1} + \dots + a_N(A + b_N I)^{-1}.$

Dzięki symetrii mamy , gdzie to wartości rozkład pojedyncza (SVD) z . Tak więc jakość aproksymacji racjonalnej macierzy jest równoważna jakości aproksymacji funkcji racjonalnej w miejscu wartości własnych. $A$

\begin{aligned} | | A^{1 / 2} - r (A) | |_{2} & = | | U (Σ^{1 / 2} - r (Σ)) U^{*} | |_{2}, \\ = max_{i} | \sqrt{σ_{i}} - r (σ_{i}) | \end{aligned}

$\begin{align} ||A^{1/2} - r(A)||_2 &= ||U\left(\Sigma^{1/2} - r(\Sigma)\right)U^*||_2, \\ &= \max_i |\sqrt{\sigma_i} - r(\sigma_i)| \end{align}$

A = U Σ U^{*}

$A = U \Sigma U^*$

A

$A$

Oznaczając numer warunku przez , możemy zastosować do dowolnej pożądanej tolerancji, wykonując dodatnio przesunięte graficznie rozwiązania Laplaciana formy, $A$ $\kappa$ $A^{1/2}b$ $O(\log \kappa)$

(A + b I) x = b .

$(A + bI)x=b.$

Rozwiązania te można wykonać za pomocą twojego ulubionego solvera Laplaciana - wolę techniki typu wielosiatkowego, ale ta w cytowanej przez ciebie pracy również powinna być w porządku. Dodatkowe pomaga tylko konwergencji solvera. $bI$

Aby uzyskać doskonały artykuł na ten temat, a także bardziej ogólne techniki analizy złożonej, które mają zastosowanie do macierzy niesymetrycznych, zobacz Obliczanie , i pokrewne funkcje macierzy przez całki konturu $A^α$ $\log(A)$ , Hale, Higham i Trefethen (2008 ).

Determinant „obliczenia”

Wyznacznik jest trudniejszy do obliczenia. O ile mi wiadomo, najlepszym sposobem jest obliczenie rozkładu Schura za pomocą algorytmu QR, a następnie odczytanie wartości własnych z przekątnej macierzy górnego trójkąta . Zajmuje to czas , gdzie jest liczbą węzłów na wykresie. $A = Q U Q^*$ $U$ $O(n^3)$ $n$

Jednak obliczanie wyznaczników jest z natury źle uwarunkowanym problemem, więc jeśli kiedykolwiek przeczytasz artykuł, który opiera się na obliczaniu wyznaczników dużej macierzy, powinieneś być bardzo sceptyczny wobec tej metody.

Na szczęście prawdopodobnie nie potrzebujesz wyznacznika. Na przykład,

Aby narysować próbki z pojedynczego rozkładu Gaussa , stała normalizacji jest taka sama we wszystkich punktach, więc nigdy nie trzeba jej obliczać. $N(0,A^{-1})$
Jeśli macierz Laplaciana reprezentuje odwrotną kowariancję lokalnego przybliżenia Gaussa w punkcie do rozkładu innego niż gaussowski, to wyznacznik rzeczywiście zmienia się z punktu do punktu. Jednak w każdym skutecznym schemacie próbkowania, jaki znam (w tym łańcuch Markowa Monte Carlo, ważność próbkowania itp.), Tak naprawdę potrzebny jest współczynnik determinujący , gdzie to bieżący punkt, a to proponowana następna próbka. $A = A_x$ $x$ $det (A_{x_{0}}^{- 1} A_{x_{p}}),$ $\det(A_{x_0}^{-1}A_{x_p}),$ $x_0$ $x_p$

Możemy postrzegać jako niskopoziomową aktualizację tożsamości, gdzie efektywna liczba rank, , aktualizacji niskiego poziomu jest lokalną miarą tego, jak prawdziwy rozkład niegaussowski; zazwyczaj jest to znacznie niższa niż pełna ranga matrycy. Rzeczywiście, jeśli jest duży, to prawdziwy rozkład jest lokalnie tak niegaussowski, że należy kwestionować całą strategię próbkowania tego rozkładu przy użyciu lokalnych przybliżeń Gaussa. $A_{x_0}^{-1}A_{x_p}$

A_{x_{0}}^{- 1} A_{x_{p}} = I + Q D Q^{*},

$A_{x_0}^{-1}A_{x_p} = I + Q D Q^*,$

r

$r$

r

$r$

Czynniki niskiego rzędu i można znaleźć w randomizowanym SVD lub Lanczos, stosując macierz do różnych wektorów, z których każde zastosowanie wymaga jednego wykresu Rozwiązanie Laplaciana. Zatem ogólna praca na rzecz uzyskania tych czynników niskiej rangi wynosi . $Q$ $D$

A_{x_{0}}^{- 1} A_{x_{p}} - I

$A_{x_0}^{-1}A_{x_p} -I$

O (r)

$O(r)$

O (r max (n, E))

$O(r \max(n,E))$

Znając , wyznacznikiem jest następnie $D = \text{diag}(d_1,d_2,\dots,d_r)$

det (A_{x_{0}}^{- 1} A_{x_{p}}) = det (I + Q D Q^{*}) = \exp (\sum_{i = 1}^{r} \log d_{i}) .

$\det(A_{x_0}^{-1}A_{x_p}) = \det(I + Q D Q^*) = \exp\left(\sum_{i=1}^r \log d_i\right).$

Te techniki obliczania racji wyznaczników niskiej rangi można znaleźć w Stochastycznej metodzie Newtona MCMC dla wielkoskalowych statystycznych odwrotnych problemów z zastosowaniem do inwersji sejsmicznej , Martin, i in. (2012). W tym artykule jest on stosowany do problemów z kontinuum, więc „wykres” jest siatką w przestrzeni 3D, a wykres Laplacian jest rzeczywistą macierzą Laplacian. Jednak wszystkie techniki mają zastosowanie do grafów Laplaciana. Prawdopodobnie są już inne artykuły stosujące tę technikę do ogólnych wykresów (rozszerzenie jest trywialne i zasadniczo to, co właśnie napisałem).

— Nick Alger
źródło