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
![{\displaystyle {\boldsymbol {Ax}}={\boldsymbol {b}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/96510817e9971a34d21f8afb23728c61424ac4eb)
donde
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[editar]
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
![{\displaystyle f(\mathbf {x} )={\frac {1}{2}}\mathbf {x} ^{\mathrm {T} }\mathbf {Ax} -\mathbf {b} ^{\mathrm {T} }\mathbf {x} {\text{.}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/abe04bb8f458386f9efebdb84e0e1989a5b0c4da)
En el método de las direcciones conjugadas se elige una serie de direcciones
que son conjugadas mutuamente, en otras palabras,
![{\displaystyle \langle \mathbf {p} _{i},\mathbf {Ap} _{i}\rangle =0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/70ea06b997c0aa5430947958e3bd8d4009a24f84)
para
. La iteración se inicia con una solución aproximada
y el correspondiente residuo
y procede con los siguientes pasos para
:
![{\displaystyle \alpha _{i}={\frac {\langle \mathbf {p} _{i},\mathbf {r} _{i}\rangle }{\langle \mathbf {p} _{i},\mathbf {Ap} _{i}\rangle }}{\text{,}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/38fac264fe173b789397d1efd6ca5fc796844850)
![{\displaystyle \mathbf {x} _{i+1}=\mathbf {x} _{i}+\alpha _{i}\mathbf {p} _{i}{\text{,}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/778575bc1a2879114afce8ae9be087f284bdf02e)
![{\displaystyle \mathbf {r} _{i+1}=\mathbf {r} _{i}+\alpha _{i}\mathbf {Ap} _{i}{\text{.}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e76e071205d922aee933a295d395feffdf656809)
Derivación de la iteración de Arnoldi/Lanczos[editar]
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[editar]
En la iteración de Arnoldi, se comienza con un vector
y construye gradualmente una base ortonormal
del subespacio de Krylov
![{\displaystyle {\mathcal {K}}({\boldsymbol {A}},{\boldsymbol {r}}_{0})=\{{\boldsymbol {r}}_{0},{\boldsymbol {Ar}}_{0},{\boldsymbol {A}}^{2}{\boldsymbol {r}}_{0},\ldots \}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/dd0ef9aabd8596fa38f7bc511481f2c90b2cc46a)
mediante la definición de
donde
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/32af2b65fd897ef205d843f88125e976b44907b3)
En otras palabras, para
,
se obtiene a través de la ortogonalización de Gram-Schmidt de
contra
seguida por la normalización.
Expresada en la forma matricial, la iteración se capta por la ecuación
![{\displaystyle {\boldsymbol {AV}}_{i}={\boldsymbol {V}}_{i+1}{\boldsymbol {\tilde {H}}}_{i}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/739bf9c2b138690032c1ac33971a680bfb8fe871)
donde
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/72ba241fd60d691aaf88e67b5e540da03559a0f8)
with
![{\displaystyle 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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/be3a6d6d30e9a396affacdab8dc72a0e28cbd0ce)
Cuando se aplica la iteración de Arnoldi a la solución de los sistemas lineals, se comienza con
, el residuo correspondiente a la aproximación inicial
. Después de cada paso de la iteración, se computa
y la nueva aproximación
.
El método de Lanczos directo[editar]
Para el resto de la discusión, se supone que
es simétrica definida positiva. Con la simetría de
, la matriz superior de Hessenberg
se hace simétrica y por lo tanto tridiagonal. Entonces se puede representar más claramente por
![{\displaystyle {\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{.}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1e9ae9b35aaba10ea324a61dc352e70d7fbc2dd9)
Esto permite una corta recurrencia de tres términos para
en la iteración, y así la iteración de Arnoldi se reduce a la iteración de Lanczos.
Ya que
es simétrica definida positiva,
también lo es. Por lo tanto,
se puede factorizar LU sin pivoteo parcial en
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9a86a3b9357c13d8b7421b05f9df123cbb89e46c)
con las recurrencias convenientes para
y
:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/07951628d38963bd1a4fa7ac46c385e9ccfb8832)
Se varía
en
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a1529f16eb967c77762e5fb17a15bae9e4906b07)
con
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/61069303d8d7feed42eef0e4ee49c58fb9fd7480)
Ahora es crucial observar que
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/121ccf23e2d5533451f07dc3709d6e92b9d5451e)
De hecho, también existen recurrencias cortas para
y
:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/818ab7e0f45243fc0d07683bf57c4a422f88978a)
Con esta formulación, llegamos a una recurrencia simple para
:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fdcb61124601d34bdcae919e02e26194d8a667a5)
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[editar]
Si se permite que
se escale y compensa el escalamiento en el factor constante, potencialmente se puede obtener recurrencias más simples de la forma:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/20d5c13a764e3d33235e9706079d92eeed11aef8)
Como las premisas de la simplificación, ahora se deriva la ortogonalidad of
y la conjugación de
, es decir, para
,
![{\displaystyle {\begin{aligned}{\boldsymbol {r}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{j}&=0{\text{,}}\\{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{j}&=0{\text{.}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/76c9dc2e992bf67cf448f0ad62268dc695db24e4)
Los residuos son mutuamente ortogonales porque
es esencialmente un múltiplo de
, dado que para
,
, y para
,
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a5603d13f9ab579a6f664a0f24b1a462bf2ee7c0)
Para ver la conjugación de
, basta mostrar que
es diagonal:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b6206b08ddfe16a744153ede5f9b1d25138e147e)
es simétria y triangular inferior simultáneamente y por lo tanto debe ser diagonal.
Ahora se puede derivar los factores constantes
and
con respecto al escalado
mediante imponer solamente la ortogonalidad de
y la conjugación de
.
Debido a la ortogonalidad de
, es necesario que
. Como resultante de eso,
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4bdc0e311f58dd4f903926cf8454e974412e314a)
De manera similar, debido a la conjugación de
, es necesario que
. Como resultante de eso,
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2d487af1395b7014a5df5d79a47dfb0211792f4c)
Esto completa la derivación.
- 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).
- Saad, Y. (2003). «Chapter 6: Krylov Subspace Methods, Part I». Iterative methods for sparse linear systems (2nd edición). SIAM. ISBN 978-0898715347.