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

Участница:Эльвира Ахиярова/Итерационный метод решения системы линейных алгебраических уравнений GMRES (обобщенный метод минимальных невязок): различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
м (Lineprinter переименовал страницу [[Итерационный метод решения системы линейных алгебраических уравнений GMRES (обобщенный метод минимальных н…)
м
Строка 1: Строка 1:
−
Автор: [[U:Эльвира Ахиярова|Эльвира Ахиярова]]
+
{{Assignment}}
 +
Итерационный метод решения системы линейных алгебраических уравнений GMRES (обобщенный метод минимальных невязок).
 +
 
 +
{{algorithm
 +
| name              = Обобщеннный метод минимальных невязок
 +
| serial_complexity = <math>O(mn+m^2+n^2)</math>
 +
| input_data        = <math>n^2+n+1</math>
 +
| output_data      = <math>n</math>
 +
}}
 +
Основные авторы описания: [[Участник:Эльвира Ахиярова|Э.А.Ахиярова]]
 +
 
 
== Свойства и структура алгоритма ==
 
== Свойства и структура алгоритма ==
 +
 
=== Общее описание алгоритма ===
 
=== Общее описание алгоритма ===
−
Обобщенный метод минимальных невязок (GMRES) является итерационным методом нахождения решения для системы линейных алгебраических уравнений с произвольной невырожденной матрицей. Метод основан на минимизации квадратичного функционала невязки на подпространствах Крылова.
+
 
−
GMRES был предложен Йозефом Саадом и Мартином Шульцем в 1986 году.<ref>Y.Saad, M.H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Scientific and Stat. Comp. 7: 856-869 (1986). </ref>
+
Обобщенный метод минимальных невязок (GMRES) - это итерационный метод для нахождения приближенного решения системы линейных алгебраических уравнений <math>Ax = b </math> с произвольной матрицей <math>A</math>. Алгоритм был преложен в 1986 году  Йозефом Саадом и Мартином Г. Шульцем
−
Наиболее популярные реализации метода основаны на модифицированном алгоритме ортогонализации Грама-Шмидта и на использовании рестартов для управления требуемым объемом памяти.
+
<ref>Y. Saad and M.H. Schultz, "GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems", SIAM J. Sci. Stat. Comput., 7:856-869, 1986. doi:10.1137/0907058.</ref>.
 +
 
 +
В то время как применение классических итерационных алгоритмов для решения СЛАУ ограничивается либо диагональным преобладанием, либо положительной определенностью матрицы, GMRES может быть использован для систем с произвольной невырожденной матрицей <math>A</math>. Данный метод опирается на итерации Арнольди для нахождения минимизирующего невязку вектора из подпространства Крылова <ref>https://en.wikipedia.org/wiki/Krylov_subspace</ref>.
 +
 
 +
Основным недостатком метода GMRES является поитерационное увеличение объема требуемой памяти. В связи с этим, для некоторых вычислительных систем выполенние программы метода может быть невозможным. Для некоторых систем используют модифицированный  обобщенный метод минимальных невязок - GMRES(m), который подразумевает использование рестартов, существенно сокращающих объемы необходимой памяти.
 +
 
 +
Ключевая идея обобщенного метода минимальных невязок основана на решении задачи наименьших квадратов на каждом итерационном шаге. На каждом <math>k</math>-м шаге точное решение <math>x^* = A^{-1}b </math> аппроксимируется вектором <math>x_k \in K_k</math>, таким образом, что:
 +
:<math>
 +
\|r_k{\|}_2 = \|Ax_k-b{\|}_2 \to \min
 +
</math>,
 +
где <math>K_k</math> - подпространство Крылова <math> k </math> - го порядка.
  
 
=== Математическое описание алгоритма ===
 
=== Математическое описание алгоритма ===
−
Основная идея метода GMRES основана на минимизации нормы невязки на каждом итерационном шаге. На каждой <math>k</math>-й итерации  точное решение <math>x^* = A^{-1}b</math> приближается к вектору <math>x_k \in K_k</math>(k-е подпространство Крылова) в том смысле, что
+
'''Исходные данные:''' Матрица <math>A \in R^{n\times n}</math>, вектор правой части <math>b \in R^n</math>.
 +
 
 +
'''Вычисляемые данные:''' Приближенное решение <math>\tilde{x} \in R^n</math> системы <math>Ax = b</math>.
 +
 
 +
На каждом <math>k</math>-м шаге метода решение <math>x_k</math> рассматривается в виде <math>x_0+Vy_k</math>, где столбцы <math>v_k</math> матрицы <math>V</math> находятся с помощью процесса Арнольди для построения ортонормированного базиса <math>\{v_1, v_2, \dots, v_k\}</math> подпространства Крылова <math>K_k(A, v_1)</math>,
 +
вектор <math>y_k</math> представляет собой решение системы <math>H_ky_k = \beta_k e_1</math>, где <math>H_k</math> - верхняя  матрица Хессенберга<ref>https://en.wikipedia.org/wiki/Hessenberg_matrix</ref> размера <math>k\times k</math>, элементы которой формируются в процессе ортогонализации Арнольди, <math>\beta_k = \|r_0{\|}_2, \ r_0 = b - Ax_0 </math>, и, наконец, <math>e_1</math> - канонический вектор <math>{\begin{pmatrix}1&0&\dots 0\end{pmatrix}}^T</math> размерности <math>k</math>
 +
 
 +
Рассмотрим подробнее процесс ортогонализации Арнольди:
 +
 
 +
==== Процесс Арнольди ====
 +
'''Исходные данные:''' начальный ортонормированный вектор <math>v_1</math> и размерность <math>k</math>
 +
 
 +
'''Вычисляемые данные:''' Ортонормированный базис <math>\{v_1, v_2, \dots, v_k\}</math> подпространства <math>K_k(A, v_1)</math>
 +
 
 +
Формулы метода:
 
:<math>
 
:<math>
 
\begin{align}
 
\begin{align}
−
\|r_k{\|}_2 &= \|Ax_k-b{\|}_2 \to \min\limits_{x_k \in K_k}
+
&\omega = Av_j \\
 +
&h_{ij} = (\omega, v_i)  \\
 +
&\omega = \omega - h_{ij}v_i \\
 +
&h_{(i+1)j} = \|\omega\| \\
 +
&v_{j+1} = \frac{\omega}{h_{(i+1)j}} \\
 +
&j  = \overline{1,k}, \ i = \overline{1,j}  
 
\end{align}
 
\end{align}
 
</math>
 
</math>
−
Приближенное решение <math>x_k</math> для системы линейных алгебраических уравнений
+
В результате выполнения процесса Арнольди получаем:
 +
:<math>V_k =\{v_1, v_2, \dots, v_k\}</math> - ортонормированный базис подпространства Крылова
 +
:<math>\overline{V_k} = \{v_1, v_2, \dots, v_k, v_{k+1}\}</math> - матрица, расширенная последним вычисленным вектором <math>v_{k+1}</math>
 +
:<math>H_k</math> - верхняя Хессенбергова матрица, элементы которой равны коэффициентам ортогонализации:
 +
:<math>H_k =
 +
\begin{pmatrix}
 +
h_{11} & h_{12} & h_{13} & h_{14} & \dots & h_{1k} \\
 +
h_{21} & h_{22} & h_{23} & h_{24} & \dots & h_{2k} \\
 +
0 & h_{32} & h_{33} & h_{34} & \dots & h_{3k} \\
 +
\vdots & &\ddots & \ddots & & \vdots \\
 +
0 & \dots & 0 & h_{k-1,k-2}& h_{k-1, k-1} & h_{k-1, k} \\
 +
0 & \dots & 0 & 0 & h_{k, k-1} & h_{k, k} \\
 +
0 & \dots & & & 0& h_{k+1, k}
 +
\end{pmatrix}
 +
</math>
 +
 
 +
==== Минимизация нормы невязки ====
 +
Выполнив несложные преобразования можно получить выражение для невязки на <math>k</math>- й итерации:
 +
:<math>
 +
r_k = \beta e_1 - H_ky_k
 +
</math>
 +
Очевидно, что для задачи
 +
:<math>
 +
\|r_k{\|}_2 \to \min
 +
</math>
 +
минимум достигается при
 +
:<math> H_ky_k = \beta e_1 </math>
 +
Для решения последней системы матрица <math>H_k</math> приводится к верхнему треугольному виду с помощью <math> k </math> вращений Гивенса:
 +
:<math> R_k = \Omega_k\Omega_{k-1}\dots \Omega_1 </math>,
 +
:<math> \Omega_i = \begin{pmatrix} 1 & &&&&&& \\
 +
&\ddots &&&&&&\\
 +
&&1&&&&& \\       
 +
&&&c_i&s_i &&& \\
 +
&&&-s_i& c_i &&&\\
 +
&&&&&1&&& \\
 +
&&&&&&\ddots&\\
 +
&&&&&&&1 \\  \end{pmatrix}, \quad 
 +
s_i = \frac{h_{i+1,i}}{\sqrt{h_{ii}^2+h_{i+1,i}^2}}, \quad c_i = \frac{h_{ii}}{\sqrt{h_{ii}^2+h_{i+1,i}^2}}</math>
 +
 
 +
=== Вычислительное ядро алгоритма ===
 +
Вычислительное ядро последовательной версии метода GMRES состоит из вычислений произведений матрицы на вектор и скалярного произведения в процессе Арнольди и произведений матрицы на вектор, обращений верхней треугольной матрицы для нахождения решения задачи минимизации невязки.
  
−
с произвольной невырожденной матрицей <math>A</math> ищется в виде:
+
=== Макроструктура алгоритма ===
 +
В методе используются следующие макрооперции:
 +
* Умножение матрицы на вектор
 +
* Вычисление скалярного произведение двух векторов
 +
* Нахождени <math>l_2</math>-нормы вектора
 +
* Умножение матриц специальной структуры
 +
* Решение системы уравнений с верхней треугольной матрицей
  
 +
=== Схема реализации последовательного алгоритма ===
 +
Для  обобщенного метода мимнимальных невязок число векторов порядка <math>n</math>, которых необходимо хранить в памяти,  растет линейно с ростом <math>k</math>, в то время как число умножений имеет квадратичную зависимость от числа итераций. Этот факт является мотивацией для того, чтобы использовать метод минимальных невязок с рестартами. Вместо того, чтобы формировать базис размерности <math>k</math>, выбирается параметр <math>m << n </math> и с помощью ортонормированного базиса размерности <math>m</math> строится приближенное решение. Посредством такой схемы необходимое количество памяти может быть оценено заранее и оставаться внутри определенных пределов. Приведем подробное описание алгоритма GMRES c рестартами:
  
−
Исходные данные: Матрица <math>A</math> (элементы <math>a_{ij}</math>), начальное приближение <math> x_0 </math>
+
1. <math>k = 0 </math>; инициализировать <math>x_0</math>; <math>r_0=b-Ax_0</math> и <math>v_1=\frac{r_0}{\|r_0 \|}</math>
  
−
Вычисляемые данные: Приближенное решение уравнения <math>x_k</math>.
+
2. Выполнить <math>m</math> шагов процессса Арнольди:
 +
 
 +
  Для <math> j  = \overline{1,m} </math> выполнять:
 +
      <math>
 +
      \omega = Av_j
 +
      </math>
 +
      Для <math> i = \overline{1,j} </math> выполнять:
 +
          <math>
 +
          \begin{align}
 +
            &h_{ij} = (\omega, v_i)  \\
 +
            &\omega = \omega - h_{ij}v_i \\
 +
            \end{align}
 +
          </math>
 +
      <math>
 +
          \begin{align}
 +
          h_{(i+1)j} = \|\omega\| \\
 +
          v_{j+1} = \frac{\omega}{h_{(i+1)j}} \\
 +
        \end{align}
 +
      </math> 
 +
 
 +
3. Вычислить <math>y_m=\arg\min_{y} \| \beta e_1 - H_m y \|,\ y \in \R^m</math> и <math>x_m=x_0+V_my_m</math>
 +
 
 +
4. Вычислить <math>r_m=b-Ax_m</math>. Если требуемая точность достигнута, остановиться
 +
 
 +
5. Рестарт:
 +
  <math>x_0 = x_m</math>
 +
  <math>r_m = \frac{r_m}{\|r_m\|}</math>
 +
  Перейти к шагу 2.
  
−
=== Вычислительное ядро алгоритма ===
 
−
=== Макроструктура алгоритма ===
 
−
=== Схема реализации последовательного алгоритма ===
 
 
=== Последовательная сложность алгоритма ===
 
=== Последовательная сложность алгоритма ===
 +
Вычислительная сложность метода GMRES складывается из сложности решения задачи минимизации невязки и сложности процесса Арнольди. Так как в 
 +
задаче минимизации невязки матрица системы имеет размер <math>m+1 \times m</math>, для нахождения решения требуется  <math>O(m^2)</math> операций. Процесс
 +
Арнольди предполагает выполнение <math>O(mn)</math> операций плюс <math>O(n^2)</math> операций для умножений матрицы на вектор.
 +
 
=== Информационный граф ===
 
=== Информационный граф ===
 +
Основным ресурсом параллелизма будут являться умножение матрицы на вектор<ref>https://algowiki-project.org/ru/Умножение_плотной_неособенной_матрицы_на_вектор_(последовательный_вещественный_вариант)</ref> и вычисление скалярных произведений<ref>https://algowiki-project.org/ru/Скалярное_произведение_векторов,_вещественная_версия,_последовательно-параллельный_вариант</ref>. Представим информационные графы для данных алгоритмов:
 +
 +
<center>
 +
<div class="thumb">
 +
<div class="thumbinner" style="width:{{#expr: (700 + 15) + 3 * (3 - 1) + 8}}px">
 +
<gallery widths=600px heights=600px>
 +
File:Matvec.png|Рисунок 1. Граф последовательного умножения плотной матрицы на вектор с отображением входных и выходных данных
 +
File:Dot.png|Рисунок 2. Последовательно-параллельное вычисление скалярного произведения
 +
</gallery>
 +
<div class="thumbcaption" style="text-align:center">
 +
</div>
 +
</div>
 +
</div>
 +
</center>
 +
 
=== Ресурс параллелизма алгоритма ===
 
=== Ресурс параллелизма алгоритма ===
 +
 +
Для '''умножения квадратной матрицы на вектор''' порядка <math>n</math> в параллельном варианте требуется последовательно выполнить следующие ярусы:
 +
 +
* по <math>n</math> ярусов умножений и сложений (в каждом из ярусов — <math>n</math> операций).
 +
 +
Для умножения матрицы размером <math>n</math> строк на <math>n</math> столбцов на вектор порядка <math>n</math> в последовательном (наиболее быстром) варианте требуется:
 +
 +
* по <math>n</math> ярусов умножений и сложений (в каждом из ярусов — <math>n</math> операций).
 +
 +
При этом использование режима накопления требует совершения умножений и сложений в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения алгоритма в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает некоторое увеличение требуемой памяти.
 +
 +
При классификации по высоте ЯПФ, таким образом, алгоритм умножения матрицы на вектор относится к алгоритмам ''с линейной сложностью''. При классификации по ширине ЯПФ его сложность также будет ''линейной''.
 +
 +
Для '''вычисления скалярного произведения векторов размерности <math>n</math>''' параллельном варианте требуется последовательно выполнить следующие ярусы:
 +
* 1 ярус вычисления произведений,
 +
* <math>k - 1</math> ярусов суммирования по частям массивов (<math>p</math> ветвей),
 +
* <math>p - 1</math> ярусов суммирования (одна последовательная ветвь),
 +
где <math>p</math> - количество процессоров, <math>k</math>- номер процессора.
 +
 +
Таким образом, в параллельном варианте критический путь алгоритма  (и соответствующая ему высота ЯПФ) будет зависеть от произведённого разбиения массива на части. В оптимальном случае (<math>p = \sqrt{n}</math>) высота ЯПФ будет равна <math> 2 \sqrt{n} - 1</math>. При классификации по высоте ЯПФ, таким образом, последовательно-параллельный метод относится к алгоритмам ''со сложностью «корень квадратный»''. При классификации по ширине ЯПФ его сложность также будет ''«корень квадратный»''.
 +
 
=== Входные и выходные данные алгоритма ===
 
=== Входные и выходные данные алгоритма ===
 +
'''Входные данные:''' Матрица <math>A \in R^{n\times n}</math>, вектор <math>b \in R^n</math>, точность <math>\varepsilon</math>. Объём: <math>n^2 + n + 1</math>
 +
 +
'''Выходные данные:''' Вектор <math>x_m \in R^n</math>. Объём: <math>n</math>
 +
 
=== Свойства алгоритма ===
 
=== Свойства алгоритма ===
−
In mathematics, the generalized minimal residual method (usually abbreviated GMRES) is an iterative method for the numerical solution of a nonsymmetric system of linear equations. The method approximates the solution by the vector in a Krylov subspace with minimal residual. The Arnoldi iteration is used to find this vector.
+
В точной арифметике GMRES всегда сходится не более чем за <math>n </math> итераций (т.к. <math>K_n = R^n = Im(A)</math>). При этом  последовательность решений сходится монотонно, в связи с тем, что <math>\|r_{m+1}\| \le \|r_m\|</math>. Этот факт справедлив из-за того, что подпространства, по которым минимизируется <math>\|r_m\|</math>, вложены друг в друга: <math>K_m \subset K_{m+1}</math>. Таким образом, минимизация по подпрастранствам большей размерности приведет нас к меньшей норме невязки.
  
−
The GMRES method was developed by Yousef Saad and Martin H. Schultz in 1986.[1] GMRES is a generalization of the MINRES method developed by Chris Paige and Michael Saunders in 1975. GMRES also is a special case of the DIIS method developed by Peter Pulay in 1980. DIIS is also applicable to non-linear systems.
+
Скорость сходимости метода зависит от расположения собственных значений матрицы <math>A</math> на комплексной плоскости. Для быстрой сходимости, собственные значения  должны быть сконцентрированны на достаточном расстоянии от начала координат. Заметим, что расположение собственных чисел намного важнее, чем число обусловленности матрицы <math>A</math>, являющимся главным критерием быстрой сходимости для метода сопряженных градиентов.
 +
 
 +
Если обобщенный метод минимальных невязок применяется к СЛАУ с матрицами, обладающими "хорошим" распределением собственных чисел, то он может работать быстрее чем стандартное <math>LU</math> разложение даже в том случае, если матрица плотна и неструктурированна.
  
 
== Программная реализация алгоритма ==
 
== Программная реализация алгоритма ==
 +
 
=== Особенности реализации последовательного алгоритма ===
 
=== Особенности реализации последовательного алгоритма ===
 +
 
=== Локальность данных и вычислений ===
 
=== Локальность данных и вычислений ===
−
==== Локальность реализации алгоритма ====
 
−
===== Структура обращений в память и качественная оценка локальности =====
 
−
===== Количественная оценка локальности =====
 
 
=== Возможные способы и особенности параллельной реализации алгоритма ===
 
=== Возможные способы и особенности параллельной реализации алгоритма ===
 
=== Масштабируемость алгоритма и его реализации ===
 
=== Масштабируемость алгоритма и его реализации ===
−
==== Масштабируемость алгоритма ====
 
−
==== Масштабируемость реализации алгоритма ====
 
 
=== Динамические характеристики и эффективность реализации алгоритма ===
 
=== Динамические характеристики и эффективность реализации алгоритма ===
 
=== Выводы для классов архитектур ===
 
=== Выводы для классов архитектур ===
 
=== Существующие реализации алгоритма ===
 
=== Существующие реализации алгоритма ===
 +
Реализации обобщенного метода минимальных невязок (GMRES) присутствуют во многих пакетах. Приведем некоторые из них:
 +
* [http://www.netlib.org/linpack/ LINPACK]
 +
* [http://math.nist.gov/tnt/download.html/ JAMA/C++]
 +
* [http://math-atlas.sourceforge.net/ ATLAS]
 +
* [http://software.intel.com/mkl‎/ Intel MLK]
 +
* [http://math.nist.gov/iml++/ IML++]
 +
* [http://netlib.org/lapack/ LAPACK]
 +
* [http://osl.iu.edu/research/mtl/ MTL]
 +
* [http://www-users.cs.umn.edu/~saad/software/SPARSKIT/ SPARSKIT]
 +
* [http://researcher.watson.ibm.com/researcher/view_group.php?id=1426 WSMP]
 +
* [https://www.mcs.anl.gov/petsc/ PETSC]
 +
* [https://www.scipy.org// SciPy]
 +
* [https://www.mathworks.com/ Matlab]
 +
* [https://fenicsproject.org/ FEniCS Project]
  
 
== Литература ==
 
== Литература ==
 +
 +
[1] A parallel implementation of the restarted GMRES iterative algorithm for nonsymmetric systems of linear equations, Rudnei Dias da Cunha, Centro de Processamento de Dados, Universidade Federal do Rio Grande do Sul, Brazil and Computing Laboratory, University of Kent, Canterbury, UK and Tim Hopkins, Computing Laboratory, University of Kent, Canterbury, UK, 1993
 +
[2] Gene H. Golub, Charles F. Van Loan, Matrix and Computations, 4th edition

Версия 21:16, 16 ноября 2016

Symbol wait.svgЭта работа ждет рассмотрения преподавателем
Дата последней правки страницы: 16.11.2016
Авторы этой статьи считают, что задание выполнено.

Итерационный метод решения системы линейных алгебраических уравнений GMRES (обобщенный метод минимальных невязок).



Обобщеннный метод минимальных невязок
Последовательный алгоритм
Последовательная сложность [math]O(mn+m^2+n^2)[/math]
Объём входных данных [math]n^2+n+1[/math]
Объём выходных данных [math]n[/math]

Основные авторы описания: Э.А.Ахиярова

Содержание

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

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

Обобщенный метод минимальных невязок (GMRES) - это итерационный метод для нахождения приближенного решения системы линейных алгебраических уравнений [math]Ax = b [/math] с произвольной матрицей [math]A[/math]. Алгоритм был преложен в 1986 году Йозефом Саадом и Мартином Г. Шульцем [1].

В то время как применение классических итерационных алгоритмов для решения СЛАУ ограничивается либо диагональным преобладанием, либо положительной определенностью матрицы, GMRES может быть использован для систем с произвольной невырожденной матрицей [math]A[/math]. Данный метод опирается на итерации Арнольди для нахождения минимизирующего невязку вектора из подпространства Крылова [2].

Основным недостатком метода GMRES является поитерационное увеличение объема требуемой памяти. В связи с этим, для некоторых вычислительных систем выполенние программы метода может быть невозможным. Для некоторых систем используют модифицированный обобщенный метод минимальных невязок - GMRES(m), который подразумевает использование рестартов, существенно сокращающих объемы необходимой памяти.

Ключевая идея обобщенного метода минимальных невязок основана на решении задачи наименьших квадратов на каждом итерационном шаге. На каждом [math]k[/math]-м шаге точное решение [math]x^* = A^{-1}b [/math] аппроксимируется вектором [math]x_k \in K_k[/math], таким образом, что:

[math] \|r_k{\|}_2 = \|Ax_k-b{\|}_2 \to \min [/math],

где [math]K_k[/math] - подпространство Крылова [math] k [/math] - го порядка.

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

Исходные данные: Матрица [math]A \in R^{n\times n}[/math], вектор правой части [math]b \in R^n[/math].

Вычисляемые данные: Приближенное решение [math]\tilde{x} \in R^n[/math] системы [math]Ax = b[/math].

На каждом [math]k[/math]-м шаге метода решение [math]x_k[/math] рассматривается в виде [math]x_0+Vy_k[/math], где столбцы [math]v_k[/math] матрицы [math]V[/math] находятся с помощью процесса Арнольди для построения ортонормированного базиса [math]\{v_1, v_2, \dots, v_k\}[/math] подпространства Крылова [math]K_k(A, v_1)[/math], вектор [math]y_k[/math] представляет собой решение системы [math]H_ky_k = \beta_k e_1[/math], где [math]H_k[/math] - верхняя матрица Хессенберга[3] размера [math]k\times k[/math], элементы которой формируются в процессе ортогонализации Арнольди, [math]\beta_k = \|r_0{\|}_2, \ r_0 = b - Ax_0 [/math], и, наконец, [math]e_1[/math] - канонический вектор [math]{\begin{pmatrix}1&0&\dots 0\end{pmatrix}}^T[/math] размерности [math]k[/math]

Рассмотрим подробнее процесс ортогонализации Арнольди:

1.2.1 Процесс Арнольди

Исходные данные: начальный ортонормированный вектор [math]v_1[/math] и размерность [math]k[/math]

Вычисляемые данные: Ортонормированный базис [math]\{v_1, v_2, \dots, v_k\}[/math] подпространства [math]K_k(A, v_1)[/math]

Формулы метода:

[math] \begin{align} &\omega = Av_j \\ &h_{ij} = (\omega, v_i) \\ &\omega = \omega - h_{ij}v_i \\ &h_{(i+1)j} = \|\omega\| \\ &v_{j+1} = \frac{\omega}{h_{(i+1)j}} \\ &j = \overline{1,k}, \ i = \overline{1,j} \end{align} [/math]

В результате выполнения процесса Арнольди получаем:

[math]V_k =\{v_1, v_2, \dots, v_k\}[/math] - ортонормированный базис подпространства Крылова
[math]\overline{V_k} = \{v_1, v_2, \dots, v_k, v_{k+1}\}[/math] - матрица, расширенная последним вычисленным вектором [math]v_{k+1}[/math]
[math]H_k[/math] - верхняя Хессенбергова матрица, элементы которой равны коэффициентам ортогонализации:
[math]H_k = \begin{pmatrix} h_{11} & h_{12} & h_{13} & h_{14} & \dots & h_{1k} \\ h_{21} & h_{22} & h_{23} & h_{24} & \dots & h_{2k} \\ 0 & h_{32} & h_{33} & h_{34} & \dots & h_{3k} \\ \vdots & &\ddots & \ddots & & \vdots \\ 0 & \dots & 0 & h_{k-1,k-2}& h_{k-1, k-1} & h_{k-1, k} \\ 0 & \dots & 0 & 0 & h_{k, k-1} & h_{k, k} \\ 0 & \dots & & & 0& h_{k+1, k} \end{pmatrix} [/math]

1.2.2 Минимизация нормы невязки

Выполнив несложные преобразования можно получить выражение для невязки на [math]k[/math]- й итерации:

[math] r_k = \beta e_1 - H_ky_k [/math]

Очевидно, что для задачи

[math] \|r_k{\|}_2 \to \min [/math]

минимум достигается при

[math] H_ky_k = \beta e_1 [/math]

Для решения последней системы матрица [math]H_k[/math] приводится к верхнему треугольному виду с помощью [math] k [/math] вращений Гивенса:

[math] R_k = \Omega_k\Omega_{k-1}\dots \Omega_1 [/math],
[math] \Omega_i = \begin{pmatrix} 1 & &&&&&& \\ &\ddots &&&&&&\\ &&1&&&&& \\ &&&c_i&s_i &&& \\ &&&-s_i& c_i &&&\\ &&&&&1&&& \\ &&&&&&\ddots&\\ &&&&&&&1 \\ \end{pmatrix}, \quad s_i = \frac{h_{i+1,i}}{\sqrt{h_{ii}^2+h_{i+1,i}^2}}, \quad c_i = \frac{h_{ii}}{\sqrt{h_{ii}^2+h_{i+1,i}^2}}[/math]

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

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

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

В методе используются следующие макрооперции:

  • Умножение матрицы на вектор
  • Вычисление скалярного произведение двух векторов
  • Нахождени [math]l_2[/math]-нормы вектора
  • Умножение матриц специальной структуры
  • Решение системы уравнений с верхней треугольной матрицей

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

Для обобщенного метода мимнимальных невязок число векторов порядка [math]n[/math], которых необходимо хранить в памяти, растет линейно с ростом [math]k[/math], в то время как число умножений имеет квадратичную зависимость от числа итераций. Этот факт является мотивацией для того, чтобы использовать метод минимальных невязок с рестартами. Вместо того, чтобы формировать базис размерности [math]k[/math], выбирается параметр [math]m \lt \lt n [/math] и с помощью ортонормированного базиса размерности [math]m[/math] строится приближенное решение. Посредством такой схемы необходимое количество памяти может быть оценено заранее и оставаться внутри определенных пределов. Приведем подробное описание алгоритма GMRES c рестартами:

1. [math]k = 0 [/math]; инициализировать [math]x_0[/math]; [math]r_0=b-Ax_0[/math] и [math]v_1=\frac{r_0}{\|r_0 \|}[/math]

2. Выполнить [math]m[/math] шагов процессса Арнольди: 

  Для [math] j  = \overline{1,m} [/math] выполнять:
     [math]
      \omega = Av_j 
      [/math]
     Для [math] i = \overline{1,j} [/math] выполнять: 
          [math]
           \begin{align}
            &h_{ij} = (\omega, v_i)   \\
            &\omega = \omega - h_{ij}v_i \\
            \end{align}
           [/math]
      [math]
          \begin{align}
           h_{(i+1)j} = \|\omega\| \\
           v_{j+1} = \frac{\omega}{h_{(i+1)j}} \\
         \end{align}
      [/math]

3. Вычислить [math]y_m=\arg\min_{y} \| \beta e_1 - H_m y \|,\ y \in \R^m[/math] и [math]x_m=x_0+V_my_m[/math]

4. Вычислить [math]r_m=b-Ax_m[/math]. Если требуемая точность достигнута, остановиться

5. Рестарт: 

  [math]x_0 = x_m[/math]
  [math]r_m = \frac{r_m}{\|r_m\|}[/math]
  Перейти к шагу 2.

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

Вычислительная сложность метода GMRES складывается из сложности решения задачи минимизации невязки и сложности процесса Арнольди. Так как в задаче минимизации невязки матрица системы имеет размер [math]m+1 \times m[/math], для нахождения решения требуется [math]O(m^2)[/math] операций. Процесс Арнольди предполагает выполнение [math]O(mn)[/math] операций плюс [math]O(n^2)[/math] операций для умножений матрицы на вектор.

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

Основным ресурсом параллелизма будут являться умножение матрицы на вектор[4] и вычисление скалярных произведений[5]. Представим информационные графы для данных алгоритмов:

  • Рисунок 1. Граф последовательного умножения плотной матрицы на вектор с отображением входных и выходных данных

  • Рисунок 2. Последовательно-параллельное вычисление скалярного произведения

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

Для умножения квадратной матрицы на вектор порядка [math]n[/math] в параллельном варианте требуется последовательно выполнить следующие ярусы:

  • по [math]n[/math] ярусов умножений и сложений (в каждом из ярусов — [math]n[/math] операций).

Для умножения матрицы размером [math]n[/math] строк на [math]n[/math] столбцов на вектор порядка [math]n[/math] в последовательном (наиболее быстром) варианте требуется:

  • по [math]n[/math] ярусов умножений и сложений (в каждом из ярусов — [math]n[/math] операций).

При этом использование режима накопления требует совершения умножений и сложений в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения алгоритма в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает некоторое увеличение требуемой памяти.

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

Для вычисления скалярного произведения векторов размерности [math]n[/math] параллельном варианте требуется последовательно выполнить следующие ярусы:

  • 1 ярус вычисления произведений,
  • [math]k - 1[/math] ярусов суммирования по частям массивов ([math]p[/math] ветвей),
  • [math]p - 1[/math] ярусов суммирования (одна последовательная ветвь),

где [math]p[/math] - количество процессоров, [math]k[/math]- номер процессора.

Таким образом, в параллельном варианте критический путь алгоритма (и соответствующая ему высота ЯПФ) будет зависеть от произведённого разбиения массива на части. В оптимальном случае ([math]p = \sqrt{n}[/math]) высота ЯПФ будет равна [math] 2 \sqrt{n} - 1[/math]. При классификации по высоте ЯПФ, таким образом, последовательно-параллельный метод относится к алгоритмам со сложностью «корень квадратный». При классификации по ширине ЯПФ его сложность также будет «корень квадратный».

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

Входные данные: Матрица [math]A \in R^{n\times n}[/math], вектор [math]b \in R^n[/math], точность [math]\varepsilon[/math]. Объём: [math]n^2 + n + 1[/math]

Выходные данные: Вектор [math]x_m \in R^n[/math]. Объём: [math]n[/math]

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

В точной арифметике GMRES всегда сходится не более чем за [math]n [/math] итераций (т.к. [math]K_n = R^n = Im(A)[/math]). При этом последовательность решений сходится монотонно, в связи с тем, что [math]\|r_{m+1}\| \le \|r_m\|[/math]. Этот факт справедлив из-за того, что подпространства, по которым минимизируется [math]\|r_m\|[/math], вложены друг в друга: [math]K_m \subset K_{m+1}[/math]. Таким образом, минимизация по подпрастранствам большей размерности приведет нас к меньшей норме невязки.

Скорость сходимости метода зависит от расположения собственных значений матрицы [math]A[/math] на комплексной плоскости. Для быстрой сходимости, собственные значения должны быть сконцентрированны на достаточном расстоянии от начала координат. Заметим, что расположение собственных чисел намного важнее, чем число обусловленности матрицы [math]A[/math], являющимся главным критерием быстрой сходимости для метода сопряженных градиентов.

Если обобщенный метод минимальных невязок применяется к СЛАУ с матрицами, обладающими "хорошим" распределением собственных чисел, то он может работать быстрее чем стандартное [math]LU[/math] разложение даже в том случае, если матрица плотна и неструктурированна.

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

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

2.2 Локальность данных и вычислений

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

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

2.5 Динамические характеристики и эффективность реализации алгоритма

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

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

Реализации обобщенного метода минимальных невязок (GMRES) присутствуют во многих пакетах. Приведем некоторые из них:

3 Литература

[1] A parallel implementation of the restarted GMRES iterative algorithm for nonsymmetric systems of linear equations, Rudnei Dias da Cunha, Centro de Processamento de Dados, Universidade Federal do Rio Grande do Sul, Brazil and Computing Laboratory, University of Kent, Canterbury, UK and Tim Hopkins, Computing Laboratory, University of Kent, Canterbury, UK, 1993

[2] Gene H. Golub, Charles F. Van Loan, Matrix and Computations, 4th edition

  1. ↑ Y. Saad and M.H. Schultz, "GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems", SIAM J. Sci. Stat. Comput., 7:856-869, 1986. doi:10.1137/0907058.
  2. ↑ https://en.wikipedia.org/wiki/Krylov_subspace
  3. ↑ https://en.wikipedia.org/wiki/Hessenberg_matrix
  4. ↑ https://algowiki-project.org/ru/Умножение_плотной_неособенной_матрицы_на_вектор_(последовательный_вещественный_вариант)
  5. ↑ https://algowiki-project.org/ru/Скалярное_произведение_векторов,_вещественная_версия,_последовательно-параллельный_вариант
Источник — https://algowiki-project.org/w/ru/index.php?title=Участница:Эльвира_Ахиярова/Итерационный_метод_решения_системы_линейных_алгебраических_уравнений_GMRES_(обобщенный_метод_минимальных_невязок)&oldid=17766