Uogólniona metoda minimalnych pozostałości - Generalized minimal residual method

Matematyki The uogólniona metoda minimalna resztkowa (GMRES) jest iteracyjna metoda do numerycznego roztworu nieokreślony niesymetrycznej układu równań liniowych . Metoda aproksymuje rozwiązanie przez wektor w podprzestrzeni Kryłowa z minimalną resztą . Arnoldi iteracja jest używany, aby znaleźć ten wektor.

Metoda GMRES została opracowana przez Yousefa Saada i Martina H. Schultza w 1986 roku. Jest uogólnieniem i udoskonaleniem metody MINRES dzięki Paige i Saunders w 1975 roku. Metoda MINRES wymaga, aby macierz była symetryczna, ale ma tę zaletę, że wymaga jedynie obsługi trzech wektorów. GMRES jest szczególnym przypadkiem metody DIIS opracowanej przez Petera Pulaya w 1980 roku. DIIS ma zastosowanie do systemów nieliniowych.

Metoda

Oznacz normę euklidesową dowolnego wektora v przez . Oznacz (kwadratowy) układ równań liniowych do rozwiązania przez ${\displaystyle \|v\|}$

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

Zakłada się, że macierz A jest odwracalna o rozmiarze m -by- m . Ponadto zakłada się, że b jest znormalizowane, to znaczy, że . ${\ Displaystyle \ | b \ | = 1}$

W n -tej Przestrzeń Kryłowa tego problemu jest

${\ Displaystyle K_ {n} = K_ {n} (A, r_ {0}) = \ operatorname {rozpiętość} \ \ {r_ {0}, Ar_ {0}, A ^ {2} r_ {0}, \ldots ,A^{n-1}r_{0}\}.\,}$

gdzie jest początkowym błędem przy początkowym zgadywaniu . Oczywiście, jeśli . ${\displaystyle r_{0}=b-Ax_{0}}$${\displaystyle x_{0}\neq 0}$${\displaystyle r_{0}=b}$${\displaystyle x_{0}=0}$

GMRES przybliża dokładne rozwiązanie przez wektor, który minimalizuje normę euklidesową reszty . ${\ Displaystyle Axe = b}$${\ Displaystyle x_ {n} \ w K_ {n}}$ ${\ Displaystyle r_ {n} = b-Ax_ {n}}$

Wektory mogą być bliskie liniowo zależne , więc zamiast tej bazy, iteracja Arnoldiego jest używana do znalezienia wektorów ortonormalnych, które tworzą bazę dla . W szczególności . ${\ 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} \ | _ ^ {-1} r_ {0}}$

Dlatego wektor można zapisać tak , jak z , gdzie jest macierzą m -by- n utworzoną przez . ${\ Displaystyle x_ {n} \ w K_ {n}}$${\ Displaystyle x_ {n} = x_ {0}+ Q_ {n} y_ {n}}$${\ Displaystyle y_ {n} \ w \ mathbb {R} ^ {n}}$${\ Displaystyle Q_ {n}}$${\displaystyle q_{1},\ldots,q_{n}}$

Proces Arnoldiego wytwarza również macierz ( )-przez- górny Hessenberg z ${\displaystyle n+1}$${\ Displaystyle n}$ ${\ Displaystyle {\ tylda {H}} _ {n}}$

${\ Displaystyle AQ_ {n} = Q_ {n + 1} {\ tylda {H}} _ {n}. \,}$

W przypadku macierzy symetrycznych faktycznie uzyskuje się symetryczną macierz trójprzekątną , co skutkuje metodą minres .

Ponieważ kolumny są ortonormalne, mamy ${\ Displaystyle Q_ {n}}$

${\ Displaystyle \ | r_ {n} \ | = \ | b-Ax_ {n} \ | = \ | bA (x_ {0}+ Q_ {n} y_ {n}) \ | = \ | r_ {0} -AQ_{n}y_{n}\|=\|\beta q_{1}-AQ_{n}y_{n}\|=\|\beta q_{1}-Q_{n+1}{\tylda {H}}_{n}y_{n}\|=\|Q_{n+1}(\beta e_{1}-{\tylda {H}}_{n}y_{n})\|= \|\beta e_{1}-{\tylda {H}}_{n}y_{n}\|,\,}$

gdzie

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

Jest to pierwszy wektor w standardowym oparciu o , i ${\ Displaystyle \ mathbb {R} ^ {n + 1}}$

${\displaystyle \beta =\|r_{0}\|\,,}$

${\displaystyle x_{0}}$będąca pierwszym wektorem próbnym (zwykle zero). W związku z tym można je znaleźć minimalizując normę euklidesową reszty ${\displaystyle x_{n}}$

${\ Displaystyle R_ {n} = {\ tylda {H}} _ {n} y_ {n} - \ beta e_ {1}.}$

Jest to liniowy problem najmniejszych kwadratów o rozmiarze n .

Daje to metodę GMRES. W -tej iteracji: ${\ Displaystyle n}$

1. obliczyć metodą Arnoldiego;${\displaystyle q_{n}}$
2. znajdź ten, który minimalizuje ;${\displaystyle y_{n}}$${\ Displaystyle \ | r_ {n} \ |}$
3. obliczyć ;${\ Displaystyle x_ {n} = x_ {0}+ Q_ {n} y_ {n}}$
4. powtórz, jeśli pozostałość nie jest jeszcze wystarczająco mała.

W każdej iteracji należy obliczyć iloczyn macierzy i wektora . To kosztuje około operacji zmiennoprzecinkowych na walnym Dense matryce wielkości , ale koszt może zmniejszyć się do rzadkich macierzy . Oprócz iloczynu macierz-wektor, operacje zmiennoprzecinkowe muszą być obliczane w n- tej iteracji. ${\ Displaystyle Aq_ {n}}$${\ Displaystyle 2m ^ {2}}$ ${\ Displaystyle m}$${\ Displaystyle O (m)}$${\ Displaystyle O (nm)}$

Konwergencja

N p iteracyjne minimalizuje resztkowej w podprzestrzeni Kryłowa . Ponieważ każda podprzestrzeń jest zawarta w następnej podprzestrzeni, reszta nie wzrasta. Po m iteracji, gdzie m oznacza wielkość macierzy A The Kryłowa przestrzeń K _m jest cała R ^m, i stąd też sposób GMRES dochodzi do dokładnego rozwiązania. Jednak pomysł jest taki, że po niewielkiej liczbie iteracji (względem m ) wektor x _n jest już dobrym przybliżeniem dokładnego rozwiązania. ${\ Displaystyle K_ {n}}$

Na ogół tak się nie dzieje. Rzeczywiście, twierdzenie Greenbauma, Ptáka i Strakoša mówi, że dla każdego nierosnącego ciągu a ₁ , …, a _{m −1} , a _m = 0, można znaleźć macierz A taką, że || r _n || = a _n dla wszystkich n , gdzie r _n jest resztą określoną powyżej. W szczególności można znaleźć macierz, dla której reszta pozostaje stała dla m  − 1 iteracji i spada do zera dopiero w ostatniej iteracji.

W praktyce jednak GMRES często działa dobrze. Można to udowodnić w określonych sytuacjach. Jeżeli symetryczna część A , czyli , jest dodatnio określona , to ${\ Displaystyle (A ^ {T} + A) / 2}$

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

gdzie i oznaczają odpowiednio najmniejszą i największą wartość własną macierzy . ${\ Displaystyle \ lambda _ {\ operatorname {min}} (M)}$${\ Displaystyle \ lambda _ {\ operatorname {max}} (M)}$${\ Displaystyle M}$

Jeśli A jest symetryczne i dodatnio określone, to mamy nawet

${\ Displaystyle \ | r_ {n} \ | \ leq \ lewo ({\ Frac {\ kappa _ {2} (A) ^ {2}-1} {\ kappa _ {2} (A) ^ {2} }}\right)^{n/2}\|r_{0}\|.}$

gdzie oznacza liczbę warunek z A w normy euklidesowej. ${\ Displaystyle \ kappa _ {2} (A)}$

W ogólnym przypadku, gdy A nie jest dodatnio określone, mamy

${\ Displaystyle {\ Frac {\ | r_ {n} \ |} \ leq \ inf _ {p \ w P_ {n}} \ | p (A) \ | \ leq \ kappa _{2}(V)\inf _{p\in P_{n}}\max _{\lambda \in \sigma (A)}|p(\lambda )|,\,}$

gdzie P _n oznacza zestaw wielomianów stopnia co najwyżej n o p (0) = 1, V jest macierzą występujące w widmowej rozkładu z A i σ ( ) jest widmem z A . Z grubsza mówiąc, mówi to, że szybka zbieżność występuje, gdy wartości własne A są zgrupowane z dala od początku i A nie jest zbyt daleko od normalności .

Wszystkie te nierówności ograniczają tylko reszty zamiast rzeczywistego błędu, czyli odległość między bieżącą iteracją x _n a dokładnym rozwiązaniem.

Rozszerzenia metody

Podobnie jak inne metody iteracyjne, GMRES jest zwykle łączony z metodą warunkowania wstępnego w celu przyspieszenia zbieżności.

Koszt iteracji rośnie jako O( n ² ), gdzie n jest liczbą iteracji. Dlatego też metoda jest czasami ponownie uruchamiana po pewnej liczbie, powiedzmy k , iteracji, z x _k jako początkowym przypuszczeniem. Wynikowa metoda nazywa się GMRES( k ) lub Restarted GMRES. W przypadku niedodatnich macierzy określonych, metoda ta może ucierpieć z powodu stagnacji w zbieżności, ponieważ zrestartowana podprzestrzeń jest często blisko wcześniejszej podprzestrzeni.

Wady GMRES i zrestartowanego GMRES są usuwane przez recykling podprzestrzeni Kryłowa w metodach typu GCRO, takich jak GCROT i GCRODR. Recykling podprzestrzeni Kryłowa w GMRES może również przyspieszyć zbieżność, gdy trzeba rozwiązać sekwencje układów liniowych.

Porównanie z innymi solverami

Iteracja Arnoldiego redukuje się do iteracji Lanczosa dla macierzy symetrycznych. Odpowiednią metodą podprzestrzeni Kryłowa jest metoda minimalnych reszt (MinRes) Paige i Saundersa. W przeciwieństwie do przypadku niesymetrycznego, metoda MinRes jest dana przez trzyokresową relację rekurencyjną . Można wykazać, że nie ma metody podprzestrzennej Kryłowa dla macierzy ogólnych, która jest określona przez krótką relację rekurencyjną, a jednak minimalizuje normy reszt, tak jak robi to GMRES.

Inna klasa metod opiera się na niesymetrycznej iteracji Lanczosa , w szczególności metoda BiCG . Wykorzystują one trzyokresową relację rekurencyjności, ale nie osiągają minimalnej reszty, a zatem reszta nie zmniejsza się monotonicznie dla tych metod. Konwergencja nie jest nawet gwarantowana.

Trzecią klasę tworzą metody takie jak CGS i BiCGSTAB . Działają one również z trójokresową relacją nawrotu (a więc bez optymalności) i mogą nawet zakończyć się przedwcześnie bez osiągnięcia zbieżności. Ideą tych metod jest odpowiedni dobór wielomianów generujących ciągu iteracji.

Żadna z tych trzech klas nie jest najlepsza dla wszystkich macierzy; zawsze są przykłady, w których jedna klasa przewyższa drugą. Dlatego w praktyce próbuje się wielu solverów, aby sprawdzić, który z nich jest najlepszy dla danego problemu.

Rozwiązywanie problemu najmniejszych kwadratów

Jedną z części metody GMRES jest znalezienie wektora, który minimalizuje ${\displaystyle y_{n}}$

${\ Displaystyle \ | {\ tylda {H}} _ {n} Y_ {n} - \ beta e_ {1} \ |. \,}$

Zauważ, że jest to macierz ( n  + 1)-by- n , stąd daje ona nadmiernie ograniczony układ liniowy n +1 równań dla n niewiadomych. ${\ Displaystyle {\ tylda {H}} _ {n}}$

Minimum można obliczyć za pomocą rozkładu QR : znajdź ( n  + 1) - przez - ( n  + 1) macierz ortogonalną Ω _n i ( n  + 1) - przez - n macierz trójkątną górną taką, że ${\ Displaystyle {\ tylda {R}} _ {n}}$

${\ Displaystyle \ Omega _ {n} {\ tylda {H}} _ {n} = {\ tylda {R}} _ {n}.}$

Macierz trójkątna ma o jeden wiersz więcej niż kolumn, więc jej dolny wiersz składa się z zera. W związku z tym można go rozłożyć jako

${\ Displaystyle {\ tylda {R}} _ {n} = {\ zacząć {bmatrix} R_ {n} \ \ 0 \ koniec {bmatrix}}}$

gdzie jest macierzą trójkątną n -by- n (a więc kwadratową). ${\ Displaystyle R_ {n}}$

Rozkład QR można tanio aktualizować z jednej iteracji do następnej, ponieważ macierze Hessenberga różnią się tylko rzędem zer i kolumną:

${\ Displaystyle {\ tylda {H}} _ {n + 1} = {\ zacząć {bmatrix} {\ tylda {H}} _ {n} i h_ {n + 1} \ \ 0 i h_ {n + 2, n + 1}\koniec{bmatrycy}},}$

gdzie h _n+1 = ( h _{1, n+1} , …, h _n+1,n+1 ) ^T . Oznacza to, że wstępne pomnożenie macierzy Hessenberga przez Ω _n , powiększoną o zera i wiersz o identyczności multiplikatywnej, daje macierz prawie trójkątną:

${\ Displaystyle {\ zacząć {bmatrix} \ Omega _ {n} i 0 \ \ 0 i 1 \ koniec {bmatrix}} {\ tylda {H}} _ {n + 1} = {\ zacząć {bmatrix} R_ {n} i r_ {n+1}\\0&\rho \\0&\sigma \end{bmatrix}}}$

Byłoby to trójkątne, gdyby σ wynosiło zero. Aby temu zaradzić, potrzebna jest rotacja Givens

${\ Displaystyle G_ {n} = {\ zacząć {bmatrix} I_ {n} i 0 i 0 \ \ 0 i c_ {n} i s_ {n} \ \ 0 i s_ {n} i c_ {n} \ koniec {bmatrix}}}$

gdzie

${\ Displaystyle c_ {n} = {\ Frac {\ rho} {\ sqrt {\ rho ^ {2} + \ sigma ^ {2}}}} \ quad {\ mbox {i}} \ quad s_ {n} ={\frac {\sigma }{\sqrt {\rho ^{2}+\sigma ^{2}}}}.}$

Dzięki tej rotacji Givens tworzymy

${\ Displaystyle \ Omega _ {n + 1} = G_ {n} {\ zacząć {bmatrix} \ Omega _ {n} i 0 \ \ 0 i 1 \ koniec {bmatrix}}.}$

W rzeczy samej,

${\ Displaystyle \ Omega _ {n + 1} {\ tylda {H}} _ {n + 1} = {\ zacznij {bmatrix} R_ {n} i r_ {n + 1} \ \ 0 i r_ {n + 1, n +1}\\0&0\end{bmatrix}}\quad {\text{z}}\quad r_{n+1,n+1}={\sqrt {\rho ^{2}+\sigma ^{2 }}}}$

jest macierzą trójkątną.

Biorąc pod uwagę rozkład QR, problem minimalizacji można łatwo rozwiązać, zauważając, że

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

Oznaczanie wektora przez ${\ Displaystyle \ beta \ Omega _ {n} e_ {1}}$

${\ Displaystyle {\ tylda {g}} _ {n} = {\ zacznij {bmatrix} g_ {n} \ \ \ gamma _ {n} \ koniec {bmatrix}}}$

z g _n ∈ R ⁿ i γ _n ∈ R , to jest

${\ Displaystyle \ | {\ tylda {H}} _ {n} r_ {n} - \ beta e_ {1} \ | = \ | {\ tylda {R}} _ {n} r_ {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\|.}$

Wektor y minimalizujący to wyrażenie jest podany przez

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

Znowu wektory można łatwo aktualizować. ${\displaystyle g_{n}}$

Przykładowy kod

Zwykłe GMRES (MATLAB / Oktawa GNU)

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;
  beta = r_norm * e1;     %Note: this is not the beta scalar in section "The method" above but the beta scalar multiplied by 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 Given 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

Zobacz też

• Metoda gradientu dwuskoniugowanego

Bibliografia

Uwagi

• A. Meister, Numerik linearer Gleichungssysteme , wydanie 2, Vieweg 2005, ISBN  978-3-528-13135-7 .
• Y. Saad, Iteracyjne metody dla rzadkich układów liniowych , wydanie 2, Towarzystwo Matematyki Przemysłowej i Stosowanej , 2003. ISBN  978-0-89871-534-7 .
• Y. Saad i MH Schultz, „GMRES: Uogólniony algorytm minimalnej reszty do rozwiązywania niesymetrycznych układów liniowych”, SIAM J. Sci. Stat. Komputer. , 7 :856-869, 1986. doi : 10.1137/0907058 .
• SC Eisenstat, HC Elman i MH Schultz, „Wariacyjne metody iteracyjne dla niesymetrycznych układów równań liniowych”, SIAM Journal on Numerical Analysis , 20(2), 345-357, 1983.
• J. Stoer i R. Bulirsch, Wprowadzenie do analizy numerycznej , wydanie 3, Springer, New York, 2002. ISBN  978-0-387-95452-3 .
• Lloyd N. Trefethen i David Bau, III, Numeryczna Algebra Liniowa , Towarzystwo Matematyki Przemysłowej i Stosowanej, 1997. ISBN  978-0-89871-361-9 .
• Dongarra i in. , Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods , wydanie 2, SIAM, Filadelfia, 1994
• Amritkar, Amit; de Sturler, Eric; Świrydowicz, Katarzyna; tafti, duński; Ahuja, Kapil (2015). „Recykling podprzestrzeni Kryłowa do zastosowań CFD i nowy hybrydowy solwer do recyklingu”. Journal of Computational Physics 303: 222. doi:10.1016/j.jcp.2015.09.040