najszybsze rozwiązanie układu liniowego dla małych macierzy kwadratowych (10x10)

9

Jestem bardzo zainteresowany optymalizacją rozwiązania liniowego rozwiązywania problemów dla małych matryc (10x10), czasami nazywanych drobnymi matrycami. Czy istnieje na to gotowe rozwiązanie? Matryca może być przyjęta jako niejedna.

Ten solver ma zostać wykonany ponad 1 000 000 razy w mikrosekundach na procesorze Intel. Mówię o poziomie optymalizacji stosowanym w grach komputerowych. Bez względu na to, czy koduję go w asemblerze i architekturze, czy badam redukcje kompromisów precyzji lub niezawodności i używam hacków zmiennoprzecinkowych (używam flagi kompilacji -ffast-matematyki, nie ma problemu). Rozwiązanie może się nawet nie powieść przez około 20% czasu!

Częściowy PivLu firmy Eigen jest najszybszy w moim bieżącym teście, wyprzedzając LAPACK, gdy jest zoptymalizowany z -O3 i dobrym kompilatorem. Ale teraz jestem w trakcie tworzenia niestandardowego solwera liniowego. Wszelkie porady będą mile widziane. Sprawię, że moje rozwiązanie stanie się open source i uzyskam kluczowe informacje w publikacjach itp.

Powiązane: Szybkość rozwiązywania układu liniowego z blokową macierzą diagonalną Jaka jest najszybsza metoda odwrócenia milionów macierzy? https://stackoverflow.com/q/50909385/1489510

rfabbri
źródło
7
To wygląda jak cel rozciągnięty. Załóżmy, że używamy najszybszej Skylake-X Xeon Platinum 8180 z teoretyczną szczytową przepustowością 4 pojedynczych precyzyjnych TFLOP i że jeden system 10x10 wymaga około 700 (około 2n ** 3/3) operacji zmiennoprzecinkowych do rozwiązania. Następnie partia 1M takich systemów mogłaby teoretycznie zostać rozwiązana w 175 mikrosekund. To nie można przekroczyć prędkości światła. Czy możesz podzielić się wydajnością, którą obecnie osiągasz, z najszybszym istniejącym kodem? BTW, czy dane są pojedynczą precyzją czy podwójną precyzją?
njuffa
@ njuffa tak, chciałem osiągnąć blisko 1ms, ale micro to inna historia. W przypadku mikro rozważałem wykorzystanie inkrementalnej struktury odwrotnej w partii poprzez wykrycie podobnych matryc, które występują często. Perf ma obecnie zasięg 10-500 ms, w zależności od procesora. Precyzja jest podwójna lub nawet złożona podwójna. Pojedyncza precyzja działa wolniej.
rfabbri
@ njuffa Mogę zmniejszyć lub zwiększyć precyzję prędkości
rfabbri
2
Wygląda na to, że precyzja / dokładność nie jest twoim priorytetem. Czy może być przydatna metoda iteracyjna skrócona przy stosunkowo niewielkiej liczbie ocen? Zwłaszcza jeśli masz początkowe domysły.
Spencer Bryngelson
1
Czy obracasz się? Czy możesz zrobić faktoryzację QR zamiast eliminacji Gaussa. Czy przeplatasz swoje systemy, abyś mógł korzystać z instrukcji SIMD i wykonywać kilka systemów jednocześnie? Czy piszesz programy liniowe bez pętli i bez adresowania pośredniego? Jakiej dokładności chcesz i jak uwarunkuję twój system? Czy mają jakąkolwiek strukturę, którą można by wykorzystać?
Carl Christian

Odpowiedzi:

7

Użycie typu macierzy Eigen, w którym liczba wierszy i kolumn jest zakodowana w typie w czasie kompilacji, daje przewagę nad LAPACK, gdzie rozmiar matrycy jest znany tylko w czasie wykonywania. Te dodatkowe informacje pozwalają kompilatorowi na pełne lub częściowe rozwijanie pętli, eliminując wiele instrukcji rozgałęzienia. Jeśli chcesz używać istniejącej biblioteki, a nie pisać własne jądra, typ danych, w którym można uwzględnić rozmiar matrycy jako parametry szablonu C ++, będzie prawdopodobnie niezbędny. Jedyną znaną mi biblioteką, która to robi, jest blask , dlatego warto porównać testy porównawcze z Eigenem.

Jeśli zdecydujesz się na wdrożenie własnej implementacji, możesz przekonać się, co robi PETSc, jeśli blokowy format CSR jest przydatnym przykładem, chociaż sam PETSc prawdopodobnie nie będzie odpowiednim narzędziem do tego, co masz na myśli. Zamiast pisać pętlę, wypisują każdą pojedynczą operację dla małych macierzy-wektorów mnożąc się jawnie (zobacz ten plik w swoim repozytorium). Gwarantuje to, że nie ma instrukcji rozgałęzienia, które można uzyskać za pomocą pętli. Wersje kodu z instrukcjami AVX są dobrym przykładem tego, jak faktycznie używać rozszerzeń wektorowych. Na przykład ta funkcja używa__m256dtyp danych do jednoczesnego działania na czterech podwójnych kartach jednocześnie. Można uzyskać znaczny wzrost wydajności, jawnie zapisując wszystkie operacje przy użyciu rozszerzeń wektorów, tylko w przypadku faktoryzacji LU zamiast mnożenia macierzy-wektora. Zamiast pisać kod C ręcznie, lepiej byłoby użyć skryptu do jego wygenerowania. Przyjemne może być także sprawdzenie, czy istnieje znacząca różnica w wydajności, gdy zmienisz kolejność niektórych operacji, aby lepiej skorzystać z potoku instrukcji.

Możesz także uzyskać przebieg dzięki narzędziu STOKE , które losowo zbada przestrzeń możliwych przekształceń programu, aby znaleźć szybszą wersję.

Daniel Shapero
źródło
tx. Już z powodzeniem używam Eigen jak Map <const Matrix <complex, 10, 10>> AA (A). sprawdzi inne rzeczy.
rfabbri
Eigen ma także AVX, a nawet nagłówek complex.h. Dlaczego PETSc do tego? W tym przypadku trudno konkurować z Eigenem. Specjalizowałem Eigen jeszcze bardziej dla mojego problemu i stosując przybliżoną strategię przestawiania, która zamiast przejmować maksimum nad kolumną, zamienia natychmiast oś, gdy znajdzie inną, która jest o 3 rzędy wielkości większa.
rfabbri
1
@rfabbri Nie sugerowałem, że używasz do tego PETSc, tylko to, co robią w tym konkretnym przypadku, może być pouczające. Zredagowałem odpowiedź, aby to wyjaśnić.
Daniel Shapero,
4

Innym pomysłem może być zastosowanie podejścia generatywnego (program piszący program). Autor (meta) programu, który wyrzuca sekwencję instrukcji C / C ++, aby wykonać niepodzielną ** LU w systemie 10x10. W zasadzie biorąc gniazdo pętli k / i / j i spłaszczając je do O (1000) lub mniej więcej linii arytmetyka skalarna. Następnie wprowadź wygenerowany program do dowolnego kompilatora optymalizującego. To, co moim zdaniem jest interesujące, polega na tym, że usunięcie pętli ujawnia każdą zależność danych i nadmiarowe podwyrażenie oraz daje kompilatorowi maksymalną możliwość zmiany kolejności instrukcji, tak aby dobrze odwzorowały rzeczywisty sprzęt (np. Liczbę jednostek wykonawczych, zagrożeń / przeciągnięć, więc na).

Jeśli znasz wszystkie macierze (lub nawet kilka z nich), możesz poprawić przepustowość, wywołując funkcje / funkcje SIMD (SSE / AVX) zamiast kodu skalarnego. W tym przypadku wykorzystujesz zawstydzający paralelizm między instancjami, zamiast gonić za paralelizmem w ramach jednej instancji. Na przykład, możesz wykonać 4 jednostki LU o podwójnej precyzji jednocześnie używając wewnętrznych elementów AVX256, pakując 4 macierze „przez” rejestr i wykonując te same operacje ** na wszystkich z nich.

** Stąd nacisk na niepodzielną LU. Obracanie psuje to podejście na dwa sposoby. Po pierwsze, wprowadza rozgałęzienia z powodu wyboru przestawnego, co oznacza, że ​​twoje zależności danych nie są tak doskonale znane. Po drugie, oznacza to, że różne „gniazda” SIMD musiałyby robić różne rzeczy, ponieważ instancja A może obracać się inaczej niż instancja B. Więc jeśli wykonasz którąkolwiek z tych czynności, sugerowałbym statyczne przestawienie macierzy przed obliczeniami (najwyższy dopuszczalny wpis każdej kolumny na przekątną).

rchilton1980
źródło
ponieważ matryce są tak małe, być może można je obrócić, jeśli zostaną wstępnie przeskalowane. Nawet przed obrotem matryc. Wszystko, czego potrzebujemy, to aby wpisy były w odległości 2-3 rzędów wielkości od siebie.
rfabbri
2

Twoje pytanie prowadzi do dwóch różnych rozważań.

Najpierw musisz wybrać odpowiedni algorytm. Dlatego należy rozważyć pytanie, czy macierze mają jakąkolwiek strukturę. Na przykład, gdy macierze są symetryczne, rozkład Cholesky'ego jest bardziej wydajny niż LU. Gdy potrzebujesz tylko ograniczonej dokładności, metoda iteracyjna może być szybsza.

Po drugie, musisz efektywnie zaimplementować algorytm. Aby to zrobić, musisz znać wąskie gardło swojego algorytmu. Czy twoja implementacja jest związana szybkością transferu pamięci lub szybkością obliczeń. Ponieważ rozważasz tylko10×10macierze, macierz powinna całkowicie zmieścić się w pamięci podręcznej procesora. Dlatego powinieneś skorzystać z jednostek SIMD (SSE, AVX itp.) I rdzeni swojego procesora, aby wykonać jak najwięcej obliczeń na cykl.

Podsumowując, odpowiedź na twoje pytanie w dużej mierze zależy od rozważanego sprzętu i matryc. Prawdopodobnie nie ma jednoznacznej odpowiedzi i musisz wypróbować kilka rzeczy, aby znaleźć optymalną metodę.

H. Rittich
źródło
Do tej pory Eigen już mocno zoptymalizował, używa SEE, AVX itp. Próbowałem metod iteracyjnych we wstępnym teście i one nie pomogły. Próbowałem Intel MKL, ale nie lepiej niż Eigen ze zoptymalizowanymi flagami GCC. Obecnie próbuję stworzyć coś lepszego i prostszego niż Eigen i przeprowadzić bardziej szczegółowe testy przy użyciu metod iteracyjnych.
rfabbri,
1

Spróbowałbym odwrócić blokowo.

https://en.wikipedia.org/wiki/Invertible_matrix#Blockwise_inversion

Eigen używa zoptymalizowanej procedury do obliczenia odwrotności macierzy 4x4, co jest prawdopodobnie najlepszą możliwą wartością. Spróbuj użyć tego jak najwięcej.

http://www.eigen.tuxfamily.org/dox/Inverse__SSE_8h_source.html

Lewy górny róg: 8x8. U góry po prawej: 8x2. U dołu po lewej: 2x8. Na dole po prawej: 2x2. Odwróć 8x8 za pomocą zoptymalizowanego kodu inwersyjnego 4x4. Reszta to produkty matrycowe.

EDYCJA: Używanie bloków 6x6, 6x4, 4x6 i 4x4 okazało się nieco szybsze niż to, co opisałem powyżej.

using namespace Eigen;

template<typename Scalar, int tl_size, int br_size>
Matrix<Scalar, tl_size + br_size, tl_size + br_size> blockwise_inversion(const Matrix<Scalar, tl_size, tl_size>& A, const Matrix<Scalar, tl_size, br_size>& B, const Matrix<Scalar, br_size, tl_size>& C, const Matrix<Scalar, br_size, br_size>& D)
{
    Matrix<Scalar, tl_size + br_size, tl_size + br_size> result;

    Matrix<Scalar, tl_size, tl_size> A_inv = A.inverse().eval();
    Matrix<Scalar, br_size, br_size> DCAB_inv = (D - C * A_inv * B).inverse();

    result.topLeftCorner<tl_size, tl_size>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<tl_size, br_size>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<br_size, tl_size>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<br_size, br_size>() = DCAB_inv;

    return result;
}

template<typename Scalar, int tl_size, int br_size>
Matrix<Scalar, tl_size + br_size, tl_size + br_size> my_inverse(const Matrix<Scalar, tl_size + br_size, tl_size + br_size>& mat)
{
    const Matrix<Scalar, tl_size, tl_size>& A = mat.topLeftCorner<tl_size, tl_size>();
    const Matrix<Scalar, tl_size, br_size>& B = mat.topRightCorner<tl_size, br_size>();
    const Matrix<Scalar, br_size, tl_size>& C = mat.bottomLeftCorner<br_size, tl_size>();
    const Matrix<Scalar, br_size, br_size>& D = mat.bottomRightCorner<br_size, br_size>();

    return blockwise_inversion<Scalar,tl_size,br_size>(A, B, C, D);
}

template<typename Scalar>
Matrix<Scalar, 10, 10> invert_10_blockwise_8_2(const Matrix<Scalar, 10, 10>& input)
{
    Matrix<Scalar, 10, 10> result;

    const Matrix<Scalar, 8, 8>& A = input.topLeftCorner<8, 8>();
    const Matrix<Scalar, 8, 2>& B = input.topRightCorner<8, 2>();
    const Matrix<Scalar, 2, 8>& C = input.bottomLeftCorner<2, 8>();
    const Matrix<Scalar, 2, 2>& D = input.bottomRightCorner<2, 2>();

    Matrix<Scalar, 8, 8> A_inv = my_inverse<Scalar, 4, 4>(A);
    Matrix<Scalar, 2, 2> DCAB_inv = (D - C * A_inv * B).inverse();

    result.topLeftCorner<8, 8>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<8, 2>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<2, 8>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<2, 2>() = DCAB_inv;

    return result;
}

template<typename Scalar>
Matrix<Scalar, 10, 10> invert_10_blockwise_6_4(const Matrix<Scalar, 10, 10>& input)
{
    Matrix<Scalar, 10, 10> result;

    const Matrix<Scalar, 6, 6>& A = input.topLeftCorner<6, 6>();
    const Matrix<Scalar, 6, 4>& B = input.topRightCorner<6, 4>();
    const Matrix<Scalar, 4, 6>& C = input.bottomLeftCorner<4, 6>();
    const Matrix<Scalar, 4, 4>& D = input.bottomRightCorner<4, 4>();

    Matrix<Scalar, 6, 6> A_inv = my_inverse<Scalar, 4, 2>(A);
    Matrix<Scalar, 4, 4> DCAB_inv = (D - C * A_inv * B).inverse().eval();

    result.topLeftCorner<6, 6>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<6, 4>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<4, 6>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<4, 4>() = DCAB_inv;

    return result;
}

Oto wyniki jednego testu porównawczego z użyciem miliona Eigen::Matrix<double,10,10>::Random()macierzy i Eigen::Matrix<double,10,1>::Random()wektorów. We wszystkich moich testach moja odwrotność jest zawsze szybsza. Moja procedura rozwiązywania polega na obliczeniu odwrotności, a następnie pomnożeniu jej przez wektor. Czasami jest szybszy niż Eigen, czasem nie. Moja metoda oznaczania na ławce może być wadliwa (nie wyłączała turbo boost itp.). Ponadto losowe funkcje Eigen mogą nie reprezentować rzeczywistych danych.

  • Odwrotność częściowa własna odwrotna: 3036 milisekund
  • Moja odwrotność z górnym blokiem 8x8: 1638 milisekund
  • Moja odwrotność z górnym blokiem 6x6: 1234 milisekundy
  • Częściowe rozwiązanie osi własnych w rozwiązaniu: 1791 milisekund
  • Moje rozwiązanie z górnym blokiem 8x8: 1739 milisekund
  • Moje rozwiązanie z górnym blokiem 6x6: 1286 milisekund

Bardzo mnie interesuje, czy ktoś może to dalej zoptymalizować, ponieważ mam aplikację elementów skończonych, która odwraca matryce gazillionów 10x10 (i tak, potrzebuję indywidualnych współczynników odwrotności, więc bezpośrednie rozwiązanie układu liniowego nie zawsze jest opcją) .

Charlie S.
źródło