Metoda gradientu sprzężonego

Porównanie zbieżności metody gradientu prostego (na zielono) i sprzężonego gradientu (na czerwono) dla minimalizacji formy kwadratowej powiązanej z zadanym układem równań liniowych. Metoda sprzężonego gradientu zbiega w co najwyżej n krokach, gdzie n jest rozmiarem macierzy zadanego układu (tutaj n=2).

Metoda gradientu sprzężonego (ang. conjugate gradient method, w skrócie CG) jest algorytmem numerycznym służącym do rozwiązywania niektórych układów równań liniowych. Pozwala rozwiązać te, których macierz jest symetryczna i dodatnio określona. Metoda gradientu sprzężonego jest metodą iteracyjną, więc może być zastosowana do układów o rzadkich macierzach, które mogą być zbyt duże dla algorytmów bezpośrednich takich jak np. rozkład Choleskiego. Takie układy pojawiają się często w trakcie numerycznego rozwiązywania równań różniczkowych cząstkowych.

Metoda gradientu sprzężonego może również zostać użyta do rozwiązania problemu optymalizacji bez ograniczeń.

Opis metody

Rozpatrzmy rozwiązania poniższego układu równań:

Ax = b,

gdzie macierz A n na n jest symetryczna, rzeczywista i dodatnio określona.

Oznaczmy rozwiązanie tego układu przez x_*.

Bezpośrednia metoda gradientu sprzężonego

Mówimy, że dwa niezerowe wektory u i v są sprzężone (względem A), jeśli

${𝐮}^{\mathrm{T}}𝐀𝐯=\mathrm{𝟎}.$

Ponieważ A jest symetryczna i dodatnio określona, lewa strona definiuje iloczyn skalarny:

$⟨𝐮,𝐯{⟩}_{𝐀}:=⟨{𝐀}^{\mathrm{T}}𝐮,𝐯⟩=⟨𝐀𝐮,𝐯⟩=⟨𝐮,𝐀𝐯⟩={𝐮}^{\mathrm{T}}𝐀𝐯.$

Więc, dwa wektory są sprzężone jeśli są ortogonalne względem tego iloczynu skalarnego. Sprzężoność jest relacją symetryczną.

Przypuśćmy, że {p_k} jest ciągiem n wzajemnie sprzężonych kierunków. Wtedy p_k tworzą bazę Rⁿ, wektor x_* będący rozwiązaniem Ax = b możemy przedstawić w postaci:

${𝐱}_{*}=\sum _{i=1}^n{\alpha }_i{𝐩}_i.$

Współczynniki otrzymujemy w następujący sposób:

$𝐀{𝐱}_{*}=\sum _{i=1}^n{\alpha }_i𝐀{𝐩}_i=𝐛,$

${𝐩}_k^{\mathrm{T}}𝐀{𝐱}_{*}=\sum _{i=1}^n{\alpha }_i{𝐩}_k^{\mathrm{T}}𝐀{𝐩}_i={𝐩}_k^{\mathrm{T}}𝐛,$

${\alpha }_k=\frac{{𝐩}_k^{\mathrm{T}}𝐛}{{𝐩}_k^{\mathrm{T}}𝐀{𝐩}_k}=\frac{⟨{𝐩}_k,𝐛⟩}{⟨{𝐩}_k,{𝐩}_k{⟩}_{𝐀}}=\frac{⟨{𝐩}_k,𝐛⟩}{\Vert {𝐩}_k{\Vert }_{𝐀}^2}.$

Co daje nam następującą metodę rozwiązywania równania Ax = b. Najpierw znajdujemy ciąg n sprzężonych kierunków, następnie obliczamy współczynniki ${\alpha }_k.$

Metoda gradientu sprzężonego jako metoda iteracyjna

Jeśli właściwie dobierzemy sprzężone wektory p_k, możemy nie potrzebować ich wszystkich do dobrej aproksymacji rozwiązania x_*. Możemy więc spojrzeć na CG jak na metodę iteracyjną. Co więcej, pozwoli nam to rozwiązać układy równań, gdzie n jest tak duże, że bezpośrednia metoda zabrałaby zbyt dużo czasu.

Oznaczmy punkt startowy przez x₀. Bez starty ogólności możemy założyć, że x₀ = 0 (w przeciwnym przypadku, rozważymy układ Az = b − Ax₀). Zauważmy, że rozwiązanie x_* minimalizuje formę kwadratową:

$f(𝐱)=\frac{1}{2}{𝐱}^{\mathrm{T}}𝐀𝐱-{𝐱}^{\mathrm{T}}𝐛,𝐱\in {𝐑}^n.$

Co sugeruje, by jako pierwszy wektor bazowy p₁ wybrać gradient f w x = x₀, który wynosi Ax₀−b, a ponieważ wybraliśmy x₀ = 0, otrzymujemy −b. Pozostałe wektory w bazie będą sprzężone do gradientu (stąd nazwa metoda gradientu sprzężonego).

Niech r_k oznacza rezyduum w k-tym kroku:

${𝐫}_k=𝐛-{𝐀𝐱}_k.$

Zauważmy, że r_k jest przeciwny do gradientu f w x = x_k, więc metoda gradientu prostego nakazywałaby ruch w kierunku r_k. Tutaj jednak założyliśmy wzajemną sprzężoność kierunków p_k, więc wybieramy kierunek najbliższy do r_k pod warunkiem sprzężoności. Co wyraża się wzorem:

${𝐩}_{k+1}={𝐫}_k-\frac{{𝐩}_k^{\mathrm{T}}𝐀{𝐫}_k}{{𝐩}_k^{\mathrm{T}}𝐀{𝐩}_k}{𝐩}_k.$

Wynikowy algorytm

Upraszczając, otrzymujemy poniższy algorytm rozwiązujący Ax = b, gdzie macierz A jest rzeczywista, symetryczna i dodatnio określona. x₀ jest punktem startowym.

$r_0:=b-A{x}_0$

$p_0:=r_0$

$k:=0$

repeat

${\alpha }_k:=\frac{r_k^{\mathrm{\top }}{r}_k}{p_k^{\mathrm{\top }}A{p}_k}$

$x_{k+1}:=x_k+{\alpha }_k{p}_k$

$r_{k+1}:=r_k-{\alpha }_{k}A{p}_k$

if r_k+1 jest "wystarczająco mały" then exit loop end if

${\beta }_k:=\frac{r_{k+1}^{\mathrm{\top }}{r}_{k+1}}{r_k^{\mathrm{\top }}{r}_k}$

$p_{k+1}:=r_{k+1}+{\beta }_k{p}_k$

$k:=k+1$

end repeat

Wynikiem jest $x_{k+1}$

Przykład metody gradientu sprzężonego w Octave/MATLAB

function [x] = conjgrad(A,b,x0)
r = b - A*x0;
w = -r;
z = A*w;
a = (r'*w)/(w'*z);
x = x0 + a*w;

for i = 1:size(A,1);
    r = r - a*z;
    if( norm(r) < 1e-10 )
        break;
    end
    B = (r'*z)/(w'*z);
    w = -r + B*w;
    z = A*w;
    a = (r'*w)/(w'*z);
    x = x + a*w;
end

end

Zobacz też

• metoda gradientu prostego
• metoda najszybszego spadku
• metoda Newtona
• metoda złotego podziału
• optymalizacja

Bibliografia

• Metoda gradientu sprzężonego została zaproponowana w:
  • Szablon:Cytuj
• Opisy meteody można znaleźć w:
  • Kendell A. Atkinson (1988), An introduction to numerical analysis (2nd ed.), Section 8.9, John Wiley and Sons. Szablon:ISBN.
  • Mordecai Avriel (2003). Nonlinear Programming: Analysis and Methods. Dover Publishing. Szablon:ISBN.
  • Gene H. Golub and Charles F. Van Loan, Matrix computations (3rd ed.), Chapter 10, Johns Hopkins University Press. Szablon:ISBN.

Linki zewnętrzne

• Conjugate Gradient Method by Nadir Soualem.
• Preconditioned conjugate gradient method by Nadir Soualem.
• An Introduction to the Conjugate Gradient Method Without the Agonizing Pain by Jonathan Richard Shewchuk.
• Iterative methods for sparse linear systems by Yousef Saad
• LSQR: Sparse Equations and Least Squares by Christopher Paige and Michael Saunders.