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

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

Метод GMRES был разработан Юсефом Саадом и Мартином Х. Шульцем в 1986 году. Это обобщение и улучшение метода MINRES , предложенного Пейдж и Сондерсом в 1975 году. Метод MINRES требует, чтобы матрица была симметричной, но имеет то преимущество, что требует обработки только трех векторов. GMRES — это частный случай метода DIIS, разработанного Питером Пулеем в 1980 году. DIIS применим к нелинейным системам.

Метод

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

n -е для подпространство Крылова этой задачи есть $$K_{n}=K_{n}(A,r_{0})=\operatorname {span} \,\{r_{0},Ar_{0},A^{2}r_{0},\ldots ,A^{n-1}r_{0}\}.\,$$ где $r_{0}=b-Ax_{0}$ это первоначальная ошибка с учетом первоначального предположения $x_{0}\neq 0$. Четко $r_{0}=b$ если $x_{0}=0$.

GMRES приближает точное решение $Ax=b$ по вектору $x_{n}\in x_{0}+K_{n}$ что минимизирует евклидову норму невязки $r_{n}=b-Ax_{n}$.

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

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

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

Поскольку столбцы $Q_{n}$ ортонормированы, мы имеем $${\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}}$$ где $$e_{1}=(1,0,0,\ldots ,0)^{T}\,$$ - первый вектор в стандартном базисе $\mathbb {R} ^{n+1}$, и $$\beta =\|r_{0}\|\,,$$ $r_{0}$ являющийся первым пробным вектором невязки (обычно $b$). Следовательно, $x_{n}$ можно найти, минимизируя евклидову норму невязки $$r_{n}={\tilde {H}}_{n}y_{n}-\beta e_{1}.$$ Это линейная задача наименьших квадратов размера n .

Это дает метод GMRES. На $n$-я итерация:

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

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

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

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

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

Однако на практике GMRES часто работает хорошо. Это можно доказать в конкретных ситуациях. Если симметричная часть A , то есть $(A^{T}+A)/2$, положительно определена , тогда $$\|r_{n}\|\leq \left(1-{\frac {\lambda _{\min }^{2}(1/2(A^{T}+A))}{\lambda _{\max }(A^{T}A)}}\right)^{n/2}\|r_{0}\|,$$ где $\lambda _{\mathrm {min} }(M)$ и $\lambda _{\mathrm {max} }(M)$ обозначаем наименьшее и наибольшее собственное значение матрицы $M$, соответственно.

Если A симметричен имеем и положительно определен, то мы даже $$\|r_{n}\|\leq \left({\frac {\kappa _{2}(A)^{2}-1}{\kappa _{2}(A)^{2}}}\right)^{n/2}\|r_{0}\|.$$ где $\kappa _{2}(A)$ обозначает число обусловленности A в евклидовой норме.

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

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

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

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

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

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

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

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

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

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

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

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

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

Минимум можно вычислить с помощью QR-разложения : найдите ( n + 1) на ( n + 1) ортогональную матрицу Ω _n и ( n + 1) на n . верхнюю треугольную матрицу ${\tilde {R}}_{n}$ такой, что $$\Omega _{n}{\tilde {H}}_{n}={\tilde {R}}_{n}.$$ Треугольная матрица имеет на одну строку больше, чем столбцов, поэтому ее нижняя строка состоит из нуля. Следовательно, его можно разложить как $${\tilde {R}}_{n}={\begin{bmatrix}R_{n}\\0\end{bmatrix}},$$ где $R_{n}$ представляет собой треугольную матрицу размера n x n (таким образом, квадратную).

QR-разложение можно дешево обновлять от одной итерации к другой, поскольку матрицы Хессенберга отличаются только строкой нулей и столбцом: $${\tilde {H}}_{n+1}={\begin{bmatrix}{\tilde {H}}_{n}&h_{n+1}\\0&h_{n+2,n+1}\end{bmatrix}},$$ где час _n+1 = ( час _{1, n +1} , ..., час _{n +1, n +1} ) ^Т . Это означает, что предварительное умножение матрицы Хессенберга на Ω _n , дополненную нулями и строкой с мультипликативной единицей, дает почти треугольную матрицу: $${\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}}$$ Это было бы треугольным, если σ равно нулю. Чтобы исправить это, необходимо вращение Гивенса. $$G_{n}={\begin{bmatrix}I_{n}&0&0\\0&c_{n}&s_{n}\\0&-s_{n}&c_{n}\end{bmatrix}}$$ где $$c_{n}={\frac {\rho }{\sqrt {\rho ^{2}+\sigma ^{2}}}}\quad {\text{and}}\quad s_{n}={\frac {\sigma }{\sqrt {\rho ^{2}+\sigma ^{2}}}}.$$ С помощью этого вращения Гивенса мы формируем $$\Omega _{n+1}=G_{n}{\begin{bmatrix}\Omega _{n}&0\\0&1\end{bmatrix}}.$$ Действительно, $$\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-разложение, проблему минимизации легко решить, заметив, что $${\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}}$$ Обозначая вектор $\beta \Omega _{n}e_{1}$ к $${\tilde {g}}_{n}={\begin{bmatrix}g_{n}\\\gamma _{n}\end{bmatrix}}$$ где g _n ∈ R ^н и γn _∈ R , это $${\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}}$$ Вектор y , который минимизирует это выражение, определяется выражением $$y_{n}=R_{n}^{-1}g_{n}.$$ И снова векторы $g_{n}$ легко обновляются.

Пример кода

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

function [x, e] = gmres(A, b, x, max_iterations, threshold)
  n = length(A);
  m = max_iterations;

  % use x as the initial vector
  r = b - A * x;

  b_norm = norm(b);
  error = norm(r) / b_norm;

  % initialize the 1D vectors
  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;
  % Note: this is not the beta scalar in section "The method" above but
  % the beta scalar multiplied by e1
  beta = r_norm * e1;
  for k = 1:m

    % run arnoldi
    [H(1:k+1, k), Q(:, k+1)] = arnoldi(A, Q, k);
    
    % eliminate the last element in H ith row and update the rotation matrix
    [H(1:k+1, k), cs(k), sn(k)] = apply_givens_rotation(H(1:k+1,k), cs, sn, k);
    
    % update the residual vector
    beta(k + 1) = -sn(k) * beta(k);
    beta(k)     = cs(k) * beta(k);
    error       = abs(beta(k + 1)) / b_norm;

    % save the error
    e = [e; error];

    if (error <= threshold)
      break;
    end
  end
  % if threshold is not reached, k = m at this point (and not m+1) 
  
  % calculate the result
  y = H(1:k, 1:k) \ beta(1:k);
  x = x + Q(:, 1:k) * y;
end

%----------------------------------------------------%
%                  Arnoldi Function                  %
%----------------------------------------------------%
function [h, q] = arnoldi(A, Q, k)
  q = A*Q(:,k);   % Krylov Vector
  for i = 1:k     % Modified Gram-Schmidt, keeping the Hessenberg matrix
    h(i) = q' * Q(:, i);
    q = q - h(i) * Q(:, i);
  end
  h(k + 1) = norm(q);
  q = q / h(k + 1);
end

%---------------------------------------------------------------------%
%                  Applying Givens Rotation to H col                  %
%---------------------------------------------------------------------%
function [h, cs_k, sn_k] = apply_givens_rotation(h, cs, sn, k)
  % apply for ith column
  for 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;
  end

  % update the next sin cos values for rotation
  [cs_k, sn_k] = givens_rotation(h(k), h(k + 1));

  % eliminate H(i + 1, i)
  h(k) = cs_k * h(k) + sn_k * h(k + 1);
  h(k + 1) = 0.0;
end

%%----Calculate the Givens rotation matrix----%%
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;  % see http://www.netlib.org/eispack/comqr.f
    sn = v2 / t;
%  end
end

См. также

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

Ссылки

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