Testowanie, czy dwie matryce 12x12 mają tę samą determinantę

Dano mi macierz , która jest symetryczna, odwracalna, dodatnia określona i gęsta. Muszę sprawdzić, czy gdzie jest macierzą wszystkich jedynek. $12 \times 12$ $Q$

det (Q) = det (12 I - Q - J) (1)

$\det(Q) = \det(12I-Q-J) \; \; (1)$

J

$J$

Obecnie robię to z biblioteką pancernika, ale okazuje się, że jest zbyt wolny. Chodzi o to, że muszę to zrobić dla tryliona macierzy i okazuje się, że obliczenie dwóch wyznaczników jest wąskim gardłem mojego programu. Dlatego mam dwa pytania

Czy jest jakiś sposób na szybsze obliczenie wyznacznika, biorąc pod uwagę jego wielkość? Czy być może niechlujne rozszerzenie dla macierzy $12 \times12$ , które mogłoby działać w tym przypadku?
Czy istnieje jakiś inny skuteczny sposób sprawdzenia równości $(1)$

Edytować. Aby odpowiedzieć na komentarze. Muszę obliczyć wszystkie połączone niekomplementarne wykresy $G$ rzędu $13$ tak aby $G$ i $\overline{G}$ miały tę samą liczbę drzew rozpinających. Motywację do tego można znaleźć w tym poście mathoverflow . Jeśli chodzi o maszynę, to równolegle pracuję na 8-rdzeniowej maszynie 3.4GHh.

Edytować. Udało mi się zmniejszyć oczekiwany czas działania o 50%, tworząc program C do specyficznego obliczania wyznacznika macierzy $12 \times 12$ . Sugestie są nadal mile widziane.

— Jernej
źródło

Jak wolne jest zbyt wolne? Jak długo trwa na jakim sprzęcie? Czy tryliony tych są niezależne, aby można było obliczyć wiele z tych determinant równolegle? Jeśli tak, na jakiej maszynie można uruchomić? Co doprowadziło do tego problemu? Czy na pewno musisz obliczyć wyznaczniki?

Q

$Q$

— Bill Barth

Jak często (dla jakiej części przypadków) wyznaczniki są takie same / różne? Jeśli są one różne przez większość czasu, może być tańszy test, aby ustalić, że mogą się różnić, a można sprawdzić, czy są one takie same tylko w przypadku niepowodzenia pierwszego testu. Odwrotnie, jeśli są one takie same przez większość czasu.

— Wolfgang Bangerth,

Jak już pytano: czy możesz podać szczegółowe informacje na temat tego, skąd pochodzi ? Może jest lepsze podejście niż ślepe obliczanie determinant.

Q

$\mathbf Q$

— JM

Pojęcie, że warunek ten należy przetestować „na trylion matryc”, sugeruje 1), że jest znane, że apriori ma jakąś specjalną strukturę (w przeciwnym razie oczekiwanie, że warunek ten jest losowy, jest niewielkie), i 2) że lepsze podejście może być charakterystyka wszystkich macierzy pomocą tej właściwości (z efektywnie sprawdzalną formułą).

Q

$Q$

Q

$Q$

— hardmath

@hardmath Tak, jest macierzą całkowitą mającą przekątne od do i jako elementy

Q

$Q$

1

$1$

12

$12$

- 1

$-1$

— poziome

Odpowiedzi:

Ponieważ już używasz C ++, a twoje macierze są symetryczne dodatnio określone, wykonałbym niepodzielną faktoryzację dla a także dla $LDL^T$ $Q$ $12I-Q-J$ . Tutaj że jest również dodatnio określony, w przeciwnym razie będzie wymagał obrotu dla stabilności numerycznej (możliwe jest również, że nawet jeśli nie jest to wartość dodatnia, nie jest konieczne, ale musisz spróbować). $12I-Q-J$ $LDL^T$

Jest to szybsze niż faktoryzacja LU, a także szybsze niż Cholesky, ponieważ unika się pierwiastków kwadratowych. Wyznacznik jest po prostu iloczynem elementów diagonalnej macierzyKod do faktoryzacji LDL jest tak prosty, że możesz napisać go w mniej niż 50 liniach C. Strona Wikipedii $D$ na to opisuje algorytm, i mam kilka prostych matrycy kod do zrobienia Cholesky tutaj . Możesz to znacznie uprościć i zmodyfikować, aby uniknąć pierwiastka kwadratowego w celu zaimplementowania faktoryzacji $LDL^T$

Ponieważ możesz także kontrolować format przechowywania, możesz dalej optymalizować procedurę, aby przechowywać tylko połowę macierzy i pakować ją w szyk liniowy, aby zmaksymalizować lokalizację pamięci. Napisałbym również prosty niestandardowy produkt kropkowy i procedury aktualizacji rangi 1, ponieważ rozmiary problemów są tak małe, że powinieneś pozwolić kompilatorowi na wprowadzenie procedur, aby zmniejszyć obciążenie połączeń. Ponieważ jest to pętla o stałym rozmiarze, kompilator powinien mieć możliwość automatycznego wstawiania i rozwijania elementów w razie potrzeby.

Unikałbym próbowania sztuczek, aby wykorzystać fakt, że zawiera w wyrażeniu. Jest prawdopodobne, że przy tak małych rozmiarach problemów sztuczki te są wolniejsze niż po prostu wykonywanie dwóch oddzielnych obliczeń determinant. Oczywiście jedynym sposobem na zweryfikowanie tych twierdzeń jest wypróbowanie ich. $12I-Q-J$ $Q$

— Victor Liu
źródło

Po drugie zalecam wdrożenie , czyli Cholesky wolnego od korzeni, ponieważ Armadillo nie wydaje się mieć możliwości skorzystania z pozytywnej definitywności / dominacji diagonalnej.

L D L^{T}

$LDL^T$

— hardmath

Bez pewnych informacji na temat budowy tych pozytywnych określonych rzeczywistych macierzy symetrycznych, sugestie, które należy poczynić, są z konieczności dość ograniczone. $12\times 12$

Pobrałem pakiet Armadillo z Sourceforge i przejrzałem dokumentację. Postaraj się poprawić wydajność osobnego obliczania i , gdzie jest macierzą pierwszego rzędu wszystkich, ustawiając np . Dokumentacja zauważa, że jest to ustawienie domyślne dla macierzy do rozmiaru , więc przez pominięcie zakładam, że opcja jest domyślna dla $\det(Q)$ $\det(12I - Q - J)$ $J$ det(Q,slow=false) $4\times 4$ slow=true $12\times 12$ .

To, co slow=true przypuszczalnie działa, to częściowe lub pełne obrócenie w celu uzyskania formy rzutu rzędowego, z której łatwo można znaleźć wyznacznik. Jakkolwiek wiesz z góry, macierz jest pozytywnie określona, więc obrót nie jest konieczny dla stabilności (przynajmniej przypuszczalnie dla większości twoich obliczeń. Nie jest jasne, czy pakiet Armadillo rzuca wyjątek, jeśli osie stają się zbyt małe, ale powinno to być cecha rozsądnego numerycznego pakietu algebry liniowej EDYCJA: Znalazłem kod Armadillo, który implementuje się w pliku nagłówkowym , używając szablonów C ++ dla znacznej funkcjonalności. $Q$ detinclude\armadillo_bits\auxlib_meat.hppslow=false nie wpływa na to, jak $12\times 12$ wyznacznik zostanie wykonany, ponieważ obliczenia zostaną „przerzucone przez ścianę” do LAPACK (lub ATLAS) w tym punkcie bez wskazania, że obrót nie jest wymagany; widziećdet_lapack i jego wywołania w tym pliku.

Inną kwestią byłoby zastosowanie się do ich zalecenia budowy pakietu Armadillo łączącego szybkie zamienniki BLAS i LAPACK, jeśli rzeczywiście ich używasz; patrz ust. 5 pliku README.TXT Armadillo, aby uzyskać szczegółowe informacje. [Użycie dedykowanej 64-bitowej wersji BLAS lub LAPACK jest również zalecane ze względu na szybkość na obecnych komputerach 64-bitowych.]

Redukcja wiersza do postaci ekhelona jest zasadniczo eliminacją Gaussa i ma złożoność arytmetyczną . W przypadku obu macierzy stanowi to dwukrotność tej pracy lub . Operacje te mogą być „wąskim gardłem” w przetwarzaniu, ale nie ma nadziei, że bez specjalnej struktury w (lub niektórych znanych zależności między trylionami przypadków testowych umożliwiających amortyzację) praca mogłaby zostać zredukowana do $\frac{2}{3} n^3 + O(n^2)$ $\frac{4}{3} n^3 + O(n^2)$ $Q$ $O(n^2)$ .

Dla porównania, ekspansja przez kofaktory ogólnej macierzy obejmujeoperacje mnożenia (i mniej więcej tyle samo dodawania / odejmowania), więc dla porównanie ( vs. ) wyraźnie sprzyja eliminacji nad kofaktorami. $n\times n$ $n!$ $n=12$ $12! = 479001600$ $\frac{2}{3} n^3 = 1152$

Innym podejściem wymagającym pracy byłoby sprowadzenie do postaci tridiagonalnej za pomocą transformacji Householdera, co również nadaje formie tridiagonal. Obliczenia i można następnie wykonać w operacjach . [Wpływ aktualizacji rangi 1 na drugą determinantę można wyrazić jako czynnik skalarny podany przez rozwiązanie jednego układu trójosiowego.] $\frac{4}{3} n^3 + O(n^2)$ $Q$ $12I - Q$ $\det(Q)$ $\det(12I - Q - J)$ $O(n)$ $-J$

Wdrożenie takiego niezależnego obliczenia może być przydatne jako sprawdzenie wyników udanych (lub nieudanych) wywołań detfunkcji Armadillo .

Przypadek specjalny: jak sugeruje Komentarz Jernej, załóżmy, że gdzie jak poprzednio, jest macierzą (rangi 1) wszystkich, a jest macierz niepołączona (dodatnia). Rzeczywiście dla proponowanego zastosowania w teorii grafów byłyby to macierze całkowite. Następnie wyraźna formuła dla $Q = D - J$ $J$ $D=\text{diag}(d_1,\ldots,d_n)$ $\det(Q)$ to:

det (Q) = (\prod_{i = 1}^{n} d_{i}) (1 - \sum_{i = 1}^{n} d_{i}^{- 1})

$\det(Q) = \left(\prod_{i=1}^n d_i \right)\left(1 - \sum_{i=1}^n d_i^{-1} \right)$

Szkic dowodu daje możliwość zilustrowania szerszego zastosowania, tj. Ilekroć ma znaną determinantę, a układ jest szybko rozwiązywany. Zacznij od faktorowania: $D$ $Dv = (1\ldots 1)^T$

det (D - J) = det (D) \cdot det (I - D^{- 1} J)

$\det(D - J) = \det(D) \cdot \det(I - D^{-1}J)$

Teraz znów zajmuje 1 pozycję, a mianowicie $D^{-1}J$ $(d_1^{-1}\ldots d_n^{-1})^T(1\ldots 1)$ . Zauważ, że drugą determinantą jest po prostu:

f (1) = det (I - D^{- 1} J)

$f(1) = \det(I- D^{-1}J)$

gdzie jest charakterystyczny wielomianem . Jako macierz rangi 1, musi mieć (co najmniej) czynniki aby uwzględnić jej pustą przestrzeń. „Brakująca” wartość własna to , jak można zobaczyć na podstawie obliczeń: $f(x)$ $D^{-1}J$ $f(x)$ $n-1$ $x$ $\sum d_i^{-1}$

D^{- 1} J (d_{1}^{- 1} \dots d_{n}^{- 1})^{T} = (\sum d_{i}^{- 1}) (d_{1}^{- 1} \dots d_{n}^{- 1})^{T}

$D^{-1}J\; (d_1^{-1}\ldots d_n^{-1})^T = \left(\sum d_i^{-1} \right) (d_1^{-1}\ldots d_n^{-1})^T$

Wynika z tego, że charakterystyczny wielomian , a jest taki jak pokazano powyżej dla , . $f(x) = x^{n-1} (x - \sum d_i^{-1})$ $f(1)$ $\det(I - D^{-1}J)$ $1 - \sum d_i^{-1}$

Zauważ też, że jeśli , to , macierz diagonalna, której wyznacznikiem jest po prostu iloczyn jej wpisów diagonalnych. $Q = D-J$ $12I - Q - J = 12I - D + J - J = 12I - D$

— matematyka
źródło

Hm .. jest w rzeczywistości gdzie jest macierzą przylegania więc myślę, że ten wynik może być niepoprawny. W szczególności oznaczałoby to, że liczba drzew spinających na wykresie jest określona przez sekwencję stopni, która się nie utrzymuje.

Q

$Q$

D - A

$D-A$

A

$A$

G

$G$

G

$G$

— Jernej

Pozuzakątowe wpisy będą wówczas ogólnie zawierać 0, a także -1. rozkładu sugeruje Victor wykorzystuje symetrii i zmniejsza prowadzącą określenie w liczbie operacji z do . Istnieje dokładne podejście do liczb całkowitych, ale prawdopodobnie nie jest to potrzebne w przypadku skromnej macierzy rozmiarów i wpisów. Jeśli rozumiem konstrukcję, jest pozytywnie określony z tego samego powodu, dla którego jest.

Q

$Q$

L D L^{T}

$LDL^T$

\frac{2}{3} n^{3}

$\frac{2}{3}n^3$

\frac{1}{3} n^{3}

$\frac{1}{3}n^3$

12 I - Q - J

$12I-Q-J$

Q

$Q$

— hardmath

@Jernej: Jeśli uważasz, że coś, co powiedziałem, jest niepoprawne, utworzyłem pokój czatu na podstawie tego pytania, w którym dyskusję można prowadzić bez niepotrzebnego komentowania tutaj.

— hardmath

Jeśli masz uporządkowany sposób wyliczania wykresów, które chcesz obliczyć wyznaczniki, być może znajdziesz aktualizacje niskiej rangi, które przenoszą cię z jednego wykresu na drugi.

Jeśli tak, to możesz użyć lematu wyznacznika macierzy, aby tanio obliczyć wyznacznik kolejnego wykresu, który ma być zliczony, korzystając z wiedzy o wyznaczniku bieżącego wykresu.

To znaczy, dla macierzy i wektorów : Można to uogólnić, jeśli U i V są macierzy, a jest : $A$ $u,v$

det (A + u v^{T}) = (1 + v^{T} A^{- 1} u) det (A)

$\det(A+uv^T)=(1+v^T A^{-1} u)\det(A)$

n \times m

$n \times m$

A

$A$

n \times n

$n \times n$

det (A + U V^{T}) = det (I_{m} + V^{T} A^{- 1} U) det (A)

$\det(A+UV^T)=\det(I_m + V^TA^{-1}U)\det(A)$

Aby efektywnie obliczyć odwrotność, możesz użyć wzoru Shermana-Morrisona, aby uzyskać odwrotność kolejnej macierzy z bieżącej:

(A + u v^{T})^{- 1} = A^{- 1} - \frac{A^{- 1} u v^{T} A^{- 1}}{1 + v^{T} A^{- 1} u}

$(A+uv^T)^{-1}=A^{-1}-\frac{A^{-1} uv^T A^{-1}}{1+v^T A^{-1} u}$

— Richard
źródło