Jak radzić sobie ze złożonością kodu numerycznego, na przykład w przypadku dużych macierzy jakobskich?

Rozwiązuję nieliniowy układ sprzężonych równań i obliczyłem jakobian układu dyskretnego. Wynik jest naprawdę skomplikowany, poniżej są (tylko!) Pierwsze 3 kolumny macierzy $3\times 9$ ,

(Złożoność wynika częściowo z faktu, że schemat numeryczny wymaga wykładniczego dopasowania dla stabilności).

Mam dość ogólne pytanie dotyczące implementacji kodów numerycznych z wykorzystaniem jakobianów.

Mogę śmiało zaimplementować tę matrycę w kodzie. Ale moja intuicja podpowiada mi, żebym oczekiwał kilku dni (może tygodni!) Żmudnego debugowania ze względu na samą złożoność i niemożność uniknięcia błędów. Jak radzić sobie ze złożonością, taką jak ta w kodzie numerycznym, wydaje się to nieuniknione ?! Czy korzystasz z automatycznego generowania kodu z pakietów symbolicznych (a następnie dostosowujesz kod ręcznie)?

Najpierw planuję debugować analityczne jakobian z przybliżoną różnicą skończoną, czy powinienem być świadomy jakichkolwiek pułapek? Jak radzisz sobie z podobnymi problemami w kodzie?

Aktualizacja

Koduję to w Pythonie i użyłem sympy do wygenerowania jakobianu. Może mogę skorzystać z funkcji generowania kodu ?

Jedno słowo: modułowość .

W twoim jakobińskim jest wiele powtarzających się wyrażeń, które można zapisać jako ich własną funkcję. Nie ma powodu, aby pisać tę samą operację więcej niż raz, a to ułatwi debugowanie; jeśli napiszesz to tylko raz, jest tylko jedno miejsce na błąd (teoretycznie).

Kod modułowy również ułatwi testowanie; Możesz pisać testy dla każdego komponentu swojego jakobiańskiego, zamiast próbować przetestować całą matrycę. Na przykład, jeśli piszesz swoją funkcję am () w sposób modułowy, możesz łatwo napisać dla niej testy poczytalności, sprawdzić, czy odpowiednio ją różnicujesz itp.

Inną sugestią byłoby przyjrzenie się automatycznym bibliotekom różnicowania do składania jakobianów. Nie ma gwarancji, że są one wolne od błędów, ale prawdopodobnie będzie mniej błędów debugowania / mniej błędów niż pisanie własnych. Oto kilka, na które warto spojrzeć:

• Sacado (Sandia Labs)
• ADIC (Argonne)

Przepraszamy, właśnie zobaczyłem, że używasz Pythona. ScientificPython obsługuje AD.

• ScientificPython

Pozwólcie mi ważyć się tutaj z kilkoma słowami ostrożności, poprzedzonymi historią. Dawno temu pracowałem z facetem, kiedy dopiero zaczynałem. Miał problem z optymalizacją do rozwiązania, z dość niechlujnym celem. Jego rozwiązaniem było wygenerowanie analitycznych pochodnych do optymalizacji.

Problem, który widziałem, był taki paskudny. Każdy z nich został wygenerowany za pomocą Macsyma i przekonwertowany na kod fortran. W rzeczywistości kompilator Fortran był tym zdenerwowany, ponieważ przekroczył maksymalną liczbę instrukcji kontynuacji. Chociaż znaleźliśmy flagę, która pozwoliła nam obejść ten problem, były też inne problemy.

• W długich wyrażeniach, które są zwykle generowane przez systemy CA, istnieje ryzyko masywnego anulowania odejmowania. Oblicz wiele dużych liczb, ale okazuje się, że wszystkie się znoszą, dając niewielką liczbę.

• Często analitycznie generowane pochodne są w rzeczywistości bardziej kosztowne do oszacowania niż pochodne generowane numerycznie z wykorzystaniem różnic skończonych. Gradient dla n zmiennych może zająć ponad n razy więcej niż koszt oceny funkcji celu. (Być może będziesz w stanie zaoszczędzić trochę czasu, ponieważ wiele terminów może być ponownie używanych w różnych pochodnych, ale to również zmusi cię do ostrożnego kodowania ręcznego, zamiast używania wyrażeń generowanych komputerowo. I za każdym razem, gdy kodujesz nieprzyjemne matematyczne wyrażeń, prawdopodobieństwo błędu nie jest trywialne. Upewnij się, że weryfikujesz pochodne pod względem dokładności).

Chodzi o to, że te wyrażenia generowane przez CA mają własne problemy. Zabawne jest to, że mój kolega był dumny ze złożoności problemu, że najwyraźniej rozwiązał naprawdę trudny problem, ponieważ algebra była tak paskudna. Nie sądzę, by zastanawiał się, czy ta algebra rzeczywiście oblicza prawidłową rzecz, czy robi to tak dokładnie i czy robi to tak skutecznie.

Gdybym był wówczas osobą starszą w tym projekcie, przeczytałbym mu akt zamieszek. Jego duma sprawiła, że ​​zastosował rozwiązanie, które prawdopodobnie było niepotrzebnie złożone, nawet nie sprawdzając, czy gradient oparty na skończonej różnicy jest odpowiedni. Założę się, że spędziliśmy może tydzień pracy na optymalizację. Przynajmniej doradziłbym mu, aby dokładnie przetestował wytworzony gradient. Czy to było dokładne? Jaka była dokładność w porównaniu do pochodnych różnic skończonych? W rzeczywistości istnieją dziś narzędzia, które również zwrócą oszacowanie błędu w ich przewidywaniu pochodnych. Jest to z pewnością prawda w przypadku adaptacyjnego kodu różnicującego (wyprowadza się) , który napisałem w MATLAB.

Przetestuj kod. Sprawdź pochodne.

Ale zanim zrobisz KAŻDĄ z tych czynności, zastanów się, czy możliwe są inne, lepsze schematy optymalizacji. Na przykład, jeśli wykonujesz dopasowanie wykładnicze, istnieje bardzo duża szansa, że ​​możesz użyć podzielonego nieliniowego najmniejszego kwadratu (czasami nazywanego separowalnym najmniejszym kwadratem. Myślę, że to był termin użyty przez Sebera i Wilda w ich książce.) Pomysł polega na rozbiciu zestawu parametrów na wewnętrznie liniowe i wewnętrznie nieliniowe zestawy. Użyj optymalizacji, która działa tylko na parametrach nieliniowych. Biorąc pod uwagę, że parametry te są „znane”, wówczas parametry wewnętrznie liniowe można oszacować za pomocą prostych liniowych najmniejszych kwadratów. Ten schemat zmniejszy przestrzeń parametrów w optymalizacji. Sprawia, że ​​problem staje się bardziej niezawodny, ponieważ nie trzeba znaleźć wartości początkowych dla parametrów liniowych. Zmniejsza wymiary przestrzeni wyszukiwania, dzięki czemu problem działa szybciej. Znowu dostarczyłemnarzędzie do tego celu , ale tylko w MATLAB.

Jeśli korzystasz z analitycznych pochodnych, koduj je, aby ponownie wykorzystać warunki. Może to być poważną oszczędnością czasu i może faktycznie zmniejszyć liczbę błędów, oszczędzając Twój czas. Ale sprawdź te liczby!

Należy rozważyć kilka strategii:

1. Znajdź pochodne w formie symbolicznej za pomocą CAS, a następnie wyeksportuj kod do obliczenia pochodnych.

2. Użyj narzędzia automatycznego różnicowania (AD), aby utworzyć kod, który oblicza pochodne z kodu w celu obliczenia funkcji.

3. Użyj przybliżonych różnic skończonych, aby przybliżyć przybliżenie jakobianów.

Automatyczne różnicowanie może wytworzyć bardziej wydajny kod do obliczania całego jakobianu, a następnie wykorzystanie obliczeń symbolicznych do utworzenia formuły dla każdego wpisu w macierzy. Skończone różnice są dobrym sposobem na podwójne sprawdzenie swoich instrumentów pochodnych.

Oto przykład zastosowania automatycznego różnicowania za pomocą Sacado w jednym kodzie: http://www.dealii.org/developer/doxygen/deal.II/step_33.html

Oprócz doskonałych sugestii BrianBorchera, innym możliwym podejściem do funkcji o wartościach rzeczywistych jest zastosowanie przybliżenia pochodnej złożonego kroku (zobacz ten artykuł (paywalled) i ten artykuł ). W niektórych przypadkach takie podejście daje dokładniejsze numeryczne pochodne kosztem zmiany wartości zmiennych w funkcji z rzeczywistej na złożoną. W drugim artykule wymieniono niektóre przypadki, w których przybliżenie funkcji kroku złożonego może się załamać.