Método de gradiente conjugado

Algoritmo de optimización matemática

Una comparación de la convergencia del descenso de gradiente con el tamaño de paso óptimo (en verde) y el vector conjugado (en rojo) para minimizar una función cuadrática asociada con un sistema lineal dado. El gradiente conjugado, suponiendo aritmética exacta, converge como máximo en n pasos, donde n es el tamaño de la matriz del sistema (aquí n  = 2).

En matemáticas , el método del gradiente conjugado es un algoritmo para la solución numérica de sistemas particulares de ecuaciones lineales , es decir aquellos cuya matriz es positiva-semidefinida . El método del gradiente conjugado a menudo se implementa como un algoritmo iterativo , aplicable a sistemas dispersos que son demasiado grandes para ser manejados mediante una implementación directa u otros métodos directos como la descomposición de Cholesky . A menudo surgen grandes sistemas dispersos al resolver numéricamente ecuaciones diferenciales parciales o problemas de optimización.

El método del gradiente conjugado también se puede utilizar para resolver problemas de optimización sin restricciones , como la minimización de energía . Se atribuye comúnmente a Magnus Hestenes y Eduard Stiefel , ^{[1] [2]} quienes lo programaron en el Z4 , ^[3] y lo investigaron exhaustivamente. ^{[4] [5]}

El método del gradiente biconjugado proporciona una generalización a matrices no simétricas. Varios métodos de gradiente conjugado no lineal buscan mínimos de problemas de optimización no lineal.

Descripción del problema abordado por gradientes conjugados.

Supongamos que queremos resolver el sistema de ecuaciones lineales.

${\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} }$

para el vector , donde la matriz conocida es simétrica (es decir, A ^T = A ), definida positiva (es decir, x ^T Ax > 0 para todos los vectores distintos de cero en R ⁿ ) y real , y también se conoce. Denotamos la solución única de este sistema por . ${\displaystyle \mathbf {x} }$ ${\displaystyle n\times n}$ ${\displaystyle \mathbf {A} }$ ${\displaystyle \mathbf {x} }$ ${\displaystyle \mathbf {b} }$ ${\displaystyle \mathbf {x} _{*}}$

Derivación como método directo

El método del gradiente conjugado se puede derivar desde varias perspectivas diferentes, incluida la especialización del método de dirección conjugada para la optimización y la variación de la iteración de Arnoldi / Lanczos para problemas de valores propios . A pesar de las diferencias en sus enfoques, estas derivaciones comparten un tema común: demostrar la ortogonalidad de los residuos y la conjugación de las direcciones de búsqueda. Estas dos propiedades son cruciales para desarrollar la conocida formulación sucinta del método.

Decimos que dos vectores distintos de cero u y v son conjugados (con respecto a ) si ${\displaystyle \mathbf {A} }$

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

Como es simétrico y definido positivo, el lado izquierdo define un producto interno ${\displaystyle \mathbf {A} }$

${\displaystyle \mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} =\langle \mathbf {u} ,\mathbf {v} \rangle _{\mathbf {A} }:=\langle \mathbf {A} \mathbf {u} ,\mathbf {v} \rangle =\langle \mathbf {u} ,\mathbf {A} ^{\mathsf {T}}\mathbf {v} \rangle =\langle \mathbf {u} ,\mathbf {A} \mathbf {v} \rangle .}$

Dos vectores son conjugados si y sólo si son ortogonales con respecto a este producto interno. Ser conjugado es una relación simétrica: si es conjugado con , entonces es conjugado con . Suponer que ${\displaystyle \mathbf {u} }$ ${\displaystyle \mathbf {v} }$ ${\displaystyle \mathbf {v} }$ ${\displaystyle \mathbf {u} }$

${\displaystyle P=\{\mathbf {p} _{1},\dots ,\mathbf {p} _{n}\}}$

es un conjunto de vectores mutuamente conjugados con respecto a , es decir, para todos . Entonces forma una base para y podemos expresar la solución de en esta base: ${\displaystyle n}$ ${\displaystyle \mathbf {A} }$ ${\displaystyle \mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{j}=0}$ ${\displaystyle i\neq j}$ ${\displaystyle P}$ ${\displaystyle \mathbb {R} ^{n}}$ ${\displaystyle \mathbf {x} _{*}}$ ${\displaystyle \mathbf {Ax} =\mathbf {b} }$

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

Multiplicar por la izquierda el problema con el vector produce ${\displaystyle \mathbf {Ax} =\mathbf {b} }$ ${\displaystyle \mathbf {p} _{k}^{\mathsf {T}}}$

${\displaystyle \mathbf {p} _{k}^{\mathsf {T}}\mathbf {b} =\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{i}=\sum _{i=1}^{n}\alpha _{i}\left\langle \mathbf {p} _{k},\mathbf {p} _{i}\right\rangle _{\mathbf {A} }=\alpha _{k}\left\langle \mathbf {p} _{k},\mathbf {p} _{k}\right\rangle _{\mathbf {A} }}$

y entonces

${\displaystyle \alpha _{k}={\frac {\langle \mathbf {p} _{k},\mathbf {b} \rangle }{\langle \mathbf {p} _{k},\mathbf {p} _{k}\rangle _{\mathbf {A} }}}.}$

Esto proporciona el siguiente método ^[4] para resolver la ecuación Ax = b : encuentre una secuencia de direcciones conjugadas y luego calcule los coeficientes . ${\displaystyle n}$ ${\displaystyle \alpha _{k}}$

Como método iterativo

Si elegimos los vectores conjugados con cuidado, es posible que no los necesitemos todos para obtener una buena aproximación a la solución . Por tanto, queremos considerar el método del gradiente conjugado como un método iterativo. Esto también nos permite resolver aproximadamente sistemas donde n es tan grande que el método directo tomaría demasiado tiempo. ${\displaystyle \mathbf {p} _{k}}$ ${\displaystyle \mathbf {x} _{*}}$

Denotamos la estimación inicial para x _∗ por x ₀ (podemos asumir sin pérdida de generalidad que x ₀ = 0 ; de lo contrario, considere el sistema Az = b − Ax ₀ ). Comenzando con x ₀ buscamos la solución y en cada iteración necesitamos una métrica que nos diga si estamos más cerca de la solución x _∗ (que no conocemos). Esta métrica proviene del hecho de que la solución x _∗ es también el único minimizador de la siguiente función cuadrática

${\displaystyle f(\mathbf {x} )={\tfrac {1}{2}}\mathbf {x} ^{\mathsf {T}}\mathbf {A} \mathbf {x} -\mathbf {x} ^{\mathsf {T}}\mathbf {b} ,\qquad \mathbf {x} \in \mathbf {R} ^{n}\,.}$

La existencia de un minimizador único es evidente ya que su matriz de Hesse de segundas derivadas es simétrica positiva-definida.

${\displaystyle \mathbf {H} (f(\mathbf {x} ))=\mathbf {A} \,,}$

y que el minimizador (use D f ( x )=0) resuelve el problema inicial se desprende de su primera derivada

${\displaystyle \nabla f(\mathbf {x} )=\mathbf {A} \mathbf {x} -\mathbf {b} \,.}$

Esto sugiere tomar el primer vector de base p ₀ como el negativo del gradiente de f en x = x ₀ . El gradiente de f es igual a Ax − b . Comenzando con una suposición inicial x ₀ , esto significa que tomamos p ₀ = b − Ax ₀ . Los otros vectores en la base se conjugarán con el gradiente, de ahí el nombre método de gradiente conjugado . Tenga en cuenta que p ₀ es también el residual proporcionado por este paso inicial del algoritmo.

Sea r _k el residual en el k ésimo paso:

${\displaystyle \mathbf {r} _{k}=\mathbf {b} -\mathbf {Ax} _{k}.}$

Como se observó anteriormente, es el gradiente negativo de at , por lo que el método de descenso _de gradiente requeriría moverse en la dirección rk . Aquí, sin embargo, insistimos en que las direcciones deben estar conjugadas entre sí. Una forma práctica de hacer cumplir esto es exigir que la siguiente dirección de búsqueda se construya a partir del residual actual y de todas las direcciones de búsqueda anteriores. La restricción de conjugación es de tipo ortonormal y, por tanto, el algoritmo puede verse como un ejemplo de ortonormalización de Gram-Schmidt . Esto da la siguiente expresión: ${\displaystyle \mathbf {r} _{k}}$ ${\displaystyle f}$ ${\displaystyle \mathbf {x} _{k}}$ ${\displaystyle \mathbf {p} _{k}}$

${\displaystyle \mathbf {p} _{k}=\mathbf {r} _{k}-\sum _{i<k}{\frac {\mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {r} _{k}}{\mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{i}}}\mathbf {p} _{i}}$

(consulte la imagen en la parte superior del artículo para conocer el efecto de la restricción de conjugación en la convergencia). Siguiendo esta dirección, la siguiente ubicación óptima está dada por

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

con

${\displaystyle \alpha _{k}={\frac {\mathbf {p} _{k}^{\mathsf {T}}(\mathbf {b} -\mathbf {Ax} _{k})}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}={\frac {\mathbf {p} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}},}$

donde la última igualdad se deriva de la definición de . La expresión para se puede derivar si se sustituye la expresión para x _{k +1} en f y se minimiza con respecto a ${\displaystyle \mathbf {r} _{k}}$ ${\displaystyle \alpha _{k}}$ ${\displaystyle \alpha _{k}}$

${\displaystyle {\begin{aligned}f(\mathbf {x} _{k+1})&=f(\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k})=:g(\alpha _{k})\\g'(\alpha _{k})&{\overset {!}{=}}0\quad \Rightarrow \quad \alpha _{k}={\frac {\mathbf {p} _{k}^{\mathsf {T}}(\mathbf {b} -\mathbf {Ax} _{k})}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}\,.\end{aligned}}}$

El algoritmo resultante

El algoritmo anterior ofrece la explicación más sencilla del método del gradiente conjugado. Aparentemente, el algoritmo tal como se indica requiere el almacenamiento de todas las direcciones de búsqueda anteriores y los vectores de residuos, así como muchas multiplicaciones de matriz-vector y, por lo tanto, puede ser computacionalmente costoso. Sin embargo, un análisis más detallado del algoritmo muestra que es ortogonal a , es decir , para i ≠ j. Y es ortogonal a , es decir , para . Esto se puede considerar que a medida que avanza el algoritmo, abarca el mismo subespacio de Krylov . Donde forman la base ortogonal con respecto al producto interno estándar, y forman la base ortogonal con respecto al producto interno inducido por . Por lo tanto, puede considerarse como una proyección del subespacio de Krylov. ${\displaystyle \mathbf {r} _{i}}$ ${\displaystyle \mathbf {r} _{j}}$ ${\displaystyle \mathbf {r} _{i}^{\mathsf {T}}\mathbf {r} _{j}=0}$ ${\displaystyle \mathbf {p} _{i}}$ ${\displaystyle \mathbf {A} }$ ${\displaystyle \mathbf {p} _{j}}$ ${\displaystyle \mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{j}=0}$ ${\displaystyle i\neq j}$ ${\displaystyle \mathbf {p} _{i}}$ ${\displaystyle \mathbf {r} _{i}}$ ${\displaystyle \mathbf {r} _{i}}$ ${\displaystyle \mathbf {p} _{i}}$ ${\displaystyle \mathbf {A} }$ ${\displaystyle \mathbf {x} _{k}}$ ${\displaystyle \mathbf {x} }$

Es decir, si el método CG comienza con , entonces ^[6] ${\displaystyle \mathbf {x} _{0}=0}$

$${\displaystyle x_{k}=\mathrm {argmin} _{y\in \mathbb {R} ^{n}}{\left\{(x-y)^{\top }A(x-y)^{T}:y\in \operatorname {span} \left\{b,Ab,\ldots ,A^{k-1}b\right\}\right\}}}$$

0 ${\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} }$ ${\displaystyle \mathbf {A} }$ ${\displaystyle \mathbf {x} _{0}}$

${\displaystyle {\begin{aligned}&\mathbf {r} _{0}:=\mathbf {b} -\mathbf {Ax} _{0}\\&{\hbox{if }}\mathbf {r} _{0}{\text{ is sufficiently small, then return }}\mathbf {x} _{0}{\text{ as the result}}\\&\mathbf {p} _{0}:=\mathbf {r} _{0}\\&k:=0\\&{\text{repeat}}\\&\qquad \alpha _{k}:={\frac {\mathbf {r} _{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 \mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}\\&\qquad {\hbox{if }}\mathbf {r} _{k+1}{\text{ is sufficiently small, then exit loop}}\\&\qquad \beta _{k}:={\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}}\\&\qquad \mathbf {p} _{k+1}:=\mathbf {r} _{k+1}+\beta _{k}\mathbf {p} _{k}\\&\qquad k:=k+1\\&{\text{end repeat}}\\&{\text{return }}\mathbf {x} _{k+1}{\text{ as the result}}\end{aligned}}}$

Este es el algoritmo más utilizado. La misma fórmula para βk también se utiliza en el método de gradiente conjugado no lineal _de Fletcher-Reeves .

Se reinicia

Observamos que se calcula mediante el método de descenso de gradiente aplicado a . La configuración también haría que se calcule mediante el método de descenso de gradiente desde , es decir, se puede utilizar como una implementación simple de un reinicio de las iteraciones de gradiente conjugado. ^[4] Los reinicios podrían ralentizar la convergencia, pero pueden mejorar la estabilidad si el método del gradiente conjugado no se comporta correctamente, por ejemplo, debido a un error de redondeo . ${\displaystyle \mathbf {x} _{1}}$ ${\displaystyle \mathbf {x} _{0}}$ ${\displaystyle \beta _{k}=0}$ ${\displaystyle \mathbf {x} _{k+1}}$ ${\displaystyle \mathbf {x} _{k}}$

Cálculo residual explícito

Las fórmulas y , que son válidas en aritmética exacta, hacen que las fórmulas y sean matemáticamente equivalentes. El primero se utiliza en el algoritmo para evitar una multiplicación adicional porque el vector ya está calculado para evaluar . Este último puede ser más preciso, sustituyendo el cálculo explícito por el implícito por la recursividad sujeta a acumulación de error de redondeo , por lo que se recomienda para una evaluación ocasional. ^[7] ${\displaystyle \mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}}$ ${\displaystyle \mathbf {r} _{k}:=\mathbf {b} -\mathbf {Ax} _{k}}$ ${\displaystyle \mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}}$ ${\displaystyle \mathbf {r} _{k+1}:=\mathbf {b} -\mathbf {Ax} _{k+1}}$ ${\displaystyle \mathbf {A} }$ ${\displaystyle \mathbf {Ap} _{k}}$ ${\displaystyle \alpha _{k}}$ ${\displaystyle \mathbf {r} _{k+1}:=\mathbf {b} -\mathbf {Ax} _{k+1}}$

Normalmente se utiliza una norma del residual como criterio de parada. La norma del residual explícito proporciona un nivel garantizado de precisión tanto en aritmética exacta como en presencia de errores de redondeo , donde la convergencia naturalmente se estanca. Por el contrario, se sabe que el residual implícito sigue disminuyendo en amplitud muy por debajo del nivel de errores de redondeo y, por tanto, no puede utilizarse para determinar el estancamiento de la convergencia. ${\displaystyle \mathbf {r} _{k+1}:=\mathbf {b} -\mathbf {Ax} _{k+1}}$ ${\displaystyle \mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}}$

Cálculo de alfa y beta.

En el algoritmo, α _k se elige de manera que sea ortogonal a . El denominador se simplifica de ${\displaystyle \mathbf {r} _{k+1}}$ ${\displaystyle \mathbf {r} _{k}}$

${\displaystyle \alpha _{k}={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{k}}}}$

desde . El β _k se elige de manera que esté conjugado con . Inicialmente, β _k es ${\displaystyle \mathbf {r} _{k+1}=\mathbf {p} _{k+1}-\mathbf {\beta } _{k}\mathbf {p} _{k}}$ ${\displaystyle \mathbf {p} _{k+1}}$ ${\displaystyle \mathbf {p} _{k}}$

${\displaystyle \beta _{k}=-{\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}}$

usando

${\displaystyle \mathbf {r} _{k+1}=\mathbf {r} _{k}-\alpha _{k}\mathbf {A} \mathbf {p} _{k}}$

y equivalentemente

${\displaystyle \mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}(\mathbf {r} _{k}-\mathbf {r} _{k+1}),}$

el numerador de β _k se reescribe como

${\displaystyle \mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}\mathbf {r} _{k+1}^{\mathsf {T}}(\mathbf {r} _{k}-\mathbf {r} _{k+1})=-{\frac {1}{\alpha _{k}}}\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}$

porque y son ortogonales por diseño. El denominador se reescribe como ${\displaystyle \mathbf {r} _{k+1}}$ ${\displaystyle \mathbf {r} _{k}}$

${\displaystyle \mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}=(\mathbf {r} _{k}+\beta _{k-1}\mathbf {p} _{k-1})^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}\mathbf {r} _{k}^{\mathsf {T}}(\mathbf {r} _{k}-\mathbf {r} _{k+1})={\frac {1}{\alpha _{k}}}\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}$

usando que las direcciones de búsqueda p _k están conjugadas y nuevamente que los residuos son ortogonales. Esto da el β en el algoritmo después de cancelar α _k .

Código de ejemplo en Julia (lenguaje de programación)

""" gradiente_conjugado!(A, b, x)Devuelva la solución a `A * x = b` usando el método de gradiente conjugado."""función conjugate_gradient! (  A :: MatrizAbstract , b :: VectorAbstract , x :: VectorAbstract ; tol = eps ( eltipo ( b ))   ) # Inicializar vector residual residual = b - A * x       # Inicializar el vector de dirección de búsqueda dirección_búsqueda = residual   # Calcular la norma residual al cuadrado inicialnorma ( x ) = raíz cuadrada ( suma ( x .^ 2 ))   old_resid_norm = norma ( residual )   # Iterar hasta la convergencia mientras que old_resid_norm > tol    A_dirección_de_búsqueda = A * dirección_de_búsqueda     tamaño_paso = norma_resid_antigua ^ 2 / ( dirección_búsqueda ' * A_dirección_búsqueda )       # Actualizar solución @. x = x + tamaño_paso * dirección_búsqueda        # Actualizar residual @. residual = residual - tamaño_paso * Una_dirección_búsqueda        new_resid_norm = norma ( residual )    # Actualizar vector de dirección de búsqueda @. dirección_búsqueda = residual +      ( nueva_resid_norm / antigua_resid_norm ) ^ 2 * dirección_búsqueda     # Actualizar la norma residual al cuadrado para la próxima iteración norma_resid_antigua = norma_resid_nueva   fin volver x fin

Ejemplo numérico

Considere el sistema lineal Ax = b dado por

${\displaystyle \mathbf {A} \mathbf {x} ={\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}={\begin{bmatrix}1\\2\end{bmatrix}},}$

Realizaremos dos pasos del método del gradiente conjugado comenzando con la suposición inicial.

${\displaystyle \mathbf {x} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}}$

para encontrar una solución aproximada al sistema.

Solución

Como referencia, la solución exacta es

${\displaystyle \mathbf {x} ={\begin{bmatrix}{\frac {1}{11}}\\\\{\frac {7}{11}}\end{bmatrix}}\approx {\begin{bmatrix}0.0909\\\\0.6364\end{bmatrix}}}$

Nuestro primer paso es calcular el vector residual r ₀ asociado con x ₀ . Este residual se calcula a partir de la fórmula r ₀ = b - Ax ₀ , y en nuestro caso es igual a

${\displaystyle \mathbf {r} _{0}={\begin{bmatrix}1\\2\end{bmatrix}}-{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}2\\1\end{bmatrix}}={\begin{bmatrix}-8\\-3\end{bmatrix}}=\mathbf {p} _{0}.}$

Dado que esta es la primera iteración, usaremos el vector residual r ₀ como nuestra dirección de búsqueda inicial p ₀ ; el método de selección de p _k cambiará en iteraciones posteriores.

Ahora calculamos el escalar α ₀ usando la relación

${\displaystyle \alpha _{0}={\frac {\mathbf {r} _{0}^{\mathsf {T}}\mathbf {r} _{0}}{\mathbf {p} _{0}^{\mathsf {T}}\mathbf {Ap} _{0}}}={\frac {{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}{{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}}={\frac {73}{331}}\approx 0.2205}$

Ahora podemos calcular x ₁ usando la fórmula

${\displaystyle \mathbf {x} _{1}=\mathbf {x} _{0}+\alpha _{0}\mathbf {p} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}+{\frac {73}{331}}{\begin{bmatrix}-8\\-3\end{bmatrix}}\approx {\begin{bmatrix}0.2356\\0.3384\end{bmatrix}}.}$

Este resultado completa la primera iteración, siendo el resultado una solución aproximada "mejorada" del sistema, x ₁ . Ahora podemos continuar y calcular el siguiente vector residual r ₁ usando la fórmula

${\displaystyle \mathbf {r} _{1}=\mathbf {r} _{0}-\alpha _{0}\mathbf {A} \mathbf {p} _{0}={\begin{bmatrix}-8\\-3\end{bmatrix}}-{\frac {73}{331}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}\approx {\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}.}$

Nuestro siguiente paso en el proceso es calcular el escalar β ₀ que eventualmente se usará para determinar la siguiente dirección de búsqueda p ₁ .

${\displaystyle \beta _{0}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {r} _{0}^{\mathsf {T}}\mathbf {r} _{0}}}\approx {\frac {{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}}{{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}}=0.0088.}$

Ahora, usando este escalar β ₀ , podemos calcular la siguiente dirección de búsqueda p ₁ usando la relación

${\displaystyle \mathbf {p} _{1}=\mathbf {r} _{1}+\beta _{0}\mathbf {p} _{0}\approx {\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}+0.0088{\begin{bmatrix}-8\\-3\end{bmatrix}}={\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}.}$

Ahora calculamos el escalar α ₁ usando nuestro p ₁ recién adquirido usando el mismo método que el usado para α ₀ .

${\displaystyle \alpha _{1}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {p} _{1}^{\mathsf {T}}\mathbf {Ap} _{1}}}\approx {\frac {{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}}{{\begin{bmatrix}-0.3511&0.7229\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}}}=0.4122.}$

Finalmente, encontramos x ₂ usando el mismo método que se usó para encontrar x ₁ .

${\displaystyle \mathbf {x} _{2}=\mathbf {x} _{1}+\alpha _{1}\mathbf {p} _{1}\approx {\begin{bmatrix}0.2356\\0.3384\end{bmatrix}}+0.4122{\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}={\begin{bmatrix}0.0909\\0.6364\end{bmatrix}}.}$

El resultado, x ₂ , es una "mejor" aproximación a la solución del sistema que x ₁ y x ₀ . Si en este ejemplo se utilizara aritmética exacta en lugar de precisión limitada, entonces, en teoría, la solución exacta se habría alcanzado después de n = 2 iteraciones ( siendo n el orden del sistema).

Propiedades de convergencia

En teoría, el método del gradiente conjugado puede verse como un método directo, ya que, en ausencia de error de redondeo, produce la solución exacta después de un número finito de iteraciones, que no es mayor que el tamaño de la matriz. En la práctica, nunca se obtiene la solución exacta ya que el método del gradiente conjugado es inestable con respecto incluso a perturbaciones pequeñas; por ejemplo, la mayoría de las direcciones no son en la práctica conjugadas, debido a la naturaleza degenerativa de la generación de los subespacios de Krylov.

Como método iterativo , el método del gradiente conjugado mejora monótonamente (en la norma energética) las aproximaciones a la solución exacta y puede alcanzar la tolerancia requerida después de un número relativamente pequeño (en comparación con el tamaño del problema) de iteraciones. La mejora suele ser lineal y su velocidad está determinada por el número de condición de la matriz del sistema : cuanto mayor es, más lenta es la mejora. ^[8] ${\displaystyle \mathbf {x} _{k}}$ ${\displaystyle \kappa (A)}$ ${\displaystyle A}$ ${\displaystyle \kappa (A)}$

Si es grande, el preacondicionamiento se usa comúnmente para reemplazar el sistema original por uno que sea más pequeño que , consulte a continuación. ${\displaystyle \kappa (A)}$ ${\displaystyle \mathbf {Ax} -\mathbf {b} =0}$ ${\displaystyle \mathbf {M} ^{-1}(\mathbf {Ax} -\mathbf {b} )=0}$ ${\displaystyle \kappa (\mathbf {M} ^{-1}\mathbf {A} )}$ ${\displaystyle \kappa (\mathbf {A} )}$

Teorema de convergencia

Definir un subconjunto de polinomios como

${\displaystyle \Pi _{k}^{*}:=\left\lbrace \ p\in \Pi _{k}\ :\ p(0)=1\ \right\rbrace \,,}$

donde es el conjunto de polinomios de máximo grado . ${\displaystyle \Pi _{k}}$ ${\displaystyle k}$

Sean las aproximaciones iterativas de la solución exacta y defina los errores como . Ahora, la tasa de convergencia se puede aproximar como ^{[4] [9]} ${\displaystyle \left(\mathbf {x} _{k}\right)_{k}}$ ${\displaystyle \mathbf {x} _{*}}$ ${\displaystyle \mathbf {e} _{k}:=\mathbf {x} _{k}-\mathbf {x} _{*}}$

${\displaystyle {\begin{aligned}\left\|\mathbf {e} _{k}\right\|_{\mathbf {A} }&=\min _{p\in \Pi _{k}^{*}}\left\|p(\mathbf {A} )\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq \min _{p\in \Pi _{k}^{*}}\,\max _{\lambda \in \sigma (\mathbf {A} )}|p(\lambda )|\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq 2\left({\frac {{\sqrt {\kappa (\mathbf {A} )}}-1}{{\sqrt {\kappa (\mathbf {A} )}}+1}}\right)^{k}\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq 2\exp \left({\frac {-2k}{\sqrt {\kappa (\mathbf {A} )}}}\right)\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\,,\end{aligned}}}$

donde denota el espectro y denota el número de condición . ${\displaystyle \sigma (\mathbf {A} )}$ ${\displaystyle \kappa (\mathbf {A} )}$

Esto muestra que las iteraciones son suficientes para reducir el error en cualquier caso . ${\displaystyle k={\tfrac {1}{2}}{\sqrt {\kappa (\mathbf {A} )}}\log \left(\left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\varepsilon ^{-1}\right)}$ ${\displaystyle 2\varepsilon }$ ${\displaystyle \varepsilon >0}$

Tenga en cuenta que el límite importante cuando tiende a ${\displaystyle \kappa (\mathbf {A} )}$ ${\displaystyle \infty }$

${\displaystyle {\frac {{\sqrt {\kappa (\mathbf {A} )}}-1}{{\sqrt {\kappa (\mathbf {A} )}}+1}}\approx 1-{\frac {2}{\sqrt {\kappa (\mathbf {A} )}}}\quad {\text{for}}\quad \kappa (\mathbf {A} )\gg 1\,.}$

Este límite muestra una tasa de convergencia más rápida en comparación con los métodos iterativos de Jacobi o Gauss-Seidel , que escalan como . ${\displaystyle \approx 1-{\frac {2}{\kappa (\mathbf {A} )}}}$

No se supone ningún error de redondeo en el teorema de convergencia, pero el límite de convergencia es comúnmente válido en la práctica, como lo explica teóricamente [ ^5] Anne Greenbaum .

Convergencia práctica

Si se inicializa aleatoriamente, la primera etapa de iteraciones suele ser la más rápida, ya que el error se elimina dentro del subespacio de Krylov que inicialmente refleja un número de condición efectivo más pequeño. La segunda etapa de convergencia suele estar bien definida por la convergencia teórica ligada a , pero puede ser superlineal, dependiendo de la distribución del espectro de la matriz y la distribución espectral del error. ^[5] En la última etapa, se alcanza la precisión más pequeña posible y la convergencia se detiene o el método puede incluso comenzar a divergir. En aplicaciones informáticas científicas típicas en formato de punto flotante de doble precisión para matrices de grandes tamaños, el método del gradiente conjugado utiliza un criterio de parada con una tolerancia que termina las iteraciones durante la primera o segunda etapa. ${\textstyle {\sqrt {\kappa (\mathbf {A} )}}}$ ${\displaystyle A}$

El método del gradiente conjugado precondicionado

En la mayoría de los casos, el preacondicionamiento es necesario para asegurar una rápida convergencia del método del gradiente conjugado. Si es simétrico positivo-definido y tiene un número de condición mejor que , se puede utilizar un método de gradiente conjugado precondicionado. Toma la siguiente forma: ^[10] ${\displaystyle \mathbf {M} ^{-1}}$ ${\displaystyle \mathbf {M} ^{-1}\mathbf {A} }$ ${\displaystyle \mathbf {A} }$

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

${\displaystyle {\textrm {Solve:}}\mathbf {M} \mathbf {z} _{0}:=\mathbf {r} _{0}}$

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

${\displaystyle k:=0\,}$

repetir

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

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

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

si r _{k +1} es suficientemente pequeño , entonces salga del final del bucle si

${\displaystyle \mathrm {Solve} \ \mathbf {M} \mathbf {z} _{k+1}:=\mathbf {r} _{k+1}}$

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

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

${\displaystyle k:=k+1\,}$

finalizar repetir

El resultado es x _{k +1}

La formulación anterior equivale a aplicar el método del gradiente conjugado regular al sistema precondicionado ^[11]

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

dónde

${\displaystyle \mathbf {EE} ^{\mathsf {T}}=\mathbf {M} ,\qquad \mathbf {\hat {x}} =\mathbf {E} ^{\mathsf {T}}\mathbf {x} .}$

Se debe utilizar la descomposición de Cholesky del precondicionador para mantener la simetría (y la precisión positiva) del sistema. Sin embargo, no es necesario calcular esta descomposición y basta con saberlo . Se puede demostrar que tiene el mismo espectro que . ${\displaystyle \mathbf {M} ^{-1}}$ ${\displaystyle \mathbf {E} ^{-1}\mathbf {A} (\mathbf {E} ^{-1})^{\mathsf {T}}}$ ${\displaystyle \mathbf {M} ^{-1}\mathbf {A} }$

La matriz precondicionadora M tiene que ser simétrica positiva-definida y fija, es decir, no puede cambiar de una iteración a otra. Si se viola cualquiera de estos supuestos sobre el precondicionador, el comportamiento del método del gradiente conjugado precondicionado puede volverse impredecible.

Un ejemplo de precondicionador de uso común es la factorización incompleta de Cholesky . ^[12]

Uso del preacondicionador en la práctica

Es importante tener en cuenta que no queremos invertir la matriz explícitamente para poder usarla en el proceso, ya que invertir tomaría más tiempo/recursos computacionales que resolver el algoritmo de gradiente conjugado en sí. Como ejemplo, digamos que estamos usando un precondicionador proveniente de una factorización de Cholesky incompleta. La matriz resultante es la matriz triangular inferior y la matriz precondicionadora es: ${\displaystyle \mathbf {M} }$ ${\displaystyle \mathbf {M} ^{-1}}$ ${\displaystyle \mathbf {M} }$ ${\displaystyle \mathbf {L} }$

${\displaystyle \mathbf {M} =\mathbf {LL} ^{\mathsf {T}}}$

Entonces tenemos que resolver:

${\displaystyle \mathbf {Mz} =\mathbf {r} }$

${\displaystyle \mathbf {z} =\mathbf {M} ^{-1}\mathbf {r} }$

Pero:

${\displaystyle \mathbf {M} ^{-1}=(\mathbf {L} ^{-1})^{\mathsf {T}}\mathbf {L} ^{-1}}$

Entonces:

${\displaystyle \mathbf {z} =(\mathbf {L} ^{-1})^{\mathsf {T}}\mathbf {L} ^{-1}\mathbf {r} }$

Tomemos un vector intermediario : ${\displaystyle \mathbf {a} }$

${\displaystyle \mathbf {a} =\mathbf {L} ^{-1}\mathbf {r} }$

${\displaystyle \mathbf {r} =\mathbf {L} \mathbf {a} }$

Dado que y y conocido, y es triangular inferior, resolverlo es fácil y computacionalmente barato mediante el uso de sustitución directa . Luego, sustituimos en la ecuación original: ${\displaystyle \mathbf {r} }$ ${\displaystyle \mathbf {L} }$ ${\displaystyle \mathbf {L} }$ ${\displaystyle \mathbf {a} }$ ${\displaystyle \mathbf {a} }$

${\displaystyle \mathbf {z} =(\mathbf {L} ^{-1})^{\mathsf {T}}\mathbf {a} }$

${\displaystyle \mathbf {a} =\mathbf {L} ^{\mathsf {T}}\mathbf {z} }$

Dado que y son conocidos, y es triangular superior, resolverlo es fácil y computacionalmente barato mediante el uso de sustitución hacia atrás . ${\displaystyle \mathbf {a} }$ ${\displaystyle \mathbf {L} ^{\mathsf {T}}}$ ${\displaystyle \mathbf {L} ^{\mathsf {T}}}$ ${\displaystyle \mathbf {z} }$

Con este método, no es necesario invertirlo ni hacerlo explícitamente, y aun así obtenemos . ${\displaystyle \mathbf {M} }$ ${\displaystyle \mathbf {L} }$ ${\displaystyle \mathbf {z} }$

El método del gradiente conjugado precondicionado flexible

En aplicaciones numéricamente desafiantes, se utilizan precondicionadores sofisticados, que pueden conducir a un precondicionamiento variable, que cambia entre iteraciones. Incluso si el precondicionador es simétrico positivo-definido en cada iteración, el hecho de que pueda cambiar hace que los argumentos anteriores sean inválidos y, en pruebas prácticas, conduce a una desaceleración significativa de la convergencia del algoritmo presentado anteriormente. Usando la fórmula de Polak-Ribière

${\displaystyle \beta _{k}:={\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\left(\mathbf {z} _{k+1}-\mathbf {z} _{k}\right)}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {z} _{k}}}}$

en lugar de la fórmula de Fletcher-Reeves

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

puede mejorar dramáticamente la convergencia en este caso. ^[13] Esta versión del método de gradiente conjugado precondicionado puede denominarse ^[14] flexible, ya que permite un precondicionamiento variable. También se ha demostrado que la versión flexible ^[15] es robusta incluso si el preacondicionador no es simétrico positivo definido (SPD).

La implementación de la versión flexible requiere almacenar un vector adicional. Para un precondicionador SPD fijo, ambas fórmulas para βk son equivalentes en aritmética exacta, es decir, sin el error _de redondeo . ${\displaystyle \mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {z} _{k}=0,}$

La explicación matemática del mejor comportamiento de convergencia del método con la fórmula de Polak-Ribière es que el método es localmente óptimo en este caso, en particular, no converge más lento que el método de descenso más pronunciado localmente óptimo. ^[dieciséis]

vs. el método de descenso más empinado localmente óptimo

Tanto en el método de gradiente conjugado original como en el precondicionado , solo es necesario configurarlo para que sean localmente óptimos, utilizando los métodos de descenso más pronunciado de búsqueda de líneas . Con esta sustitución, los vectores p son siempre los mismos que los vectores z , por lo que no es necesario almacenar los vectores p . Por lo tanto, cada iteración de estos métodos de descenso más pronunciado es un poco más barata en comparación con la de los métodos de gradiente conjugado. Sin embargo, estos últimos convergen más rápido, a menos que se utilice un precondicionador (altamente) variable y/o no SPD , ver arriba. ${\displaystyle \beta _{k}:=0}$

Método de gradiente conjugado como controlador de retroalimentación óptimo para doble integrador

El método del gradiente conjugado también se puede derivar utilizando la teoría de control óptimo . ^[17] En este enfoque, el método del gradiente conjugado resulta ser un controlador de retroalimentación óptimo .

$${\displaystyle u=k(x,v):=-\gamma _{a}\nabla f(x)-\gamma _{b}v}$$

sistema doble integrador

$${\displaystyle {\dot {x}}=v,\quad {\dot {v}}=u}$$

^[17] ${\displaystyle \gamma _{a}}$ ${\displaystyle \gamma _{b}}$

gradiente conjugado en las ecuaciones normales

El método del gradiente conjugado se puede aplicar a una matriz arbitraria de n por m aplicándola a las ecuaciones normales A ^T A y al vector del lado derecho A ^T b , ya que A ^T A es una matriz semidefinida positiva simétrica para cualquier A . El resultado es un gradiente conjugado en las ecuaciones normales ( CGN o CGNR ).

A ^T Ax = A ^T b

Como método iterativo, no es necesario formar A ^T A explícitamente en la memoria, sino solo realizar las multiplicaciones matriz-vector y transponer matriz-vector. Por lo tanto, CGNR es particularmente útil cuando A es una matriz dispersa ya que estas operaciones suelen ser extremadamente eficientes. Sin embargo, la desventaja de formar ecuaciones normales es que el número de condición κ( A ^T A ) es igual a κ ² ( A ), por lo que la tasa de convergencia de CGNR puede ser lenta y la calidad de la solución aproximada puede ser sensible al redondeo. errores. Encontrar un buen preacondicionador suele ser una parte importante del uso del método CGNR.

Se han propuesto varios algoritmos (por ejemplo, CGLS, LSQR). El algoritmo LSQR supuestamente tiene la mejor estabilidad numérica cuando A está mal condicionado, es decir, A tiene un número de condición grande .

Método de gradiente conjugado para matrices hermitianas complejas

El método del gradiente conjugado con una modificación trivial se puede ampliar para resolver, dada la matriz A de valores complejos y el vector b, el sistema de ecuaciones lineales para el vector x de valores complejos, donde A es hermitiano (es decir, A' = A) y positivo -matriz definida , y el símbolo 'denota la transpuesta conjugada . La modificación trivial es simplemente sustituir la transpuesta conjugada por la transpuesta real en todas partes. ${\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} }$

Ventajas y desventajas

Las ventajas y desventajas de los métodos de gradiente conjugado se resumen en las notas de clase de Nemirovsky y BenTal. ^{[18] : Sección 7.3}

Un ejemplo patológico

Este ejemplo es de ^[19] Let y define ${\textstyle t\in (0,1)}$

$${\displaystyle W={\begin{bmatrix}t&{\sqrt {t}}&&&&\\{\sqrt {t}}&1+t&{\sqrt {t}}&&&\\&{\sqrt {t}}&1+t&{\sqrt {t}}&&\\&&{\sqrt {t}}&\ddots &\ddots &\\&&&\ddots &&\\&&&&&{\sqrt {t}}\\&&&&{\sqrt {t}}&1+t\end{bmatrix}},\quad b={\begin{bmatrix}1\\0\\\vdots \\0\end{bmatrix}}}$$

${\displaystyle W}$ ${\textstyle Wx=b}$

$${\displaystyle \|b-Wx_{k}\|^{2}=(1/t)^{k},\quad \|b-Wx_{n}\|^{2}=0}$$

Ver también

• Método de gradiente biconjugado (BiCG)
• Método de gradiente cuadrado conjugado (CGS)
• Método residual conjugado
• Propagación de la creencia gaussiana
• Método iterativo: sistemas lineales
• subespacio de Krylov
• Método de gradiente conjugado no lineal
• Preacondicionamiento
• Multiplicación de vectores y matrices dispersas

Referencias

1. Hestenes, Magnus R .; Stiefel, Eduard (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. Straeter, TA (1971). Sobre la extensión de la clase Davidon-Broyden de rango uno, métodos de minimización cuasi-Newton a un espacio de Hilbert de dimensión infinita con aplicaciones a problemas de control óptimo (tesis doctoral). Universidad Estatal de Carolina del Norte. hdl : 2060/19710026200 – a través del servidor de informes técnicos de la NASA.
3. Speiser, Ambros (2004). "Konrad Zuse und die ERMETH: Ein weltweiter Architektur-Vergleich" [Konrad Zuse y ERMETH: una comparación mundial de arquitecturas]. En Hellige, Hans Dieter (ed.). Geschichten der Informatik. Visionen, Paradigmen, Leitmotive (en alemán). Berlín: Springer. pag. 185.ISBN​ 3-540-00217-0.
4. Polyak, Boris (1987). Introducción a la Optimización.
5. Greenbaum, Anne (1997). Métodos iterativos para resolver sistemas lineales . doi :10.1137/1.9781611970937. ISBN 978-0-89871-396-1.
6. Paqueta, Elliot; Trogdon, Thomas (marzo de 2023). "Universalidad de los algoritmos de gradiente conjugado y MINRES en matrices de covarianza de muestra". Comunicaciones sobre Matemática Pura y Aplicada . 76 (5): 1085-1136. arXiv : 2007.00640 . doi :10.1002/cpa.22081. ISSN  0010-3640.
7. Shewchuk, Jonathan R (1994). Introducción al método del gradiente conjugado sin dolor agonizante (PDF) .
8. Saad, Yousef (2003). Métodos iterativos para sistemas lineales dispersos (2ª ed.). Filadelfia, Pensilvania: Sociedad de Matemáticas Industriales y Aplicadas. págs.195. ISBN 978-0-89871-534-7.
9. Hackbusch, W. (21 de junio de 2016). Solución iterativa de grandes sistemas dispersos de ecuaciones (2ª ed.). Suiza: Springer. ISBN 978-3-319-28483-5. OCLC  952572240.
10. Barrett, Richard; Baya, Michael; Chan, Tony F.; Demmel, James; Donato, junio; Dongarra, Jack; Eijkhout, Víctor; Pozo, Roldán; Romina, Carlos; van der Vorst, Henk. Plantillas para la solución de sistemas lineales: componentes básicos de métodos iterativos (PDF) (2ª ed.). Filadelfia, PA: SIAM. pag. 13 . Consultado el 31 de marzo de 2020 .
11. Golub, gen H.; Préstamo de Van, Charles F. (2013). Cálculos matriciales (4ª ed.). Prensa de la Universidad Johns Hopkins. segundo. 11.5.2. ISBN 978-1-4214-0794-4.
12. Conco, P.; Golub, GH; Meurant, G. (1985). "Precondicionamiento de bloques para el método de gradiente conjugado". Revista SIAM de Computación Científica y Estadística . 6 (1): 220–252. doi :10.1137/0906018.
13. Golub, gen H.; Sí, Qiang (1999). "Método de gradiente conjugado precondicionado inexacto con iteración interior-exterior". Revista SIAM de Computación Científica . 21 (4): 1305. CiteSeerX 10.1.1.56.1755 . doi :10.1137/S1064827597323415. 
14. Notay, Yvan (2000). "Gradientes conjugados flexibles". Revista SIAM de Computación Científica . 22 (4): 1444-1460. CiteSeerX 10.1.1.35.7473 . doi :10.1137/S1064827599362314. 
15. Bouwmeester, Henricus; Dougherty, Andrés; Knyazev, Andrés V. (2015). "Precondicionamiento no simétrico para métodos de gradiente conjugado y descenso más pronunciado 1". Procedia Ciencias de la Computación . 51 : 276–285. arXiv : 1212.6680 . doi : 10.1016/j.procs.2015.05.241 . S2CID  51978658.
16. Knyazev, Andrés V.; Lashuk, Ilya (2008). "Métodos de gradiente conjugado y descenso más pronunciado con preacondicionamiento variable". Revista SIAM sobre Análisis y Aplicaciones de Matrices . 29 (4): 1267. arXiv : matemáticas/0605767 . doi :10.1137/060675290. S2CID  17614913.
17. Ross, IM , "Una teoría de control óptimo para la optimización acelerada", arXiv :1902.09004, 2019.
18. Nemirovsky y Ben-Tal (2023). "Optimización III: Optimización convexa" (PDF) .
19. Pennington, Fabián Pedregosa, Courtney Paquette, Tom Trogdon, Jeffrey. "Tutorial de teoría de matrices aleatorias y aprendizaje automático". aprendizaje-matrix-aleatorio.github.io . Consultado el 5 de diciembre de 2023 .{{cite web}}: CS1 maint: multiple names: authors list (link)

Otras lecturas

• Atkinson, Kendell A. (1988). "Sección 8.9". Una introducción al análisis numérico (2ª ed.). John Wiley e hijos. ISBN 978-0-471-50023-0.
• Avriel, Mardoqueo (2003). Programación no lineal: análisis y métodos . Publicación de Dover. ISBN 978-0-486-43227-4.
• Golub, Gene H.; Préstamo de Van, Charles F. (2013). "Capítulo 11". Cálculos matriciales (4ª ed.). Prensa de la Universidad Johns Hopkins. ISBN 978-1-4214-0794-4.
• Saad, Yousef (1 de abril de 2003). "Capítulo 6" . Métodos iterativos para sistemas lineales dispersos (2ª ed.). SIAM. ISBN 978-0-89871-534-7.
• Gérard Meurant: "Detección y corrección de errores silenciosos en el algoritmo de gradiente conjugado", Algoritmos numéricos, vol.92 (2023), pp.869-891. URL=https://doi.org/10.1007/s11075-022-01380-1
• Meurant, Gerard; Tichy, Petr (2024). "Estimación de la norma de error en el algoritmo de gradiente conjugado" . SIAM. ISBN 978-1-61197-785-1.

enlaces externos

• "Gradientes conjugados, método de", Enciclopedia de Matemáticas , EMS Press , 2001 [1994]