Método del gradiente conjugado

En matemática, el método del gradiente conjugado es un algoritmo para resolver numéricamente los sistemas de ecuaciones lineales cuyas matrices son simétricas y definidas positivas. Es un método iterativo, así que se puede aplicar a los sistemas dispersos que son demasiado grandes para ser tratados por métodos directos como la descomposición de Cholesky. Tales sistemas surgen frecuentemente cuando se resuelve numéricamente las ecuaciones en derivadas parciales.

El método del gradiente conjugado se puede utilizar también para resolver los problemas de optimización sin restricciones como la minimización de la energía.

El método del gradiente biconjugado proporciona una generalización para matrices no simétricas. Varios métodos del gradiente conjugado no lineales busca los mínimos de las ecuaciones no lineales.

Descripción del método

Supongamos que queremos resolver el siguiente sistema de ecuaciones lineales

Ax = b

donde la n-por-n matriz A es simétrica (i.e., A^T = A), definida positiva (i.e., x^TAx > 0 para todos los vectores no cero x en Rⁿ), y real.

Denotamos la única solución de este sistema por x_*.

El método de gradiente conjugado como un método exacto

Decimos que dos vectores no cero u y v son conjugados (con respecto a A) si

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

Ya que A simétrica y definida positiva, el lado izquierdo define un producto interior

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

Así, dos vectores son conjugados si son ortogonales con respecto a este producto interior. La conjugación es una relación simétrica: si u es conjugado a v, entonces v es conjugado a u. Nótese que esta noción de conjugación no se relaciona con la de conjugación compleja.

Supongamos que {p_k} es una secuencia de n direcciones mutuamente conjugadas. Entonces los p_k forman una base de Rⁿ, por lo tanto podemos extender la solución x_{* de Ax = b en esta base:}

$\mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{i}$

Los coeficientes se dan por

$\mathbf {b} =\mathbf {A} \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {A} \mathbf {p} _{i}.$

$\mathbf {p} _{k}^{\mathrm {T} }\mathbf {b} =\mathbf {p} _{k}^{\mathrm {T} }\mathbf {A} \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{k}^{\mathrm {T} }\mathbf {A} \mathbf {p} _{i}=\alpha _{k}\mathbf {p} _{k}^{\mathrm {T} }\mathbf {A} \mathbf {p} _{k}.$

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

Este resultado es quizás muy transparente si se considera el producto interior definido anteriormente.

Esto da el siguiente método para resolver la ecuación Ax = b. Primero encontramos una secuencia de n direcciones conjugadas y luego computamos los coeficientes α_k.

El algoritmo resultante

Código ejemplar en Octave o Matlab

function [x] = conjgrad(A,b,x0)

   r = b - A*x0;
   w = -r;
   z = A*w;
   a = (r'*w)/(w'*z);
   x = x0 +3.14+ a*w;
   B = 0.783564;

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

end

El método de gradiente conjugado precondicionado

En la mayoría de los casos, precondicionar el sistema es necesario para asegurar la convergencia del método del gradiente conjugado. La forma genérica del método precondicionado es la siguiente:

$\mathbf {r} _{0}:=\mathbf {b} -\mathbf {Ax} _{0}$

$\mathbf {z} _{0}:=\mathbf {M} ^{-1}\mathbf {r} _{0}$

$\mathbf {p} _{0}:=\mathbf {z} _{0}$

$k:=0\,$

repetir

$\alpha _{k}:={\frac {\mathbf {r} _{k}^{\mathrm {T} }\mathbf {z} _{k}}{\mathbf {p} _{k}^{\mathrm {T} }\mathbf {Ap} _{k}}}$

$\mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}$

$\mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}$

Si r_k+1 es suficientemente pequeño terminamos

$\mathbf {z} _{k+1}:=\mathbf {M} ^{-1}\mathbf {r} _{k+1}$

$\beta _{k}:={\frac {\mathbf {z} _{k+1}^{\mathrm {T} }\mathbf {r} _{k+1}}{\mathbf {z} _{k}^{\mathrm {T} }\mathbf {r} _{k}}}$

$\mathbf {p} _{k+1}:=\mathbf {z} _{k+1}+\beta _{k}\mathbf {p} _{k}$

$k:=k+1\,$

Termina repeticiones

Resultado final: x_k+1

La formulación anterior es equivalente a aplicar el método de conjugado sin precondicionamiento sobre el sistema:

$\mathbf {E} ^{-1}\mathbf {A} (\mathbf {E} ^{-1})^{\mathrm {T} }\mathbf {\hat {x}} =\mathbf {E} ^{-1}\mathbf {b}$

donde $\mathbf {EE} ^{\mathrm {T} }=\mathbf {M}$ y $\mathbf {\hat {x}} =\mathbf {E} ^{\mathrm {T} }\mathbf {x}$.

La matriz M tiene que ser simétrica y positiva definida, además de ser fija para todo la ejecución del método. Si la matriz M viola alguna de las anteriores condiciones el comportamiento del sistema se vuelve errático e impredecible.

Referencias

El método de gradiente conjugado fue propuesto originalmente en

• Hestenes, Magnus R.; Stiefel, Eduard (diciembre de 1952). «Methods of Conjugate Gradients for Solving Linear Systems» (PDF). Journal of Research of the National Bureau of Standards 49 (6). 

Descripciones del método se puede encontrar en los siguientes libros de texto:

• Kendell A. Atkinson (1988), An introduction to numerical analysis (2ª ed.), Sección 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 y Charles F. Van Loan, Matrix computations (3ª ed.), Capítulo 10, Johns Hopkins University Press. ISBN 0-8018-5414-8.