Metoda gradientu sprzężonego

Z Wikipedii, wolnej encyklopedii
Skocz do: nawigacja, szukaj
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 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[edytuj | edytuj kod]

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[edytuj | edytuj kod]

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

 \mathbf{u}^{\mathrm{T}} \mathbf{A} \mathbf{v} = \mathbf{0}.

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

 \langle \mathbf{u},\mathbf{v} \rangle_\mathbf{A} := \langle \mathbf{A}^{\mathrm{T}} \mathbf{u}, \mathbf{v}\rangle = \langle \mathbf{A} \mathbf{u}, \mathbf{v}\rangle = \langle \mathbf{u}, \mathbf{A}\mathbf{v} \rangle = \mathbf{u}^{\mathrm{T}} \mathbf{A} \mathbf{v}.

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

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

 \mathbf{x}_* = \sum^{n}_{i=1} \alpha_i \mathbf{p}_i

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

 \mathbf{A}\mathbf{x}_* = \sum^{n}_{i=1} \alpha_i  \mathbf{A} \mathbf{p}_i = \mathbf{b}.
 \mathbf{p}_k^{\mathrm{T}} \mathbf{A}\mathbf{x}_* = \sum^{n}_{i=1} \alpha_i\mathbf{p}_k^{\mathrm{T}} \mathbf{A} \mathbf{p}_i= \mathbf{p}_k^{\mathrm{T}} \mathbf{b}.
 \alpha_k = \frac{\mathbf{p}_k^{\mathrm{T}} \mathbf{b}}{\mathbf{p}_k^{\mathrm{T}} \mathbf{A} \mathbf{p}_k} = \frac{\langle \mathbf{p}_k, \mathbf{b}\rangle}{\,\,\,\langle \mathbf{p}_k,  \mathbf{p}_k\rangle_\mathbf{A}} = \frac{\langle \mathbf{p}_k, \mathbf{b}\rangle}{\,\,\,\|\mathbf{p}_k\|_\mathbf{A}^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 αk.

Metoda gradientu sprzężonego jako metoda iteracyjna[edytuj | edytuj kod]

Jeśli właściwie dobierzemy sprzężone wektory pk, możemy nie potrzebować ich wszystkich do dobrej aproksymacji rozwiązania x*Więc chcielibyśmy 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 x0. Bez starty ogólności możemy założyć, że x0 = 0 (w przeciwnym przypadku, rozważymy układ Az = bAx0). Zauważmy, że rozwiązanie x* minimalizuje formę kwadratową:

 f(\mathbf{x}) = \frac12 \mathbf{x}^{\mathrm{T}} \mathbf{A}\mathbf{x} - \mathbf{b}^{\mathrm{T}} \mathbf{x} , \quad \mathbf{x}\in\mathbf{R}^n.

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

Niech rk oznacza rezyduum w k-tym kroku:

 \mathbf{r}_k = \mathbf{b} - \mathbf{Ax}_k.  \,

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

 \mathbf{p}_{k+1} = \mathbf{r}_k - \frac{\mathbf{p}_k^{\mathrm{T}} \mathbf{A} \mathbf{r}_k}{\mathbf{p}_k^{\mathrm{T}}\mathbf{A} \mathbf{p}_k} \mathbf{p}_k

Wynikowy algorytm[edytuj | edytuj kod]

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

r_0 := b - A x_0 \,
p_0 := r_0 \,
k := 0 \,
repeat
\alpha_k := \frac{r_k^\top r_k}{p_k^\top A p_k}  \,
x_{k+1} := x_k + \alpha_k p_k \,
r_{k+1} := r_k - \alpha_k A p_k \,
if rk+1 jest "wystarczająco mały" then exit loop end if
\beta_k := \frac{r_{k+1}^\top r_{k+1}}{r_k^\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[edytuj | edytuj kod]

function [x] = conjgrad(A,b,x0)
  
    r = b - A*x0;
    w = -r;
    z = A*w;
    a = (r'*w)/(w'*z);
    x = x0 + a*w;
    B = 0;
  
    for i = 1:size(A)(1);
       r = r - a*z;
       if( norm(r) < 1e-10 )
            break;
       endif
       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ż[edytuj | edytuj kod]

Bibliografia[edytuj | edytuj kod]

Metoda gradientu sprzężonego została zaproponowana w:

Opisy meteody można znaleźć w:

  • Kendell A. Atkinson (1988), An introduction to numerical analysis (2nd ed.), Section 8.9, John Wiley and Sons. ISBN 0-471-50023-2.
  • Mordecai Avriel (2003). Nonlinear Programming: Analysis and Methods. Dover Publishing. ISBN 0-486-43227-0.
  • Gene H. Golub and Charles F. Van Loan, Matrix computations (3rd ed.), Chapter 10, Johns Hopkins University Press. ISBN 0-8018-5414-8.

Linki zewnętrzne[edytuj | edytuj kod]