Об автоматическом дифференцировании, методе Ньютона и решении СЛАУ на Delphi. Часть 1
Об автоматическом дифференцировании (АД) на Хабре уже писалось здесь и здесь. В данной статье предлагается реализация АД для Delphi (протестировано в Embarcadero XE2, XE6) вместе с удобными классами методов Ньютона для решения нелинейных уравнений f(x) = 0 и систем F(X) = 0. Любые ссылки на готовые аналогичные библиотеки приветствуются, сам же я подобного не нашел, не считая отличного решателя СЛАУ с разреженной матрицей (см. под катом). Давайте в самом начале отметим для себя, что выбор Delphi обусловлен legacy-кодом, тем не менее на C++ задачу можно решать следующим образом. Во-первых, для методов Ньютоновского (базовый метод Ньютона, метод Бройдена) типа имеются Numerical Recipes in C++. Ранее «Рецепты» были только на чистом C и приходилось делать свои классовые обертки. Во-вторых, можно взять одну из АД-библиотек из списка Autodiff.org. По моему опыту использования CPPAD быстрее FADBAD и Trilinos::Sacado примерно на 30%, самая же быстрая, судя по описанию, библиотека это новая ADEPT. В-третьих, для матрично-векторных операций можно взять проверенный временем uBlas , либо новые и быстрые конкуренты Armadillo и blaze-lib — это если для решения СЛАУ использовать отдельные библиотеки (например, SuiteSparce или Pardiso для прямых и ITL для итерационных методов). Привлекательным является также использование интегрированных библиотек линейной алгебры и решателей СЛАУ Eigen, MTL, PETSc (имеются также Ньютоновские решатели на C). Вся триада из заголовка полностью реализована, насколько мне известно, только в Trilinos.
О применении численного дифференцирования
Производные можно вычислять аналитически и численно. К аналитическим методам относятся ручное дифференцирование, символьное (Maple, Wolfram и т.п.) и непосредственно автоматическое дифференцирование, выраженное в средствах выбранного языка программирования.
Современный тренд на использование АД оправдан одной простой причиной — с помощью этой техники устраняется избыточность кода и его дублирование. Другой аргумент состоит в том, что, например, при решении нелинейных дифференциальных уравнений (систем) сеточными методами способ вычисления F(X) сам по себе является нетривиальной задачей. В реальных задачах невязка F(X) представлена суперпозицией вызовов функций с разных слоев программы и ручное дифференцирование теряет свою наглядность. Следует также отметить, что при моделировании нестационарных процессов F(X) меняется на каждом шаге по времени, также может меняться и сам вектор неизвестных X. Использование АД позволяет нам сконцентрироваться непосредственно на формировании F(X), однако не снимает вопрос о верификации получаемой матрицы Якобиана dF(X)/dX. Дело в том, что при вычислении невязок мы можем забыть изменить тип промежуточной переменной со стандартного double на тип АД (а многие библиотеки имеют неявное приведение типа АД к double), тем самым некоторые производные будут вычислены некорректно. Проблема в этом случае состоит в том, что даже при наличии ошибок в формулах для производных метод Ньютона может сходиться, хоть и за возросшее число итераций. Это может быть незаметно при одних начальных данных и приводить к расходимости процесса при других.
Таким образом, какой бы аналитической способ дифференцирования df/dx не был выбран, его крайне желательно дополнить сравнением с численным дифференцированием (f(x + h) — f(x)) / h, иначе всегда будут оставаться сомнения в правильности кода. Естественно, численные производные никогда не совпадут с точностью с правильными аналитическими, тем не менее можно порекомендовать следующий прием юнит-тестирования. Можно выгрузить в текстовые файлы вектора X, F(X) и матрицу dF(X)/dX и выложить на SVN. Затем получить dF(X)/dX численно и сохранить файл поверх существующего, после чего визуально сравнивать файлы между собой. Здесь Вы всегда увидите насколько поменялись значения и сможете локализовать координаты элементов матрицы с большими отклонениями (не в долях) — в этом месте находится ошибка аналитического дифференцирования.
Реализация АД
В Embarcadero Delphi до версии XE5 отсутствует возможность перегрузки арифметических операций для классов, но есть такая возможность для структур record (ссылка):
Обычно в АД на C++ размерность вектора производных является переменной величиной и инициализируется в конструкторе. В Delphi нельзя (есть попытки обойти) перегрузить оператор присваивания для структур и в связи с побитовым копированием вектор производных приходится делать фиксированной длины. Соответствующая константа TAutoDiff.n_all должна изначально задаваться вручную.
Пример 1На данный момент в библиотеке перегружены почти все бинарные и унарные операторы, за исключением тригонометрических функций и булевых операций сдвига. Недостающие функции легко доработать самостоятельно.
Реализация АД для разреженных матриц
С одной стороны фиксированное значение n_all это существенное ограничение, ведь размерность вектора может поступать извне. С другой стороны для некоторых задач его можно ослабить. Дело в том, что если говорить о численных методах решения уравнений механики сплошных сред, то возникающие в них матрицы имеют разреженную структуру — классический пример это «схема крест» для оператора Лапласа, когда в уравнении для узла с координатами (i, j) (ограничимся 2D) задействованы только 5 узлов: (i, j), (i-1, j), (i+1, j), (i, j-1), (i, j+1). Обобщая идею мы должны заложить следующее для данной конкретной задачи:
Пусть общее число узлов в решаемой задаче N. Тогда в матрице Якобиана N_all = N * n_junc_vars столбцов, из них ненулевых только n_all. Если завести теперь внутри структуры TAutoDiff целочисленный вектор n_juncs, каждый элемент которого определяет глобальный индекс узла от 0 до N-1, то функцию dx с локальным индексом из диапазона [0, n_all-1] можно дополнить функцией dx_global с уже глобальным индексом из диапазона [0, N_all-1]. Пример 2 иллюстрирует детали использования такого интерфейса, плюсы которого будут наиболее полно видны при реализации метода Ньютона ниже.