Método del gradiente conjugado

En matemáticas , el método del gradiente conjugado es un algoritmo para la solución numérica de sistemas de ecuaciones lineales específicos , concretamente aquellos cuya matriz es semidefinida positiva . Este método se suele implementar como un algoritmo iterativo , aplicable a sistemas dispersos demasiado grandes para ser resueltos mediante una implementación directa u otros métodos directos como la descomposición de Cholesky . Los sistemas dispersos de gran tamaño suelen surgir al resolver numéricamente ecuaciones diferenciales parciales o problemas de optimización.

El método del gradiente conjugado también puede utilizarse 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. Diversos métodos de gradiente conjugado no lineal buscan mínimos en problemas de optimización no lineal.

Supongamos que queremos resolver el sistema de ecuaciones lineales

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

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

El método del gradiente conjugado puede derivarse desde diversas perspectivas, incluyendo la especialización del método de la 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 concisa del método.

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

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

Desde${\displaystyle \mathbf {A} }$es simétrica y definida positiva, el lado izquierdo define un producto interno

${\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} ^{\mathbf {T}}\mathbf {v} \rangle =\langle \mathbf {u} ,\mathbf {A} \mathbf {v} \rangle .}$

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

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

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

${\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 el problema por la izquierda${\displaystyle \mathbf {Hacha} =\mathbf {b} }$con el vector${\displaystyle \mathbf {p} _{k}^{\mathsf {T}}}$rendimientos

${\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 da como resultado el siguiente método ^{[ 4 ]} para resolver la ecuación.${\displaystyle \mathbf {Hacha} =\mathbf {b} }$: encontrar una secuencia de${\displaystyle n}$direcciones conjugadas y luego calcular los coeficientes${\displaystyle \alpha _{k}}$.

Si elegimos los vectores conjugados${\displaystyle \mathbf {p} _{k}}$Si lo hacemos con cuidado, es posible que no necesitemos todos ellos para obtener una buena aproximación a la solución.${\displaystyle \mathbf {x} _{*}}$. 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${\displaystyle n}$es tan grande que el método directo llevaría demasiado tiempo.

Denotamos la estimación inicial por${\displaystyle \mathbf {x} _{*}}$por${\displaystyle \mathbf {x} _{0}}$(podemos asumir sin pérdida de generalidad que${\displaystyle \mathbf {x} _ {0}=\mathbf {0} }$, de lo contrario considere el sistema${\displaystyle \mathbf {Az} =\mathbf {b} -\mathbf {Ax} _{0}}$en su lugar). Empezando por${\displaystyle \mathbf {x} _{0}}$Buscamos la solución y en cada iteración necesitamos una métrica que nos indique si estamos más cerca de la solución.${\displaystyle \mathbf {x} _{*}}$(que desconocemos). Esta métrica proviene del hecho de que la solución${\displaystyle \mathbf {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 \mathbb {R} ^{n}\,.}$

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

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

y que el minimizador (usar${\displaystyle Df(\mathbf {x} )=0}$) resuelve el problema inicial que se deduce de su primera derivada

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

Esto sugiere tomar el primer vector base.${\displaystyle \mathbf {p} _{0}}$ser el negativo del gradiente de${\displaystyle f}$en${\displaystyle \mathbf {x} =\mathbf {x} _{0}}$. El gradiente de${\displaystyle f}$igual${\displaystyle \mathbf {Ax} -\mathbf {b} }$Comenzando con una suposición inicial${\displaystyle \mathbf {x} _{0}}$, esto significa que tomamos${\displaystyle \mathbf {p} _{0}=\mathbf {b} -\mathbf {Ax} _{0}}$. Los demás vectores de la base serán conjugados al gradiente, de ahí el nombre de método del gradiente conjugado . Tenga en cuenta que${\displaystyle \mathbf {p} _{0}}$es también el residuo proporcionado por este paso inicial del algoritmo.

Dejar${\displaystyle \mathbf {r} _{k}}$ser el residuo en el${\displaystyle k}$º paso:

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

Como se observó anteriormente,${\displaystyle \mathbf {r} _{k}}$es el gradiente negativo de${\displaystyle f}$en${\displaystyle \mathbf {x} _{k}}$, 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${\displaystyle \mathbf {p} _{k}}$Deben ser conjugadas entre sí. Una forma práctica de imponer esto es exigiendo 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 considerarse un ejemplo de ortonormalización de Gram-Schmidt . Esto da como resultado la siguiente expresión:

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

(véase la imagen en la parte superior del artículo para observar el efecto de la restricción de conjugación en la convergencia). Siguiendo esta dirección, la siguiente ubicación óptima viene 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 deduce de la definición de${\displaystyle \mathbf {r} _{k}}$. La expresión para${\displaystyle \alpha _{k}}$se puede derivar si se sustituye la expresión para x _{k +1} en f y se minimiza con respecto a${\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}}}$

utilizando álgebra lineal""" x = gradiente_conjugado(A, b, x0 = cero(b); atol=longitud(b)*eps(norma(b))Devuelve la solución de `A * x = b` utilizando el método del gradiente conjugado.`A` debe ser una matriz definida positiva u otro operador lineal.`x0` es la estimación inicial para la solución (por defecto es el vector cero).`atol` es la tolerancia absoluta en la magnitud del residuo `b - A * x`.para convergencia (el valor predeterminado es épsilon de máquina).Devuelve el vector de solución aproximada `x`."""función gradiente_conjugado (A , b :: AbstractVector , x0 :: AbstractVector = zero ( b ); atol = length ( b ) * eps ( norm ( b )))x = copiar ( x0 ) # inicializar la soluciónr = b - A * x0 # residuo inicialp = copiar ( r ) # dirección de búsqueda inicialr²old = r ' * r # norma al cuadrado del residuok = 0mientras r²old > atol ^ 2 # iterar hasta la convergenciaAp = A * p # dirección de búsquedaα = r²old / ( p ' * Ap ) # tamaño del paso@. x += α * p # actualizar solución# Actualizar residuo:Si ( k + 1 ) % 16 == 0 # cada 16 iteraciones, recalcular el residuo desde ceror .= b .- A * x # para evitar la acumulación de errores numéricosdemás@. r -= α * Ap # utiliza la fórmula de actualización que ahorra un producto matriz-vectorfinr²nuevo = r ' * r@. p = r + ( r²nuevo / r²antiguo ) * p # actualizar la dirección de búsquedar²old = r²new # actualiza la norma residual al cuadradok += 1findevolver xfin

función x = gradiente_conjugado ( A, b, x0, tol )% Devuelve la solución de `A * x = b` utilizando el método del gradiente conjugado.% Recordatorio: A debe ser simétrica y definida positiva.si nargin < 4tol = eps ;finr = b - A * x0 ;p = r ;rsold = r ' * r ;x = x0 ;mientras sqrt ( rsold ) > tolAp = A * p ;alfa = rsold / ( p ' * Ap );x = x + alfa * p ;r = r - alfa * Ap ;rsnew = r ' * r ;p = r + ( rsnew / rsold ) * p ;rsold = rsnew ;finfin

Con aritmética exacta, el número de iteraciones necesarias no supera el orden de la matriz. Este comportamiento se conoce como la propiedad de terminación finita del método del gradiente conjugado. Se refiere a la capacidad del método para alcanzar la solución exacta de un sistema lineal en un número finito de pasos —como máximo igual a la dimensión del sistema— cuando se utiliza aritmética exacta. Esta propiedad surge del hecho de que, en cada iteración, el método genera un vector residual ortogonal a todos los residuos anteriores. Estos residuos forman un conjunto mutuamente ortogonal.

En un espacio n -dimensional, es imposible construir más de n vectores linealmente independientes y mutuamente ortogonales, a menos que uno de ellos sea el vector cero. Por lo tanto, una vez que aparece un residuo cero, el método ha alcanzado la solución y debe finalizar. Esto garantiza que el método del gradiente conjugado converja en un máximo de n pasos.

Para demostrar esto, consideremos el sistema:

${\displaystyle A={\begin{bmatrix}3&-2\\-2&4\end{bmatrix}},\quad \mathbf {b} ={\begin{bmatrix}1\\1\end{bmatrix}}}$

Partimos de una suposición inicial.${\displaystyle \mathbf {x} _{0}={\begin{bmatrix}1\\2\end{bmatrix}}}$. Desde${\displaystyle A}$Si la ecuación es simétrica definida positiva y el sistema es bidimensional, el método del gradiente conjugado debería encontrar la solución exacta en no más de dos pasos. El siguiente código de MATLAB demuestra este comportamiento:

A = [ 3 , - 2 ; - 2 , 4 ]; x_true = [ 1 ; 1 ]; b = A * x_true ;x = [ 1 ; 2 ]; % estimación inicial r = b - A * x ; p = r ;para k = 1 : 2 Ap = A * p ; alpha = ( r ' * r ) / ( p ' * Ap ); x = x + alpha * p ; r_new = r - alpha * Ap ; beta = ( r_new ' * r_new ) / ( r ' * r ); p = r_new + beta * p ; r = r_new ; findisp ( 'Solución exacta:' ); disp ( x );

La salida confirma que el método alcanza${\displaystyle {\begin{bmatrix}1\\1\end{bmatrix}}}$Tras dos iteraciones, se obtuvo un resultado consistente con la predicción teórica. Este ejemplo ilustra cómo el método del gradiente conjugado se comporta como un método directo en condiciones idealizadas.

En teoría, el método del gradiente conjugado puede considerarse un método directo, ya que, en ausencia de errores de redondeo, produce la solución exacta tras un número finito de iteraciones, que no supera el tamaño de la matriz. En la práctica, nunca se obtiene la solución exacta, puesto que el método del gradiente conjugado es inestable incluso ante pequeñas perturbaciones; por ejemplo, la mayoría de las direcciones no son conjugadas en la práctica, 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 de energía) las aproximaciones.${\displaystyle \mathbf {x} _{k}}$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.${\displaystyle \kappa (A)}$de la matriz del sistema${\displaystyle A}$: el más grande${\displaystyle \kappa (A)}$es, cuanto más lenta sea la mejora. ^{[ 9 ]}

Sin embargo, surge un caso interesante cuando los valores propios están espaciados logarítmicamente para una matriz simétrica grande. Por ejemplo, sea${\displaystyle A=QDQ^{T}}$dónde${\displaystyle Q}$es una matriz ortogonal aleatoria y${\displaystyle D}$es una matriz diagonal con valores propios que van desde${\displaystyle \lambda _{n}=1}$a${\displaystyle \lambda _{1}=10^{6}}$, espaciados logarítmicamente. A pesar de la propiedad de terminación finita de CGM, donde teóricamente se debería alcanzar la solución exacta en como máximo${\displaystyle n}$pasos, el método puede mostrar estancamiento en la convergencia. En tal escenario, incluso después de muchas más iteraciones, por ejemplo, diez veces el tamaño de la matriz, el error puede disminuir solo modestamente (por ejemplo, a${\displaystyle 10^{-5}}$). Además, el error iterativo puede oscilar significativamente, lo que lo hace poco fiable como condición de parada. Esta mala convergencia no se explica solo por el número de condición (por ejemplo,${\displaystyle \kappa _{2}(A)=10^{6}}$), sino más bien por la propia distribución de los valores propios. Cuando los valores propios están más uniformemente espaciados o distribuidos aleatoriamente, estos problemas de convergencia suelen estar ausentes, lo que pone de relieve que el rendimiento de CGM no solo depende de${\displaystyle \kappa (A)}$pero también sobre cómo se distribuyen los valores propios. ^{[ 10 ]}

Si${\displaystyle \kappa (A)}$es grande, el preacondicionamiento se usa comúnmente para reemplazar el sistema original${\displaystyle \mathbf {Ax} -\mathbf {b} =0}$con${\displaystyle \mathbf {M} ^{-1}(\mathbf {Ax} -\mathbf {b} )=0}$de tal manera que${\displaystyle \kappa (\mathbf {M} ^{-1}\mathbf {A} )}$es más pequeño que${\displaystyle \kappa (\mathbf {A} )}$, vea abajo.

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

${\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 salir del bucle. Fin del bucle.

${\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 ^{[ 13 ].}

${\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 positividad definida) del sistema. Sin embargo, no es necesario calcular esta descomposición, y basta con saber${\displaystyle \mathbf {M} ^{-1}}$Se puede demostrar que${\displaystyle \mathbf {E} ^{-1}\mathbf {A} (\mathbf {E} ^{-1})^{\mathsf {T}}}$tiene el mismo espectro que${\displaystyle \mathbf {M} ^{-1}\mathbf {A} }$.

La matriz de preacondicionamiento${\displaystyle \mathbf {M} }$Debe ser simétrica, definida positiva y fija, es decir, no puede cambiar de una iteración a otra. Si se incumple alguna de estas condiciones sobre el precondicionador, el comportamiento del método del gradiente conjugado precondicionado puede volverse impredecible.

Un ejemplo de precondicionador comúnmente utilizado es la factorización de Cholesky incompleta . ^{[ 14 ]}

En aplicaciones numéricamente complejas, se utilizan precondicionadores sofisticados, lo que puede dar lugar a un precondicionamiento variable, que cambia entre iteraciones. Incluso si el precondicionador es simétrico definido positivo en cada iteración, el hecho de que pueda cambiar invalida los argumentos anteriores y, en pruebas prácticas, provoca una ralentización significativa de la convergencia del algoritmo presentado anteriormente. Utilizando 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 drásticamente la convergencia en este caso. ^{[ 15 ]} Esta versión del método del gradiente conjugado precondicionado puede llamarse ^{[ 16 ]} flexible , ya que permite un precondicionamiento variable. También se ha demostrado ^{[ 17 ]} que la versión flexible es robusta incluso si el precondicionador no es simétrico definido positivo (SPD).

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

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 lentamente que el método de descenso más pronunciado localmente óptimo. ^{[ 18 ]}

En ambos métodos de gradiente conjugado, el original y el precondicionado, solo es necesario establecer ${\displaystyle \beta _{k}:=0}$Para optimizarlos localmente, se utilizan métodos de búsqueda lineal y descenso más pronunciado . Con esta sustitución, los vectores p son siempre iguales a 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 económica que 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 que no sea SPD ( véase más arriba).

El método del gradiente conjugado también se puede derivar utilizando la teoría de control óptimo . ^{[ 19 ]} 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}$$para el sistema de doble integrador ,$${\displaystyle {\dot {x}}=v,\quad {\dot {v}}=u}$$Las cantidades${\displaystyle \gamma _{a}}$y${\displaystyle \gamma _{b}}$son ganancias de retroalimentación variables. ^{[ 19 ]}

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 simétrica semidefinida positiva para cualquier A . El resultado es el 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 memoria, sino solo realizar las multiplicaciones matriz-vector y transpuesta 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 las 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 a errores de redondeo. Encontrar un buen precondicionador suele ser una parte importante del uso del método CGNR.

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

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

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