Mnożenie małej macierzy rzadkiej

Question

dla klasy Muszę napisać własną liniową procedurę rozwiązywania równań dla macierzy rzadkich. Mam swobodę korzystania z dowolnej struktury danych dla rzadkich macierzy i muszę zaimplementować kilka rozwiązań, w tym koniugować gradient.

Zastanawiam się, czy istnieje znany sposób przechowywania rzadkich matryc, tak że mnożenie za pomocą wektora jest stosunkowo szybkie.

W tej chwili moje rzadkie macierze są w zasadzie zaimplementowane jako zapakowane std::map< std::pair<int, int>, double>, które przechowuje dane, jeśli takie istnieją. Przekształca to mnożenie macierzy od wektora do złożoności O (n²) do O (n²log (n)), ponieważ muszę wykonać wyszukiwanie dla każdego elementu macierzy. Zajrzałem do formatu macierzy Yale Sparse i wydaje się, że pobieranie elementu jest również w O (log (n)), więc nie jestem pewien, czy byłoby to znacznie szybsze.

Dla odniesienia mam macierz 800x800, która jest wypełniona 5000 wpisów. Zajmuje to około 450 sekund, aby rozwiązać taki system metodą gradientu sprzężonego.

Czy myślisz, że można to zrobić szybciej z inną strukturą danych?

dziękuję!

Answer 1

Najczęściej wybieranymi opcjami są CSC or CSR storage. Oba są skuteczne w mnożeniu macierzy-wektor. Kodowanie tych procedur mnożenia jest bardzo proste, jeśli musisz to zrobić sam.

To powiedziawszy, przechowywanie w Yale daje także bardzo wydajne mnożenie wektorów macierzy. Jeśli przeprowadzasz wyszukiwanie w elementach macierzy, błędnie rozumiesz, jak używać tego formatu. Proponuję zapoznać się ze standardowymi bibliotekami rzadkimi, aby dowiedzieć się, w jaki sposób zaimplementowano multiplikację macierzy-wektor.

Nawet z bieżącym miejscem przechowywania można wykonywać mnożenie macierzy w złożoności O (n). Wszystkie rzadkie algorytmy mnożenia wektorów, które kiedykolwiek widziałem, sprowadzają się do tych samych kroków. Na przykład rozważ y = Siekiera.

1. Zeruj wektor wyników, y.
2. Inicjalizacja iteratora dla niezerowych elementów macierzy, A.
3. Uzyskaj następny niezerowy element macierzy, A [i, j] say. Zwróć uwagę, że wzorzec i, j nie ma regularnego wzorca. To po prostu odzwierciedla kolejność, w której przechowywane są niezerowe elementy A.
4. y [i] + = A [i, j] * x [j]
5. Jeśli istnieje więcej elementów, goto 3.

Podejrzewam piszesz klasyczny podwójny pętli for gęsty kod mnożenia:

for (i=0; i<N; i++) 
    for (j=0; j<N; j++) 
     y[i] += A[i,j]*x[j] 

i to właśnie prowadzi do wykonywania wyszukiwań.

Ale nie sugeruję, abyś trzymał się swojej pamięci masowej std::map. To nie będzie super wydajne. Polecam CSC głównie dlatego, że jest najczęściej używany.