Участник:N Zakharov/Метод Ньютона для решения систем нелинейных уравнений

Основные авторы описания: Александр Чернышев, Никита Захаров. Вклад во все разделы равнозначный.

Метод Ньютона
Последовательный алгоритм
Последовательная сложность	[math]O(kn^3)[/math]
Объём входных данных	[math]2n+1[/math]
Объём выходных данных	[math]n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]000[/math]
Ширина ярусно-параллельной формы	[math]000[/math]

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

Метод Ньютона, алгоритм Ньютона — это итерационный численный метод нахождения корня (нуля) заданной функции. Метод был впервые предложен английским физиком, математиком и астрономом Исааком Ньютоном (1643—1727). Поиск решения осуществляется путём построения последовательных приближений и основан на принципах простой итерации. Метод обладает квадратичной сходимостью. Модификацией метода является метод хорд и касательных. Также метод Ньютона может быть использован для решения задач оптимизации, в которых требуется определить ноль первой производной либо градиента в случае многомерного пространства.

1.1.1 Историческая справка

Метод был описан Исааком Ньютоном в рукописи «Об анализе уравнениями бесконечных рядов» (лат. «De analysi per aequationes numero terminorum infinitas»), адресованной в 1669 году Барроу, и в работе «Метод флюксий и бесконечные ряды» (лат. «De metodis fluxionum et serierum infinitarum») или «Аналитическая геометрия» (лат. «Geometria analytica») в собраниях трудов Ньютона, которая была написана в 1671 году. В своих работах Ньютон вводит такие понятия, как разложение функции в ряд, бесконечно малые и флюксии (производные в нынешнем понимании). Указанные работы были изданы значительно позднее: первая вышла в свет в 1711 году благодаря Уильяму Джонсону, вторая была издана Джоном Кользоном в 1736 году уже после смерти создателя. Однако описание метода существенно отличалось от его нынешнего изложения: Ньютон применял свой метод исключительно к полиномам. Он вычислял не последовательные приближения [math]x_{n}[/math], а последовательность полиномов и в результате получал приближённое решение [math]x[/math].

Впервые метод был опубликован в трактате «Алгебра» Джона Валлиса в 1685 году, по просьбе которого он был кратко описан самим Ньютоном. В 1690 году Джозеф Рафсон опубликовал упрощённое описание в работе «Общий анализ уравнений» (лат. «Analysis aequationum universalis»). Рафсон рассматривал метод Ньютона как чисто алгебраический и ограничил его применение полиномами, однако при этом он описал метод на основе последовательных приближений [math]x_{n}[/math] вместо более трудной для понимания последовательности полиномов, использованной Ньютоном. Наконец, в 1740 году метод Ньютона был описан Томасом Симпсоном как итеративный метод первого порядка решения нелинейных уравнений. В той же публикации Симпсон обобщил метод на случай системы из двух уравнений и отметил, что метод Ньютона также может быть применён для решения задач оптимизации путём нахождения нуля производной или градиента.

В 1879 году Артур Кэли в работе «Проблема комплексных чисел Ньютона — Фурье» (англ. «The Newton-Fourier imaginary problem») был первым, кто отметил трудности в обобщении метода Ньютона на случай мнимых корней полиномов степени выше второй и комплексных начальных приближений. Эта работа открыла путь к изучению теории фракталов.

1.2 Математическое описание алгоритма

В этой статье мы сразу сформулируем задачу для многомерного случая и не будем останавливаться на методе Ньютона для одномерных задач .

Пусть необходимо найти решение системы:

[math] \left\{\begin{matrix}f_1(x_1, ..., x_n) = 0 \\ ... \\ f_n(x_1, ..., x_n) = 0, \end{matrix}\right. [/math]

где [math]f = (f_1, ..., f_n) : \mathbb{R}^n \rightarrow \mathbb{R}^n[/math] - непрерывно дифференцируемое отображение в окрестности решения.

При начальном приближении [math]x_0[/math] и сделанных предположениях [math](k+1)[/math]-я итерация метода будет выглядеть следующим образом: [math] x_{k+1} = x_k - [f'(x_k)]^{-1}f(x_k) [/math]

Основная сложность на каждой итерации заключается в решении системы уравнений. Один из самых эффективных методов их решения базируется на идеях Крылова, а именно проекционных итерационных методах. Они основываются на тех же принципах, что и методы, которые используются для задач поиска собственных значений. На самом деле, многие методы Крылова напрямую вытекают из алгоритмов Арнольди. Приближенное решение можно найти по формуле [math]\tilde{x} = x_0 + V_my[/math], а после подставления условий Петрова-Гарелкина оно принимает вид [math]\tilde{x} = x_0 + V (W^TAV)^{-1}W^Tr_0[/math], где [math]r_0 = b - Ax_0[/math], [math]V = \{v_1, v_2, ..., v_m\}[/math] - матрица, чьи столбцы являются ортонормированным базисом пространства, а матрица [math]W[/math] находится из соотношения [math]W^HV = I[/math].
[math]V^T_mAV_m = H_m[/math]
[math]G(y) = \left \| b-Ax \right \| = \left \| b-A(x_0+V_my) \right \| = \left \| r_0 -AV_my \right \|[/math]
Решение задачи минимизации
Наша цель состоит в нахождении вектора [math]y_n[/math], минимизирующего задачу [math] \| \tilde{H}_n y_n - \beta e_1 \|, \, [/math] где [math]\tilde{H}_n[/math] матрица размера (n+1) на n.
В результате применения QR-декомпозиции к матрице [math]\tilde{H}_n[/math] находим ортогональную матрицу [math]\Omega_n[/math] размера (n+1) на (n+1) и верхнюю треугольную матрицу [math]\tilde{R}_n[/math] размера (n+1) на n.
При этом в матрице [math]\tilde{R}_n[/math] последняя строка будет нулевой, то есть ее можно представить в следующем виде

[math] \tilde{R}_n = \begin{bmatrix} R_n \\ 0 \end{bmatrix}. [/math]

Перепишем нашу исходную задачу [math] \| \tilde{H}_n y_n - \beta e_1 \| = \| \Omega_n (\tilde{H}_n y_n - \beta e_1) \| = \| \tilde{R}_n y_n - \beta \Omega_n e_1 \|. [/math]

Положим [math]\beta\Omega_ne_1=\tilde{g}_n[/math],

где [math] \tilde{g}_n = \begin{bmatrix} g_n \\ \gamma_n \end{bmatrix} [/math]

Вектор [math]y_n[/math] находится из следующего соотношения [math] y_n = R_n^{-1} g_n. [/math]

1.3 Вычислительное ядро алгоритма

Описание алгоритма Арнольди
В результате мы получаем матрицу [math]H^{{m+1}\times{m}}[/math]
Пусть [math]r_0[/math] - произвольный входной вектор [math]v_1[/math] = [math]r_0/\begin{Vmatrix}r_0\end{Vmatrix}[/math]

 Для [math]m = 1,2,...,k [/math]
     [math]\omega_m[/math] = [math]Av_k[/math]
     Для [math]j = 1,2, ..,m [/math]
         [math]h_{jm}=(\omega_m, v_j)[/math]
         [math]\omega_m=\omega_m-h_{jm}v_j[/math]
     Конец
     [math]h_{m+1,m}=\begin{Vmatrix}\omega_m\end{Vmatrix}[/math]
     [math]v_{m+1}=\omega_m/h_{m+1,m}[/math]
 Конец

Не будем подробно останавливаться на описании решения задачи оптимизации.
Однако, важно отметить, что расчет матриц [math]\Omega_m, R_m[/math] и вектора [math]g_m[/math] можно упростить, используя формулы перехода от m к m+1 .
Модифицируя алгоритм Арнольди соответствующим образом, мы можем ощутимо снизить вычислительные затраты, необходимые нахождения QR-разложения матрицы [math]H_m[/math] и вектора [math]g_m[/math].

В результате, основной интерес для нас представляет решение СЛАУ [math]R_my_m=g_m[/math]. Так как матрица [math]R_m[/math] является верхнетреугольной, то вычисление вектора [math]y_m[/math] совпадает с обратным ходом метода Гаусса.

1.4 Макроструктура алгоритма

В данном алгоритме можно выделить 2 базовых макрооперации:

• Умножение матрицы на вектор

• Скалярное произведение векторов

1.5 Схема реализации последовательного алгоритма

 Выбор начального приближения x_k
 Вычисление вектора F(x_k)
 WHILE (условие остановки не выполнено)
   Вычисление якобиана J(x_k) = F'(x_k)
   Решение линейной системы уравнений алгоритмом GMRES J(x_k)s_k = -F(x_k)
   x_k = x_k + s_k
   Вычисление вектора F(x_k)

1.6 Последовательная сложность алгоритма

Будем оценивать сложность одной итерации метода Ньютона.
Для начала оценим сложность решения СЛАУ
Будем считать, что наша исходная матрица была порядка n. Предположим, что мы нашли решение СЛАУ за [math]k[/math] итераций. В результате выполнения алгоритма Арнольди мы получаем матрицу [math]H_k[/math] размера [math](k+1)*k[/math] и матрицу [math]Q_k[/math] размера [math]n*k[/math]. QR-факторизация матрицы [math]H_k[/math] будет стоить нам порядка [math]O(k^2)[/math] арифметических операций. Сложность алгоритма Арнольди при этом составляет [math]O(k*n)[/math] арифметических операций плюс сложность [math]k[/math] раз умножения матрицы на вектор - [math]O(k*n^2)[/math].
При решении задачи минимизации нам необходимо:
1) Решить СЛАУ с верхнетреугольной матрицей порядка [math]k[/math] для нахождения [math]y_k[/math] - сложность [math]O(k^2)[/math] арифметических операций.
2) Умножить матрицу [math]Q_k[/math] на найденный вектор - сложность [math]O(k * n)[/math] арифметических операций.
3) Сложить 2 вектора длины n - [math]O(n).[/math] Сложность оставшихся операций из алгоритма Ньютона (сложение векторов + подсчет нормы) составляет [math]O(2*n).[/math]
В результате мы получаем, что на одну итерацию метода Ньютона нужно [math]O(k*n^2)[/math] операций.
Отметим, что данная оценка справедлива при небольших [math]k[/math] относительно [math]n[/math].

1.7 Информационный граф

Рис. 1. Общая схема

Рис. 2. Итерация

1.8 Ресурс параллелизма алгоритма

Вообще говоря, метод Ньютона является последовательным методом, но на определенных шагах, таких как решение СЛАУ, приближенное вычисление Якобиана, можно применить идеи параллельных вычислений.

Основным ресурсом параллелизма в представленном алгоритме являются: умножение матрицы на вектор, скалярное произведение и элементарные операции над векторами.

1.9 Входные и выходные данные алгоритма

Входные данные:
[math]n^2 + n[/math] функций от [math]n[/math] переменных
[math]n[/math] вещественных чисел
Выходные данные:
[math]n[/math] вещественных чисел

1.10 Свойства алгоритма

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

2 Программная реализация алгоритма

2.7 Существующие реализации алгоритма

Из существующих реализаций важно отметить пакет PETSc и его компонент SNES, позволяющий решать системы нелинейных уравнений методом Ньютона
https://www.mcs.anl.gov/petsc/

3 Литература

<references \>