Faktoryzacja jako klucz do optymalizacji obliczeń matematycznych

Faktoryzacja nie przyspiesza obliczeń dlatego, że zapis wygląda krócej. Jej przewaga polega na tym, że zmienia strukturę problemu: usuwa powtarzające się działania, ujawnia miejsca zerowe, pozwala zastosować szybszy algorytm albo umożliwia wielokrotne wykorzystanie raz wykonanego rozkładu.

Różnicę widać już w prostym wyrażeniu:

(6x^3 + 9x^2 = 3x^2(2x+3))

Obie postacie są matematycznie równoważne, ale odpowiadają na inne pytania. Postać rozwinięta ułatwia odczytanie współczynników, natomiast postać iloczynowa natychmiast pokazuje wspólny czynnik, miejsca zerowe i możliwość ograniczenia liczby operacji. W dużych obliczeniach numerycznych taka zmiana reprezentacji potrafi zmniejszyć czas pracy z godzin do minut. W źle dobranym przypadku może natomiast nie dać żadnego zysku, a nawet zwiększyć koszt obliczeń i zużycie pamięci.

Największy zysk daje usunięcie działań wykonywanych wielokrotnie

Najprostsza forma optymalizacji polega na wyłączeniu wspólnego czynnika przed nawias. Zamiast obliczać osobno:

(ax + ay + az)

można zapisać:

(a(x+y+z))

W pierwszej wersji potrzebne są trzy mnożenia i dwa dodawania. W drugiej — jedno mnożenie i dwa dodawania. Przy pojedynczym rachunku oszczędność jest niewielka. Jeżeli jednak wyrażenie jest liczone dla milionów rekordów, dwóch kanałów obrazu o rozdzielczości 4K albo każdego punktu symulacji, różnica zaczyna być mierzalna.

Ten mechanizm wykorzystują kompilatory oraz systemy algebry komputerowej, takie jak Mathematica, Maple, MATLAB, SymPy czy SageMath. Jedną z typowych optymalizacji jest eliminacja wspólnych podwyrażeń. Jeżeli program wielokrotnie oblicza fragment:

(x^2 + 2xy + y^2)

można zastąpić go przez:

((x+y)^2)

Pierwsza postać wymaga między innymi wyznaczenia dwóch kwadratów, iloczynu (xy), dodatkowego mnożenia przez 2 i dwóch dodawań. Druga sprowadza rachunek do jednego dodawania i jednego podniesienia do kwadratu. Trzeba jednak uwzględnić arytmetykę zmiennoprzecinkową. Matematyczna równoważność nie oznacza identycznego wyniku bit po bicie, ponieważ zmiana kolejności działań może zmienić błąd zaokrąglenia.

Faktoryzacja szczególnie dobrze działa w trzech sytuacjach:

• gdy ten sam fragment wyrażenia pojawia się wielokrotnie,
• gdy zapis iloczynowy umożliwia pominięcie dalszych obliczeń, na przykład po wykryciu czynnika równego zero,
• gdy wynik rozkładu można wykorzystać dla wielu zestawów danych.

Dobrym przykładem jest rozwiązywanie układu równań liniowych:

(Ax=b)

Jeżeli macierz (A) pozostaje taka sama, ale zmienia się wektor (b), nie należy za każdym razem rozpoczynać obliczeń od początku. Dla gęstej macierzy kwadratowej o wymiarze (n \times n) rozkład LU wymaga w przybliżeniu (\frac{2}{3}n^3) operacji arytmetycznych. Późniejsze rozwiązanie układu dla jednego nowego wektora prawej strony kosztuje już około (2n^2) operacji.

Dla macierzy (1000 \times 1000) oznacza to orientacyjnie:

• około 667 mln operacji na wykonanie rozkładu LU,
• około 2 mln operacji na rozwiązanie układu po wykonaniu rozkładu.

Jeżeli trzeba rozwiązać sto układów z tą samą macierzą (A), ponowne liczenie faktoryzacji sto razy jest kosztownym błędem. Rozkład wykonuje się raz, zapisuje jego wynik, a następnie wykorzystuje go dla kolejnych wektorów (b).

W praktyce pierwszym pytaniem nie powinno więc być: „Czy da się rozłożyć to wyrażenie?”, lecz: „Który fragment obliczeń powtarzam i czy mogę policzyć go tylko raz?”

Rodzaj rozkładu trzeba dopasować do właściwości problemu

Słowo „faktoryzacja” obejmuje kilka różnych operacji. Rozkład liczby na czynniki pierwsze, rozkład wielomianu i rozkład macierzy służą innym celom. Łączenie ich w jedną metodę prowadzi do błędnych oczekiwań.

W przypadku liczb całkowitych najprostszy algorytm polega na sprawdzaniu kolejnych dzielników do (\sqrt{n}). Dla liczby 10 000 wystarczy przetestować dzielniki nie większe niż 100. Dla liczby mającej kilkaset cyfr takie podejście staje się bezużyteczne.

W zależności od wielkości liczby i spodziewanych czynników wykorzystuje się między innymi:

• dzielenie próbne dla małych liczb i małych czynników,
• metodę Pollarda rho do wyszukiwania stosunkowo niewielkich czynników,
• metodę krzywych eliptycznych ECM, gdy liczba ma duży rozmiar, ale może zawierać czynnik liczący kilkadziesiąt cyfr,
• sito kwadratowe dla większych liczb złożonych,
• ogólne sito ciała liczbowego, czyli GNFS, dla bardzo dużych liczb bez znanej szczególnej struktury.

Nie istnieje szybki klasyczny algorytm, który w rozsądnym czasie rozłoży każdą dowolnie dużą liczbę. To ograniczenie ma praktyczne znaczenie w kryptografii. W systemie RSA moduł publiczny jest iloczynem dwóch dużych liczb pierwszych. Znajomość tego iloczynu nie daje efektywnej klasycznej metody odzyskania czynników, jeżeli parametry zostały wygenerowane poprawnie i mają odpowiedni rozmiar.

Nie należy jednak wyciągać z tego wniosku, że każda faktoryzacja jest kosztowna. Liczba:

(2^{20} \cdot 3^{10})

może mieć dużą wartość, ale jej struktura jest oczywista. Trudność zależy nie tylko od liczby cyfr, lecz także od wielkości najmniejszego czynnika, budowy liczby oraz dostępnej wiedzy o jej pochodzeniu.

W przypadku wielomianów rozkład pozwala szybko odczytać pierwiastki. Dla:

(x^2-5x+6=(x-2)(x-3))

równanie (x^2-5x+6=0) można rozwiązać bez dalszych przekształceń: (x=2) lub (x=3).

Nie każdy wielomian rozkłada się jednak na czynniki w wybranym zbiorze liczb. Wyrażenie:

(x^2+1)

jest nierozkładalne na czynniki liniowe nad liczbami rzeczywistymi, ale nad liczbami zespolonymi przyjmuje postać:

((x-i)(x+i))

Przed rozpoczęciem faktoryzacji trzeba więc ustalić, czy obliczenia są prowadzone nad liczbami całkowitymi, wymiernymi, rzeczywistymi, zespolonymi czy w ciele skończonym. Wynik może być inny w każdym z tych przypadków.

Dla macierzy wybór metody jest jeszcze bardziej praktyczny:

• LU sprawdza się w ogólnych kwadratowych układach równań,
• Cholesky jest właściwym wyborem dla macierzy symetrycznych dodatnio określonych,
• QR jest zwykle bezpieczniejszy numerycznie przy rozwiązywaniu problemów najmniejszych kwadratów,
• SVD pomaga w analizie macierzy bliskich osobliwości, redukcji wymiaru i wyznaczaniu pseudoodwrotności.

Rozkład Cholesky’ego dla gęstej macierzy wymaga około (\frac{1}{3}n^3) operacji. To mniej więcej połowa głównego kosztu rozkładu LU. Zysk jest konkretny, lecz warunkowy: macierz musi być symetryczna i dodatnio określona. Próba wymuszenia tej metody na niewłaściwych danych kończy się błędem albo wynikiem pozbawionym sensu.

QR kosztuje więcej niż LU, ale lepiej radzi sobie z częścią problemów źle uwarunkowanych. SVD jest jeszcze droższy, za to dostarcza informacji o wartościach osobliwych, randze numerycznej i kierunkach, w których dane tracą informację. Stosowanie SVD do każdego zwykłego układu równań jest przesadą. Rezygnowanie z niego przy macierzy niemal osobliwej bywa natomiast fałszywą oszczędnością.

Optymalizacja kończy się tam, gdzie zaczynają się błędy numeryczne

Krótszy wzór nie zawsze oznacza lepszy algorytm. Najbardziej problematyczny przypadek pojawia się wtedy, gdy przekształcenie zwiększa ryzyko przepełnienia, utraty cyfr znaczących albo niestabilności numerycznej.

Rozważmy wyrażenie:

(x^2-y^2=(x-y)(x+y))

Dla liczb całkowitych obie formy są równoważne bez zastrzeżeń, o ile nie dojdzie do przepełnienia zakresu typu danych. W arytmetyce zmiennoprzecinkowej wyniki mogą się jednak minimalnie różnić. Jeżeli (x) i (y) są bardzo duże i prawie równe, obliczenie (x-y) prowadzi do katastrofalnej utraty cyfr znaczących. Z drugiej strony postać rozwinięta może najpierw wyznaczyć dwa ogromne kwadraty i dopiero później je odjąć, co również jest niebezpieczne.

Nie ma tu uniwersalnej odpowiedzi. Trzeba sprawdzić:

• zakres wartości wejściowych,
• typ danych, na przykład `float32`, `float64` lub liczby całkowite 64-bitowe,
• dopuszczalny błąd względny i bezwzględny,
• możliwość wystąpienia wartości bardzo bliskich zeru,
• uwarunkowanie całego problemu.

W praktyce najpierw optymalizuje się algorytm, a dopiero później pojedyncze działania. Zmniejszenie liczby mnożeń z dziesięciu do ośmiu nie pomoże, jeżeli wcześniej wybrano metodę o złożoności (O(n^3)) tam, gdzie można było wykorzystać strukturę rzadką i zejść do kosztu zależnego od liczby elementów niezerowych.

To częsty błąd przy macierzach. Gęsta macierz (10,000 \times 10,000) zawiera 100 mln elementów. Przy zapisie w formacie `float64` same wartości zajmują około 800 MB, bez uwzględnienia kopii roboczych i struktur pomocniczych. Jeżeli tylko 0,1% elementów jest niezerowych, przechowywanie jej jako macierzy gęstej marnuje pamięć i może uniemożliwić wykonanie obliczeń na zwykłej stacji roboczej.

Faktoryzacja macierzy rzadkiej wymaga jednak uwagi na zjawisko fill-in. Podczas rozkładu w miejscach pierwotnie zerowych mogą pojawić się wartości niezerowe. Źle dobrana kolejność wierszy i kolumn potrafi wielokrotnie zwiększyć zużycie pamięci. Dlatego biblioteki takie jak SuiteSparse, Intel oneMKL, SciPy czy MATLAB stosują permutacje ograniczające fill-in. Sam wybór formatu rzadkiego nie gwarantuje więc szybkiego wyniku.

Podobny problem występuje w obliczeniach symbolicznych. Rozłożenie wyrażenia na czynniki może ułatwić analizę pierwiastków, ale późniejsze mnożenie, różniczkowanie albo podstawianie danych nie zawsze będzie szybsze. Czasem bardziej użyteczna jest postać rozwinięta, czasem iloczynowa, a czasem schemat Hornera.

Przykładowo wielomian:

(a_0+a_1x+a_2x^2+a_3x^3)

można obliczać jako:

(((a_3x+a_2)x+a_1)x+a_0)

Schemat Hornera ogranicza liczbę mnożeń i nie wymaga osobnego wyznaczania kolejnych potęg (x). Dla wielomianu stopnia (n) potrzebuje zasadniczo (n) mnożeń i (n) dodawań. To ważne rozróżnienie: najlepsza postać do analizy nie musi być najlepszą postacią do obliczeń.

Praktyczna kolejność działań powinna wyglądać następująco:

1. Zidentyfikuj fragment pochłaniający najwięcej czasu lub pamięci za pomocą profilera.
2. Sprawdź, czy te same dane lub podwyrażenia są przetwarzane wielokrotnie.
3. Wybierz rozkład odpowiadający właściwościom problemu, a nie najpopularniejszą metodę.
4. Porównaj czas, pamięć i błąd numeryczny przed oraz po zmianie.
5. Dopiero na końcu upraszczaj pojedyncze wzory i mikrooptymalizuj kod.

Bez pomiarów łatwo zoptymalizować fragment odpowiadający za 1% czasu pracy programu, jednocześnie pozostawiając nietkniętą operację, która zajmuje pozostałe 99%.

FAQ

Czy faktoryzacja zawsze zmniejsza liczbę obliczeń?
Nie. Zysk pojawia się wtedy, gdy rozkład usuwa powtarzające się działania, pozwala zastosować szybszy algorytm albo jest wykorzystywany wielokrotnie. Jednorazowe rozłożenie prostego wyrażenia może kosztować więcej, niż później zaoszczędzi.

Czym różni się faktoryzacja liczby od faktoryzacji macierzy?
Liczbę całkowitą rozkłada się na iloczyn liczb pierwszych. Macierz rozkłada się na iloczyn macierzy o wygodniejszych właściwościach, na przykład trójkątnych, ortogonalnych lub diagonalnych. Cele i algorytmy są inne.

Kiedy wybrać rozkład LU, a kiedy Cholesky’ego?
Cholesky jest priorytetem dla macierzy symetrycznej dodatnio określonej, ponieważ zużywa mniej operacji i pamięci. W ogólnym układzie kwadratowym wybiera się LU z częściowym wyborem elementu głównego.

Czy do rozwiązania układu równań trzeba obliczać macierz odwrotną?
Zwykle nie. Wyznaczenie (A^{-1}), a następnie obliczenie (x=A^{-1}b), jest przeważnie wolniejsze i mniej stabilne niż bezpośrednie rozwiązanie układu na podstawie rozkładu LU, Cholesky’ego lub QR.

Czy rozkład na czynniki może pogorszyć dokładność wyniku?
Tak. Zmiana kolejności działań wpływa na błędy zaokrągleń w arytmetyce zmiennoprzecinkowej. Każdą optymalizację numeryczną trzeba sprawdzić na danych typowych, granicznych i źle uwarunkowanych.

Czy komputer kwantowy może szybko rozkładać duże liczby?
Algorytm Shora daje teoretyczną możliwość znacznego przyspieszenia faktoryzacji na odpowiednio dużym, odpornym na błędy komputerze kwantowym. Nie oznacza to, że współczesny użytkownik może dziś rozłożyć dowolny klucz RSA za pomocą dostępnej usługi chmurowej.

Od czego zacząć optymalizację obliczeń?
Nie od ręcznego upraszczania wzorów. Najpierw zmierz czas i pamięć, znajdź powtarzany koszt oraz sprawdź strukturę danych. Pierwszym błędem do usunięcia jest ponowne wykonywanie tej samej kosztownej faktoryzacji — zwłaszcza rozkładu macierzy — dla danych, które się nie zmieniły.