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

Эта работа прошла предварительную проверку Дата последней правки страницы: 22.11.2016 Данная работа соответствует формальным критериям. Проверено VadimVV.

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

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]\frac{r_0}{\|r_0\|}[/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}=\|\omega_m\|[/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]nk[/math]. QR-факторизация матрицы [math]H_k[/math] будет стоить нам порядка [math]O(k^3)[/math] арифметических операций. Сложность алгоритма Арнольди при этом составляет [math]O(k^2n)[/math] арифметических операций плюс сложность [math]k[/math] раз умножения матрицы на вектор - [math]O(kn^2)[/math].
При решении задачи минимизации нам необходимо:
1) Решить СЛАУ с верхнетреугольной матрицей порядка [math]k[/math] для нахождения [math]y_k[/math] - сложность [math]O(k^2)[/math] арифметических операций.
2) Умножить матрицу [math]Q_k[/math] на найденный вектор - сложность [math]O(kn)[/math] арифметических операций.
3) Сложить 2 вектора длины n - [math]O(n).[/math] Предположим, что сложность расчета одного элемента вектор-функции F и матрицы якобиана J будет [math]O(1)[/math].

Тогда сложность расчета вектор-функции будет [math]O(n)[/math], матрицы якобиана - [math]O(n^2)[/math].

Сложность оставшихся операций из алгоритма Ньютона (сложение векторов + подсчет нормы) составляет [math]O(2n).[/math]
В результате мы получаем, что на одну итерацию метода Ньютона нужно [math]O(kn^2)[/math] операций.
Отметим, что данная оценка справедлива при небольших [math]k[/math] относительно [math]n[/math].

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

Рис. 1. Информационный граф

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

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

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

Пусть нам будет дана матрица порядка n. Будем предполагать, что доступно p процессоров.

Как и ранее будем предполагать, что сложность расчета одного элемента вектор-функции F и матрицы якобиана J будет [math]O(1)[/math].

При получении оценки последовательной сложности алгоритма было показано, что для одной итерации метода Ньютона необходимо выполнить:

На итерациях Арнольди:

  [math]k[/math] операций умножения матрицы на вектор - [math]O(\frac{kn^2}{p})[/math]
  [math]k^2[/math] операций расчета нормы, умножения вектора на число, сложения векторов - [math]O(\frac{k^2n}{p})[/math]
  QR-факторизация матрицы - [math]O(k^3)[/math]

На решении задачи минимизации:

  Обратный ход метода Гаусса - [math]O(k^2)[/math]
  Умножение матрицы на вектор - [math]O(\frac{kn}{p})[/math]
  Сложение векторов - [math]O(\frac{n}{p})[/math]

На оставшиеся операции одной итерации метода Ньютона:

  Сложение векторов - [math]O(\frac{n}{p})[/math]
  Вычисление вектор-функции F в точке - [math]O(\frac{n}{p})[/math]
  Вычисление матрицы якобиана J в точке - [math]O(\frac{n^2}{p})[/math]
  Расчет нормы - [math]O(\frac{n}{p})[/math]

В результате, при небольших значениях параметра [math]k[/math] получим суммарную оценку [math]O(\frac{kn^2}{p})[/math]. Выполнение данного ограничения можно достичь следующими способами:

1) Увеличение скорости сходимости решения СЛАУ при использовании предобуславливания

2) Используя модификацию данного алгоритма - GMRES(m)

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

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

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

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

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

2.4 Масштабируемость алгоритма и его реализации

В качестве модельной задачи рассмотрим один из примеров , поставляемых вместе с модулем SNES пакета PETSc.

Тестирование алгоритма проводилось на суперкомпьютере "Ломоносов"Суперкомпьютерного комплекса Московского университета.

Версия MPI - impi/4.1.0

Реализация BLAS - mkl/11.2.0

Запуски проводились на regulal4

Модель используемого CPU - Intel Xeon X5570 2.93GHz

Версия PETSc - 3.7.4

Флаги компиляции: -fPIC -Wall -Wwrite-strings -Wno-strict-aliasing -Wno-unknown-pragmas -fvisibility=hidden -g3

Кроме того, использовались следующие параметры:

1) Pестарт алгоритма GMRES через 300 итераций

2) Максимальное число итераций для нахождения решения СЛАУ - 1500

3) Максимальное число итераций метода Ньютона на данном этапе решения - 20

Набор начальных параметров:

• Число процессоров [1 2 4 8 16 32 48 64 80 96 112 128];
• Порядок матрицы [82 122 162 202 242].

Рисунок 1. Изменение производительности в зависимости от числа процессоров и размера матрицы.

Рисунок 2. Изменение ускорения в зависимости от числа процессоров и размера матрицы.

Рисунок 3. Изменение эффективности в зависимости от числа процессоров и размера матрицы

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

При фиксированном числе процессоров наблюдается рост ускорения при увеличении вычислительной сложности задачи

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

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

3 Литература

<references \>