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 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

$\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 {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:

$\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

Jeśli właściwie dobierzemy sprzężone wektory p_k, 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 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(\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 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:

$\mathbf{r}_k = \mathbf{b} - \mathbf{Ax}_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:

$\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

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^\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 r_k+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

 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ż

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

Bibliografia

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

• Magnus R. Hestenes, Eduard Stiefel. Methods of Conjugate Gradients for Solving Linear Systems. „Journal of Research of the National Bureau of Standards”. 49 (6), grudzień 1952. 

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

• 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.