Обобщенный метод минимальной невязки

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

Метод GMRES был разработан Юсефом Саадом и Мартином Х. Шульцем в 1986 году. ^[1] GMRES является обобщением метода MINRES ^[2], разработанного Крисом Пейджем и Майклом Сондерсом в 1975 году. GMRES также является частным случаем метода DIIS. разработан Питером Пулаем в 1980 году. DIIS также применим к нелинейным системам.

Метод [ править ]

Обозначим евклидову норму любого вектора v через . Обозначим (квадратную) систему линейных уравнений, которую необходимо решить, через${\ Displaystyle \ | v \ |}$

${\ Displaystyle Ax = b. \,}$

Матрица считаются обратимым размером м матрица с размерностью м . Кроме того, предполагается, что b нормализовано, т. Е. Что .${\ Displaystyle \ | б \ | = 1}$

В н -го Крылов подпространство для этой задачи

${\ displaystyle K_ {n} = K_ {n} (A, r_ {0}) = \ operatorname {span} \, \ {r_ {0}, Ar_ {0}, A ^ {2} r_ {0}, \ ldots, A ^ {n-1} r_ {0} \}. \,}$

где - начальная ошибка с учетом первоначального предположения . Ясно, если .${\ displaystyle r_ {0} = b-Ax_ {0}}$${\ displaystyle x_ {0} \ neq 0}$${\ displaystyle r_ {0} = b}$${\ displaystyle x_ {0} = 0}$

GMRES аппроксимирует точное решение вектором, который минимизирует евклидову норму невязки .${\ displaystyle Ax = b}$${\ displaystyle x_ {n} \ in K_ {n}}$ ${\ displaystyle r_ {n} = Ax_ {n} -b}$

Векторы могут быть близки к линейно зависимым , поэтому вместо этого базиса используется итерация Арнольди для поиска ортонормированных векторов, которые образуют основу для . В частности, . ${\ displaystyle r_ {0}, Ar_ {0}, \ ldots A ^ {n-1} r_ {0}}$${\ displaystyle q_ {1}, q_ {2}, \ ldots, q_ {n} \,}$${\ displaystyle K_ {n}}$${\ displaystyle q_ {1} = \ | r_ {0} \ | _ {2} ^ {- 1} r_ {0}}$

Следовательно, вектор может быть записан как с , где - матрица размером m на n, образованная посредством .${\ displaystyle x_ {n} \ in K_ {n}}$${\ displaystyle x_ {n} = x_ {0} + Q_ {n} y_ {n}}$${\ Displaystyle у_ {п} \ в \ mathbb {R} ^ {п}}$${\ displaystyle Q_ {n}}$${\ displaystyle q_ {1}, \ ldots, q_ {n}}$

Процесс Арнольди также дает матрицу Хессенберга ( ) по верхнему краю с${\displaystyle n+1}$${\displaystyle n}$ ${\displaystyle {\tilde {H}}_{n}}$

${\displaystyle AQ_{n}=Q_{n+1}{\tilde {H}}_{n}.\,}$

Поскольку столбцы ортонормированы, мы имеем${\displaystyle Q_{n}}$

${\displaystyle \|r_{n}\|=\|b-Ax_{n}\|=\|b-A(x_{0}+z_{n})\|=\|r_{0}-Az_{n}\|=\|\beta q_{1}-AQ_{n}y_{n}\|=\|\beta q_{1}-Q_{n+1}{\tilde {H}}_{n}y_{n}\|=\|Q_{n+1}(\beta e_{1}-{\tilde {H}}_{n}y_{n})\|=\|\beta e_{1}-{\tilde {H}}_{n}y_{n}\|,\,}$

куда

${\displaystyle e_{1}=(1,0,0,\ldots ,0)^{T}\,}$

первый вектор в стандартном базисе из и ${\displaystyle \mathbb {R} ^{n+1}}$

${\displaystyle \beta =\|b-Ax_{0}\|\,,}$

${\displaystyle x_{0}}$являющийся первым пробным вектором (обычно нулевым). Следовательно, можно найти, минимизируя евклидову норму невязки${\displaystyle x_{n}}$

${\displaystyle r_{n}={\tilde {H}}_{n}y_{n}-\beta e_{1}.}$

Это линейная задача наименьших квадратов размера n .

Это дает метод GMRES. На -й итерации:${\displaystyle n}$

1. рассчитать по методу Арнольди;${\displaystyle q_{n}}$
2. найти то, что минимизирует ;${\displaystyle y_{n}}$${\displaystyle \|r_{n}\|}$
3. вычислить ;${\displaystyle x_{n}=x_{0}+Q_{n}y_{n}}$
4. повторите, если остаток еще недостаточно мал.

На каждой итерации должно быть вычислено произведение матрицы и вектора . Это требует операций с плавающей запятой для плотных матриц общего размера , но для разреженных матриц стоимость может снизиться до . В дополнение к произведению матрицы и вектора на n -й итерации должны быть вычислены операции с плавающей запятой .${\displaystyle Aq_{n}}$${\displaystyle 2m^{2}}$ ${\displaystyle m}$${\displaystyle O(m)}$${\displaystyle O(nm)}$

Конвергенция [ править ]

П - й итерации минимизирует остаточную в подпространстве Крылова К _п . Поскольку каждое подпространство содержится в следующем подпространстве, невязка не увеличивается. После m итераций, где m - размер матрицы A , пространство Крылова K _m - это все R ^m, и, следовательно, метод GMRES приходит к точному решению. Однако идея состоит в том, что после небольшого количества итераций (относительно m ) вектор x _n уже является хорошим приближением к точному решению.

Так не бывает вообще. Действительно, теорема Гринбаума, Патака и Стракоша утверждает, что для любой невозрастающей последовательности a ₁ ,…, a _{m −1} , a _m = 0 можно найти такую ​​матрицу A , что || r _n || = a _n для всех n , где r _n - невязка, определенная выше. В частности, можно найти матрицу, для которой невязка остается постоянной в течение m  - 1 итераций и падает до нуля только на последней итерации.

Однако на практике GMRES часто работает хорошо. Это можно доказать в конкретных ситуациях. Если симметричная часть А , то есть , является положительно определенным , то${\displaystyle (A^{T}+A)/2}$

${\displaystyle \|r_{n}\|\leq \left(1-{\frac {\lambda _{\min }^{2}(1/2(A^{T}+A))}{\lambda _{\max }(A^{T}A)}}\right)^{n/2}\|r_{0}\|,}$

где и обозначают наименьшее и наибольшее собственное значение матрицы соответственно. ^[3]${\displaystyle \lambda _{\mathrm {min} }(M)}$${\displaystyle \lambda _{\mathrm {max} }(M)}$${\displaystyle M}$

Если является симметричной и положительно определена, то у нас даже есть

${\displaystyle \|r_{n}\|\leq \left({\frac {\kappa _{2}(A)^{2}-1}{\kappa _{2}(A)^{2}}}\right)^{n/2}\|r_{0}\|.}$

где обозначает число обусловленности из А в евклидовой норме.${\displaystyle \kappa _{2}(A)}$

В общем случае, когда A не является положительно определенным, имеем

${\displaystyle {\frac {\|r_{n}\|}{\|b\|}}\leq \inf _{p\in P_{n}}\|p(A)\|\leq \kappa _{2}(V)\inf _{p\in P_{n}}\max _{\lambda \in \sigma (A)}|p(\lambda )|,\,}$

где Р _п обозначает множество многочленов степени не выше п с р (0) = 1, V является матрицей , появляющаяся в спектральном разложении из A , и σ ( ) является спектр из A . Грубо говоря, это означает, что быстрая сходимость происходит, когда собственные значения A сгруппированы далеко от начала координат, а A не слишком далеко от нормальности . ^[4]

Все эти неравенства ограничивают только остатки, а не фактическую ошибку, то есть расстояние между текущей итерацией x _n и точным решением.

Расширения метода [ править ]

Как и другие итерационные методы, GMRES обычно сочетается с методом предварительной обработки , чтобы ускорить сходимость.

Стоимость итераций растет как O ( n ² ), где n - номер итерации. Поэтому метод иногда перезапускается после некоторого количества, скажем k , итераций с x _{k в} качестве начального предположения. Полученный метод называется GMRES ( k ) или Restarted GMRES. Для неположительно определенных матриц этот метод может страдать от стагнации сходимости, поскольку перезапущенное подпространство часто близко к предыдущему подпространству.

Недостатки GMRES и перезапущенного GMRES устраняются повторным использованием подпространства Крылова в методах типа GCRO, таких как GCROT и GCRODR. ^[5] Повторное использование подпространств Крылова в GMRES также может ускорить сходимость, когда необходимо решить последовательности линейных систем. ^[6]

Сравнение с другими решателями [ править ]

Итерация Арнольди сводится к итерации Ланцоша для симметричных матриц. Соответствующий метод подпространств Крылова является методом минимальных невязок (MinRes) Пейдж и Сондерс. В отличие от несимметричного случая, метод MinRes задается трехчленным рекуррентным соотношением . Можно показать, что не существует метода подпространств Крылова для общих матриц, который задается коротким рекуррентным соотношением и все же минимизирует нормы остатков, как это делает GMRES.

Другой класс методов основан на несимметричной итерации Ланцоша , в частности метод BiCG . Они используют трехчленное рекуррентное соотношение, но они не достигают минимальной невязки, и, следовательно, невязка не уменьшается монотонно для этих методов. Сходимость даже не гарантируется.

Третий класс образован такими методами, как CGS и BiCGSTAB . Они также работают с трехчленным рекуррентным отношением (следовательно, без оптимальности) и могут даже преждевременно завершиться без достижения сходимости. Идея этих методов заключается в подходящем выборе порождающих полиномов итерационной последовательности.

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

Решение задачи наименьших квадратов [ править ]

Одна часть метода GMRES - найти вектор, который минимизирует${\displaystyle y_{n}}$

${\displaystyle \|{\tilde {H}}_{n}y_{n}-\beta e_{1}\|.\,}$

Обратите внимание, что это матрица  размером ( n + 1) на n , поэтому она дает чрезмерно ограниченную линейную систему из n +1 уравнений для n неизвестных.${\displaystyle {\tilde {H}}_{n}}$

Минимум может быть вычислен с помощью QR-разложения : найти ( n  + 1) -by- ( n  + 1) ортогональную матрицу Ω _n и ( n  + 1) -by- n матрицу верхнего треугольника такие, что${\displaystyle {\tilde {R}}_{n}}$

${\displaystyle \Omega _{n}{\tilde {H}}_{n}={\tilde {R}}_{n}.}$

Треугольная матрица имеет на одну строку больше, чем столбцов, поэтому ее нижняя строка состоит из нуля. Следовательно, его можно разложить как

${\displaystyle {\tilde {R}}_{n}={\begin{bmatrix}R_{n}\\0\end{bmatrix}},}$

где это н матрицу с размерностью п (таким образом , квадрат) треугольная матрица.${\displaystyle R_{n}}$

QR-разложение можно дешево обновлять от одной итерации к следующей, потому что матрицы Хессенберга отличаются только строкой нулей и столбцом:

${\displaystyle {\tilde {H}}_{n+1}={\begin{bmatrix}{\tilde {H}}_{n}&h_{n+1}\\0&h_{n+2,n+1}\end{bmatrix}},}$

где ч _{п + 1} = ( ч _{1, п + 1} , ..., ч _{п + 1, п + 1} ) ^Т . Это означает, что предварительное умножение матрицы Хессенберга на Ω _n , дополненное нулями и строкой с мультипликативной единицей, дает почти треугольную матрицу:

${\displaystyle {\begin{bmatrix}\Omega _{n}&0\\0&1\end{bmatrix}}{\tilde {H}}_{n+1}={\begin{bmatrix}R_{n}&r_{n+1}\\0&\rho \\0&\sigma \end{bmatrix}}}$

Это было бы треугольником, если σ равно нулю. Чтобы исправить это, нужно вращение Гивенса

${\displaystyle G_{n}={\begin{bmatrix}I_{n}&0&0\\0&c_{n}&s_{n}\\0&-s_{n}&c_{n}\end{bmatrix}}}$

куда

${\displaystyle c_{n}={\frac {\rho }{\sqrt {\rho ^{2}+\sigma ^{2}}}}\quad {\mbox{and}}\quad s_{n}={\frac {\sigma }{\sqrt {\rho ^{2}+\sigma ^{2}}}}.}$

С помощью этого вращения Гивенса мы формируем

${\displaystyle \Omega _{n+1}=G_{n}{\begin{bmatrix}\Omega _{n}&0\\0&1\end{bmatrix}}.}$

В самом деле,

${\displaystyle \Omega _{n+1}{\tilde {H}}_{n+1}={\begin{bmatrix}R_{n}&r_{n+1}\\0&r_{n+1,n+1}\\0&0\end{bmatrix}}\quad {\text{with}}\quad r_{n+1,n+1}={\sqrt {\rho ^{2}+\sigma ^{2}}}}$

- треугольная матрица.

Учитывая разложение QR, проблема минимизации легко решается, если заметить, что

${\displaystyle \|{\tilde {H}}_{n}y_{n}-\beta e_{1}\|=\|\Omega _{n}({\tilde {H}}_{n}y_{n}-\beta e_{1})\|=\|{\tilde {R}}_{n}y_{n}-\beta \Omega _{n}e_{1}\|.}$

Обозначая вектор через${\displaystyle \beta \Omega _{n}e_{1}}$

${\displaystyle {\tilde {g}}_{n}={\begin{bmatrix}g_{n}\\\gamma _{n}\end{bmatrix}}}$

при g _n ∈ R ⁿ и γ _n ∈ R это

${\displaystyle \|{\tilde {H}}_{n}y_{n}-\beta e_{1}\|=\|{\tilde {R}}_{n}y_{n}-\beta \Omega _{n}e_{1}\|=\left\|{\begin{bmatrix}R_{n}\\0\end{bmatrix}}y_{n}-{\begin{bmatrix}g_{n}\\\gamma _{n}\end{bmatrix}}\right\|.}$

Вектор y, который минимизирует это выражение, задается следующим образом:

${\displaystyle y_{n}=R_{n}^{-1}g_{n}.}$

Опять же, векторы легко обновить. ^[7]${\displaystyle g_{n}}$

Пример кода [ править ]

Обычный GMRES (MATLAB / GNU Octave) [ править ]

function  [x, e] = gmres ( A, b, x, max_iterations, порог )  n  =  длина ( A );  m  =  max_iterations ; % использовать x как начальный вектор  r  =  b  -  A  *  x ; b_norm  =  норма ( б );  ошибка  =  норма ( r )  /  b_norm ; % инициализировать одномерные векторы  sn  =  zeros ( m ,  1 );  cs  =  нули ( m ,  1 );  % e1 = нули (n, 1);  e1  =  нули ( m + 1 ,  1 );  е1 ( 1 )  =  1 ;  е  =  [ ошибка ];  r_norm  =  норма ( г );  Q (:, 1 )  =  r  /  r_norm ; бета  =  r_norm  *  e1 ;  % Примечание: это не бета-скаляр из раздела «Метод» выше, а бета-скаляр, умноженный на e1  для  k  =  1 : m % run arnoldi  [ H ( 1 : k + 1 ,  k )  Q (:,  k + 1 )]  =  arnoldi ( A ,  Q ,  k );  % удалить последний элемент в строке H i и обновить матрицу вращения  [ H ( 1 : k + 1 ,  k )  cs ( k )  sn ( k )]  =  apply_givens_rotation ( H ( 1 : k + 1 , k ),  cs ,  сн ,  к );  % обновить остаточный вектор  beta ( k  +  1 )  =  - sn ( k )  *  beta ( k );  бета ( к )  =  сс ( к )  *  бета ( к );  error  =  abs ( beta ( k  +  1 ))  /  b_norm ; % сохранить ошибку  e  =  [ e ;  ошибка ]; если  ( ошибка  <=  порог )  перерыв ;  end  end  %, если порог не достигнут, k = m в этой точке (а не m + 1)  % вычислить результат  y  =  H ( 1 : k ,  1 : k )  \  beta ( 1 : k );  х  =  х  +  Q (:,  1 : к )  *  у ; конец% ------------------------------------------------- --- %% Функция Арнольди%% ----------------------------------------- -----------% function  [h, q] = arnoldi ( A, Q, k )  q  =  A * Q (:, k );  % Вектор Крылова  для  i  =  1 : k  % Модифицированный по Граму-Шмидту, сохраняя матрицу Хессенберга  h ( i )  =  q '  *  Q (:,  i );  д  =  д  -  ч ( я )  *  Q (:,  я );  конец  h (k  +  1 )  =  norm ( q );  q  =  q  / h ( k  +  1 ); конец% ------------------------------------------------- --------------------%% Применение вращения Гивенса к H col%% -------------------- -------------------------------------------------% функция  [h, cs_k, sn_k] = apply_givens_rotation ( h, cs, sn, k )  % применить для i-го столбца  для  i  =  1 : k - 1  temp  =  cs ( i )  *  h ( i )  +  sn ( i )  *  h ( i  +  1 );  h ( i + 1 )  =  - sn ( i )  *  h ( i )  +  cs ( i )  *  h ( i  + 1 );  h ( i )  =  темп ;  конец %  Обновление  следующий грех сов значений для вращения [ cs_k sn_k ] = givens_rotation ( ч ( к ), ч ( к + 1 ));              % удалить H (i + 1, i)  h ( k )  =  cs_k  *  h ( k )  +  sn_k  *  h ( k  +  1 );  h ( k  +  1 )  =  0,0 ; конец%% ---- Вычислить  Учитывая вращение матрицы ---- %% функции [CS, Sn] = givens_rotation ( v1, v2 )     % если (v1 == 0) % cs = 0; % sn = 1; % else  t  =  sqrt ( v1 ^ 2  +  v2 ^ 2 ); % cs = abs (v1) / t; % sn = cs * v2 / v1;  cs  =  v1  /  t ;  % см. http://www.netlib.org/eispack/comqr.f  sn  =  v2  /  t ; % конец конец

См. Также [ править ]

• Метод двусопряженного градиента

Ссылки [ править ]

Примечания [ править ]

• А. Мейстер, Numerik linearer Gleichungssysteme , 2-е издание, Vieweg 2005, ISBN 978-3-528-13135-7 . 
• Ю. Саад, Итерационные методы для разреженных линейных систем , 2-е издание, Общество промышленной и прикладной математики , 2003. ISBN 978-0-89871-534-7 . 
• Ю. Саад и М. Х. Шульц, "GMRES: обобщенный алгоритм минимальной невязки для решения несимметричных линейных систем", SIAM J. Sci. Стат. Comput. , 7 : 856-869, 1986. DOI : 10,1137 / 0907058 .
• SC Eisenstat, HC Elman и MH Schultz, "Вариационные итерационные методы для несимметричных систем линейных уравнений", журнал SIAM по численному анализу , 20 (2), 345–357, 1983.
• Дж. Стоер и Р. Булирш, Введение в численный анализ , 3-е издание, Спрингер, Нью-Йорк, 2002. ISBN 978-0-387-95452-3 . 
• Ллойд Н. Трефетен и Дэвид Бау, III, Численная линейная алгебра , Общество промышленной и прикладной математики, 1997. ISBN 978-0-89871-361-9 . 
• Dongarra et al. , Шаблоны для решения линейных систем: строительные блоки для итерационных методов , 2-е издание, SIAM, Филадельфия, 1994
• Амриткар, Амит; де Стерлер, Эрик; Свиридович, Катажина; Тафти, Данеш; Ахуджа, Капил (2015). «Переработка подпространств Крылова для приложений CFD и новый гибридный решатель рециклинга». Журнал вычислительной физики 303: 222. DOI: 10.1016 / j.jcp.2015.09.040