Usuario:Kxx/Sandbox/Derivación del método del gradiente conjugado

En la álgebra lineal numérica, el método del gradiente conjugado es un método iterativo para la resolución numérica del sistema lineal

${\boldsymbol {Ax}}={\boldsymbol {b}}$

donde ${\boldsymbol {A}}$ es simétrica definida positiva. El método del gradiente conjugado se puede derivar de varias perspectivas diferentes, incluidas la especialización del método de las direcciones conjugadas para la optimización, y la variación de la iteración de Arnoldi/Lanczos para los problemas de los valores propios.

El intento de este artículo es documentar los paso importantes en estas derivaciones.

Derivación del método de las direcciones conjugadas

Hestenes y Stiefel, los inventores del método del gradiente conjugado, derivaron el método del más general método de las direcciones conjugadas en el contexto de la minimización de la función cuadrática

$f(\mathbf {x} )={\frac {1}{2}}\mathbf {x} ^{\mathrm {T} }\mathbf {Ax} -\mathbf {b} ^{\mathrm {T} }\mathbf {x} {\text{.}}$

En el método de las direcciones conjugadas se elige una serie de direcciones $\mathbf {p} _{0},\mathbf {p} _{1},\mathbf {p} _{2},\ldots$ que son conjugadas mutuamente, en otras palabras,

$\langle \mathbf {p} _{i},\mathbf {Ap} _{i}\rangle =0$

para $i\neq j$. La iteración se inicia con una solución aproximada $\mathbf {x} _{0}$ y el correspondiente residuo $\mathbf {r} _{0}$ y procede con los siguientes pasos para $i=0,1,2,\ldots$:

$\alpha _{i}={\frac {\langle \mathbf {p} _{i},\mathbf {r} _{i}\rangle }{\langle \mathbf {p} _{i},\mathbf {Ap} _{i}\rangle }}{\text{,}}$

$\mathbf {x} _{i+1}=\mathbf {x} _{i}+\alpha _{i}\mathbf {p} _{i}{\text{,}}$

$\mathbf {r} _{i+1}=\mathbf {r} _{i}+\alpha _{i}\mathbf {Ap} _{i}{\text{.}}$

Derivación de la iteración de Arnoldi/Lanczos

Véanse también: Iteración de Arnoldi e Iteración de Lanczos.

El método del gradiente conjugado también se puede ver como una variante de la iteración de Arnoldi/Lanczos aplicada a la solución de los sistemas lineales.

El método de Arnoldi general

En la iteración de Arnoldi, se comienza con un vector ${\boldsymbol {r}}_{0}$ y construye gradualmente una base ortonormal $\{{\boldsymbol {v}}_{1},{\boldsymbol {v}}_{2},{\boldsymbol {v}}_{3},\ldots \}$ del subespacio de Krylov

${\mathcal {K}}({\boldsymbol {A}},{\boldsymbol {r}}_{0})=\{{\boldsymbol {r}}_{0},{\boldsymbol {Ar}}_{0},{\boldsymbol {A}}^{2}{\boldsymbol {r}}_{0},\ldots \}$

mediante la definición de ${\boldsymbol {v}}_{i}={\boldsymbol {w}}_{i}/\lVert {\boldsymbol {w}}_{i}\rVert _{2}$ donde

${\boldsymbol {w}}_{i}={\begin{cases}{\boldsymbol {r}}_{0}&{\text{si }}i=1{\text{,}}\\{\boldsymbol {Av}}_{i-1}-\sum _{j=1}^{i-1}({\boldsymbol {v}}_{j}^{\mathrm {T} }{\boldsymbol {Av}}_{i-1}){\boldsymbol {v}}_{j}&{\text{si }}i>1{\text{.}}\end{cases}}$

En otras palabras, para $i>1$, ${\boldsymbol {v}}_{i}$ se obtiene a través de la ortogonalización de Gram-Schmidt de ${\boldsymbol {Av}}_{i-1}$ contra $\{{\boldsymbol {v}}_{1},{\boldsymbol {v}}_{2},\ldots ,{\boldsymbol {v}}_{i-1}\}$ seguida por la normalización.

Expresada en la forma matricial, la iteración se capta por la ecuación

${\boldsymbol {AV}}_{i}={\boldsymbol {V}}_{i+1}{\boldsymbol {\tilde {H}}}_{i}$

donde

${\begin{aligned}{\boldsymbol {V}}_{i}&={\begin{bmatrix}{\boldsymbol {v}}_{1}&{\boldsymbol {v}}_{2}&\cdots &{\boldsymbol {v}}_{i}\end{bmatrix}}{\text{,}}\\{\boldsymbol {\tilde {H}}}_{i}&={\begin{bmatrix}h_{11}&h_{12}&h_{13}&\cdots &h_{1,i}\\h_{21}&h_{22}&h_{23}&\cdots &h_{2,i}\\&h_{32}&h_{33}&\cdots &h_{3,i}\\&&\ddots &\ddots &\vdots \\&&&h_{i,i-1}&h_{i,i}\\&&&&h_{i+1,i}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {H}}_{i}\\h_{i+1,i}{\boldsymbol {e}}_{i}^{\mathrm {T} }\end{bmatrix}}\end{aligned}}$

with

$h_{ji}={\begin{cases}{\boldsymbol {v}}_{j}^{\mathrm {T} }{\boldsymbol {Av}}_{i}&{\text{si }}j\leq i{\text{,}}\\\lVert {\boldsymbol {w}}_{i+1}\rVert _{2}&{\text{si }}j=i+1{\text{,}}\\0&{\text{si }}j>i+1{\text{.}}\end{cases}}$

Cuando se aplica la iteración de Arnoldi a la solución de los sistemas lineals, se comienza con ${\boldsymbol {r}}_{0}={\boldsymbol {b}}-{\boldsymbol {Ax}}_{0}$, el residuo correspondiente a la aproximación inicial ${\boldsymbol {x}}_{0}$. Después de cada paso de la iteración, se computa ${\boldsymbol {y}}_{i}={\boldsymbol {H}}_{i}^{-1}(\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {e}}_{1})$ y la nueva aproximación ${\boldsymbol {x}}_{i}={\boldsymbol {x}}_{0}+{\boldsymbol {V}}_{i}{\boldsymbol {y}}_{i}$.

El método de Lanczos directo

Para el resto de la discusión, se supone que ${\boldsymbol {A}}$ es simétrica definida positiva. Con la simetría de ${\boldsymbol {A}}$, la matriz superior de Hessenberg ${\boldsymbol {H}}_{i}={\boldsymbol {V}}_{i}^{\mathrm {T} }{\boldsymbol {AV}}_{i}$ se hace simétrica y por lo tanto tridiagonal. Entonces se puede representar más claramente por

${\boldsymbol {H}}_{i}={\begin{bmatrix}a_{1}&b_{2}\\b_{2}&a_{2}&b_{3}\\&\ddots &\ddots &\ddots \\&&b_{i-1}&a_{i-1}&b_{i}\\&&&b_{i}&a_{i}\end{bmatrix}}{\text{.}}$

Esto permite una corta recurrencia de tres términos para ${\boldsymbol {v}}_{i}$ en la iteración, y así la iteración de Arnoldi se reduce a la iteración de Lanczos.

Ya que ${\boldsymbol {A}}$ es simétrica definida positiva, ${\boldsymbol {H}}_{i}$ también lo es. Por lo tanto, ${\boldsymbol {H}}_{i}$ se puede factorizar LU sin pivoteo parcial en

${\boldsymbol {H}}_{i}={\boldsymbol {L}}_{i}{\boldsymbol {U}}_{i}={\begin{bmatrix}1\\c_{2}&1\\&\ddots &\ddots \\&&c_{i-1}&1\\&&&c_{i}&1\end{bmatrix}}{\begin{bmatrix}d_{1}&b_{2}\\&d_{2}&b_{3}\\&&\ddots &\ddots \\&&&d_{i-1}&b_{i}\\&&&&d_{i}\end{bmatrix}}$

con las recurrencias convenientes para $c_{i}$ y $d_{i}$:

${\begin{aligned}c_{i}&=b_{i}/d_{i-1}{\text{,}}\\d_{i}&={\begin{cases}a_{1}&{\text{si }}i=1{\text{,}}\\a_{i}-c_{i}b_{i}&{\text{si }}i>1{\text{.}}\end{cases}}\end{aligned}}$

Se varía ${\boldsymbol {x}}_{i}={\boldsymbol {x}}_{0}+{\boldsymbol {V}}_{i}{\boldsymbol {y}}_{i}$ en

${\begin{aligned}{\boldsymbol {x}}_{i}&={\boldsymbol {x}}_{0}+{\boldsymbol {V}}_{i}{\boldsymbol {H}}_{i}^{-1}(\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {e}}_{1})\\&={\boldsymbol {x}}_{0}+{\boldsymbol {V}}_{i}{\boldsymbol {U}}_{i}^{-1}{\boldsymbol {L}}_{i}^{-1}(\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {e}}_{1})\\&={\boldsymbol {x}}_{0}+{\boldsymbol {P}}_{i}{\boldsymbol {z}}_{i}\end{aligned}}$

con

${\begin{aligned}{\boldsymbol {P}}_{i}&={\boldsymbol {V}}_{i}{\boldsymbol {U}}_{i}^{-1}{\text{,}}\\{\boldsymbol {z}}_{i}&={\boldsymbol {L}}_{i}^{-1}(\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {e}}_{1}){\text{.}}\end{aligned}}$

Ahora es crucial observar que

${\begin{aligned}{\boldsymbol {P}}_{i}&={\begin{bmatrix}{\boldsymbol {P}}_{i-1}&{\boldsymbol {p}}_{i}\end{bmatrix}}{\text{,}}\\{\boldsymbol {z}}_{i}&={\begin{bmatrix}{\boldsymbol {z}}_{i-1}\\\zeta _{i}\end{bmatrix}}{\text{.}}\end{aligned}}$

De hecho, también existen recurrencias cortas para ${\boldsymbol {p}}_{i}$ y $\zeta _{i}$:

${\begin{aligned}{\boldsymbol {p}}_{i}&={\frac {1}{d_{i}}}({\boldsymbol {v}}_{i}-b_{i}{\boldsymbol {p}}_{i-1}){\text{,}}\\\alpha _{i}&=-c_{i}\zeta _{i-1}{\text{.}}\end{aligned}}$

Con esta formulación, llegamos a una recurrencia simple para ${\boldsymbol {x}}_{i}$:

${\begin{aligned}{\boldsymbol {x}}_{i}&={\boldsymbol {x}}_{0}+{\boldsymbol {P}}_{i}{\boldsymbol {z}}_{i}\\&={\boldsymbol {x}}_{0}+{\boldsymbol {P}}_{i-1}{\boldsymbol {z}}_{i-1}+\zeta _{i}{\boldsymbol {p}}_{i}\\&={\boldsymbol {x}}_{i-1}+\zeta _{i}{\boldsymbol {p}}_{i}{\text{.}}\end{aligned}}$

Las anteriores relaciones resultan sencillamente en el método de Lanczos directo, que es un poco más complejo.

El método del gradiente conjugado de imponer la ortogonalidad y la conjugación

Si se permite que ${\boldsymbol {p}}_{i}$ se escale y compensa el escalamiento en el factor constante, potencialmente se puede obtener recurrencias más simples de la forma:

${\begin{aligned}{\boldsymbol {x}}_{i}&={\boldsymbol {x}}_{i-1}+\alpha _{i-1}{\boldsymbol {p}}_{i-1}{\text{,}}\\{\boldsymbol {r}}_{i}&={\boldsymbol {r}}_{i-1}-\alpha _{i-1}{\boldsymbol {Ap}}_{i-1}{\text{,}}\\{\boldsymbol {p}}_{i}&={\boldsymbol {r}}_{i}+\beta _{i-1}{\boldsymbol {p}}_{i-1}{\text{.}}\end{aligned}}$

Como las premisas de la simplificación, ahora se deriva la ortogonalidad of ${\boldsymbol {r}}_{i}$ y la conjugación de ${\boldsymbol {p}}_{i}$, es decir, para $i\neq j$,

${\begin{aligned}{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{j}&=0{\text{,}}\\{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{j}&=0{\text{.}}\end{aligned}}$

Los residuos son mutuamente ortogonales porque ${\boldsymbol {r}}_{i}$ es esencialmente un múltiplo de ${\boldsymbol {v}}_{i+1}$, dado que para $i=0$, ${\boldsymbol {r}}_{0}=\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {v}}_{1}$, y para $i>0$,

${\begin{aligned}{\boldsymbol {r}}_{i}&={\boldsymbol {b}}-{\boldsymbol {Ax}}_{i}\\&={\boldsymbol {b}}-{\boldsymbol {A}}({\boldsymbol {x}}_{0}+{\boldsymbol {V}}_{i}{\boldsymbol {y}}_{i})\\&={\boldsymbol {r}}_{0}-{\boldsymbol {AV}}_{i}{\boldsymbol {y}}_{i}\\&={\boldsymbol {r}}_{0}-{\boldsymbol {V}}_{i+1}{\boldsymbol {\tilde {H}}}_{i}{\boldsymbol {y}}_{i}\\&={\boldsymbol {r}}_{0}-{\boldsymbol {V}}_{i}{\boldsymbol {H}}_{i}{\boldsymbol {y}}_{i}-h_{i+1,i}({\boldsymbol {e}}_{i}^{\mathrm {T} }{\boldsymbol {y}}_{i}){\boldsymbol {v}}_{i+1}\\&=\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {v}}_{1}-{\boldsymbol {V}}_{i}(\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {e}}_{1})-h_{i+1,i}({\boldsymbol {e}}_{i}^{\mathrm {T} }{\boldsymbol {y}}_{i}){\boldsymbol {v}}_{i+1}\\&=-h_{i+1,i}({\boldsymbol {e}}_{i}^{\mathrm {T} }{\boldsymbol {y}}_{i}){\boldsymbol {v}}_{i+1}{\text{.}}\end{aligned}}$

Para ver la conjugación de ${\boldsymbol {p}}_{i}$, basta mostrar que ${\boldsymbol {P}}_{i}^{\mathrm {T} }{\boldsymbol {AP}}_{i}$ es diagonal:

${\begin{aligned}{\boldsymbol {P}}_{i}^{\mathrm {T} }{\boldsymbol {AP}}_{i}&={\boldsymbol {U}}_{i}^{-\mathrm {T} }{\boldsymbol {V}}_{i}^{\mathrm {T} }{\boldsymbol {AV}}_{i}{\boldsymbol {U}}_{i}^{-1}\\&={\boldsymbol {U}}_{i}^{-\mathrm {T} }{\boldsymbol {H}}_{i}{\boldsymbol {U}}_{i}^{-1}\\&={\boldsymbol {U}}_{i}^{-\mathrm {T} }{\boldsymbol {L}}_{i}{\boldsymbol {U}}_{i}{\boldsymbol {U}}_{i}^{-1}\\&={\boldsymbol {U}}_{i}^{-\mathrm {T} }{\boldsymbol {L}}_{i}\end{aligned}}$

es simétria y triangular inferior simultáneamente y por lo tanto debe ser diagonal.

Ahora se puede derivar los factores constantes $\alpha _{i}$ and $\beta _{i}$ con respecto al escalado ${\boldsymbol {p}}_{i}$ mediante imponer solamente la ortogonalidad de ${\boldsymbol {r}}_{i}$ y la conjugación de ${\boldsymbol {p}}_{i}$.

Debido a la ortogonalidad de ${\boldsymbol {r}}_{i}$, es necesario que ${\boldsymbol {r}}_{i+1}^{\mathrm {T} }{\boldsymbol {r}}_{i}=({\boldsymbol {r}}_{i}-\alpha _{i}{\boldsymbol {Ap}}_{i})^{\mathrm {T} }{\boldsymbol {r}}_{i}=0$. Como resultante de eso,

${\begin{aligned}\alpha _{i}&={\frac {{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{i}}{{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}\\&={\frac {{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{i}}{({\boldsymbol {p}}_{i}-\beta _{i-1}{\boldsymbol {p}}_{i-1})^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}\\&={\frac {{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{i}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}{\text{.}}\end{aligned}}$

De manera similar, debido a la conjugación de ${\boldsymbol {p}}_{i}$, es necesario que ${\boldsymbol {p}}_{i+1}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}=({\boldsymbol {r}}_{i+1}+\beta _{i}{\boldsymbol {p}}_{i})^{\mathrm {T} }{\boldsymbol {Ap}}_{i}=0$. Como resultante de eso,

${\begin{aligned}\beta _{i}&=-{\frac {{\boldsymbol {r}}_{i+1}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}\\&=-{\frac {{\boldsymbol {r}}_{i+1}^{\mathrm {T} }({\boldsymbol {r}}_{i}-{\boldsymbol {r}}_{i+1})}{\alpha _{i}{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}\\&={\frac {{\boldsymbol {r}}_{i+1}^{\mathrm {T} }{\boldsymbol {r}}_{i+1}}{{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{i}}}{\text{.}}\end{aligned}}$

Esto completa la derivación.

References

1. Hestenes, M. R.; Stiefel, E. (December de 1952). «Methods of conjugate gradients for solving linear systems» (PDF). Journal of Research of the National Bureau of Standards 49 (6). 
2. Saad, Y. (2003). «Chapter 6: Krylov Subspace Methods, Part I». Iterative methods for sparse linear systems (2nd edición). SIAM. ISBN 978-0898715347.