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

Метод численного решения некоторых систем уравнений

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

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

Метод

Обозначим евклидову норму любого вектора v через . Обозначим (квадратную) систему линейных уравнений, которую нужно решить, через Матрица A предполагается обратимой размером m на m . Кроме того, предполагается, что b нормализовано , т.е. ${\displaystyle \|v\|}$ $${\displaystyle Ax=b.}$$ ${\displaystyle \|b\|=1}$

n -ое подпространство Крылова для этой задачи равно , где — начальная ошибка при заданном начальном предположении . Очевидно, что если . $${\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 x_{0}+K_{n}}$ ${\displaystyle r_{n}=b-Ax_{n}}$

Векторы могут быть близки к линейно зависимым , поэтому вместо этого базиса используется итерация Арнольди для нахождения ортонормальных векторов , которые образуют базис для . В частности, . ${\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, образованная . Другими словами, нахождение n -го приближения решения (т.е. ) сводится к нахождению вектора , который определяется посредством минимизации остатка, как описано ниже. ${\displaystyle x_{n}\in x_{0}+K_{n}}$ ${\displaystyle x_{n}=x_{0}+Q_{n}y_{n}}$ ${\displaystyle y_{n}\in \mathbb {R} ^{n}}$ ${\displaystyle Q_{n}}$ ${\displaystyle q_{1},\ldots ,q_{n}}$ ${\displaystyle x_{n}}$ ${\displaystyle y_{n}}$

Процесс Арнольди также создает , матрицу Хессенберга ( ) , которая удовлетворяет равенству, которое используется для упрощения вычисления (см. § Решение задачи наименьших квадратов). Обратите внимание, что для симметричных матриц фактически достигается симметричная трехдиагональная матрица, что приводит к методу MINRES . ${\displaystyle {\tilde {H}}_{n}}$ ${\displaystyle n+1}$ ${\displaystyle n}$ $${\displaystyle AQ_{n}=Q_{n+1}{\tilde {H}}_{n}\,}$$ ${\displaystyle y_{n}}$

Поскольку столбцы ортонормальны, то имеем , где — первый вектор в стандартном базисе , а — первый пробный вектор остатка (обычно ). Следовательно, можно найти, минимизируя евклидову норму остатка Это линейная задача наименьших квадратов размера n . ${\displaystyle Q_{n}}$ $${\displaystyle {\begin{aligned}\left\|r_{n}\right\|&=\left\|b-Ax_{n}\right\|\\&=\left\|b-A(x_{0}+Q_{n}y_{n})\right\|\\&=\left\|r_{0}-AQ_{n}y_{n}\right\|\\&=\left\|\beta q_{1}-AQ_{n}y_{n}\right\|\\&=\left\|\beta q_{1}-Q_{n+1}{\tilde {H}}_{n}y_{n}\right\|\\&=\left\|Q_{n+1}(\beta e_{1}-{\tilde {H}}_{n}y_{n})\right\|\\&=\left\|\beta e_{1}-{\tilde {H}}_{n}y_{n}\right\|\end{aligned}}}$$ $${\displaystyle e_{1}=(1,0,0,\ldots ,0)^{T}\,}$$ ${\displaystyle \mathbb {R} ^{n+1}}$ $${\displaystyle \beta =\|r_{0}\|\,,}$$ ${\displaystyle r_{0}}$ ${\displaystyle b}$ ${\displaystyle x_{n}}$ $${\displaystyle r_{n}={\tilde {H}}_{n}y_{n}-\beta e_{1}.}$$

Это дает метод 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)}$

Конвергенция

Итерация n минимизирует остаток в подпространстве Крылова . Поскольку каждое подпространство содержится в следующем подпространстве, остаток не увеличивается. После m итераций, где m — размер матрицы A , пространство Крылова K _m равно всему R ^m и, следовательно, метод GMRES приходит к точному решению. Однако идея заключается в том, что после небольшого числа итераций (относительно m ) вектор x _n уже является хорошим приближением к точному решению. ${\displaystyle K_{n}}$

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

На практике, однако, GMRES часто работает хорошо. Это можно доказать в конкретных ситуациях. Если симметричная часть A , то есть , положительно определена , то где и обозначают наименьшее и наибольшее собственное значение матрицы , соответственно. ^[4] ${\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}\|,}$$ ${\displaystyle \lambda _{\mathrm {min} }(M)}$ ${\displaystyle \lambda _{\mathrm {max} }(M)}$ ${\displaystyle M}$

Если A симметричен и положительно определен, то мы даже имеем , где обозначает число обусловленности A в евклидовой норме. $${\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 не является положительно определенным, мы имеем где P _n обозначает множество полиномов степени не выше n с p (0) = 1, V — матрица, появляющаяся в спектральном разложении A , а σ ( A ) — спектр A. Грубо говоря, это говорит о том, что быстрая сходимость происходит, когда собственные значения A сгруппированы вдали от начала координат и A не слишком далеко от нормальности . ^[5] $${\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 )|,\,}$$

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

Расширения метода

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

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

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

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

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

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

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

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

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

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

Минимум можно вычислить с помощью QR-разложения : найти ортогональную матрицу Ω _n размером ( n  + 1) на ( n  + 1) и верхнюю треугольную матрицу ( n  + 1) на n, такую, что Треугольная матрица имеет на одну строку больше, чем столбцов, поэтому ее нижняя строка состоит из нуля. Следовательно, ее можно разложить следующим образом, где — треугольная матрица размером n на 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-разложение можно легко обновить от одной итерации к другой, поскольку матрицы Хессенберга отличаются только строкой нулей и столбцом: где h _n+1 = ( h _{1, n +1} , ..., h _{n +1, n +1} ) ^T . Это означает, что предварительное умножение матрицы Хессенберга на Ω _n , дополненное нулями и строкой с мультипликативным тождеством, дает почти треугольную матрицу: Она была бы треугольной, если бы σ было равно нулю. Чтобы исправить это, нужно вращение Гивенса , где С этим вращением Гивенса мы формируем Действительно, является треугольной матрицей с . $${\displaystyle {\tilde {H}}_{n+1}={\begin{bmatrix}{\tilde {H}}_{n}&h_{n+1}\\0&h_{n+2,n+1}\end{bmatrix}},}$$ $${\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 {\text{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}}}$$ ${\textstyle r_{n+1,n+1}={\sqrt {\rho ^{2}+\sigma ^{2}}}}$

Учитывая QR-разложение, задача минимизации легко решается, если заметить, что Обозначим вектор через с g _n ∈ R ⁿ и γ _n ∈ R , тогда Вектор y , который минимизирует это выражение, задается как Опять же, векторы легко обновлять. ^[8] $${\displaystyle {\begin{aligned}\left\|{\tilde {H}}_{n}y_{n}-\beta e_{1}\right\|&=\left\|\Omega _{n}({\tilde {H}}_{n}y_{n}-\beta e_{1})\right\|\\&=\left\|{\tilde {R}}_{n}y_{n}-\beta \Omega _{n}e_{1}\right\|.\end{aligned}}}$$ ${\displaystyle \beta \Omega _{n}e_{1}}$ $${\displaystyle {\tilde {g}}_{n}={\begin{bmatrix}g_{n}\\\gamma _{n}\end{bmatrix}}}$$ $${\displaystyle {\begin{aligned}\left\|{\tilde {H}}_{n}y_{n}-\beta e_{1}\right\|&=\left\|{\tilde {R}}_{n}y_{n}-\beta \Omega _{n}e_{1}\right\|\\&=\left\|{\begin{bmatrix}R_{n}\\0\end{bmatrix}}y_{n}-{\begin{bmatrix}g_{n}\\\gamma _{n}\end{bmatrix}}\right\|.\end{aligned}}}$$ $${\displaystyle y_{n}=R_{n}^{-1}g_{n}.}$$ ${\displaystyle g_{n}}$

Пример кода

Регулярный GMRES (MATLAB / GNU Octave)

функция [x, e] = gmres ( A, b, x, max_iterations, threshold ) n = длина ( A ); m = max_iterations ;        % использовать x как начальный вектор r = b - A * x ;        b_норма = норма ( b ); ошибка = норма ( r ) / b_норма ;        % инициализируем одномерные векторы sn = zeros ( m , 1 ); cs = zeros ( m , 1 ); %e1 = zeros(n, 1); e1 = zeros ( m + 1 , 1 ); e1 ( 1 ) = 1 ; e = [ error ]; r_norm = norm ( r ); Q (:, 1 ) = r / r_norm ; % Примечание: это не бета-скаляр в разделе "Метод" выше, а % бета-скаляр, умноженный на e1 beta = r_norm * e1 ; для k = 1 : m                                       % запустить arnoldi [ H ( 1 : k + 1 , k ), Q (:, k + 1 )] = arnoldi ( A , Q , k ); % исключить последний элемент в H -й строке и обновить матрицу вращения [ H ( 1 : k + 1 , k ), cs ( k ), sn ( k )] = apply_givens_rotation ( H ( 1 : k + 1 , k ), cs , sn , k ); % обновить вектор остатков beta ( k + 1 ) = - sn ( k ) * beta ( k ); beta ( k ) = cs ( k ) * beta ( k ); error = abs ( beta ( k + 1 )) / b_norm ;                                         % сохранить ошибку e = [ e ; error ];     если ( ошибка <= порог ) break ; конец конец % если порог не достигнут, k = m в этой точке (а не m+1) % вычислить результат y = H ( 1 : k , 1 : k ) \ beta ( 1 : k ); x = x + Q (:, 1 : k ) * y ; конец                       %----------------------------------------------------% % Функция Арнольди % %----------------------------------------------------% function [h, q] = arnoldi ( A, Q, k ) q = A * Q (:, k ); % Вектор Крылова для i = 1 : k % Модифицированный Грам-Шмидт, сохраняющий матрицу Хессенберга h ( i ) = q ' * Q (:, i ); q = q - h ( i ) * Q (:, i ); end h ( k + 1 ) = norm ( q ); q = q / h ( k + 1 ); end                                     %---------------------------------------------------------------------% % Применение вращения Гивенса к столбцу H % %---------------------------------------------------------------------% function [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 ) = temp ; конец                                 % обновить следующие значения sin cos для вращения [ cs_k , sn_k ] = givens_rotation ( h ( k ), h ( k + 1 ));        % исключить H(i + 1, i) h ( k ) = cs_k * h ( k ) + sn_k * h ( k + 1 ); h ( k + 1 ) = 0.0 ; конец                %%----Вычислить матрицу вращения Гивенса----%% function [cs, sn] = givens_rotation ( v1, v2 ) % if (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 ; % end end                 

Смотрите также

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

Ссылки

1. Саад, Юсеф; Шульц, Мартин Х. (1986). «GMRES: обобщенный алгоритм минимального остатка для решения несимметричных линейных систем». Журнал SIAM по научным и статистическим вычислениям . 7 (3): 856–869. doi :10.1137/0907058. ISSN  0196-5204.
2. Пейдж и Сондерс, «Решение разреженных неопределенных систем линейных уравнений», SIAM J. Numer. Anal., т. 12, стр. 617 (1975) https://doi.org/10.1137/0712047
3. Нифа, Науфал (2017). Solveurs Performanceants pour l'optimisation sous contraintes en идентифицирование параметров [ Эффективные решатели для ограниченной оптимизации в задачах идентификации параметров ] (Диссертация) (на французском языке).
4. Eisenstat, Elman & Schultz 1983, Thm 3.3. Примечание: все результаты для GCR справедливы также для GMRES, см. Saad & Schultz 1986
5. Трефетен, Ллойд Н.; Бау, Дэвид, III. (1997). Численная линейная алгебра . Филадельфия: Общество промышленной и прикладной математики. Теорема 35.2. ISBN 978-0-89871-361-9.{{cite book}}: CS1 maint: multiple names: authors list (link)
6. Амриткар, Амит; де Стурлер, Эрик; Свиридович, Катажина; Тафти, Данеш; Ахуджа, Капил (2015). «Переработка подпространств Крылова для приложений вычислительной гидродинамики и новый гибридный решатель задач переработки». Журнал вычислительной физики . 303 : 222. arXiv : 1501.03358 . Bibcode :2015JCoPh.303..222A. doi :10.1016/j.jcp.2015.09.040. S2CID  2933274.
7. Гауль, Андре (2014). Повторное использование методов подпространства Крылова для последовательностей линейных систем (Ph.D.). Технический университет Берлина. doi :10.14279/depositonce-4147.
8. Стоер, Йозеф; Булирш, Роланд (2002). Введение в численный анализ . Тексты по прикладной математике (3-е изд.). Нью-Йорк: Springer. §8.7.2. ISBN 978-0-387-95452-3.

• Мейстер, Андреас; Фёмель, Кристоф (2005). Числовая линейная система Gleichungssysteme . Висбаден: Просмотрег. ISBN 978-3-528-13135-7.
• Саад, И. (2003). Итерационные методы для разреженных линейных систем (2-е изд.). Филадельфия: SIAM. ISBN 978-0-89871-534-7.
• Эйзенштат, Стэнли К.; Элман, Говард К.; Шульц, Мартин Х. (1983). «Вариационные итерационные методы для несимметричных систем линейных уравнений». Журнал SIAM по численному анализу . 20 (2): 345–357. doi :10.1137/0720023. ISSN  0036-1429.
• Донгарра и др., Шаблоны для решения линейных систем: строительные блоки для итерационных методов, 2-е издание, SIAM, Филадельфия, 1994 г.
• Иманкулов, Тимур; Лебедев, Данил; Маткерим, Базаргуль; Дарибаев, Беймбет; Касымбек, Нурислам (2021-10-08). "Численное моделирование многофазного многокомпонентного течения в пористых средах: анализ эффективности метода Ньютона". Fluids . 6 (10): 355. doi : 10.3390/fluids6100355 . ISSN  2311-5521.