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

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

Материал из Алговики
Перейти к навигации Перейти к поиску
[непроверенная версия][досмотренная версия]
м (Add \Delta x^{(k)} in two equations.)
 
(не показано 19 промежуточных версий 2 участников)
Строка 1: Строка 1:
...пока в процессе работы (Игорь Коньшин)... (используются описания студентов и 4 их реализации)
+
<!--------------------------------------------------------------------
 
 
 
??? существует также чья-то страница (октябрь 2017) без реализаций :::
 
??? существует также чья-то страница (октябрь 2017) без реализаций :::
 
https://algowiki-project.org/ru/Метод_Ньютона
 
https://algowiki-project.org/ru/Метод_Ньютона
 
+
--------------------------------------------------------------------->
<!---------------------------------------------------------------->
 
  
 
Основные авторы статьи:  
 
Основные авторы статьи:  
 
[https://algowiki-project.org/ru/Участник:SKirill<b>К.О.Шохин</b>] и
 
[https://algowiki-project.org/ru/Участник:SKirill<b>К.О.Шохин</b>] и
[https://algowiki-project.org/ru/Участник:Лебедев Артём<b>А.А.Лебедев</b>] (разделы 1, 2.4.1, 2.7)<br/>
+
[https://algowiki-project.org/ru/Участник:Лебедев_Артём<b>А.А.Лебедев</b>],<br>
Авторы отдельных разделов:
 
[https://algowiki-project.org/ru/Участник:Александр Чернышев<b>А.Чернышев</b>] и
 
[https://algowiki-project.org/ru/Участник:N Zakharov<b>Н.Захаров</b>] (раздел 2.4.2),
 
 
[https://algowiki-project.org/ru/Участник:Oleggium<b>О.Р.Гирняк</b>] и
 
[https://algowiki-project.org/ru/Участник:Oleggium<b>О.Р.Гирняк</b>] и
[https://algowiki-project.org/ru/Участник:Dimx19<b>Д.А.Васильков</b>] (разделы 1.7, 2.4.3),
+
[https://algowiki-project.org/ru/Участник:Dimx19<b>Д.А.Васильков</b>] (раздел 1.7).
[https://algowiki-project.org/ru/Участник:<b></b>] и
 
[https://algowiki-project.org/ru/Участник:<b></b>] (раздел 2.4.4)
 
  
 
{{algorithm
 
{{algorithm
Строка 75: Строка 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>
Строка 83: Строка 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>
  
 
Увеличение числа итераций, сопровождающее такую модификацию, компенсируется «дешевизной» одного шага итерации.
 
Увеличение числа итераций, сопровождающее такую модификацию, компенсируется «дешевизной» одного шага итерации.
Строка 103: Строка 96:
 
Формулы метода описаны выше. На рис 1. представлена блок-схема алгоритма, в которой:
 
Формулы метода описаны выше. На рис 1. представлена блок-схема алгоритма, в которой:
  
# n - число уравнений в СЛАУ
+
# <math>n</math> - число уравнений в СЛАУ
 
# Max - предельное число итераций
 
# Max - предельное число итераций
 
# <math>\varepsilon</math> - точность вычислений
 
# <math>\varepsilon</math> - точность вычислений
Строка 112: Строка 105:
 
=== Последовательная сложность алгоритма ===
 
=== Последовательная сложность алгоритма ===
  
Пусть вычисление значения первой производной функции <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>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>
Строка 128: Строка 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>.
Строка 148: Строка 141:
 
=== Ресурс параллелизма алгоритма ===
 
=== Ресурс параллелизма алгоритма ===
  
Хоть сам по себе алгоритм и является строго последовательным, в нем все же присутствует возможность ускорить выполнение за счет распараллеливания.
+
Хотя сам по себе алгоритм и является строго последовательным, в нем все же присутствует возможность ускорить выполнение за счет распараллеливания отдельных этапов.
 
Основной ресурс параллелизма заключен в решении СЛАУ <math>(1)</math>. Применять можно любой из известных алгоритмов, ориентируясь на известные особенности структуры матрицы Якоби исходной системы.
 
Основной ресурс параллелизма заключен в решении СЛАУ <math>(1)</math>. Применять можно любой из известных алгоритмов, ориентируясь на известные особенности структуры матрицы Якоби исходной системы.
  
Строка 159: Строка 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>.
  
 
=== Входные и выходные данные алгоритма ===
 
=== Входные и выходные данные алгоритма ===
Строка 189: Строка 182:
 
Существует теорема о достаточном условии сходимости метода Ньютона:
 
Существует теорема о достаточном условии сходимости метода Ньютона:
  
Пусть функция <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>.
  
 
Использовались следующие обозначения:
 
Использовались следующие обозначения:
Строка 204: Строка 197:
  
 
=== Особенности реализации последовательного алгоритма ===
 
=== Особенности реализации последовательного алгоритма ===
 
=== Локальность данных и вычислений ===
 
  
 
=== Возможные способы и особенности параллельной реализации алгоритма ===
 
=== Возможные способы и особенности параллельной реализации алгоритма ===
  
===Масштабируемость алгоритма и его реализации===
+
=== Результаты прогонов ===
 
 
Масштабируемость алгоритма и его реализаций определяется главным образом масштабируемостью реализации алгоритма решения СЛАУ, используемого для нахождения изменения текущего решения на очередной итерации.
 
 
 
В настоящем разделе проведено исследование масштабируемости различных параллельных реализаций Метод Ньютона для систем нелинейных уравнений согласно [[Scalability methodology|методике]] AlgoWiki. В основном исследование проводилось на суперкомпьютере "Ломоносов" [http://parallel.ru/cluster Суперкомпьютерного комплекса Московского университета].
 
 
 
==== Реализация 1 ====
 
 
 
В исследуемой реализации для решения СЛАУ использовался метод Гаусса. Исследование проводилось для систем размером от 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 (2)
 
</math>
 
 
 
Данная система выбрана так как для нее известно решение, а также начальное приближение <math> x_0 = (0.87, 0.87, ..., 0.87) </math>,  при котором метод Ньютона успешно сходится.
 
 
 
Количество ядер процессоров, на которых производились вычисления, изменялось от 16 до 128 по степеням 2. Эксперименты проводились на суперкомпьютере Ломоносов в разделе test. Исследование проводилось на узлах, обладающих следующими характеристиками:
 
 
 
*Количество ядер: 8 ядер архитектуры x86
 
*Количество памяти: 12Гб
 
*Количество GPU: 0
 
 
 
На рис. 5 показаны результаты измерений времени решения системы <math>(2)</math> в зависимости от размера системы и количества процессоров. Можно утверждать, что увеличение числа процессоров дает выигрыш во времени, однако, из-за особенностей реализации параллельного решения СЛАУ, при большом количестве узлов и большом размере системы будет осуществляться пересылка между узлами больших объемов данных и возникающие при этом задержки нивелируют выигрыш от увеличения числа процессоров. Поэтому масштабируемость исследуемой реализации весьма ограничена.
 
 
 
[[file:Снимок экрана 2016-11-15 в 22.36.33.png|thumb|center|700px|Рис.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 <ref name="Программная реализация метода Ньютона">https://github.com/ArtyLebedev/newton</ref>
 
 
 
==== Реализация 2 ====
 
 
 
В качестве модельной задачи рассматривался один из примеров<ref>http://www.mcs.anl.gov/petsc/petsc-3.5/src/snes/examples/tutorials/ex30.c.html</ref>, поставляемых вместе с модулем SNES пакета PETSc.
 
 
 
Тестирование алгоритма проводилось на суперкомпьютере "Ломоносов"[http://parallel.ru/cluster Суперкомпьютерного комплекса Московского университета].
 
 
 
* Версия 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
 
 
 
Кроме того, использовались следующие параметры:
 
 
 
# Pестарт алгоритма GMRES через 300 итераций
 
# Максимальное число итераций для нахождения решения СЛАУ - 1500
 
# Максимальное число итераций метода Ньютона на данном этапе решения - 20
 
 
 
Набор начальных параметров:
 
 
 
* Число процессоров [1 2 4 8 16 32 48 64 80 96 112 128];
 
* Порядок матрицы [82 122 162 202 242].
 
 
 
[[Файл:NewtonFlops.jpg|thumb|center|700px|Рис. 6. Изменение производительности в зависимости от числа процессоров и размера матрицы.]]
 
[[Файл:NewtonAcc.jpg|thumb|center|700px|Рис. 7. Изменение ускорения в зависимости от числа процессоров и размера матрицы.]]
 
[[Файл:NewtonEff.jpg|thumb|center|700px|Рис. 8. Изменение эффективности в зависимости от числа процессоров и размера матрицы]]
 
 
 
Как видно из приведенных графиков, эффективность распараллеливания алгоритма довольно быстро убывает при увеличении числа процессоров.
 
 
 
При фиксированном числе процессоров наблюдается рост ускорения при увеличении вычислительной сложности задачи.
 
 
 
==== Реализация 3 ====
 
 
 
Для исследования масштабируемости алгоритм был реализован нами самостоятельно на библиотеке CUDA 8.0 для языка C++(компилятор nvcc) для Windows на домашнем компьютере(OS Windows 10 x64, CPU Intel Core i7-2600 3.4 GHz, 8 GB RAM, GPU NVIDIA GTX 550 Ti). GPU поддерживает спецификацию  CUDA 2.1 и максимальное число потоков в блоке [https://ru.wikipedia.org/wiki/CUDA 1024]. Однако физически у нас есть всего 192 CUDA ядра(см [http://www.geforce.com/hardware/desktop-gpus/geforce-gtx-550ti/specifications спецификации] видеокарты). Все потоки должны быть элементами двумерного квадратного блока. Поэтому исследования масштабируемости были проведены на 4, 16, 36, 64, 100, 144, 196, 256, 324, 400 потоках(на большем количестве потоков проводить исследования бессмысленно ввиду ухудшения производительности). Параметры компиляции:
 
<source>"nvcc.exe" -gencode=arch=compute_20,code=\"sm_20,compute_20\" --use-local-env --cl-version 2015 -ccbin "C:\Program Files (x86)\Microsoft Visual Studio 14.0\VC\bin"  -I"C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v8.0\include" -I"C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v8.0\include" -G --keep-dir Debug -maxrregcount=0 --machine 32 --compile -cudart static  -g  -DWIN32 -D_DEBUG -D_CONSOLE -D_MBCS -Xcompiler "/EHsc /W3 /nologo /Od /FS /Zi /RTC1 /MDd " -o Debug\kernel.cu.obj "kernel.cu" </source>
 
Решение СЛАУ было реализовано в 2 этапа:
 
 
 
1) Нахождение обратной матрицы с помощью метода Гаусса-Жордана. Для этого было написано ядро для CUDA:
 
<source lang="c++">__global__ void gaussjordan(float *A, float *I, int n, int i)</source>
 
 
 
2) Умножение обратной матрицы на вектор правых частей. Для этого было написано ядро для CUDA:
 
<source lang="c++">__global__ void multiply_by_vec(float *mat, float *vec, float *out, const int N)</source>
 
 
 
Необходимо отметить, что использование своих ядер для CUDA непременно ведет к снижению производительности. Для улучшения производительности предпочтительно использовать матричные операции из библиотеки cuBLAS. Однако, при таком подходе невозможно явно задать количество потоков - для конкретной видеокарты сама библиотека выбирает, какое количество потоков ей использовать. В связи с этим и были написаны свои ядра, позволяющие задавать количество потоков.
 
 
 
Считалось, что на вход алгоритму были заданы нелинейные функции и их производные. С данной реализацией можно ознакомиться [https://dl.dropboxusercontent.com/u/43317547/ParallelingCuda.zip здесь].
 
 
 
Говоря о масштабируемости данного алгоритма, важно отметить следующее:
 
 
 
1. Время вычисления набора функций и набора их производных очень сильно зависит от их сложности и не связано с мощностью алгоритма.
 
 
 
2. Скорость отработки алгоритма от начала до конца может не быть связана с его мощностью: заранее неизвестно начальное приближение, а значит и количество итераций, за которое он сойдется.
 
 
 
Вывод: для оценки масштабируемости необходимо сравнивать одну итерацию алгоритма, и притом, без учета времени вычисления значения функций и их производных. В итоге, приходим к выводу, что целесообразно измерить время выполнения одного решения СЛАУ. Также стоит отметить, что в нашей реализации при хорошем выборе начального приближения алгоритм сходился за 10-20 итераций.
 
 
 
Можно рассмотреть зависимость масштабируемости от одного из двух критериев: размерности задачи или же количества потоков. На графике ниже представлена зависимость времени решения СЛАУ от количества потоков и размерности матрицы.
 
 
 
[[Файл:Newton_plot.png|1100px|Рис. 9. Время решения СЛАУ в зависимости от размерности системы и количества потоков|мини|центр]]
 
 
 
Как уже было сказано, из-за наличия только 192 физических процессоров, данный алгоритм показывает лучшее время работы при использовании 144 либо 256 потоков. 144 потока (размер блока 12х12) - наиболее близкое количество потоков к 192 физическим процессорам. Хорошее время работы на 256 потоках объясняется тем, что мы используем стандартный размер блока 16х16, а также многие операции деления или умножения на количество потоков или размер блока в этом случае компилятор может заменить на операции побитового сдвига, которые работают быстрее обычных.
 
 
 
==== Реализация 4 ====
 
 
 
 
 
=== Динамические характеристики и эффективность реализации алгоритма ===
 
 
 
 
=== Выводы для классов архитектур ===
 
=== Выводы для классов архитектур ===
  
=== Существующие реализации алгоритма ===
+
== Литература ==
  
*Последовательные реализации
+
<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