Comment gérer la complexité du code numérique, par exemple, lorsqu'il s'agit de grandes matrices jacobiennes?

Je suis en train de résoudre un système non linéaire d'équations couplées et j'ai calculé le jacobien du système discrétisé. Le résultat est vraiment compliqué, voici (seulement!) Les 3 premières colonnes d'une matrice ,$3\times 9$

(La complexité se pose, en partie, parce que le schéma numérique nécessite un ajustement exponentiel pour la stabilité.)

J'ai une question assez générale sur l'implémentation de codes numériques utilisant des jacobiens.

Je peux aller de l'avant et implémenter cette matrice dans le code. Mais mon intuition me dit de m'attendre à quelques jours (peut-être des semaines!) De débogage fastidieux en raison de la complexité et de l'inévitabilité de l'introduction d'erreurs. Comment faire face à une complexité comme celle-ci en code numérique, cela semble incontournable?! Utilisez-vous la génération automatique de code à partir de packages symboliques (puis ajustez le code à la main)?

Je prévois d'abord de déboguer le jacobien analytique avec une approximation de différence finie, dois-je être au courant des pièges? Comment traitez-vous des problèmes similaires dans votre code?

Mise à jour

Je code cela en Python et j'ai utilisé sympy pour générer le jacobien. Peut-être que je peux utiliser la fonction de génération de code ?

Un mot: Modularité .

Il y a beaucoup d'expressions répétées dans votre jacobien qui pourraient être écrites comme leur propre fonction. Il n'y a aucune raison pour que vous écriviez la même opération plus d'une fois, et cela facilitera le débogage; si vous ne l'écrivez qu'une seule fois, il n'y a qu'un seul endroit pour une erreur (en théorie).

Le code modulaire facilitera également les tests; Vous pouvez écrire des tests pour chaque composant de votre jacobien au lieu d'essayer de tester la matrice entière. Par exemple, si vous écrivez votre fonction am () de manière modulaire, vous pouvez facilement écrire des tests d'intégrité pour elle, vérifier si vous la différenciez correctement, etc.

Une autre suggestion serait d'examiner les bibliothèques de différenciation automatique pour assembler le jacobien. Il n'y a aucune garantie qu'ils sont exempts d'erreurs, mais il y aura probablement moins de débogage / moins d'erreurs que d'écrire les vôtres. Voici quelques-uns que vous voudrez peut-être regarder:

• Sacado (Sandia Labs)
• ADIC (Argonne)

Désolé, je viens de voir que vous utilisez python. ScientificPython prend en charge AD.

• ScientificPython

Permettez-moi de peser ici avec quelques mots de prudence, précédés d'une histoire. Il y a longtemps, je travaillais avec un gars quand je débutais. Il avait un problème d'optimisation à résoudre, avec un objectif plutôt salissant. Sa solution a été de générer les dérivés analytiques pour une optimisation.

Le problème que j'ai vu était que ces dérivés étaient désagréables. Générés à l'aide de Macsyma, convertis en code fortran, ils comportaient chacun des dizaines d'instructions de continuation. En fait, le compilateur Fortran s'est énervé car il dépassait le nombre maximal d'instructions de continuation. Alors que nous avons trouvé un drapeau qui nous a permis de contourner ce problème, il y avait d'autres problèmes.

• Dans les expressions longues, comme celles généralement générées par les systèmes CA, il existe un risque d'annulation massive soustractive. Calculez beaucoup de grands nombres, seulement pour découvrir qu'ils s'annulent tous pour donner un petit nombre.

• Souvent, les dérivés générés analytiquement sont en fait plus coûteux à évaluer que les dérivés générés numériquement utilisant des différences finies. Un gradient pour n variables peut prendre plus de n fois le coût de l'évaluation de votre fonction objectif. (Vous pourrez peut-être gagner du temps, car de nombreux termes peuvent être réutilisés dans les divers dérivés, mais cela vous obligera également à effectuer un codage manuel minutieux, au lieu d'utiliser des expressions générées par ordinateur. Et chaque fois que vous codez de mauvaises mathématiques expressions, la probabilité d'une erreur n'est pas anodine. Assurez-vous de vérifier l'exactitude de ces dérivés.)

Le point de mon histoire est que ces expressions générées par CA ont leurs propres problèmes. Le plus drôle, c'est que mon collègue était en fait fier de la complexité du problème, qu'il résolvait clairement un problème vraiment difficile parce que l'algèbre était si méchante. Ce que je ne pense pas qu'il ait considéré, c'est si cette algèbre calculait réellement la bonne chose, si elle le faisait avec précision et si elle le faisait efficacement.

Si j'avais été la personne âgée à l'époque sur ce projet, je lui aurais lu l'acte anti-émeute. Sa fierté l'a amené à utiliser une solution qui était probablement inutilement complexe, sans même vérifier qu'un gradient basé sur les différences finies était adéquat. Je parie que nous avions peut-être passé une semaine-homme pour lancer cette optimisation. À tout le moins, je lui aurais conseillé de tester soigneusement le gradient produit. Était-ce exact? Quelle était sa précision par rapport aux dérivés à différences finies? En fait, il existe aujourd'hui des outils qui renverront également une estimation de l'erreur dans leur prédiction dérivée. C'est certainement vrai pour le code de différenciation adaptative (dérivé) que j'ai écrit dans MATLAB.

Testez le code. Vérifiez les dérivés.

Mais avant de faire TOUT cela, demandez-vous si d'autres, de meilleurs schémas d'optimisation sont une option. Par exemple, si vous faites un ajustement exponentiel, il y a de très bonnes chances que vous puissiez utiliser des moindres carrés non linéaires partitionnés (parfois appelés moindres carrés séparables. Je pense que c'était le terme utilisé par Seber et Wild dans leur livre.) L'idée consiste à décomposer l'ensemble des paramètres en ensembles intrinsèquement linéaires et intrinsèquement non linéaires. Utilisez une optimisation qui fonctionne uniquement sur les paramètres non linéaires. Étant donné que ces paramètres sont "connus", les paramètres intrinsèquement linéaires peuvent être estimés à l'aide de moindres carrés linéaires simples. Ce schéma réduira l'espace des paramètres dans l'optimisation. Cela rend le problème plus robuste, car vous n'avez pas besoin de trouver les valeurs de départ pour les paramètres linéaires. Il réduit la dimensionnalité de votre espace de recherche, ce qui accélère le problème. Encore une fois, j'ai fourniun outil à cet effet , mais uniquement dans MATLAB.

Si vous utilisez les dérivés analytiques, codez-les pour réutiliser les termes. Cela peut être un gain de temps considérable et peut en fait réduire les bugs, ce qui vous fait gagner du temps. Mais alors vérifiez ces chiffres!

Il y a plusieurs stratégies à considérer:

1. Trouvez les dérivés sous forme symbolique à l'aide d'un CAS, puis exportez le code pour calculer les dérivés.

2. Utilisez un outil de différenciation automatique (AD) pour produire du code qui calcule les dérivées du code pour calculer les fonctions.

3. Utilisez des approximations aux différences finies pour approximer le jacobien.

La différenciation automatique pourrait produire un code plus efficace pour calculer le jacobien entier puis utiliser le calcul symbolique pour produire une formule pour chaque entrée de la matrice. Les différences finies sont un bon moyen de vérifier vos dérivés.

Voici un exemple où nous avons utilisé la différenciation automatique en utilisant Sacado dans un code: http://www.dealii.org/developer/doxygen/deal.II/step_33.html

En plus des excellentes suggestions de BrianBorcher, une autre approche possible pour les fonctions à valeur réelle consiste à utiliser l'approximation dérivée à étapes complexes (voir cet article (paywalled) et cet article ). Dans certains cas, cette approche donne des dérivées numériques plus précises au prix de changer les valeurs des variables de votre fonction de réelles à complexes. Le deuxième article répertorie certains cas où l'approximation de la fonction d'étape complexe peut échouer.