Как бороться со сложностью в числовом коде, например, при работе с большими якобиевыми матрицами?

Я решаю нелинейную систему связанных уравнений и вычисляю якобиан дискретной системы. Результат действительно сложный, ниже приведены (только!) Первые 3 столбца матрицы ,$3\times 9$

(Сложность возникает отчасти потому, что численная схема требует экспоненциального подбора для устойчивости.)

У меня довольно общий вопрос о реализации числовых кодов с использованием якобианов.

Я могу пойти дальше и реализовать эту матрицу в коде. Но моя интуиция подсказывает мне ожидать утомительной отладки в течение нескольких дней (возможно, недель!) Из-за явной сложности и неизбежности появления ошибок. Как справиться со сложностью, такой как это в числовом коде, это кажется неизбежным ?! Используете ли вы автоматическую генерацию кода из символических пакетов (затем настройте код вручную)?

Сначала я планирую отладить аналитический якобиан с конечно-разностным приближением. Должен ли я знать о каких-либо подводных камнях? Как вы решаете подобные проблемы в вашем коде?

Обновить

Я кодирую это на Python и использовал sympy для генерации якобиана. Может быть, я могу использовать функцию генерации кода ?

Одним словом: модульность .

В вашем якобиане много повторяющихся выражений, которые можно записать как их собственную функцию. У вас нет причин писать одну и ту же операцию более одного раза, и это облегчит отладку; если вы пишете только один раз, есть только одно место для ошибки (теоретически).

Модульный код также облегчит тестирование; Вы можете написать тесты для каждого компонента вашего якобиана, а не пытаться проверить всю матрицу. Например, если вы пишете свою функцию am () по модульному принципу, вы можете легко написать для нее тесты работоспособности, проверить, правильно ли вы ее дифференцируете, и т. Д.

Другим предложением было бы взглянуть на библиотеки автоматического дифференцирования для сборки якобиана. Там нет никакой гарантии, что они без ошибок, но, вероятно, будет меньше отладки / меньше ошибок, чем писать свои собственные. Вот некоторые из них, которые вы можете посмотреть на это:

• Sacado (Сандия Лабс)
• ADIC (Аргонн)

Извините, только что увидел, что вы используете Python. ScientificPython имеет поддержку AD.

• ScientificPython

Позвольте мне сказать несколько слов предостережения, предваряя рассказ. Давным-давно я работал с парнем, когда только начинал. У него была проблема с оптимизацией, с довольно грязной целью. Его решением было генерировать аналитические производные для оптимизации.

Проблема, которую я увидел, заключалась в том, что эти производные были неприятными. Созданные с использованием Macsyma, преобразованные в код на языке Fortran, каждый из них представлял собой десятки операторов продолжения. Фактически, компилятор Фортрана расстроился из-за того, что он превысил максимальное количество операторов продолжения. Хотя мы нашли флаг, который позволил нам обойти эту проблему, были и другие проблемы.

• В длинных выражениях, которые обычно генерируются системами CA, существует риск массивного вычитания. Вычислите много больших чисел, только чтобы найти, что все они взаимно компенсируют друг друга, чтобы получить небольшое число.

• Часто аналитически сгенерированные производные на самом деле более дороги для оценки, чем численно сгенерированные производные с использованием конечных разностей. Градиент для n переменных может занять более чем в n раз дороже оценки вашей целевой функции. (Вы можете сэкономить некоторое время, потому что многие термины могут быть повторно использованы в различных производных, но это также заставит вас делать осторожное ручное кодирование вместо использования сгенерированных компьютером выражений. И каждый раз, когда вы кодируете неприятные математические В выражениях вероятность ошибки не тривиальна. Убедитесь, что вы проверяете эти производные на точность.)

Суть моей истории в том, что эти сгенерированные CA выражения имеют свои собственные проблемы. Самое смешное, что мой коллега на самом деле гордился сложностью проблемы, потому что он явно решал действительно сложную задачу, потому что алгебра была настолько противной. Я не думаю, что он думал о том, что эта алгебра действительно вычисляет правильную вещь, делает ли она это точно, и делает ли она это эффективно.

Если бы я был старшим человеком в то время в этом проекте, я бы прочитал его акт бунта. Его гордость заставила его использовать решение, которое, вероятно, было излишне сложным, даже не проверяя, что градиент на основе конечных разностей был адекватным. Могу поспорить, что мы потратили, возможно, человеческую неделю, чтобы запустить эту оптимизацию. По крайней мере, я бы посоветовал ему тщательно проверить полученный градиент. Было ли это точно? Насколько это было точно по сравнению с производными с конечной разностью? Фактически, сегодня существуют инструменты, которые также будут возвращать оценку ошибки в их производном прогнозе. Это, безусловно, верно для адаптивного дифференцирующего кода (производного), который я написал в MATLAB.

Протестируйте код. Проверьте производные.

Но прежде чем делать ЛЮБОЕ из этого, подумайте, есть ли другие, более эффективные схемы оптимизации. Например, если вы делаете экспоненциальную подгонку, то есть очень хороший шанс, что вы можете использовать разделенные нелинейные наименьшие квадраты (иногда называемые сепарабельными наименьшими квадратами. Я думаю, что именно этот термин использовали Себер и Уайлд в своей книге.) Идея состоит в том, чтобы разбить набор параметров на внутренне линейные и внутренне нелинейные множества. Используйте оптимизацию, которая работает только с нелинейными параметрами. Если эти параметры «известны», то линейно линейные параметры могут быть оценены с использованием простых линейных наименьших квадратов. Эта схема уменьшит пространство параметров при оптимизации. Это делает проблему более устойчивой, поскольку вам не нужно находить начальные значения для линейных параметров. Это уменьшает размерность вашего пространства поиска, поэтому проблема запускается быстрее. Я снова поставилинструмент для этой цели , но только в MATLAB.

Если вы используете аналитические производные, используйте их для повторного использования терминов. Это может быть серьезной экономией времени и может фактически уменьшить количество ошибок, экономя ваше собственное время. Но тогда проверь эти цифры!

Есть несколько стратегий для рассмотрения:

1. Найдите производные в символической форме, используя CAS, затем экспортируйте код для вычисления производных.

2. Используйте инструмент автоматического дифференцирования (AD) для создания кода, который вычисляет производные от кода для вычисления функций.

3. Используйте конечно-разностные аппроксимации для аппроксимации якобиана.

Автоматическое дифференцирование может создать более эффективный код для вычисления всего якобиана, чем использование символьных вычислений для получения формулы для каждой записи в матрице. Конечные различия - хороший способ перепроверить ваши производные.

Вот пример, где мы использовали автоматическое дифференцирование с использованием Sacado в одном коде: http://www.dealii.org/developer/doxygen/deal.II/step_33.html

В дополнение к превосходным предложениям BrianBorcher, другой возможный подход для вещественных функций заключается в использовании аппроксимации производного сложного шага (см. Эту статью (paywalled) и эту статью ). В некоторых случаях этот подход дает более точные числовые производные за счет изменения значений переменных в вашей функции с реальных на сложные. Во второй статье перечислены некоторые случаи, когда приближение комплексной ступенчатой ​​функции может нарушиться.