Участник:Alexander34396/Обобщенный метод минимальных невязок: различия между версиями
| Строка 134: | Строка 134: | ||
'''Входные данные:''' | '''Входные данные:''' | ||
* <math> A </math> – квадратная матрица размера <math> n </math>-на-<math> n </math>; | * <math> A </math> – квадратная матрица размера <math> n </math>-на-<math> n </math>; | ||
* <math> b </math> – вектор правой части размера <math> n </math> | * <math> b </math> – вектор правой части размера <math> n </math>; | ||
* <math> x_0 </math> – начальное приближение размера <math> n </math> | * <math> x_0 </math> – начальное приближение размера <math> n </math>; | ||
* <math> \epsilon </math> – погрешность решения. | * <math> \epsilon </math> – погрешность решения. | ||
| Строка 144: | Строка 144: | ||
'''Выходные данные''': | '''Выходные данные''': | ||
* <math> x_k </math> - вектор приближённого решения размера <math> n </math> | * <math> x_k </math> - вектор приближённого решения размера <math> n </math>. | ||
'''Объём выходных данных''': <math> n </math>. | '''Объём выходных данных''': <math> n </math>. | ||
Версия от 15:54, 17 октября 2016
| Файл:Symbol wait.svg | Эта работа ждет рассмотрения преподавателем Дата последней правки страницы: 17.10.2016 Авторы этой статьи считают, что задание выполнено. |
Свойства и структура алгоритма
Общее описание алгоритма
Обобщённый метод минимальных невязок (англ. Generalized minimal residual method, GMRES) - итерационный метод численного решения системы линейных алгебраических уравнений с невырожденной матрицей. Метод основан на минимизации квадратичного функционала невязки на подпространствах Крылова. Разработан Юсефом Саадом и Мартином Шульцем[1] в 1986 году как обобщение метода MINRES[2] на случай систем с несимметричными матрицами.
Математическое описание алгоритма
Исходные данные:
- система линейных алгебраических уравнений вида [math]\displaystyle{ Ax = b }[/math], где [math]\displaystyle{ A }[/math] — невырожденная матрица размера [math]\displaystyle{ n }[/math]-на-[math]\displaystyle{ n }[/math];
- [math]\displaystyle{ x_0 }[/math] - начальное приближение.
Вычисляемые данные:
- [math]\displaystyle{ x_k }[/math] - приближённое решение исходной системы.
Метод GMRES приближает точное решение исходной системы [math]\displaystyle{ Ax = b }[/math] вектором [math]\displaystyle{ x_k }[/math], который минимизирует [math]\displaystyle{ l_2 }[/math]-норму невязки [math]\displaystyle{ r_k = Ax_k - b }[/math] на подпространстве Крылова:
- [math]\displaystyle{ K_k = K_k(A,b) = \operatorname{span} \, \{ b, Ab, A^2b, \ldots, A^{k-1}b \} }[/math].
На каждой итерации метода решение уточняется поправкой, представленной в виде разложения по [math]\displaystyle{ l_2 }[/math]-ортонормированному базису пространства [math]\displaystyle{ K_k }[/math]:
- [math]\displaystyle{ x_k = x_0 + z_k }[/math], где [math]\displaystyle{ x_0 }[/math] - некоторое начальное приближение, [math]\displaystyle{ z_k \in K_k }[/math] - поправка решения.
Ортогонализация Арнольди
Для построения ортонормированного базиса [math]\displaystyle{ K_k }[/math] метод использует процесс Арнольди.
Исходные данные:
- [math]\displaystyle{ v_1 }[/math], такой что [math]\displaystyle{ \|v_1\|_2 = 1 }[/math];
- матрица [math]\displaystyle{ A }[/math] размером [math]\displaystyle{ n }[/math]-на-[math]\displaystyle{ n }[/math].
Формула процесса:
- [math]\displaystyle{ h_{j+1j}v_{j+1} = Av_j - \sum_{i=1}^j h_{ij}v_{j} }[/math], где [math]\displaystyle{ h_{ij} = (Av_j, v_i) }[/math], [math]\displaystyle{ \quad j=1,\ldots,k }[/math].
Алгоритм процесса:
Выполнять для [math]\displaystyle{ \quad j=1,\ldots,k }[/math]:
- [math]\displaystyle{ h_{ij} = (Av_j, v_i), \quad i=1,\ldots,j }[/math];
- [math]\displaystyle{ \hat{v}_{j+1} = Av_j - \sum_{i=1}^j h_{ij}v_{i} }[/math];
- [math]\displaystyle{ h_{j+1j} = \|\hat{v}_{j+1}\|_2 }[/math];
- [math]\displaystyle{ v_{j+1} = \frac{\hat{v}_{j+1}}{h_{j+1j}} }[/math].
При введении матричных обозначений можно записать:
- [math]\displaystyle{ z_k = V_ky_k }[/math], где [math]\displaystyle{ y_k \in \mathbb{R}^k }[/math] - вектор коэффициентов;
- [math]\displaystyle{ AV_k = V_kH_k + h_{k+1k}v_{k+1} e_k^T = V_{k+1}\overline{H}_k }[/math], где [math]\displaystyle{ V_k = [v_1|v_2|...|v_k] }[/math], а [math]\displaystyle{ \overline{H}_k }[/math] - верхняя матрица Хессенберга размерности [math]\displaystyle{ (j+1) }[/math] на [math]\displaystyle{ j }[/math].
Минимизация функционала невязки
Для решения исходной системы GMRES вычисляет приближённое решение [math]\displaystyle{ x_k }[/math], которое минимизизует функционал невязки:
- [math]\displaystyle{ J(y) = \|b- A x_k\|_2 = \|b - A(x_0 + V_ky)\|_2 = \|\beta e_1 - \overline{H}_ky_k \|_2 }[/math], где [math]\displaystyle{ \beta = \|r_0\|_2 }[/math].
Вектор [math]\displaystyle{ y_k }[/math] может найден как решение линейной задачи наименьших квадратов размера [math]\displaystyle{ k+1 }[/math]-на-[math]\displaystyle{ k }[/math]:
- [math]\displaystyle{ y_k = \arg\min_{y}\| \beta e_1 - \overline{H}_ky\|_2 }[/math].
Для решения задачи матрица [math]\displaystyle{ \overline{H}_k }[/math] приводится к верхнему треугольному виду с помощью [math]\displaystyle{ k }[/math] последовательных вращений Гивенса.
Общая схема метода
В общем виде k-aя итерация алгоритма GMRES может быть записан следующим образом:
- найти ортонормированный базис [math]\displaystyle{ V_k }[/math] подпространства [math]\displaystyle{ K_k }[/math] с помощью ортогонализации Арнольди;
- найти [math]\displaystyle{ y_k }[/math], минимизирующий [math]\displaystyle{ \|r_k\|_2 }[/math];
- вычислить [math]\displaystyle{ x_k = x_0 + V_ky_k }[/math];
- вычислить [math]\displaystyle{ r_k }[/math] и остановиться,если требуемая точность была достигнута, иначе повторить.
Вычислительное ядро алгоритма
Вычислительное ядро последовательной версии метода GMRES состоит из двух частей:
- Вычисление ортонормированного базиса [math]\displaystyle{ K_k }[/math];
- Формирование приближенного решения [math]\displaystyle{ x_k }[/math].
Для вычисления ортонормированного базиса [math]\displaystyle{ K_k }[/math] метод использует процесс Арнольди:
- [math]\displaystyle{ h_{j+1j}v_{j+1} = Av_j - \sum_{i=1}^j h_{ij}v_{j} }[/math], где [math]\displaystyle{ h_{ij} = (Av_j, v_i) }[/math], [math]\displaystyle{ \quad j=1,\ldots,k }[/math].
Для нахождения приближённого решения [math]\displaystyle{ x_k }[/math] метод использует формулу:
- [math]\displaystyle{ x_k = x_0 + V_ky_k }[/math].
Макроструктура алгоритма
В алгоритме можно выделить следующие макрооперации:
- Умножение матрицы на вектор;
- Вычисление скалярного произведения векторов;
- Вычисление [math]\displaystyle{ l_2 }[/math]-нормы вектора;
- Умножение вектора на скаляр;
- Деление вектора на скаляр.
Схема реализации последовательного алгоритма
На практике используется перезапускаемая версия метода - GMRES(m):
1 Подготовка перед итерационным процессом:
- 1.1 Выбрать начальное приближение [math]\displaystyle{ x_0 }[/math];
- 1.2 Посчитать невязку [math]\displaystyle{ r_0 = b - Ax_0 }[/math];
- 1.3 Вычислить [math]\displaystyle{ v_1 = \frac{r_0}{\|r_0\|_2} }[/math].
2 Построение ортонормированного базиса [math]\displaystyle{ K_m }[/math]:
- Для всех [math]\displaystyle{ j }[/math] от 1 до m по нарастанию выполнять:
- 2.1 [math]\displaystyle{ h_{ij} := (Av_j, v_i), \quad i=1,\ldots,j }[/math];
- 2.2 [math]\displaystyle{ \hat{v}_{j+1} := Av_j - \sum_{i=1}^j h_{ij}v_{i} }[/math];
- 2.3 [math]\displaystyle{ h_{j+1j} = \|\hat{v}_{j+1}\|_2 }[/math];
- 2.4 [math]\displaystyle{ v_{j+1} = \frac{\hat{v}_{j+1}}{h_{j+1j}} }[/math].
3 Вычисление [math]\displaystyle{ x_m }[/math]:
- 3.1 [math]\displaystyle{ x_m = x_0 + V_my_m }[/math], где [math]\displaystyle{ y_m }[/math] минимизирует [math]\displaystyle{ \|r_0 - AV_my_m\|_2 }[/math];
- 3.2 Вычислить [math]\displaystyle{ r_m }[/math];
- 3.3 Если требуемая точность достигнута, остановиться.
4 Перезапуск:
- 4.1 [math]\displaystyle{ x_0 = x_m }[/math];
- 4.2 [math]\displaystyle{ v_1 = \frac{r_m}{\|r_m\|_2} }[/math];
- 4.3 Перейти к шагу 2.
Последовательная сложность алгоритма
Если пренебречь сложностью вычисления [math]\displaystyle{ y_k }[/math], то общую сложность алгоритма GMRES можно разделить на две части:
1. Сложность вычисления ортонормированного базиса пространства [math]\displaystyle{ K_k }[/math]:
- a) Для вычисления [math]\displaystyle{ j }[/math]-го вектора базиса [math]\displaystyle{ K_k, j \lt k }[/math] требуется:
- [math]\displaystyle{ NZ + (2j + 1)n }[/math] мультипликативных операций, где [math]\displaystyle{ NZ }[/math] - количество ненулевых элементов матрицы [math]\displaystyle{ A }[/math].
- b) Вычисление последнего вектора базиса требует:
- [math]\displaystyle{ n(k + 1) }[/math] мультипликативных операций.
- c) Общая мультипликативная сложность вычисления ортонормированного базиса [math]\displaystyle{ K_k }[/math]:
- [math]\displaystyle{ nk(k + 1) + kNZ }[/math].
2. Сложность вычисления приближённого решения [math]\displaystyle{ x_k }[/math]:
- a) Вычисление этой формулы требует [math]\displaystyle{ nk }[/math] мультипликативных операций.
Информационный граф
Ресурс параллелизма алгоритма
Входные и выходные данные алгоритма
Входные данные:
- [math]\displaystyle{ A }[/math] – квадратная матрица размера [math]\displaystyle{ n }[/math]-на-[math]\displaystyle{ n }[/math];
- [math]\displaystyle{ b }[/math] – вектор правой части размера [math]\displaystyle{ n }[/math];
- [math]\displaystyle{ x_0 }[/math] – начальное приближение размера [math]\displaystyle{ n }[/math];
- [math]\displaystyle{ \epsilon }[/math] – погрешность решения.
Для перезапускаемой версии алгоритма дополнительно:
- [math]\displaystyle{ m }[/math] – число итераций.
Объём входных данных: [math]\displaystyle{ n^2 + 2n + 1 }[/math].
Выходные данные:
- [math]\displaystyle{ x_k }[/math] - вектор приближённого решения размера [math]\displaystyle{ n }[/math].
Объём выходных данных: [math]\displaystyle{ n }[/math].
Свойства алгоритма
Программная реализация алгоритма
Особенности реализации последовательного алгоритма
Локальность данных и вычислений
Возможные способы и особенности параллельной реализации алгоритма
Масштабируемость алгоритма и его реализации
Динамические характеристики и эффективность реализации алгоритма
Выводы для классов архитектур
Существующие реализации алгоритма
Реализации алгоритма GMRES присутствуют во многих программных пакетах:
Литература
<references \>
- ↑ 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)
- ↑ C. C. Paige and M. A. Saunders. Solution of sparse indefinite systems of linear equations, SIAM J. Numerical Analysis 12, 617-629 (1975)