Обобщенный метод минимальной невязки - Generalized minimal residual method

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

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

Метод

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

${displaystyle Ax=b.,}$

Матрица А предполагается обратимый размера м-от-м. Кроме того, предполагается, что б нормализовано, т. е. что ${displaystyle |b|=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}eq 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}}$.

Следовательно, вектор ${displaystyle x_{n}in 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 n+1}$)-от-${displaystyle n}$ верхний Матрица Гессенберга ${displaystyle { ilde {H}}_{n}}$ с участием

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

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

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

где

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

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

${displaystyle eta =|b-Ax_{0}|,,}$

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

${displaystyle r_{n}={ ilde {H}}_{n}y_{n}-eta 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. повторить, если остаток еще недостаточно мал.

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

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

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

Так не бывает вообще. Действительно, теорема Гринбаума, Птака и Стракоша утверждает, что для любой невозрастающей последовательности а₁, …, а_м−1, а_м = 0, можно найти матрицу А такой, что ||р_п|| = а_п для всех п, где р_п - невязка, определенная выше. В частности, можно найти матрицу, для которой невязка остается постоянной при м - 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)}}ight)^{n/2}|r_{0}|,}$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

${displaystyle |{ ilde {H}}_{n}y_{n}-eta e_{1}|.,}$

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

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

${displaystyle Omega _{n}{ ilde {H}}_{n}={ ilde {R}}_{n}.}$

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

${displaystyle { ilde {R}}_{n}={egin{bmatrix}R_{n}�end{bmatrix}},}$

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

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

${displaystyle { ilde {H}}_{n+1}={egin{bmatrix}{ ilde {H}}_{n}&h_{n+1}�&h_{n+2,n+1}end{bmatrix}},}$

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

${displaystyle {egin{bmatrix}Omega _{n}&0�&1end{bmatrix}}{ ilde {H}}_{n+1}={egin{bmatrix}R_{n}&r_{n+1}�&ho �&sigma end{bmatrix}}}$

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

${displaystyle G_{n}={egin{bmatrix}I_{n}&0&0�&c_{n}&s_{n}�&-s_{n}&c_{n}end{bmatrix}}}$

где

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

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

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

Действительно,

${displaystyle Omega _{n+1}{ ilde {H}}_{n+1}={egin{bmatrix}R_{n}&r_{n+1}�&r_{n+1,n+1}�&0end{bmatrix}}quad { ext{with}}quad r_{n+1,n+1}={sqrt {ho ^{2}+sigma ^{2}}}}$

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

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

${displaystyle |{ ilde {H}}_{n}y_{n}-eta e_{1}|=|Omega _{n}({ ilde {H}}_{n}y_{n}-eta e_{1})|=|{ ilde {R}}_{n}y_{n}-eta Omega _{n}e_{1}|.}$

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

${displaystyle { ilde {g}}_{n}={egin{bmatrix}g_{n}gamma _{n}end{bmatrix}}}$

с участием г_п ∈ р^п и γ_п ∈ р, это

${displaystyle |{ ilde {H}}_{n}y_{n}-eta e_{1}|=|{ ilde {R}}_{n}y_{n}-eta Omega _{n}e_{1}|=left|{egin{bmatrix}R_{n}�end{bmatrix}}y_{n}-{egin{bmatrix}g_{n}gamma _{n}end{bmatrix}}ight|.}$

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

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

И снова векторы ${displaystyle g_{n}}$ легко обновляются.^[7]

Пример кода

Обычный GMRES (MATLAB / GNU Octave)

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

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

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

использованная литература

1. Ю. Саад и М.Х. Шульц
2. Пейдж и Сондерс, "Решение разреженных неопределенных систем линейных уравнений", SIAM J. Numer. Anal., Том 12, стр. 617 (1975) https://doi.org/10.1137/0712047
3. Eisenstat, Elman & Schultz, Thm 3.3. NB: все результаты для GCR справедливы и для GMRES, см. Саад и Шульц
4. Trefethen & Bau, Thm 35.2
5. Амриткар, Амит; де Стерлер, Эрик; Свиридович, Катажина; Тафти, Данеш; Ахуджа, Капил (2015). «Переработка подпространств Крылова для приложений CFD и новый гибридный решатель рециклинга». Журнал вычислительной физики. 303: 222. arXiv:1501.03358. Bibcode:2015JCoPh.303..222A. Дои:10.1016 / j.jcp.2015.09.040.
6. Галлия, Андре (2014). Переработка методов подпространства Крылова для последовательностей линейных систем (Кандидат наук.). ТУ Берлин. Дои:10.14279 / depositonce-4147.
7. Стоер и Булирш, §8.7.2

Заметки

• А. Мейстер, 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. Дои:10.1137/0907058.
• С.К. Эйзенстат, Х.К. Эльман и М. Шульц, "Вариационные итерационные методы для несимметричных систем линейных уравнений", Журнал 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