Участник:Alexander34396/Обобщенный метод минимальных невязок: различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
Строка 54: Строка 54:
 
* Вычисление Евклидовой нормы вектора;
 
* Вычисление Евклидовой нормы вектора;
 
* Умножение вектора на скаляр;
 
* Умножение вектора на скаляр;
* Деление вектора на скаляр.
+
* Деление вектора на скаляр;
 +
* Поиск <math> underset{y \in K_k}{\operatorname{argmin}}\|r_0 - AV_ky_k\|_2 </math>
 +
 
  
 
=== Схема реализации последовательного алгоритма ===
 
=== Схема реализации последовательного алгоритма ===
Строка 60: Строка 62:
 
Для решения исходной системы методом GMRES можно воспользоваться следующим алгоритмом:
 
Для решения исходной системы методом GMRES можно воспользоваться следующим алгоритмом:
  
1 Подготовка перед итерационным процессом:
+
'''Подготовка перед итерационным процессом:'''
:1.1 Выбрать начальное приближение <math> x_0 </math>;
+
# Выбрать начальное приближение <math> x_0 </math>;
:1.2 Посчитать невязку <math> r_0 = b - Ax_0 </math>;
+
# Посчитать невязку <math> r_0 = b - Ax_0 </math>;
:1.3 Вычислить <math> v_1 = \frac{r_0}{\|r_0\|_2} </math>.
+
# Вычислить <math> v_1 = \frac{r_0}{\|r_0\|_2} </math>.
 
 
2 Построение ортонормального базиса <math> K_m </math>:
 
 
 
:Для всех <math> j </math> от 1 до m по нарастанию выполнять:
 
:2.1 <math> h_{ij} := (Av_j, v_i), \quad i=1,\ldots,j </math>;
 
:2.2 <math> \hat{v}_{j+1} := Av_j - \sum_{i=1}^j h_{ij}v_{i} </math>;
 
:2.3 <math> h_{j+1j} = \|\hat{v}_{j+1}\|_2 </math>;
 
:2.4 <math> v_{j+1} = \frac{\hat{v}_{j+1}}{h_{j+1j}} </math>.
 
 
 
3 Вычисление приближённого решения <math> x_m </math>:
 
:3.1 <math> x_m = x_0 + V_my_m </math>, где <math>y_m</math> минимизирует <math>\|r_0 - AV_my_m\|_2</math>;
 
:3.2 Вычислить <math> r_m </math>;
 
:3.3 Если требуемая точность достигнута, остановиться.
 
  
4 Рестарт:
+
'''k-ая итерация:'''
:4.1 <math>x_0 = x_m</math>;
+
# <math> h_{ik} := (Av_k, v_i), \quad i=1,\ldots,k </math>;
:4.2 <math>v_1 = \frac{r_m}{\|r_m\|_2}</math>;
+
# <math> \hat{v}_{k+1} := Av_k - \sum_{i=1}^k h_{ik}v_{i} </math>;
:4.3 Перейти к шагу 2.
+
# <math> h_{k+1k} = \|\hat{v}_{k+1}\|_2 </math>;
 +
# <math> v_{k+1} = \frac{\hat{v}_{k+1}}{h_{k+1k}} </math>.
 +
# <math> x_k = x_0 + V_ky_k </math>, где <math>y_k</math> минимизирует <math>\|r_0 - AV_ky_k\|_2</math>

Версия 22:57, 15 октября 2016

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

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

Обобщённый метод минимальных невязок (англ. Generalized minimal residual method, GMRES) - итерационный метод численного решения системы линейных алгебраических уравнений с невырожденной матрицей. Метод основан на минимизации квадратичного функционала невязки на подпространствах Крылова. Разработан Юсефом Саадом и Мартином Шульцем в 1986 году как обобщение метода MINRES на случай систем с несимметричными матрицами.

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

Исходные данные:

  • система линейных алгебраических уравнений вида [math] Ax = b [/math], где [math] A [/math] — невырожденная матрица размера [math]n[/math]-на-[math] n [/math].

Вычисляемые данные:

  • [math] x_k [/math] - приближённое решение исходной системы.

Подпространство Крылова размерности [math] k, k \leq n [/math] для решения исходной системы:

[math] K_k = K_k(A,b) = \operatorname{span} \, \{ b, Ab, A^2b, \ldots, A^{k-1}b \}. \, [/math]

Метод GMRES приближает точное решение исходной системы [math] Ax = b [/math] вектором [math] x_k \in K_k [/math], минимизирующим Евклидову норму невязки [math]r_k = Ax_k-b[/math].

Для решения исходной системы GMRES, используя [math] l_2 [/math]-ортонормальный базис пространства [math] K_k [/math], выполняет поиск приближённого решения [math] x_k [/math] в виде:

[math] x_k = x_0 + z_k [/math],

где [math] x_0 [/math] - некоторое начальное приближение, [math] z_k \in K_k [/math] - поправка решения.

Для построения ортонормального базиса [math] K_k [/math] метод использует ортогонализацию Арнольди. При введении для базиса [math] K_k [/math] матричного обозначения [math] V_k [/math] можно записать:

[math] z_k = V_ky_k [/math],

где [math] y_k \in \mathbb{R}^k [/math] - вектор коэффициентов.

В общем виде k-aя итерация алгоритма GMRES может быть записан следующим образом:

  1. найти ортонормальный базис [math] V_k [/math] подпространства [math] K_k [/math] с помощью ортогонализации Арнольди;
  2. найти [math] y_k [/math], минимизирующий [math] \|r_k\|_2 [/math];
  3. вычислить [math] x_k = x_0 + V_ky_k [/math];
  4. вычислить [math] r_k [/math] и остановиться,если требуемая точность была достигнута, иначе повторить для k + 1.

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

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

  • Вычисление ортонормального базиса [math] K_k [/math];
  • Формирование приближенного решения [math] x_k [/math].

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

[math] \hat{v}_{k+1} := Av_k - \sum_{i=1}^k h_{ik}v_{k} [/math], где [math] h_{ik} := (Av_k, v_i) [/math].

Этот процесс требует [math] k(k + 1) n + kNZ [/math] мультипликативных операций, где NZ - количество ненулевых элементов матрицы [math] A [/math].

Для нахождения на каждой итерации приближённого решения метод использует формулу:

[math] x_k = x_0 + V_ky_k [/math].

Вычисление этой формулы требует [math]nk[/math] мультипликативных операций.

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

В алгоритме можно выделить следующие макрооперации:

  • Умножение матрицы на вектор;
  • Вычисление скалярного произведения векторов;
  • Вычисление Евклидовой нормы вектора;
  • Умножение вектора на скаляр;
  • Деление вектора на скаляр;
  • Поиск [math] underset{y \in K_k}{\operatorname{argmin}}\|r_0 - AV_ky_k\|_2 [/math]


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

Для решения исходной системы методом GMRES можно воспользоваться следующим алгоритмом:

Подготовка перед итерационным процессом:

  1. Выбрать начальное приближение [math] x_0 [/math];
  2. Посчитать невязку [math] r_0 = b - Ax_0 [/math];
  3. Вычислить [math] v_1 = \frac{r_0}{\|r_0\|_2} [/math].

k-ая итерация:

  1. [math] h_{ik} := (Av_k, v_i), \quad i=1,\ldots,k [/math];
  2. [math] \hat{v}_{k+1} := Av_k - \sum_{i=1}^k h_{ik}v_{i} [/math];
  3. [math] h_{k+1k} = \|\hat{v}_{k+1}\|_2 [/math];
  4. [math] v_{k+1} = \frac{\hat{v}_{k+1}}{h_{k+1k}} [/math].
  5. [math] x_k = x_0 + V_ky_k [/math], где [math]y_k[/math] минимизирует [math]\|r_0 - AV_ky_k\|_2[/math]