Método del gradiente conjugado

Algoritmo de optimización matemática

Comparación de la convergencia del descenso de gradiente con tamaño de paso óptimo (en verde) y vector conjugado (en rojo) para minimizar una función cuadrática asociada con un sistema lineal dado. El gradiente conjugado, suponiendo una aritmética exacta, converge en un máximo de 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 se implementa a menudo como un algoritmo iterativo , aplicable a sistemas dispersos que son demasiado grandes para ser manejados por una implementación directa u otros métodos directos como la descomposición de Cholesky . Los sistemas dispersos grandes a menudo surgen 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 para 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} _{*}}$

La 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: probar 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 respecto de este producto interno. Ser conjugado es una relación simétrica: si es conjugado a , entonces es conjugado a . Supóngase 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 todo . 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}.}$

Multiplicando por la izquierda el problema con el vector obtenemos ${\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 : encontrar una secuencia de direcciones conjugadas y luego calcular 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 lo 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 suponer sin pérdida de generalidad que x ₀ = 0 , de lo contrario, consideremos el sistema Az = b − Ax ₀ en su lugar). 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 nos es desconocida). Esta métrica proviene del hecho de que la solución x _∗ es también el minimizador único 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 hessiana de derivadas segundas es simétrica positiva-definida.

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

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

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

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

Sea r _k el residuo en el k -ésimo paso:

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

Como se observó anteriormente, es el gradiente negativo de en , por lo que el método de descenso de gradiente requeriría moverse en la dirección r _k . Aquí, sin embargo, insistimos en que las direcciones deben ser conjugadas entre sí. Una forma práctica de hacer cumplir esto es requerir que la siguiente dirección de búsqueda se construya a partir del residuo actual y todas las direcciones de búsqueda anteriores. La restricción de conjugación es una restricción de tipo ortonormal y, por lo 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}}$

(ver la imagen en la parte superior del artículo para ver 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 sigue 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 proporciona 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 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, y abarcan 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, se puede considerar como la proyección de en el 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] A continuación se detalla el algoritmo para resolver donde es una matriz real, simétrica y definida positiva. El vector de entrada puede ser una solución inicial aproximada o 0 . Es una formulación diferente del procedimiento exacto descrito anteriormente. ${\displaystyle \mathbf {x} _{0}=0}$ $${\displaystyle x_{k}=\mathrm {argmin} _{y\in \mathbb {R} ^{n}}{\left\{(x-y)^{\top }A(x-y):y\in \operatorname {span} \left\{b,Ab,\ldots ,A^{k-1}b\right\}\right\}}}$$ ${\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 .

Reinicios

Observamos que se calcula mediante el método de descenso de gradiente aplicado a . La configuración haría que se calcule mediante el método de descenso de gradiente de manera similar a partir de , 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 de gradiente conjugado se comporta mal, 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 de residuos explícitos

Las fórmulas y , que se cumplen en aritmética exacta, hacen que las fórmulas y sean matemáticamente equivalentes. La primera se utiliza en el algoritmo para evitar una multiplicación adicional por ya que el vector ya está calculado para evaluar . La última puede ser más precisa, sustituyendo el cálculo explícito por el implícito por la recursión sujeta a la acumulación de error de redondeo , y por lo tanto 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 residuo como criterio de parada. La norma del residuo explícito proporciona un nivel garantizado de precisión tanto en aritmética exacta como en presencia de errores de redondeo , donde la convergencia se estanca naturalmente. Por el contrario, se sabe que el residuo implícito se vuelve cada vez más pequeño en amplitud muy por debajo del nivel de errores de redondeo y, por lo tanto, no se puede utilizar 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, se elige α _k de modo que sea ortogonal a . El denominador se simplifica a partir 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}}}}$

ya que . El β _k se elige de manera que sea conjugado a . 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}}$

Utilizando 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 enJulia (lenguaje de programación)

""" gradiente_conjugado!(A, b, x)Devuelve la solución a `A * x = b` usando el método de gradiente conjugado."""función conjugate_gradient! (  A :: MatrizAbstracta , b :: VectorAbstracto , x :: VectorAbstracto ; tol = eps ( eltype ( b ))   ) # Inicializar vector residual residuo = b - A * x       # Inicializar el vector de dirección de búsqueda direccion_busqueda = copiar ( residual )   # Calcular la norma residual cuadrática inicialnorma ( x ) = sqrt ( suma ( x .^ 2 ))   old_resid_norm = norma ( residual )   # Iterar hasta la convergencia mientras que old_resid_norm > tol    A_direccion_de_busqueda = A * direccion_de_busqueda     tamaño_de_paso = antigua_norma_resid ^ 2 / ( dirección_de_búsqueda ' * A_dirección_de_búsqueda )       #Actualizar solución @. x = x + tamaño_de_paso * dirección_de_búsqueda        # Actualizar residual @. residual = residual - tamaño_de_paso * A_dirección_de_búsqueda        new_resid_norm = norma ( residual )    # Actualizar vector de dirección de búsqueda @ .direccion_de_busqueda = residual +      ( nueva_norma_residual / antigua_norma_residual ) ^ 2 * dirección_de_búsqueda     # Actualizar la norma residual al cuadrado para la siguiente iteración antigua_norma_resid = nueva_norma_resid   fin devolver x fin

Ejemplo numérico

Consideremos 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 estimación inicial

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

con el fin de 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 a 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, utilizaremos el vector residual r ₀ como nuestra dirección de búsqueda inicial p ₀ ; el método de selección 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, y el resultado es una solución aproximada "mejorada" del sistema, x ₁ . Ahora podemos continuar y calcular el siguiente vector residual r ₁ utilizando 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 utilizará 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, utilizando este escalar β ₀ , podemos calcular la siguiente dirección de búsqueda p ₁ utilizando 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 α ₁ utilizando nuestro p ₁ recién adquirido y utilizando el mismo método que utilizamos 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 usamos 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 aproximación "mejor" 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 teóricamente se habría alcanzado la solución exacta después de n = 2 iteraciones ( siendo n el orden del sistema).

Propiedades de convergencia

El método del gradiente conjugado puede considerarse teóricamente 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, la solución exacta nunca se obtiene ya que el método del gradiente conjugado es inestable con respecto incluso a pequeñas perturbaciones, por ejemplo, la mayoría de las direcciones no son conjugadas en la práctica, debido a una naturaleza degenerativa de la generación de los subespacios de Krylov.

Como método iterativo , el método del gradiente conjugado mejora de manera monótona (en la norma de energía) 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 es típicamente lineal y su velocidad está determinada por el número de condición de la matriz del sistema : cuanto mayor sea, más lenta será la mejora. ^[8] ${\displaystyle \mathbf {x} _{k}}$ ${\displaystyle \kappa (A)}$ ${\displaystyle A}$ ${\displaystyle \kappa (A)}$

Si es grande, el preacondicionamiento se utiliza comúnmente para reemplazar el sistema original por uno que sea más pequeño que , vea 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 grado máximo . ${\displaystyle \Pi _{k}}$ ${\displaystyle k}$

Sean las aproximaciones iterativas de la solución exacta y definamos 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 demuestra que las iteraciones son suficientes para reducir el error a cualquier . ${\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}$

Nótese que el límite importante 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} )}}}$

En el teorema de convergencia no se supone ningún error de redondeo, pero el límite de convergencia es comúnmente válido en la práctica, como lo explicó 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 efectiva más pequeño. La segunda etapa de convergencia suele estar bien definida por el límite de convergencia teórico con , pero puede ser superlineal, dependiendo de una 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 alcanzable y la convergencia se detiene o el método puede incluso comenzar a divergir. En aplicaciones típicas de computación científica en formato de punto flotante de doble precisión para matrices de grandes tamaños, el método de gradiente conjugado utiliza un criterio de detención 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 preacondicionado

En la mayoría de los casos, es necesario un preacondicionamiento para garantizar una convergencia rápida 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 del gradiente conjugado preacondicionado. Tiene 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 finaliza el bucle de salida 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\,}$

Fin de repetición

El resultado es x _{k +1}

La formulación anterior es equivalente a aplicar el método del gradiente conjugado regular al sistema preacondicionado ^[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} .}$

La descomposición de Cholesky del precondicionador debe utilizarse para mantener la simetría (y la definitividad positiva) del sistema. Sin embargo, no es necesario calcular esta descomposición y es suficiente saber que . 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 preacondicionadora M debe ser simétrica, definida positiva y fija, es decir, no puede cambiar de iteración a iteración. Si se viola alguna de estas suposiciones sobre el preacondicionador, el comportamiento del método del gradiente conjugado preacondicionado puede volverse impredecible.

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

Utilizando el 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 invertirla demandaría más tiempo y recursos computacionales que resolver el algoritmo de gradiente conjugado en sí. Como ejemplo, digamos que estamos usando un precondicionador que proviene de una factorización de Cholesky incompleta. La matriz resultante es la matriz triangular inferior y la matriz del precondicionador es: ${\displaystyle \mathbf {M} }$ ${\displaystyle \mathbf {M} ^{-1}}$ ${\displaystyle \mathbf {M} }$ ${\displaystyle \mathbf {L} }$

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

Luego 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 intermedio : ${\displaystyle \mathbf {a} }$

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

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

Como y y se conocen, y es triangular inferior, resolver para es fácil y computacionalmente económico usando sustitución hacia adelante . 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, resolver es fácil y computacionalmente económico mediante sustitución hacia atrás . ${\displaystyle \mathbf {a} }$ ${\displaystyle \mathbf {L} ^{\mathsf {T}}}$ ${\displaystyle \mathbf {L} ^{\mathsf {T}}}$ ${\displaystyle \mathbf {z} }$

Usando este método, no hay necesidad de invertir ni de forma explícita, y aún así obtenemos . ${\displaystyle \mathbf {M} }$ ${\displaystyle \mathbf {L} }$ ${\displaystyle \mathbf {z} }$

El método del gradiente conjugado preacondicionado flexible

En aplicaciones numéricamente desafiantes, se utilizan preacondicionadores sofisticados, que pueden llevar a un preacondicionamiento variable, que cambia entre iteraciones. Incluso si el preacondicionador 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. Uso de 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 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 drásticamente la convergencia en este caso. ^[13] Esta versión del método de gradiente conjugado preacondicionado se puede llamar ^[14] flexible , ya que permite un preacondicionamiento variable. También se muestra ^[15] que la versión flexible 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 preacondicionador 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. ^[16]

Vs. el método de descenso más pronunciado óptimo a nivel local

En los métodos de gradiente conjugado originales y preacondicionados, solo es necesario establecer para que sean localmente óptimos, utilizando la búsqueda de línea , métodos de descenso más pronunciado . 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 preacondicionador (altamente) variable y/o no SPD , consulte más arriba. ${\displaystyle \beta _{k}:=0}$

Método de gradiente conjugado como controlador de retroalimentación óptimo para integradores dobles

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 para el sistema de doble integrador . Las cantidades y son ganancias de retroalimentación variables. ^[17] $${\displaystyle u=k(x,v):=-\gamma _{a}\nabla f(x)-\gamma _{b}v}$$ $${\displaystyle {\dot {x}}=v,\quad {\dot {v}}=u}$$ ${\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ándolo 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 matriz transpuesta-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 las ecuaciones normales es que el número de condición κ( A ^T A ) es igual a κ ² ( A ) y, por lo tanto, la tasa de convergencia de CGNR puede ser lenta y la calidad de la solución aproximada puede ser sensible a errores de redondeo. Encontrar un buen preacondicionador es a menudo 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 hermíticas complejas

El método del gradiente conjugado con una modificación trivial es extensible para resolver, dada una matriz A y un vector b de valor complejo, el sistema de ecuaciones lineales para el vector x de valor complejo, donde A es hermítica (es decir, A' = A) y la matriz definida positiva , 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] : Sec.7.3}

Un ejemplo patológico

Este ejemplo es de ^[19] Sea , y defina Dado que es invertible, existe una solución única para . Resolverlo por descenso de gradiente conjugado nos da una convergencia bastante mala: En palabras, durante el proceso de CG, el error crece exponencialmente, hasta que de repente se vuelve cero cuando se encuentra la solución única. ${\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}$$

Véase también

• Método del gradiente biconjugado (BiCG)
• Método del gradiente cuadrado conjugado (CGS)
• Método del residuo conjugado
• Propagación de creencias gaussianas
• Método iterativo: sistemas lineales
• Subespacio de Krylov
• Método de gradiente conjugado no lineal
• Preacondicionamiento
• Multiplicación de matrices dispersas por vectores

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 – vía NASA Technical Reports Server.
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. Paquette, Elliot; Trogdon, Thomas (marzo de 2023). "Universalidad de los algoritmos de gradiente conjugado y MINRES en matrices de covarianza de muestra". Communications on Pure and Applied Mathematics . 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 el dolor agonizante (PDF) .
8. Saad, Yousef (2003). Métodos iterativos para sistemas lineales dispersos (2.ª ed.). Filadelfia, Pensilvania: Society for Industrial and Applied Mathematics. pp. 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; Berry, Michael; Chan, Tony F.; Demmel, James; Donato, June; Dongarra, Jack; Eijkhout, Victor; Pozo, Roldan; Romine, Charles; van der Vorst, Henk. Plantillas para la solución de sistemas lineales: bloques de construcción para métodos iterativos (PDF) (2.ª ed.). Filadelfia, PA: SIAM. p. 13. Consultado el 31 de marzo de 2020 .
11. Golub, Gene H.; Van Loan, Charles F. (2013). Cálculos matriciales (4.ª ed.). Johns Hopkins University Press. Sec. 11.5.2. ISBN 978-1-4214-0794-4.
12. Concus, P.; Golub, GH; Meurant, G. (1985). "Preacondicionamiento de bloques para el método del gradiente conjugado". Revista SIAM sobre computación científica y estadística . 6 (1): 220–252. doi :10.1137/0906018.
13. Golub, Gene H.; Ye, Qiang (1999). "Método de gradiente conjugado preacondicionado inexacto con iteración interna-externa". Revista SIAM de informática 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 informática científica . 22 (4): 1444–1460. CiteSeerX 10.1.1.35.7473 . doi :10.1137/S1064827599362314. 
15. Bouwmeester, Henricus; Dougherty, Andrew; Knyazev, Andrew V. (2015). "Preacondicionamiento no simétrico para métodos de gradiente conjugado y descenso más pronunciado 1". Procedia Computer Science . 51 : 276–285. arXiv : 1212.6680 . doi : 10.1016/j.procs.2015.05.241 . S2CID  51978658.
16. Knyazev, Andrew V.; Lashuk, Ilya (2008). "Métodos de descenso más pronunciado y gradiente conjugado con preacondicionamiento variable". Revista SIAM sobre análisis de matrices y aplicaciones . 29 (4): 1267. arXiv : math/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, Fabian Pedregosa, Courtney Paquette, Tom Trogdon, Jeffrey. "Tutorial de teoría de matrices aleatorias y aprendizaje automático". random-matrix-learning.github.io . Consultado el 5 de diciembre de 2023 .{{cite web}}: CS1 maint: multiple names: authors list (link)

Lectura adicional

• Atkinson, Kendell A. (1988). "Sección 8.9". Introducción al análisis numérico (2.ª ed.). John Wiley and Sons. ISBN 978-0-471-50023-0.
• Avriel, Mordecai (2003). Programación no lineal: análisis y métodos . Dover Publishing. ISBN 978-0-486-43227-4.
• Golub, Gene H.; Van Loan, Charles F. (2013). "Capítulo 11". Cálculos matriciales (4.ª ed.). Johns Hopkins University Press. 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", Numerical Algorithms, 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]