Wie geht man mit Komplexität im numerischen Code um, wenn man zum Beispiel mit großen Jacobi-Matrizen umgeht?

Ich löse ein nichtlineares System gekoppelter Gleichungen und habe den Jacobi des diskretisierten Systems berechnet. Das Ergebnis ist wirklich kompliziert, unten sind (nur!) Die ersten 3 Spalten einer $3\times 9$ Matrix,

(Die Komplexität entsteht teilweise, weil das numerische Schema eine exponentielle Anpassung für die Stabilität erfordert.)

Ich habe eine ganz allgemeine Frage zur Implementierung numerischer Codes mit Jacobianern.

Ich kann fortfahren und diese Matrix in Code implementieren. Aber meine Intuition sagt mir, dass ich aufgrund der Komplexität und der Unvermeidbarkeit der Einführung von Fehlern einige Tage (vielleicht Wochen!) Langwieriges Debuggen erwarten muss. Wie geht man mit einer solchen Komplexität im numerischen Code um, sie scheint unvermeidlich?! Verwenden Sie die automatische Codegenerierung aus symbolischen Paketen (und optimieren Sie den Code dann von Hand)?

Zuerst habe ich vor, den analytischen Jacobian mit einer Finite-Differenzen-Näherung zu debuggen. Sollte ich mir irgendwelche Fallstricke bewusst sein? Wie gehen Sie mit ähnlichen Problemen in Ihrem Code um?

Aktualisieren

Ich codiere dies in Python und habe Sympy verwendet , um den Jacobian zu generieren. Vielleicht kann ich die Codegenerierungsfunktion verwenden ?

Ein Wort: Modularität .

Es gibt viele wiederholte Ausdrücke in Ihrem Jacobian, die als ihre eigene Funktion geschrieben werden könnten. Es gibt keinen Grund für Sie, denselben Vorgang mehrmals zu schreiben, und das erleichtert das Debuggen. Wenn Sie es nur einmal schreiben, gibt es (theoretisch) nur einen Ort für einen Fehler.

Modularer Code erleichtert auch das Testen. Sie können Tests für jede Komponente Ihres Jacobian schreiben, anstatt zu versuchen, die gesamte Matrix zu testen. Wenn Sie beispielsweise Ihre Funktion am () modular schreiben, können Sie problemlos Vernunftstests dafür schreiben, prüfen, ob Sie sie richtig differenzieren usw.

Ein weiterer Vorschlag wäre, sich automatische Differenzierungsbibliotheken für die Zusammenstellung des Jacobian anzusehen. Es gibt keine Garantie, dass sie fehlerfrei sind, aber es wird wahrscheinlich weniger Debugging / weniger Fehler geben als beim Schreiben Ihrer eigenen. Hier sind einige, die Sie sich ansehen möchten:

• Sacado (Sandia Labs)
• ADIC (Argonne)

Entschuldigung, ich habe gerade gesehen, dass Sie Python verwenden. ScientificPython unterstützt AD.

• ScientificPython

Lassen Sie mich hier mit ein paar vorsichtigen Worten abwägen, denen eine Geschichte vorangestellt ist. Vor langer Zeit habe ich mit einem Kollegen gearbeitet, als ich gerade anfing. Er hatte ein Optimierungsproblem zu lösen, mit einem ziemlich chaotischen Ziel. Seine Lösung bestand darin, die analytischen Derivate für eine Optimierung zu generieren.

Das Problem, das ich sah, war, dass diese Derivate böse waren. Mit Macsyma generiert und in Fortran-Code konvertiert, waren sie jeweils Dutzende von Fortsetzungsanweisungen lang. Tatsächlich war der Fortran-Compiler darüber verärgert, da er die maximale Anzahl von Fortsetzungsanweisungen überschritten hatte. Während wir eine Flagge fanden, die es uns ermöglichte, dieses Problem zu umgehen, gab es andere Probleme.

• In langen Ausdrücken, wie sie üblicherweise von CA-Systemen erzeugt werden, besteht das Risiko einer massiven subtraktiven Stornierung. Berechnen Sie viele große Zahlen, nur um festzustellen, dass sich alle gegenseitig aufheben, um eine kleine Zahl zu erhalten.

• Oft sind analytisch erzeugte Derivate tatsächlich teurer zu bewerten als numerisch erzeugte Derivate unter Verwendung endlicher Differenzen. Ein Gradient für n Variablen kann mehr als das N-fache der Kosten für die Bewertung Ihrer Zielfunktion in Anspruch nehmen. (Möglicherweise können Sie Zeit sparen, da viele der Begriffe in den verschiedenen Ableitungen wiederverwendet werden können. Dies zwingt Sie jedoch auch dazu, eine sorgfältige Handcodierung durchzuführen, anstatt computergenerierte Ausdrücke zu verwenden. Und jedes Mal, wenn Sie böse mathematische Begriffe codieren Ausdrücke, die Wahrscheinlichkeit eines Fehlers ist nicht trivial. Stellen Sie sicher, dass Sie diese Ableitungen auf Richtigkeit überprüfen.)

Der Punkt meiner Geschichte ist, dass diese CA-generierten Ausdrücke ihre eigenen Probleme haben. Das Lustige ist, dass mein Kollege tatsächlich stolz auf die Komplexität des Problems war, dass er eindeutig ein wirklich schwieriges Problem löste, weil die Algebra so böse war. Ich glaube nicht, dass er darüber nachgedacht hat, ob diese Algebra tatsächlich das Richtige berechnet hat, ob sie so genau funktioniert und ob sie so effizient funktioniert.

Wäre ich zu diesem Zeitpunkt die leitende Person bei diesem Projekt gewesen, hätte ich ihm den Aufruhr vorgelesen. Sein Stolz veranlasste ihn, eine Lösung zu verwenden, die wahrscheinlich unnötig komplex war, ohne zu überprüfen, ob ein Gradient auf der Basis endlicher Differenzen angemessen war. Ich wette, wir hatten vielleicht eine Mannwoche Zeit damit verbracht, diese Optimierung zum Laufen zu bringen. Zumindest hätte ich ihm geraten, den erzeugten Gradienten sorgfältig zu testen. War es genau? Wie genau war es im Vergleich zu Finite-Differenzen-Derivaten? Tatsächlich gibt es heutzutage Tools, die auch eine Schätzung des Fehlers in ihrer abgeleiteten Vorhersage zurückgeben. Dies gilt sicherlich für den adaptiven Differenzierungscode (Derivat), den ich in MATLAB geschrieben habe.

Testen Sie den Code. Überprüfen Sie die Derivate.

Bevor Sie dies jedoch tun, sollten Sie überlegen, ob andere, bessere Optimierungsschemata in Frage kommen. Wenn Sie beispielsweise eine Exponentialanpassung durchführen, besteht eine sehr gute Chance, dass Sie partitionierte nichtlineare kleinste Quadrate verwenden (manchmal als trennbare kleinste Quadrate bezeichnet. Ich denke, das war der Begriff, den Seber und Wild in ihrem Buch verwendet haben.) Die Idee besteht darin, den Parametersatz in intrinsisch lineare und intrinsisch nichtlineare Mengen zu unterteilen. Verwenden Sie eine Optimierung, die nur für die nichtlinearen Parameter funktioniert. Wenn diese Parameter "bekannt" sind, können die intrinsisch linearen Parameter unter Verwendung einfacher linearer kleinster Quadrate geschätzt werden. Dieses Schema reduziert den Parameterraum in der Optimierung. Dies macht das Problem robuster, da Sie keine Startwerte für die linearen Parameter finden müssen. Dadurch wird die Dimensionalität Ihres Suchraums verringert, sodass das Problem schneller ausgeführt wird. Wieder habe ich geliefertein Werkzeug für diesen Zweck , aber nur in MATLAB.

Wenn Sie die analytischen Derivate verwenden, codieren Sie sie, um Begriffe wiederzuverwenden. Dies kann eine erhebliche Zeitersparnis bedeuten und die Fehler tatsächlich reduzieren, wodurch Sie Ihre eigene Zeit sparen. Aber dann überprüfen Sie diese Zahlen!

Es sind verschiedene Strategien zu berücksichtigen:

1. Suchen Sie die Ableitungen in symbolischer Form mithilfe eines CAS und exportieren Sie den Code zur Berechnung der Ableitungen.

2. Verwenden Sie ein automatisches Differenzierungswerkzeug (AD), um Code zu erstellen, der die Ableitungen vom Code berechnet, um die Funktionen zu berechnen.

3. Verwenden Sie endliche Differenznäherungen, um den Jacobi zu approximieren.

Die automatische Differenzierung könnte einen effizienteren Code für die Berechnung des gesamten Jacobi erzeugen als die symbolische Berechnung, um eine Formel für jeden Eintrag in der Matrix zu erstellen. Endliche Unterschiede sind eine gute Möglichkeit, Ihre Derivate zu überprüfen.

Hier ist ein Beispiel dafür, wo wir die automatische Differenzierung mit Sacado in einem Code verwendet haben: http://www.dealii.org/developer/doxygen/deal.II/step_33.html

Neben den hervorragenden Vorschlägen von BrianBorcher besteht ein weiterer möglicher Ansatz für realwertige Funktionen darin, die Ableitung der Näherung in komplexen Schritten zu verwenden (siehe diesen Artikel (paywalled) und diesen Artikel ). In einigen Fällen liefert dieser Ansatz genauere numerische Ableitungen auf Kosten der Änderung der Werte von Variablen in Ihrer Funktion von real zu komplex. Der zweite Artikel listet einige Fälle auf, in denen die komplexe Schrittfunktionsnäherung zusammenbrechen könnte.