Уровень алгоритма

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

Материал из Алговики
Перейти к навигации Перейти к поиску
м (Add \Delta x^{(k)} in two equations.)
 
(не показано 30 промежуточных версий 2 участников)
Строка 1: Строка 1:
В процессе работы (Игорь Коньшин)... (описания студентов и 4 реализации)
+
<!--------------------------------------------------------------------
 
+
??? существует также чья-то страница (октябрь 2017) без реализаций :::
??? существует чья-то среднего качества страница без реализаций :::
 
 
https://algowiki-project.org/ru/Метод_Ньютона
 
https://algowiki-project.org/ru/Метод_Ньютона
 +
--------------------------------------------------------------------->
  
 
+
Основные авторы статьи:  
 
+
[https://algowiki-project.org/ru/Участник:SKirill<b>К.О.Шохин</b>] и
Основные авторы описания: [[Участник:SKirill|Шохин К.О.]](1.1,1.2,1.4,1.6,1.8,1.10,2.7), [[Участник:Лебедев Артём|Лебедев А.А.]](1.1,1.3,1.5,1.7,1.9,2.7)
+
[https://algowiki-project.org/ru/Участник:Лебедев_Артём<b>А.А.Лебедев</b>],<br>
 +
[https://algowiki-project.org/ru/Участник:Oleggium<b>О.Р.Гирняк</b>] и
 +
[https://algowiki-project.org/ru/Участник:Dimx19<b>Д.А.Васильков</b>] (раздел 1.7).
  
 
{{algorithm
 
{{algorithm
 
| name              = Метод Ньютона решения систем нелинейных уравнений
 
| name              = Метод Ньютона решения систем нелинейных уравнений
 
| serial_complexity = <math>O(L \cdot n^3)</math>, если для решения СЛАУ использовать метод Гаусса, алгоритм сойдется за <math>L</math> итераций, вычисление значения каждой функции системы и их производных в точке имеет оценку сложности <math>O(n)</math>
 
| serial_complexity = <math>O(L \cdot n^3)</math>, если для решения СЛАУ использовать метод Гаусса, алгоритм сойдется за <math>L</math> итераций, вычисление значения каждой функции системы и их производных в точке имеет оценку сложности <math>O(n)</math>
| input_data        = <math>n</math> функций от <math>n</math> переменных, <math>n^2</math> частных производных этих функций по каждой переменной, <math>n</math>-мерный вектор(начальное приближение решения), <math>\varepsilon</math> - требуемая точность, <math>max\_iter</math> - максимальное число итераций
+
| input_data        = <math>n</math> функций от <math>n</math> переменных, <math>n^2</math> частных производных этих функций по каждой переменной, <math>n</math>-мерный вектор(начальное приближение решения), <math>\varepsilon</math> - требуемая точность, <math>{\rm max\_iter}</math> - максимальное число итераций
 
| output_data      = <math>n</math>-мерный вектор
 
| output_data      = <math>n</math>-мерный вектор
 
}}
 
}}
Строка 19: Строка 21:
 
=== Общее описание алгоритма ===
 
=== Общее описание алгоритма ===
  
Метод Ньютона<ref name="Tyrtysh">Тыртышников Е. Е. "Методы численного анализа" — М., Академия, 2007. - 320 c.</ref><ref name="Bahvalov">Бахвалов Н. С., Жидков Н. П., Кобельков. Г. М. "Численные Методы" — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.</ref> решения систем нелинейных уравнений является обобщением метода Ньютона решения нелинейных уравнений, который основан на идее линеаризации. Пусть <math> F(x) : \R^1 \to \R^1</math> - дифференцируемая функция и необходимо решить уравнение  
+
Метод Ньютона<ref name="Tyrtysh">Тыртышников Е. Е. "Методы численного анализа" — М., Академия, 2007. - 320 c.</ref><ref name="Bahvalov">Бахвалов Н. С., Жидков Н. П., Кобельков. Г. М. "Численные Методы" — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.</ref> решения систем нелинейных уравнений является обобщением метода Ньютона решения нелинейных уравнений, который основан на идее линеаризации. Пусть <math> F(x) : {\mathbb R}^1 \to {\mathbb R}^1</math> - дифференцируемая функция и необходимо решить уравнение  
 
<math>F(x) = 0</math>.
 
<math>F(x) = 0</math>.
  
 
Взяв некоторое <math>x_0</math> в качестве начального приближения решения, мы можем построить линейную аппроксимацию  
 
Взяв некоторое <math>x_0</math> в качестве начального приближения решения, мы можем построить линейную аппроксимацию  
  
<math>F(x)</math> в окрестности <math>x_0 : F(x_0+h) \approx F(x_0)+F^'(x_0)h</math> и решить получающееся линейное уравнение <math>F(x_0 )+F^' (x_0 )h =0</math>.
+
<math>F(x)</math> в окрестности <math>x_0 : F(x_0+h) \approx F(x_0)+F^\prime(x_0)h</math> и решить получающееся линейное уравнение <math>F(x_0 )+F^\prime (x_0 )h =0</math>.
  
Таким образом получаем итеративный метод : <math> x_{k+1} = x_k - {F^'(x_k)}^{-1}F(x_k) , k = 0,1,\ldots </math> .
+
Таким образом получаем итеративный метод : <math> x_{k+1} = x_k - {F^\prime(x_k)}^{-1}F(x_k), \ k = 0,1,\ldots </math> .
  
 
Данный метод был предложен Ньютоном в 1669 году<ref name = Polyak>Поляк Б. Т. "Метод Ньютона и его роль в оптимизации и вычислительной математике" - Труды ИСА РАН 2006. Т. 28</ref>. Более точно, Ньютон оперировал только с полиномами; в выражении для <math>F(x+h)</math>  он отбрасывал члены более высокого порядка по h , чем линейные.  Ученик Ньютона Рафсон в 1690 г. предложил общую форму метода (т. е. не предполагалось что <math>F(x)</math> обязательно полином и использовалось понятие производной),  поэтому часто говорят о методе Ньютона—Рафсона.
 
Данный метод был предложен Ньютоном в 1669 году<ref name = Polyak>Поляк Б. Т. "Метод Ньютона и его роль в оптимизации и вычислительной математике" - Труды ИСА РАН 2006. Т. 28</ref>. Более точно, Ньютон оперировал только с полиномами; в выражении для <math>F(x+h)</math>  он отбрасывал члены более высокого порядка по h , чем линейные.  Ученик Ньютона Рафсон в 1690 г. предложил общую форму метода (т. е. не предполагалось что <math>F(x)</math> обязательно полином и использовалось понятие производной),  поэтому часто говорят о методе Ньютона—Рафсона.
Строка 43: Строка 45:
 
\end{matrix}\right.
 
\end{matrix}\right.
 
</math>
 
</math>
, где <math>f_i(x_1, \ldots,x_n) : \R^n \to \R, i = 1, \ldots ,n</math> - нелинейные функции, определенные и непрерывно дифференцируемые в некоторой области <math>G \subset \R^n</math>.
+
, где <math>f_i(x_1, \ldots,x_n) : {\mathbb R}^n \to {\mathbb R}, \ i = 1,\ldots,n</math> - нелинейные функции, определенные и непрерывно дифференцируемые в некоторой области <math>G \subset {\mathbb R}^n</math>.
  
 
Запишем ее в векторном виде:  
 
Запишем ее в векторном виде:  
Строка 66: Строка 68:
  
 
Основная сложность метода Ньютона заключается в обращении матрицы Якоби. Вводя обозначение <math>\Delta x^{(k)} = x^{(k+1)} - x^{(k)}</math> получаем СЛАУ для вычисления <math>\Delta x^{(k)}:</math>
 
Основная сложность метода Ньютона заключается в обращении матрицы Якоби. Вводя обозначение <math>\Delta x^{(k)} = x^{(k+1)} - x^{(k)}</math> получаем СЛАУ для вычисления <math>\Delta x^{(k)}:</math>
<math>\frac{\partial{F(x^{(k)})}}{\partial{x}} = -F(x^{(k)})</math>
+
<math>\frac{\partial{F(x^{(k)})}}{\partial{x}} \Delta x^{(k)} = -F(x^{(k)})</math>
  
 
Тогда <math>x^{(k+1)} = x^{(k)}+ \Delta x^{(k)}.</math>
 
Тогда <math>x^{(k+1)} = x^{(k)}+ \Delta x^{(k)}.</math>
Строка 74: Строка 76:
 
При <math>n_i \le k < n_{i+1}</math> вычисление <math>\Delta x^{(k)}</math> осуществляют по следующей формуле:
 
При <math>n_i \le k < n_{i+1}</math> вычисление <math>\Delta x^{(k)}</math> осуществляют по следующей формуле:
  
<math>\frac{\partial{F(x^{n_i})}}{\partial{x}} = -F(x^{(k)}).</math>
+
<math>\frac{\partial{F(x^{n_i})}}{\partial{x}} \Delta x^{(k)} = -F(x^{(k)}).</math>
  
 
Увеличение числа итераций, сопровождающее такую модификацию, компенсируется «дешевизной» одного шага итерации.
 
Увеличение числа итераций, сопровождающее такую модификацию, компенсируется «дешевизной» одного шага итерации.
Строка 82: Строка 84:
 
Основная вычислительная нагрузка алгоритма заключается в решении СЛАУ:
 
Основная вычислительная нагрузка алгоритма заключается в решении СЛАУ:
  
<math>\frac{\partial F(x^{(k)})}{\partial x}\Delta x^{(k)} = -F(x^{(k)}) \qquad(*)</math>
+
<math>\frac{\partial F(x^{(k)})}{\partial x}\Delta x^{(k)} = -F(x^{(k)}) \qquad (1)</math>
  
 
Для нахождения <math>\Delta x^{(k)}</math>, по которому вычисляется значение вектора <math>\overline{x}</math> на очередной итерации: <math>x^{(k+1)} = x^{(k)} + \Delta x^{(k)}.</math>
 
Для нахождения <math>\Delta x^{(k)}</math>, по которому вычисляется значение вектора <math>\overline{x}</math> на очередной итерации: <math>x^{(k+1)} = x^{(k)} + \Delta x^{(k)}.</math>
Строка 88: Строка 90:
 
=== Макроструктура алгоритма ===
 
=== Макроструктура алгоритма ===
  
Как уже было сказано, основную часть каждой итерации метода составляет нахождение матрицы Якоби и решение СЛАУ <math>(*)</math> для нахождения <math>\Delta x^{(k)}.</math>
+
Как уже было сказано, основную часть каждой итерации метода составляет нахождение матрицы Якоби и решение СЛАУ <math>(1)</math> для нахождения <math>\Delta x^{(k)}.</math>
  
 
=== Схема реализации последовательного алгоритма ===
 
=== Схема реализации последовательного алгоритма ===
Строка 94: Строка 96:
 
Формулы метода описаны выше. На рис 1. представлена блок-схема алгоритма, в которой:
 
Формулы метода описаны выше. На рис 1. представлена блок-схема алгоритма, в которой:
  
# n - число уравнений в СЛАУ
+
# <math>n</math> - число уравнений в СЛАУ
 
# Max - предельное число итераций
 
# Max - предельное число итераций
 
# <math>\varepsilon</math> - точность вычислений
 
# <math>\varepsilon</math> - точность вычислений
Строка 103: Строка 105:
 
=== Последовательная сложность алгоритма ===
 
=== Последовательная сложность алгоритма ===
  
Пусть вычисление значения первой производной функции <math>f_i(x_1,x_2,\ldots,x_n)</math> в точке <math>x^{(k)}</math> имеет оценку сложности <math>D_i</math>, решение СЛАУ <math>(*)</math> имеет оценку сложности <math>S</math>(в эту оценку также включаются операции по нахождению элементов матрицы Якоби), а оценка сложности сложения векторов – <math>V</math>. Тогда оценка сложности одной итерации представляется в виде:
+
Пусть вычисление значения первой производной функции <math>f_i(x_1,x_2,\ldots,x_n)</math> в точке <math>x^{(k)}</math> имеет оценку сложности <math>D_i</math>, решение СЛАУ <math>(1)</math> имеет оценку сложности <math>S</math> (в эту оценку также включаются операции по нахождению элементов матрицы Якоби), а оценка сложности сложения векторов – <math>V</math>. Тогда оценка сложности одной итерации представляется в виде:
  
 
<math>\sum^{n}_{i=1} {D_i} + S + V</math>
 
<math>\sum^{n}_{i=1} {D_i} + S + V</math>
Строка 119: Строка 121:
 
<math>O(n^2) + O(n^3) + O(n) = O(n^3)</math>
 
<math>O(n^2) + O(n^3) + O(n) = O(n^3)</math>
  
Предположим также, что выполнены все условия для квадратичной сходимости метода, и для решения системы потребуется небольшое(<100) количество итераций.
+
Предположим также, что выполнены все условия для квадратичной сходимости метода, и для решения системы потребуется небольшое (<100) количество итераций.
  
 
Пусть алгоритм сойдется за <math>L</math> итераций, тогда итоговая сложность будет равна  <math>O(L \cdot n^3)</math>.
 
Пусть алгоритм сойдется за <math>L</math> итераций, тогда итоговая сложность будет равна  <math>O(L \cdot n^3)</math>.
Строка 132: Строка 134:
  
 
[[file:NewtoneIteration.png|thumb|center|700px|Рис.3 Блок-схема итерации алгоритма]]
 
[[file:NewtoneIteration.png|thumb|center|700px|Рис.3 Блок-схема итерации алгоритма]]
 +
 +
Информационный граф метода Ньютона может быть также представлен в виде следующего рисунка:
 +
 +
[[Файл:Graph_newton004.png|1000px|Рис.4 Информационный граф алгоритма.|мини|центр]]
  
 
=== Ресурс параллелизма алгоритма ===
 
=== Ресурс параллелизма алгоритма ===
  
Хоть сам по себе алгоритм и является строго последовательным, в нем все же присутствует возможность ускорить выполнение за счет распараллеливания.
+
Хотя сам по себе алгоритм и является строго последовательным, в нем все же присутствует возможность ускорить выполнение за счет распараллеливания отдельных этапов.
Основной ресурс параллелизма заключен в решении СЛАУ <math>(*)</math>. Применять можно любой из известных алгоритмов, ориентируясь на известные особенности структуры матрицы Якоби исходной системы.
+
Основной ресурс параллелизма заключен в решении СЛАУ <math>(1)</math>. Применять можно любой из известных алгоритмов, ориентируясь на известные особенности структуры матрицы Якоби исходной системы.
  
 
В качестве примера рассмотрим метод Ньютона, использующий для решения СЛАУ параллельный вариант метода Гаусса, в котором между процессорами распределяются строки исходной матрицы. Оставаясь в тех же предположениях, что были сделаны при вычислении оценки последовательной сложности, получим оценку параллельной сложности алгоритма.
 
В качестве примера рассмотрим метод Ньютона, использующий для решения СЛАУ параллельный вариант метода Гаусса, в котором между процессорами распределяются строки исходной матрицы. Оставаясь в тех же предположениях, что были сделаны при вычислении оценки последовательной сложности, получим оценку параллельной сложности алгоритма.
Строка 146: Строка 152:
 
Таким образом, получаем следующую оценку параллельной сложности одной итерации алгоритма: <math>O(\frac{n^2}{p}) + O(\frac{n^3}{p}) + O(n) = O(\frac{n^3}{p})</math>.
 
Таким образом, получаем следующую оценку параллельной сложности одной итерации алгоритма: <math>O(\frac{n^2}{p}) + O(\frac{n^3}{p}) + O(n) = O(\frac{n^3}{p})</math>.
  
Пусть алгоритм сойдется за <math>L</math> итераций, тогда итоговая сложность будет равна <math>O(L \cdot \frac{n^3}{p})</math>.
+
Пусть алгоритм сходится за <math>L</math> итераций, тогда итоговая сложность будет равна <math>O(L \cdot \frac{n^3}{p})</math>.
  
 
=== Входные и выходные данные алгоритма ===
 
=== Входные и выходные данные алгоритма ===
Строка 171: Строка 177:
 
=== Свойства алгоритма ===
 
=== Свойства алгоритма ===
  
Соотношение последовательной и параллельной сложностей зависит от того, какие алгоритмы применяются для решения СЛАУ <math>(*)</math>.
+
Соотношение последовательной и параллельной сложностей зависит от того, какие алгоритмы применяются для решения СЛАУ <math>(1)</math>.
 
Результат работы алгоритма сильно зависит от выбора начального приближения решения. При удачном выборе алгоритм достаточно быстро сходится к искомому решению. Глобальная сходимость для многих задач отсутствует.  
 
Результат работы алгоритма сильно зависит от выбора начального приближения решения. При удачном выборе алгоритм достаточно быстро сходится к искомому решению. Глобальная сходимость для многих задач отсутствует.  
  
 
Существует теорема о достаточном условии сходимости метода Ньютона:
 
Существует теорема о достаточном условии сходимости метода Ньютона:
  
Пусть функция <math>F(x)</math> непрерывно дифференцируема в открытом выпуклом множестве <math>G \subset \R^n</math>. Предположим, что существуют <math>\overline{x^*} \in \R^n</math> и <math>r,\beta > 0 </math>, такие что <math>N(\overline{x^*},r) \subset G , F(\overline{x^*}) = 0</math>  , и существует <math>W^{-1}(\overline{x^*})</math>, причем <math>\left \| W^{-1}(\overline{x^*})\right \| \le \beta </math> и <math>W(x) \in Lip_{\gamma} (N(\overline{x^*},r))</math>. Тогда существует <math>\varepsilon > 0</math> такое, что для всех <math>x^{(0)} \in N(\overline{x^*}, \varepsilon)</math> последовательность <math>x^{(1)},x^{(2)},\ldots</math> порождаемая итерационным процессом сходится к <math>\overline{x^*}</math> и удовлетворяет неравенству <math>\left \|x^{(k+1)} - \overline{x^*}\right \| \le \beta\cdot\gamma\cdot{\left \| x^{(k)}- \overline{x^*}\right \|}^2</math>.
+
Пусть функция <math>F(x)</math> непрерывно дифференцируема в открытом выпуклом множестве <math>G \subset {\mathbb R}^n</math>. Предположим, что существуют <math>\overline{x^*} \in {\mathbb R}^n</math> и <math>r,\beta > 0 </math>, такие что <math>N(\overline{x^*},r) \subset G , F(\overline{x^*}) = 0</math>  , и существует <math>W^{-1}(\overline{x^*})</math>, причем <math>\left \| W^{-1}(\overline{x^*})\right \| \le \beta </math> и <math>W(x) \in Lip_{\gamma} (N(\overline{x^*},r))</math>. Тогда существует <math>\varepsilon > 0</math> такое, что для всех <math>x^{(0)} \in N(\overline{x^*}, \varepsilon)</math> последовательность <math>x^{(1)},x^{(2)},\ldots</math> порождаемая итерационным процессом сходится к <math>\overline{x^*}</math> и удовлетворяет неравенству <math>\left \|x^{(k+1)} - \overline{x^*}\right \| \le \beta\cdot\gamma\cdot{\left \| x^{(k)}- \overline{x^*}\right \|}^2</math>.
  
 
Использовались следующие обозначения:
 
Использовались следующие обозначения:
Строка 191: Строка 197:
  
 
=== Особенности реализации последовательного алгоритма ===
 
=== Особенности реализации последовательного алгоритма ===
 
=== Локальность данных и вычислений ===
 
  
 
=== Возможные способы и особенности параллельной реализации алгоритма ===
 
=== Возможные способы и особенности параллельной реализации алгоритма ===
  
=== Масштабируемость алгоритма и его реализации ===
+
=== Результаты прогонов ===
 
 
==== Масштабируемость алгоритма ====
 
 
 
Так как алгоритм по сути своей является последовательным, то масштабируемость конкретной реализации алгоритма определяется масштабируемостью реализации алгоритма решения СЛАУ, используемого для нахождения изменения текущего решения на очередной итерации.
 
 
 
==== Масштабируемость реализации алгоритма ====
 
 
 
В исследуемой реализации для решения СЛАУ использовался метод Гаусса. Исследование проводилось для систем размером от 1024 до 4096 уравнений, с шагом 512. Исследуемая система имеет следующий вид:
 
 
 
<math>
 
\left\{\begin{matrix}
 
cos(x_1) - 1 = 0,
 
\\ cos(x_2) - 1 = 0,
 
\\ \vdots
 
\\ cos(x_n) - 1 = 0.
 
\end{matrix}\right. \qquad(**)
 
</math>
 
 
 
Данная система выбрана так как для нее известно решение, а также известно начальное приближение <math> x_0 = (0.87, 0.87, ..., 0.87) </math>,  при котором метод Ньютона успешно сходится.
 
 
 
Количество ядер процессоров, на которых производились вычисления, изменялось от 16 до 128 по степеням 2. Эксперименты проводились на суперкомпьютере Ломоносов в разделе test. Исследование проводилось на узлах, обладающих следующими характеристиками:
 
 
 
*Количество ядер: 8 ядер архитектуры x86
 
*Количество памяти: 12Гб
 
*Количество GPU: 0
 
 
 
На рис. 3 показаны результаты измерений времени решения системы <math>(**)</math> в зависимости от размера системы и количества процессоров. Можно утверждать, что увеличение числа процессоров дает выигрыш во времени, однако, из-за особенностей реализации параллельного решения СЛАУ, при большом количестве узлов и большом размере системы будет осуществляться пересылка между узлами больших объемов данных и возникающие при этом задержки нивелируют выигрыш от увеличения числа процессоров. Поэтому масштабируемость исследуемой реализации весьма ограничена.
 
 
 
[[file:Снимок экрана 2016-11-15 в 22.36.33.png|thumb|center|700px|Рис.3 Время решения системы уравнений в зависимости от числа процессоров и размера задачи]]
 
 
 
Использовался компилятор g++, входящий в состав gcc version 4.4.7 20120313 (Red Hat 4.4.7-4). При компиляции указывался оптимизационный ключ -O2, использовалась библиотека Intel MPI 5.0.1. Ссылка на программную реализацию метода Ньютона для решения систем нелинейных уравнений: https://github.com/ArtyLebedev/newton <ref name="Программная реализация метода Ньютона">https://github.com/ArtyLebedev/newton</ref>
 
 
 
=== Динамические характеристики и эффективность реализации алгоритма ===
 
 
 
 
=== Выводы для классов архитектур ===
 
=== Выводы для классов архитектур ===
  
=== Существующие реализации алгоритма ===
+
== Литература ==
  
*Последовательные реализации
+
<references />
**ALIAS C++. Язык реализации - C++. Распространяется бесплатно, исходные коды и примеры использования можно скачать на сайте<ref name="ALIAS">https://www-sop.inria.fr/coprin/logiciels/ALIAS/ALIAS-C++/ALIAS-C++.html</ref>.
 
**Numerical Recipes. Язык реализации - C++. Исходный код можно найти в секции 9.7 книги Numerical recipes(third edition)<ref name="RECIPES">http://numerical.recipes</ref>. Бесплатно доступны<ref name="RECIPES_old">http://numerical.recipes/oldverswitcher.html</ref> предыдущие издания для языков C, Fortran77, Fortran90.
 
**Sundials. Язык реализации - C, также есть интерфейс для использования в Fortran-программах. Распространяется<ref name="SUNDIALS_seq">http://computation.llnl.gov/projects/sundials/kinsol</ref> по лицензии BSD.
 
**Numerical Mathematics- NewtonLib. Язык реализации - C, Fortran. Исходные коды доступны<ref name="LIB_book">http://elib.zib.de/pub/elib/codelib/NewtonLib/index.html</ref> в качестве приложения к книге<ref name="LIB">https://www.amazon.com/Newton-Methods-Nonlinear-Problems-Computational/dp/364223898X/</ref> Peter Deuflhards "Newton Methods for Nonlinear Problems -- Affine Invariance and Adaptive Algorithms".
 
**PETSc. Язык реализации - C, есть интерфейс для Java и Python. Распространяется<ref name =  PETSC>https://www.mcs.anl.gov/petsc/</ref> по лицензии BSD.
 
  
*Параллельные реализации
+
[[Категория:Статьи в работе]]
**Sundials. Язык реализации - C, также есть интерфейс для использования в Fortran-программах. Распространяется<ref name="SUNDIALS_seq">http://computation.llnl.gov/projects/sundials/kinsol</ref> по лицензии BSD.
 
**PETSc. Язык реализации - C, есть интерфейс для Java и Python. Распространяется<ref name =  PETSC>https://www.mcs.anl.gov/petsc/</ref> по лицензии 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" <ref name="IMPL_OPT">http://www.ccas.ru/personal/evtush/p/CMMP1303.pdf</ref>.
 
**Построение параллельной реализации метода Ньютона, а также результаты некоторых экспериментов, описываются в работе 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" <ref name="SUPG">http://www.nacad.ufrj.br/~rnelias/papers/CIL14-020.pdf</ref>.
 
**Описание использования и результатов экспериментов с параллельной реализацией метода Ньютона в задаче фильтрации вязкой сжимаемой многофазной многокомпонентной смеси в пористой среде описано в работе К.Ю.Богачева "Эффективное решение задачи фильтрации вязкой сжимаемой многофазной многокомпонентной смеси на параллельных ЭВМ"<ref name="BOGACHEV">http://academy.hpc-russia.ru/files/reservoirdynamics_bogachev.pdf</ref>.
 
 
 
== Литература ==
 
  
<references \>
+
[[en:Newton's method for systems of nonlinear equations]]

Текущая версия на 00:17, 4 декабря 2022


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



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


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

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

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

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

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

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

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

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

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

[math] \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]f_i(x_1, \ldots,x_n) : {\mathbb R}^n \to {\mathbb R}, \ i = 1,\ldots,n[/math] - нелинейные функции, определенные и непрерывно дифференцируемые в некоторой области [math]G \subset {\mathbb R}^n[/math].

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

[math]\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]\overline{x^*} = {(x^*_1,x^*_2,\ldots,x^*_n)}^T[/math], который, при подстановке в исходную систему, превращает каждое уравнение в верное числовое равенство.

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

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

[math] 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]F(\cdot)[/math] при выборе начального приближения [math]x^{(0)}[/math] из достаточно малой окрестности решения [math]\overline{x^*}[/math] имеет место сходимость последовательности [math]\{x^{(k)}\}[/math]. При дополнительном предположении [math]F(\cdot) \in C^2[/math] имеет место квадратичная сходимость метода.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[math]O(n^2) + O(n^3) + O(n) = O(n^3)[/math]

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1 Особенности реализации последовательного алгоритма

2.2 Возможные способы и особенности параллельной реализации алгоритма

2.3 Результаты прогонов

2.4 Выводы для классов архитектур

3 Литература

  1. Тыртышников Е. Е. "Методы численного анализа" — М., Академия, 2007. - 320 c.
  2. Бахвалов Н. С., Жидков Н. П., Кобельков. Г. М. "Численные Методы" — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.
  3. Поляк Б. Т. "Метод Ньютона и его роль в оптимизации и вычислительной математике" - Труды ИСА РАН 2006. Т. 28