Derivación del método del gradiente conjugado

En álgebra lineal numérica , el método del gradiente conjugado es un método iterativo para resolver numéricamente el sistema lineal.

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

donde es simétrico positivo-definido , sin calcular explícitamente. El método del gradiente conjugado se puede derivar desde varias perspectivas diferentes, incluida la especialización del método de dirección conjugada ^[1] para la optimización y la variación de la iteración de Arnoldi / Lanczos para problemas de valores propios . ${\displaystyle {\boldsymbol {A}}}$ ${\displaystyle {\boldsymbol {A}}^{-1}}$

La intención de este artículo es documentar los pasos importantes en estas derivaciones.

El método del gradiente conjugado puede verse como un caso especial del método de dirección conjugada aplicado a la minimización de la función cuadrática.

${\displaystyle f({\boldsymbol {x}})={\boldsymbol {x}}^{\mathrm {T} }{\boldsymbol {A}}{\boldsymbol {x}}-2{\boldsymbol {b}}^{\mathrm {T} }{\boldsymbol {x}}{\text{.}}}$

lo que nos permite aplicar la intuición geométrica.

Geométricamente, la función cuadrática se puede representar de forma equivalente escribiendo su valor en cada punto del espacio. Los puntos de igual valor forman sus superficies de contorno, que son elipsoides concéntricos con la ecuación para la variación de . A medida que disminuye, los elipsoides se hacen cada vez más pequeños, hasta que en su valor mínimo, el elipsoide se encoge hasta su centro compartido. $${\displaystyle {\boldsymbol {x}}^{\mathrm {T}} {\boldsymbol {A}} {\boldsymbol {x}}-2{\boldsymbol {b}}^{\mathrm {T}} {\boldsymbol {x}}=C}$$${\estilo de visualización C}$${\estilo de visualización C}$

Minimizar la función cuadrática es entonces un problema de desplazamiento por el plano, buscando ese centro compartido de todos esos elipsoides. El centro se puede encontrar mediante un cálculo explícito, pero esto es precisamente lo que estamos tratando de evitar. ${\displaystyle {\boldsymbol {A}}^{-1}}$

El método más simple es la búsqueda de línea voraz , donde comenzamos en algún punto , elegimos una dirección de alguna manera y luego minimizamos . Esto tiene una solución simple de forma cerrada que no involucra inversión de matriz: Geométricamente, comenzamos en algún punto en algún elipsoide, luego elegimos una dirección y viajamos a lo largo de esa dirección, hasta que llegamos al punto donde el elipsoide se minimiza en esa dirección. Este no es necesariamente el mínimo, pero es un progreso hacia él. Visualmente, se mueve a lo largo de una línea y se detiene tan pronto como alcanzamos un punto tangente al elipsoide de contorno. ${\displaystyle {\boldsymbol {x}}_{0}}$${\displaystyle {\boldsymbol {p}}_{0}}$${\displaystyle f({\boldsymbol {x}}_{0}+{\boldsymbol {p}}_{0}\alpha _{0})}$$${\displaystyle \alpha _{0}={\frac {{\boldsymbol {p}}_{0}^{\mathrm {T} }({\boldsymbol {b}}-{\boldsymbol {Ax}}_{0})}{{\boldsymbol {p}}_{0}^{\mathrm {T} }{\boldsymbol {A}}{\boldsymbol {p}}_{0}}}}$$${\displaystyle {\boldsymbol {x}}_{0}}$

Ahora podemos repetir este procedimiento, comenzando en nuestro nuevo punto , elegir una nueva dirección , calcular , etc. ${\displaystyle {\boldsymbol {x}}_{1}={\boldsymbol {x}}_{0}+\alpha _{0}{\boldsymbol {p}}_{0}}$${\displaystyle {\boldsymbol {p}}_{1}}$${\displaystyle \alpha _{1}}$

Podemos resumir esto como el siguiente algoritmo:

Comience eligiendo una estimación inicial y calcule el residuo inicial , luego itere: ${\displaystyle {\boldsymbol {x}}_{0}}$${\displaystyle {\boldsymbol {r}}_{0}={\boldsymbol {b}}-{\boldsymbol {Ax}}_{0}}$

${\displaystyle {\begin{aligned}\alpha _{i}&={\frac {{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {r}}_{i}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}{\text{,}}\\{\boldsymbol {x}}_{i+1}&={\boldsymbol {x}}_{i}+\alpha _{i}{\boldsymbol {p}}_{i}{\text{,}}\\{\boldsymbol {r}}_{i+1}&={\boldsymbol {r}}_{i}-\alpha _{i}{\boldsymbol {Ap}}_{i}\end{aligned}}}$

donde se deben elegir. Observe en particular cómo el residuo se calcula iterativamente paso a paso, en lugar de hacerlo de nuevo cada vez: Es posible que sea cierto que prematuramente, lo que traería problemas numéricos. Sin embargo, para elecciones particulares de , esto no ocurrirá antes de la convergencia, como probaremos a continuación. ${\displaystyle {\boldsymbol {p}}_{0},{\boldsymbol {p}}_{1},{\boldsymbol {p}}_{2},\ldots }$$${\displaystyle {\boldsymbol {r}}_{i+1}={\boldsymbol {b}}-{\boldsymbol {Ax}}_{i+1}={\boldsymbol {b}}-{\boldsymbol {A}}({\boldsymbol {x}}_{i}+\alpha _{i}{\boldsymbol {p}}_{i})={\boldsymbol {r}}_{i}-\alpha _{i}{\boldsymbol {A}}{\boldsymbol {p}}_{i}}$$${\displaystyle \alpha _{i}=0}$${\displaystyle {\boldsymbol {p}}_{0},{\boldsymbol {p}}_{1},{\boldsymbol {p}}_{2},\ldots }$

Si las direcciones no se eligen bien, el progreso será lento. En particular, el método de descenso de gradiente sería lento. Esto se puede ver en el diagrama, donde la línea verde es el resultado de elegir siempre la dirección del gradiente local. Zigzaguea hacia el mínimo, pero se pasa repetidamente. Por el contrario, si elegimos que las direcciones sean un conjunto de direcciones mutuamente conjugadas , entonces no habrá ningún sobreimpulso y obtendremos el mínimo global después de los pasos, donde es el número de dimensiones. ${\displaystyle {\boldsymbol {p}}_{0},{\boldsymbol {p}}_{1},{\boldsymbol {p}}_{2},\ldots }$${\displaystyle n}$${\displaystyle n}$

El concepto de direcciones conjugadas proviene de la geometría clásica de la elipse. Para una elipse, dos semiejes con centro son mutuamente conjugados con respecto a la elipse si y solo si las líneas son paralelas al paralelogramo tangente que limita, como se muestra en la imagen. El concepto se generaliza a elipsoides n -dimensionales, donde n semiejes son mutuamente conjugados con respecto al elipsoide si y solo si cada eje es paralelo al paralelepípedo tangente que limita . En otras palabras, para cualquier , el plano tangente al elipsoide en es un hiperplano abarcado por los vectores , donde es el centro del elipsoide. ${\displaystyle t_{0}{\boldsymbol {p}}_{0},\dots ,t_{n-1}{\boldsymbol {p}}_{n-1}}$${\displaystyle i}$${\displaystyle {\boldsymbol {c}}+t_{i}{\boldsymbol {p}}_{i}}$${\displaystyle \{{\boldsymbol {p}}_{j}:j\neq i\}}$${\displaystyle {\boldsymbol {c}}}$

Tenga en cuenta que necesitamos escalar cada vector direccional mediante un escalar , de modo que caiga exactamente en el elipsoide. ${\displaystyle {\boldsymbol {p}}_{i}}$${\displaystyle t_{i}}$${\displaystyle {\boldsymbol {c}}+t_{i}{\boldsymbol {p}}_{i}}$

Dado un elipsoide con ecuación para alguna constante , podemos traducirlo de modo que su centro esté en el origen. Esto cambia la ecuación a para alguna otra constante . La condición de tangencia es entonces: es decir, para cualquier . $${\displaystyle {\boldsymbol {x}}^{\mathrm {T} }{\boldsymbol {A}}{\boldsymbol {x}}-2{\boldsymbol {b}}^{\mathrm {T} }{\boldsymbol {x}}=C}$$${\displaystyle C}$$${\displaystyle {\boldsymbol {x}}^{\mathrm {T} }{\boldsymbol {A}}{\boldsymbol {x}}=C'}$$${\displaystyle C'}$$${\displaystyle (t_{i}{\boldsymbol {p}}_{i}+{\boldsymbol {p}}_{j}dt_{j})^{\mathrm {T} }{\boldsymbol {A}}(t_{i}{\boldsymbol {p}}_{i}+{\boldsymbol {p}}_{j}dt_{j})=C'+O(dt_{j}^{2}),\quad \forall i\neq j}$$${\displaystyle {\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{j}=0}$${\displaystyle i\neq j}$

El método de dirección conjugada es impreciso en el sentido de que no se dan fórmulas para la selección de las direcciones . Las elecciones específicas conducen a varios métodos, incluidos el método del gradiente conjugado y la eliminación gaussiana . ${\displaystyle {\boldsymbol {p}}_{0},{\boldsymbol {p}}_{1},{\boldsymbol {p}}_{2},\ldots }$

Podemos tabular las ecuaciones que necesitamos poner a cero:

Esto se parece al problema de ortogonalización, que requiere para cualquier , y para cualquier . Por lo tanto, el problema de encontrar ejes conjugados es menos restringido que el problema de ortogonalización, por lo que el proceso de Gram-Schmidt funciona, con grados de libertad adicionales que podemos usar más adelante para elegir los que simplificarían el cálculo: ${\displaystyle {\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {p}}_{j}=0}$${\displaystyle i\neq j}$${\displaystyle {\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {p}}_{j}=1}$${\displaystyle i=j}$

• Establecer arbitrariamente .${\displaystyle {\boldsymbol {p}}_{0}}$
• Establezca arbitrariamente y luego modifíquelo a .${\displaystyle {\boldsymbol {p}}_{10}}$${\displaystyle {\boldsymbol {p}}_{1}={\boldsymbol {p}}_{10}-{\frac {{\boldsymbol {p}}_{0}^{\mathrm {T} }{\boldsymbol {Ap}}_{10}}{{\boldsymbol {p}}_{0}^{\mathrm {T} }{\boldsymbol {Ap}}_{0}}}{\boldsymbol {p}}_{0}}$
• Establezca arbitrariamente y luego modifíquelo a .${\displaystyle {\boldsymbol {p}}_{20}}$${\displaystyle {\boldsymbol {p}}_{2}={\boldsymbol {p}}_{20}-\sum _{i=0}^{1}{\frac {{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{20}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}{\boldsymbol {p}}_{i}}$
• ...
• Establezca arbitrariamente y luego modifíquelo a .${\displaystyle {\boldsymbol {p}}_{n-1,0}}$${\displaystyle {\boldsymbol {p}}_{n-1}={\boldsymbol {p}}_{n-1,0}-\sum _{i=0}^{n-2}{\frac {{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{n-1,0}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}{\boldsymbol {p}}_{i}}$

La opción más natural de es el gradiente. Es decir, . Dado que las direcciones conjugadas se pueden escalar por un valor distinto de cero, lo escalamos por para limpieza de notación, obteniendo Por lo tanto, tenemos . Sustituyéndolo, tenemos el algoritmo del gradiente conjugado: Proposición. Si en algún punto, , entonces el algoritmo ha convergido, es decir, . ${\displaystyle {\boldsymbol {p}}_{k,0}}$${\displaystyle {\boldsymbol {p}}_{k,0}=\nabla f({\boldsymbol {x}}_{k})}$${\displaystyle -1/2}$$${\displaystyle {\boldsymbol {p}}_{k,0}=\mathbf {r} _{k}=\mathbf {b} -\mathbf {Ax} _{k}}$$${\displaystyle {\boldsymbol {p}}_{k}={\boldsymbol {r}}_{k}-\sum _{i=0}^{k-1}{\frac {{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ar}}_{k}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}{\boldsymbol {p}}_{i}}$$${\displaystyle {\begin{aligned}&\mathbf {r} _{0}:=\mathbf {b} -\mathbf {Ax} _{0}\\&\mathbf {p} _{0}:=\mathbf {r} _{0}\\&k:=0\\&{\text{do while }}k<n\\&\qquad \alpha _{k}:={\frac {\mathbf {p} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{k}}}\\&\qquad \mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}\\&\qquad {\text{if }}|\alpha _{k}|{\text{ is sufficiently small, then exit loop}}\\&\qquad \mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}\\&\qquad \mathbf {p} _{k+1}:={\boldsymbol {r}}_{k+1}-\sum _{i=0}^{k}{\frac {{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ar}}_{k+1}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}{\boldsymbol {p}}_{i}\\&\qquad k:=k+1\\&{\text{return }}\mathbf {x} _{k+1}{\text{ as the result}}\end{aligned}}}$$${\displaystyle \alpha _{k}=0}$${\displaystyle \nabla f(\mathrm {x} _{k+1})=0}$

Demostración. Por construcción, significaría que , es decir, si damos un paso de gradiente conjugado, volvemos exactamente al punto en el que estábamos. Esto solo es posible si el gradiente local ya es cero. ${\displaystyle \mathbf {x} _{k+1}=\mathbf {x} _{k}}$

Este algoritmo se puede simplificar significativamente mediante algunos lemas, dando como resultado el algoritmo de gradiente conjugado.

Lema 1. y . ${\displaystyle \mathbf {p} _{i}^{T}\mathbf {r} _{j}=0,\;\forall i<j}$${\displaystyle \mathbf {r} _{i}^{T}\mathbf {r} _{j}=0,\;\forall i<j}$

Demostración. Por la construcción geométrica, el plano tangente al elipsoide en contiene cada uno de los vectores de dirección conjugados anteriores . Además, es perpendicular a la tangente, por lo tanto . La segunda ecuación es verdadera ya que por construcción, es una transformada lineal de . ${\displaystyle \mathbf {x} _{j}}$${\displaystyle \mathbf {p} _{0},\mathbf {p} _{1},\dots ,\mathbf {p} _{j-1}}$${\displaystyle \mathbf {r} _{j}}$${\displaystyle \mathbf {p} _{i}^{T}\mathbf {r} _{j}=0,\;\forall i<j}$${\displaystyle \mathbf {r} _{0},\mathbf {r} _{1},\dots ,\mathbf {r} _{j-1}}$${\displaystyle \mathbf {p} _{0},\mathbf {p} _{1},\dots ,\mathbf {p} _{j-1}}$

Lema 2. . ${\displaystyle \mathbf {p} _{k}^{T}\mathbf {r} _{k}=\mathbf {r} _{k}^{T}\mathbf {r} _{k}}$

Demostración. Por construcción, , ahora aplicamos el lema 1. ${\displaystyle \mathbf {p} _{k}:={\boldsymbol {r}}_{k}-\sum _{i=0}^{k-1}{\frac {{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ar}}_{k-1}}{{\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ap}}_{i}}}{\boldsymbol {p}}_{i}}$

Lema 3. . ${\displaystyle {\boldsymbol {p}}_{i}^{\mathrm {T} }{\boldsymbol {Ar}}_{k+1}={\begin{cases}0,\;i<k\\-{\boldsymbol {r}}_{k+1}^{T}{\boldsymbol {r}}_{k+1}/\alpha _{k},\;i=k\end{cases}}}$

Demostración. Por construcción, tenemos , por lo tanto Ahora apliquemos el lema 1. ${\displaystyle \mathbf {r} _{i+1}=\mathbf {r} _{i}-\alpha _{k}\mathbf {Ap} _{i}}$$${\displaystyle {\boldsymbol {r}}_{k+1}^{T}{\boldsymbol {A}}{\boldsymbol {p}}_{i}={\boldsymbol {r}}_{k+1}^{T}{\frac {{\boldsymbol {r}}_{i}-{\boldsymbol {r}}_{i+1}}{\alpha _{i}}}}$$

Sustituyendo los lemas 1 a 3, tenemos y , que es el algoritmo de gradiente conjugado adecuado. ${\displaystyle \alpha _{k}={\frac {\mathbf {r} _{k}^{\top }\mathbf {r} _{k}}{\mathbf {p} _{k}^{\top }\mathbf {A} \mathbf {p} _{k}}}}$${\displaystyle \mathbf {p} _{k+1}:={\boldsymbol {r}}_{k+1}+{\frac {\mathbf {r} _{k+1}^{\top }\mathbf {r} _{k+1}}{\mathbf {r} _{k}^{\top }\mathbf {r} _{k}}}\mathbf {p} _{k}}$

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

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

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

definiendo dónde ${\displaystyle {\boldsymbol {v}}_{i}={\boldsymbol {w}}_{i}/\lVert {\boldsymbol {w}}_{i}\rVert _{2}}$

${\displaystyle {\boldsymbol {v}}_{i}={\begin{cases}{\boldsymbol {r}}_{0}&{\text{if }}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{if }}i>1{\text{.}}\end{cases}}}$

En otras palabras, para , se encuentra mediante la ortogonalización de Gram-Schmidt contra seguida de normalización. ${\displaystyle i>1}$${\displaystyle {\boldsymbol {v}}_{i}}$ ${\displaystyle {\boldsymbol {Av}}_{i-1}}$${\displaystyle \{{\boldsymbol {v}}_{1},{\boldsymbol {v}}_{2},\ldots ,{\boldsymbol {v}}_{i-1}\}}$

Expresada en forma matricial, la iteración queda plasmada en la ecuación

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

dónde

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

con

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

Al aplicar la iteración de Arnoldi para resolver sistemas lineales, se comienza con , el residuo corresponde a una estimación inicial . Después de cada paso de la iteración, se calcula y la nueva iteración . ${\displaystyle {\boldsymbol {r}}_{0}={\boldsymbol {b}}-{\boldsymbol {Ax}}_{0}}$${\displaystyle {\boldsymbol {x}}_{0}}$${\displaystyle {\boldsymbol {y}}_{i}={\boldsymbol {H}}_{i}^{-1}(\lVert {\boldsymbol {r}}_{0}\rVert _{2}{\boldsymbol {e}}_{1})}$${\displaystyle {\boldsymbol {x}}_{i}={\boldsymbol {x}}_{0}+{\boldsymbol {V}}_{i}{\boldsymbol {y}}_{i}}$

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

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

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

Dado que es simétrico positivo-definido, también lo es . Por lo tanto, se puede factorizar LU sin pivotar parcialmente en ${\displaystyle {\boldsymbol {A}}}$${\displaystyle {\boldsymbol {H}}_{i}}$${\displaystyle {\boldsymbol {H}}_{i}}$

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

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

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

Reescribir como ${\displaystyle {\boldsymbol {x}}_{i}={\boldsymbol {x}}_{0}+{\boldsymbol {V}}_{i}{\boldsymbol {y}}_{i}}$

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

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

Ahora es importante 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}}}$

De hecho, también hay recurrencias breves para y :${\displaystyle {\boldsymbol {p}}_{i}}$${\displaystyle \zeta _{i}}$

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

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

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

Las relaciones anteriores conducen directamente al método directo de Lanczos, que resulta ser ligeramente más complejo.

Si permitimos escalar y compensar el escalamiento en el factor constante, potencialmente podemos tener recurrencias más simples de la forma: ${\displaystyle {\boldsymbol {p}}_{i}}$

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

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

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

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

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

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

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

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

Ahora podemos derivar los factores constantes y con respecto a la escala imponiendo únicamente la ortogonalidad de y la conjugación de . ${\displaystyle \alpha _{i}}$${\displaystyle \beta _{i}}$${\displaystyle {\boldsymbol {p}}_{i}}$${\displaystyle {\boldsymbol {r}}_{i}}$${\displaystyle {\boldsymbol {p}}_{i}}$

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

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

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

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

Esto completa la derivación.

1. Hestenes, MR ; Stiefel, E. (diciembre de 1952). "Métodos de gradientes conjugados para resolver sistemas lineales" (PDF) . Revista de investigación de la Oficina Nacional de Normas . 49 (6): 409. doi :10.6028/jres.049.044.
2. Shewchuk, Jonathan Richard. " Una introducción al método del gradiente conjugado sin el dolor agonizante ". (1994)
3. Saad, Y. (2003). "Capítulo 6: Métodos del subespacio de Krylov, Parte I". Métodos iterativos para sistemas lineales dispersos (2.ª ed.). SIAM. ISBN 978-0-89871-534-7.