Метод Ньютона для систем нелинейных уравнений

...пока в процессе работы (Игорь Коньшин)... (используются описания студентов и 4 их реализации)

??? существует также чья-то страница (октябрь 2017) без реализаций ::: https://algowiki-project.org/ru/Метод_Ньютона

Основные авторы статьи: К.О.Шохин и АртёмА.А.Лебедев (разделы 1, 2.4.1, 2.7)
Авторы отдельных разделов: ЧернышевА.Чернышев и ZakharovН.Захаров (раздел 2.4.2),

Авторы: Гирняк О.Р., Васильков Д.А. (разделы 1.7, 2.4.Х)

Метод Ньютона решения систем нелинейных уравнений
Последовательный алгоритм
Последовательная сложность	[math]\displaystyle{ O(L \cdot n^3) }[/math], если для решения СЛАУ использовать метод Гаусса, алгоритм сойдется за [math]\displaystyle{ L }[/math] итераций, вычисление значения каждой функции системы и их производных в точке имеет оценку сложности [math]\displaystyle{ O(n) }[/math]
Объём входных данных	[math]\displaystyle{ n }[/math] функций от [math]\displaystyle{ n }[/math] переменных, [math]\displaystyle{ n^2 }[/math] частных производных этих функций по каждой переменной, [math]\displaystyle{ n }[/math]-мерный вектор(начальное приближение решения), [math]\displaystyle{ \varepsilon }[/math] - требуемая точность, [math]\displaystyle{ {\rm max\_iter} }[/math] - максимальное число итераций
Объём выходных данных	[math]\displaystyle{ n }[/math]-мерный вектор

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

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

Метод Ньютона решения систем нелинейных уравнений является обобщением метода Ньютона решения нелинейных уравнений, который основан на идее линеаризации. Пусть [math]\displaystyle{ F(x) : {\mathbb R}^1 \to {\mathbb R}^1 }[/math] - дифференцируемая функция и необходимо решить уравнение [math]\displaystyle{ F(x) = 0 }[/math].

Взяв некоторое [math]\displaystyle{ x_0 }[/math] в качестве начального приближения решения, мы можем построить линейную аппроксимацию

[math]\displaystyle{ F(x) }[/math] в окрестности [math]\displaystyle{ x_0 : F(x_0+h) \approx F(x_0)+F^\prime(x_0)h }[/math] и решить получающееся линейное уравнение [math]\displaystyle{ F(x_0 )+F^\prime (x_0 )h =0 }[/math].

Таким образом получаем итеративный метод : [math]\displaystyle{ x_{k+1} = x_k - {F^\prime(x_k)}^{-1}F(x_k), \ k = 0,1,\ldots }[/math] .

Данный метод был предложен Ньютоном в 1669 году. Более точно, Ньютон оперировал только с полиномами; в выражении для [math]\displaystyle{ F(x+h) }[/math] он отбрасывал члены более высокого порядка по h , чем линейные. Ученик Ньютона Рафсон в 1690 г. предложил общую форму метода (т. е. не предполагалось что [math]\displaystyle{ F(x) }[/math] обязательно полином и использовалось понятие производной), поэтому часто говорят о методе Ньютона—Рафсона. Дальнейшее развитие исследований связано с именами таких известных математиков, как Фурье, Коши и другие. Например, Фурье доказал в 1818 г., что метод сходится квадратично в окрестности корня, а Коши (1829, 1847) предложил многомерное обобщение метода и использовал метод для доказательства существования решения уравнения.

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

Пусть дана система из [math]\displaystyle{ n }[/math] нелинейных уравнений с [math]\displaystyle{ n }[/math] неизвестными.

[math]\displaystyle{ \left\{\begin{matrix} f_1(x_1, \ldots, x_n) = 0, \\ f_2(x_1, \ldots, x_n) = 0, \\ \vdots \\ f_n(x_1, \ldots, x_n) = 0. \end{matrix}\right. }[/math] , где [math]\displaystyle{ f_i(x_1, \ldots,x_n) : {\mathbb R}^n \to {\mathbb R}, \ i = 1,\ldots,n }[/math] - нелинейные функции, определенные и непрерывно дифференцируемые в некоторой области [math]\displaystyle{ G \subset {\mathbb R}^n }[/math].

Запишем ее в векторном виде:

[math]\displaystyle{ \overline{x} = {(x_1,x_2,\ldots,x_n)}^T, F(x) ={[f_1(x),f_2(x),\ldots,f_n(x)]}^T, F(x)=0 }[/math]

Требуется найти такой вектор [math]\displaystyle{ \overline{x^*} = {(x^*_1,x^*_2,\ldots,x^*_n)}^T }[/math], который, при подстановке в исходную систему, превращает каждое уравнение в верное числовое равенство.

При таком подходе формула для нахождения решения является естественным обобщением формулы одномерного итеративного метода:

[math]\displaystyle{ x^{(k+1)} = x^{(k)} - W^{-1}(x^{(k)})\cdot F(x^{(k)}) , k=0,1,2,\ldots }[/math], где

[math]\displaystyle{ W = \begin{pmatrix} \frac{\partial{f_1(x_1)}}{\partial{x_1}} & \cdots & \frac{\partial{f_1(x_n)}}{\partial{x_n}} \\ \vdots & \ddots & \vdots \\ \frac{\partial{f_n(x_1)}}{\partial{x_1}} & \cdots & \frac{\partial{f_n(x_n)}}{\partial{x_n}} \end{pmatrix} }[/math] – матрица Якоби.

В рассмотренных предположениях относительно функции [math]\displaystyle{ F(\cdot) }[/math] при выборе начального приближения [math]\displaystyle{ x^{(0)} }[/math] из достаточно малой окрестности решения [math]\displaystyle{ \overline{x^*} }[/math] имеет место сходимость последовательности [math]\displaystyle{ \{x^{(k)}\} }[/math]. При дополнительном предположении [math]\displaystyle{ F(\cdot) \in C^2 }[/math] имеет место квадратичная сходимость метода.

В качестве критерия окончания процесса итераций обычно берут условие [math]\displaystyle{ \left \| x^{(k+1)} - x^{(k)} \right \| \lt \varepsilon }[/math], где [math]\displaystyle{ \varepsilon }[/math] - требуемая точность решения.

Основная сложность метода Ньютона заключается в обращении матрицы Якоби. Вводя обозначение [math]\displaystyle{ \Delta x^{(k)} = x^{(k+1)} - x^{(k)} }[/math] получаем СЛАУ для вычисления [math]\displaystyle{ \Delta x^{(k)}: }[/math] [math]\displaystyle{ \frac{\partial{F(x^{(k)})}}{\partial{x}} = -F(x^{(k)}) }[/math]

Тогда [math]\displaystyle{ x^{(k+1)} = x^{(k)}+ \Delta x^{(k)}. }[/math]

Часто метод Ньютона модифицируют следующим образом. По ходу вычислений или заранее выбирают возрастающую последовательность чисел [math]\displaystyle{ n_0=0, n_1,\ldots }[/math]

При [math]\displaystyle{ n_i \le k \lt n_{i+1} }[/math] вычисление [math]\displaystyle{ \Delta x^{(k)} }[/math] осуществляют по следующей формуле:

[math]\displaystyle{ \frac{\partial{F(x^{n_i})}}{\partial{x}} = -F(x^{(k)}). }[/math]

Увеличение числа итераций, сопровождающее такую модификацию, компенсируется «дешевизной» одного шага итерации.

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

Основная вычислительная нагрузка алгоритма заключается в решении СЛАУ:

[math]\displaystyle{ \frac{\partial F(x^{(k)})}{\partial x}\Delta x^{(k)} = -F(x^{(k)}) \qquad (1) }[/math]

Для нахождения [math]\displaystyle{ \Delta x^{(k)} }[/math], по которому вычисляется значение вектора [math]\displaystyle{ \overline{x} }[/math] на очередной итерации: [math]\displaystyle{ x^{(k+1)} = x^{(k)} + \Delta x^{(k)}. }[/math]

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

Как уже было сказано, основную часть каждой итерации метода составляет нахождение матрицы Якоби и решение СЛАУ [math]\displaystyle{ (1) }[/math] для нахождения [math]\displaystyle{ \Delta x^{(k)}. }[/math]

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

Формулы метода описаны выше. На рис 1. представлена блок-схема алгоритма, в которой:

1. n - число уравнений в СЛАУ
2. Max - предельное число итераций
3. [math]\displaystyle{ \varepsilon }[/math] - точность вычислений
4. [math]\displaystyle{ x^{(0)} }[/math] - начальное приближение

Рис.1 Блок-схема последовательного алгоритма

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

Пусть вычисление значения первой производной функции [math]\displaystyle{ f_i(x_1,x_2,\ldots,x_n) }[/math] в точке [math]\displaystyle{ x^{(k)} }[/math] имеет оценку сложности [math]\displaystyle{ D_i }[/math], решение СЛАУ [math]\displaystyle{ (1) }[/math] имеет оценку сложности [math]\displaystyle{ S }[/math](в эту оценку также включаются операции по нахождению элементов матрицы Якоби), а оценка сложности сложения векторов – [math]\displaystyle{ V }[/math]. Тогда оценка сложности одной итерации представляется в виде:

[math]\displaystyle{ \sum^{n}_{i=1} {D_i} + S + V }[/math]

Количество итераций определяется скоростью сходимости алгоритма. Скорость сходимости метода Ньютона в общем случае является квадратичной. Однако успех или неуспех применения алгоритма во многом определяется выбором начального приближения решения. Если начальное приближение выбрано достаточно хорошо и матрица системы линейных уравнений на каждой итерации хорошо обусловлена и имеет обратную матрицу, то метод Ньютона сходится к единственному в данной окрестности решению. На практике критерием работоспособности метода является число итераций: если оно оказывается большим (для большинства задач >100), то начальное приближение выбрано плохо.

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

• Вычисление значения каждой функции в заданной точке требует [math]\displaystyle{ O(n) }[/math] операций.
• Вычисление значения производной каждой функции в заданной точке требует [math]\displaystyle{ O(n) }[/math] операций.

В таких предположениях указанное выражение для вычисления сложности примет следующий вид:

[math]\displaystyle{ O(n^2) + O(n^3) + O(n) = O(n^3) }[/math]

Предположим также, что выполнены все условия для квадратичной сходимости метода, и для решения системы потребуется небольшое(<100) количество итераций.

Пусть алгоритм сойдется за [math]\displaystyle{ L }[/math] итераций, тогда итоговая сложность будет равна [math]\displaystyle{ O(L \cdot n^3) }[/math].

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

Рис.2 Информационный граф алгоритма

На рис. 2 изображена структура информационного графа алгоритма. Блок IN представляет собой выбор начального приближения, необходимой точности. Каждый блок IT соответствует одной итерации алгоритма, блок-схема итерации алгоритма представлена на рис. 3. Блок OUT соответствует успешному окончанию алгоритма.

Как видно, алгоритм является строго последовательным, каждая итерация зависит от результата предыдущей.

Рис.3 Блок-схема итерации алгоритма

Информационный граф метода Ньютона может быть также представлен в виде следующего рисунка:

Рис.4 Информационный граф алгоритма.

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

Хоть сам по себе алгоритм и является строго последовательным, в нем все же присутствует возможность ускорить выполнение за счет распараллеливания. Основной ресурс параллелизма заключен в решении СЛАУ [math]\displaystyle{ (1) }[/math]. Применять можно любой из известных алгоритмов, ориентируясь на известные особенности структуры матрицы Якоби исходной системы.

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

Пусть задача решается на [math]\displaystyle{ p }[/math] процессорах, тогда сложность решения СЛАУ имеет оценку [math]\displaystyle{ O(\frac{n^3}{p}) }[/math].

Теперь каждому процессору нужно вычислить не [math]\displaystyle{ n }[/math], а [math]\displaystyle{ \frac{n}{p} }[/math] функций, поэтому для вычисления элементов матрицы Якоби и правой части СЛАУ потребуется [math]\displaystyle{ O(\frac{n^2}{p}) }[/math] операций.

Таким образом, получаем следующую оценку параллельной сложности одной итерации алгоритма: [math]\displaystyle{ O(\frac{n^2}{p}) + O(\frac{n^3}{p}) + O(n) = O(\frac{n^3}{p}) }[/math].

Пусть алгоритм сойдется за [math]\displaystyle{ L }[/math] итераций, тогда итоговая сложность будет равна [math]\displaystyle{ O(L \cdot \frac{n^3}{p}) }[/math].

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

При определении объема входных и выходных данных алгоритма будем придерживаться следующих предположений:

• Рассматривается 64-битная система, таким образом каждый адрес занимает 64 бита.
• Все числа, с которыми оперирует алгоритм также занимают 64 бита.
• Каждая функция будет представлена своим адресом, таким образом вклад функций в объем входных данных будет равен [math]\displaystyle{ n }[/math] чисел.
• Производная каждой функции также представлена своим адресом, и вклад производных в объем входных данных будет равен [math]\displaystyle{ n^2 }[/math] чисел.

Входными данными алгоритма являются:

1. Функции, образующие систему уравнений.
2. Производные функций, образующих систему уравнений, если их можно задать аналитически.
3. Точность, с которой требуется найти решение.
4. Максимальное число итераций, которое может быть проделано.
5. Начальное приближение [math]\displaystyle{ x^{(0)} }[/math].

Объем входных данных алгоритма равен [math]\displaystyle{ n ^2 + 2\cdot n + 2 }[/math] в случае, когда на вход алгоритму передаются частные производные ([math]\displaystyle{ n^2 }[/math] адресов частных производных [math]\displaystyle{ n }[/math] функций по [math]\displaystyle{ n }[/math] переменным, [math]\displaystyle{ n }[/math] адресов функций, образующих систему уравнений, [math]\displaystyle{ n }[/math]-мерный вектор начального приближения, требуемая точность и максимальное количество итераций), когда частные производные приближенно вычисляются непосредственно в программе, объем входных данных равен [math]\displaystyle{ 2\cdot n + 2. }[/math]

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

Объем выходных данных: [math]\displaystyle{ n. }[/math]

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

Соотношение последовательной и параллельной сложностей зависит от того, какие алгоритмы применяются для решения СЛАУ [math]\displaystyle{ (1) }[/math]. Результат работы алгоритма сильно зависит от выбора начального приближения решения. При удачном выборе алгоритм достаточно быстро сходится к искомому решению. Глобальная сходимость для многих задач отсутствует.

Существует теорема о достаточном условии сходимости метода Ньютона:

Пусть функция [math]\displaystyle{ F(x) }[/math] непрерывно дифференцируема в открытом выпуклом множестве [math]\displaystyle{ G \subset \R^n }[/math]. Предположим, что существуют [math]\displaystyle{ \overline{x^*} \in \R^n }[/math] и [math]\displaystyle{ r,\beta \gt 0 }[/math], такие что [math]\displaystyle{ N(\overline{x^*},r) \subset G , F(\overline{x^*}) = 0 }[/math] , и существует [math]\displaystyle{ W^{-1}(\overline{x^*}) }[/math], причем [math]\displaystyle{ \left \| W^{-1}(\overline{x^*})\right \| \le \beta }[/math] и [math]\displaystyle{ W(x) \in Lip_{\gamma} (N(\overline{x^*},r)) }[/math]. Тогда существует [math]\displaystyle{ \varepsilon \gt 0 }[/math] такое, что для всех [math]\displaystyle{ x^{(0)} \in N(\overline{x^*}, \varepsilon) }[/math] последовательность [math]\displaystyle{ x^{(1)},x^{(2)},\ldots }[/math] порождаемая итерационным процессом сходится к [math]\displaystyle{ \overline{x^*} }[/math] и удовлетворяет неравенству [math]\displaystyle{ \left \|x^{(k+1)} - \overline{x^*}\right \| \le \beta\cdot\gamma\cdot{\left \| x^{(k)}- \overline{x^*}\right \|}^2 }[/math].

Использовались следующие обозначения:

• [math]\displaystyle{ N(x,r) }[/math] - открытая окрестность радиуса [math]\displaystyle{ r }[/math] с центром в точке [math]\displaystyle{ x }[/math];

• [math]\displaystyle{ W(x) \in Lip_{\gamma}(N(x,r)) }[/math] означает, что [math]\displaystyle{ W(x) }[/math] непрерывна по Липшицу с константой [math]\displaystyle{ \gamma }[/math] в области [math]\displaystyle{ N(x,r) }[/math].

Оставаясь в предположениях, описанных при вычислении объема входных данных, имеем объем входных данных равный [math]\displaystyle{ n^2 + 2\cdot n + 2 = O(n^2) }[/math] чисел.

Оценку числа операций возьмем для реализации метода Ньютона с помощью метода Гаусса: [math]\displaystyle{ O(n^3) }[/math]. Тогда вычислительная мощность есть [math]\displaystyle{ O(n) }[/math].

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

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

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

В настоящем разделе проведено исследование масштабируемости различных параллельных реализаций Метод Ньютона для систем нелинейных уравнений согласно методике AlgoWiki. В основном исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета.

Реализация 1

В исследуемой реализации для решения СЛАУ использовался метод Гаусса. Исследование проводилось для систем размером от 1024 до 4096 уравнений, с шагом 512. Исследуемая система имеет следующий вид:

[math]\displaystyle{ \left\{\begin{matrix} cos(x_1) - 1 = 0, \\ cos(x_2) - 1 = 0, \\ \vdots \\ cos(x_n) - 1 = 0. \end{matrix}\right. \qquad (2) }[/math]

Данная система выбрана так как для нее известно решение, а также начальное приближение [math]\displaystyle{ x_0 = (0.87, 0.87, ..., 0.87) }[/math], при котором метод Ньютона успешно сходится.

Количество ядер процессоров, на которых производились вычисления, изменялось от 16 до 128 по степеням 2. Эксперименты проводились на суперкомпьютере Ломоносов в разделе test. Исследование проводилось на узлах, обладающих следующими характеристиками:

• Количество ядер: 8 ядер архитектуры x86
• Количество памяти: 12Гб
• Количество GPU: 0

На рис. 5 показаны результаты измерений времени решения системы [math]\displaystyle{ (2) }[/math] в зависимости от размера системы и количества процессоров. Можно утверждать, что увеличение числа процессоров дает выигрыш во времени, однако, из-за особенностей реализации параллельного решения СЛАУ, при большом количестве узлов и большом размере системы будет осуществляться пересылка между узлами больших объемов данных и возникающие при этом задержки нивелируют выигрыш от увеличения числа процессоров. Поэтому масштабируемость исследуемой реализации весьма ограничена.

Рис.5 Время решения системы уравнений в зависимости от числа процессоров и размера задачи

Использовался компилятор g++, входящий в состав gcc version 4.4.7 20120313 (Red Hat 4.4.7-4). При компиляции указывался оптимизационный ключ -O2, использовалась библиотека Intel MPI 5.0.1. Ссылка на программную реализацию метода Ньютона для решения систем нелинейных уравнений: https://github.com/ArtyLebedev/newton

Реализация 2

В качестве модельной задачи рассматривался один из примеров, поставляемых вместе с модулем 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].

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

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

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

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

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

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

• Последовательные реализации
  • ALIAS C++. Язык реализации - C++. Распространяется бесплатно, исходные коды и примеры использования можно скачать на сайте.
  • Numerical Recipes. Язык реализации - C++. Исходный код можно найти в секции 9.7 книги Numerical recipes(third edition). Бесплатно доступны предыдущие издания для языков C, Fortran77, Fortran90.
  • Sundials. Язык реализации - C, также есть интерфейс для использования в Fortran-программах. Распространяется по лицензии BSD.
  • Numerical Mathematics- NewtonLib. Язык реализации - C, Fortran. Исходные коды доступны в качестве приложения к книге Peter Deuflhards "Newton Methods for Nonlinear Problems -- Affine Invariance and Adaptive Algorithms".
  • PETSc. Язык реализации - C, есть интерфейс для Java и Python. Распространяется по лицензии BSD.

• Параллельные реализации
  • Sundials. Язык реализации - C, также есть интерфейс для использования в Fortran-программах. Распространяется по лицензии BSD.
  • PETSc. Язык реализации - C, есть интерфейс для Java и Python. Распространяется по лицензии BSD.
  • Построение параллельной реализации метода Ньютона для решения задач оптимизации, описываются в работе V. A. Garanzha, A. I. Golikov, Yu. G. Evtushenko, and M. Kh. Nguen "Parallel Implementation of Newton’s Method for Solving Large-Scale Linear Programs" .
  • Построение параллельной реализации метода Ньютона, а также результаты некоторых экспериментов, описываются в работе Renato N. Elias Alvaro L. G. A. Coutinho Marcos A. D. Martins Rubens M. Sydenstricker "PARALLEL INEXACT NEWTON-TYPE METHODS FOR THE SUPG/PSPG SOLUTION OF STEADY INCOMPRESSIBLE 3D NAVIER-STOKES EQUATIONS IN PC CLUSTERS" .
  • Описание использования и результатов экспериментов с параллельной реализацией метода Ньютона в задаче фильтрации вязкой сжимаемой многофазной многокомпонентной смеси в пористой среде описано в работе К.Ю.Богачева "Эффективное решение задачи фильтрации вязкой сжимаемой многофазной многокомпонентной смеси на параллельных ЭВМ".

Литература

<references \>