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

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

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ę.

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ą).

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 tylko$10 \times 10$macierze, 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ę.

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ą) .