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-definida . El método de gradiente conjugado a menudo se implementa 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 grandes sistemas dispersos a menudo surgen cuando se resuelven numéricamente ecuaciones diferenciales parciales o problemas de optimización.
El método de gradiente conjugado también se puede utilizar para resolver problemas de optimización sin restricciones , como la minimización de energía . Es comúnmente atribuido a Magnus Hestenes y Eduard Stiefel , [1] [2] quienes lo programaron en el Z4 , [3] y lo investigaron extensamente. [4] [5]
El método de gradiente biconjugado proporciona una generalización a matrices no simétricas. Varios métodos de gradiente conjugado no lineal buscan mínimos de ecuaciones no lineales y funciones objetivas de caja negra.
Descripción del problema abordado por gradientes conjugados
Supongamos que queremos resolver el sistema de ecuaciones lineales.
para el vector , donde lo conocido matriz es simétrico (es decir, A T = A ), positivo-definido (es decir, x T Ax > 0 para todos los vectores distintos de ceroen R n ), y real , ytambién es conocido. Denotamos la solución única de este sistema por.
Derivación como método directo
El método de gradiente conjugado puede derivarse de 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 y sucinta formulación del método.
Decimos que los no-cero dos vectores u y v son conjugado (con respecto a) Si
Desde es simétrico y positivo-definido, el lado izquierdo define un producto interno
Dos vectores se conjugan si y solo si son ortogonales con respecto a este producto interno. Ser conjugado es una relación simétrica: si es conjugado a , luego es conjugado a . Suponer que
es un conjunto de vectores mutuamente conjugados con respecto a , es decir . Luegoforma una base para, y podemos expresar la solución de en esta base:
Multiplicar por la izquierda por rendimientos
Esto da el siguiente método [4] para resolver la ecuación Ax = b : encuentre una secuencia deconjugar direcciones y luego calcular los coeficientes α k .
Como método iterativo
Si elegimos los vectores conjugados cuidadosamente, entonces es posible que no los necesitemos todos para obtener una buena aproximación a la solución . Entonces, queremos considerar el método de 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.
Denotamos la suposición inicial para x ∗ por x 0 (podemos suponer sin pérdida de generalidad que x 0 = 0 ; de lo contrario, considere el sistema Az = b - Ax 0 en su lugar). Comenzando con x 0 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 desconocemos). Esta métrica proviene del hecho de que la solución x ∗ también es el minimizador único de la siguiente función cuadrática
La existencia de un minimizador único es evidente ya que su segunda derivada está dada por una matriz simétrica positiva definida.
y que el minimizador (use D f ( x ) = 0) resuelve el problema inicial es obvio a partir de su primera derivada
Esto sugiere tomar el primer vector base p 0 como el negativo del gradiente de f en x = x 0 . El gradiente de f es igual a Ax - b . Comenzando con una suposición inicial x 0 , esto significa que tomamos p 0 = b - Ax 0 . Los otros vectores de la base se conjugarán con el gradiente, de ahí el nombre método de gradiente conjugado . Tenga en cuenta que p 0 es también el residuo proporcionado por este paso inicial del algoritmo.
Sea r k el residual en el k- ésimo paso:
Como se observó anteriormente, es el gradiente negativo de a , 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 direccionesser conjugados 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 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:
(vea 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 viene dada por
con
donde la última igualdad se sigue de la definición de . La expresión parase puede derivar si se sustituye la expresión de x k +1 en f y se minimiza wrt
El algoritmo resultante
El algoritmo anterior ofrece la explicación más sencilla del método de gradiente conjugado. Aparentemente, el algoritmo 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. Yes -ortogonal a , es decir , por . Esto puede considerarse que a medida que avanza el algoritmo, y abarcan el mismo subespacio de Krylov . Dónde forman la base ortogonal con respecto al producto interior estándar, y forman la base ortogonal con respecto al producto interno inducido por . Por lo tanto, puede considerarse como la proyección de en el subespacio de Krylov.
El algoritmo se detalla a continuación para resolver Ax = b dondees una matriz real, simétrica y definida positiva. El vector de entradapuede ser una solución inicial aproximada o 0 . Es una formulación diferente del procedimiento exacto descrito anteriormente.
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 .
Reinicia
Notamos eso se calcula mediante el método de descenso de gradiente aplicado a. Configuración haría de manera similar calculado por el método de descenso de gradiente 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 .
Cálculo residual explícito
Las fórmulas y , que ambos se mantienen en aritmética exacta, hacen que las fórmulas y matemáticamente equivalente. El primero se utiliza en el algoritmo para evitar una multiplicación extra por desde el vector ya está calculado para evaluar . Este último puede ser más preciso, sustituyendo el cálculo explícitopara el implícito por la recursividad sujeta a acumulación de errores de redondeo , por lo que se recomienda para una evaluación ocasional. [6]
Cálculo de alfa y beta
En el algoritmo, α k se elige de manera que es ortogonal a . El denominador se simplifica de
desde . El β k se elige de manera que es conjugado a . Inicialmente, β k es
utilizando
y equivalentemente
el numerador de β k se reescribe como
porque y son ortogonales por diseño. El denominador se reescribe como
usando que las direcciones de búsqueda p k están conjugadas y nuevamente que los residuos son ortogonales. Esto da la β en el algoritmo después de cancelar α k .
Código de ejemplo en MATLAB / GNU Octave
función x = conjgrad ( A, b, x ) r = b - A * x ; p = r ; rsold = r ' * r ; para i = 1 : longitud ( b ) Ap = A * p ; alfa = rsvendido / ( p ' * Ap ); x = x + alfa * p ; r = r - alfa * Ap ; rsnew = r ' * r ; si sqrt ( rsnew ) < 1e-10 rotura final p = r + ( rsnew / rsold ) * p ; rsold = rsnew ; finalfinal
Ejemplo numérico
Considere el sistema lineal Ax = b dado por
Realizaremos dos pasos del método de gradiente conjugado comenzando con la suposición inicial
para encontrar una solución aproximada al sistema.
Solución
Como referencia, la solución exacta es
Nuestro primer paso es calcular el vector residual r 0 asociado con x 0 . Este residual se calcula a partir de la fórmula r 0 = b - Ax 0 , y en nuestro caso es igual a
Dado que esta es la primera iteración, usaremos el vector residual r 0 como nuestra dirección de búsqueda inicial p 0 ; el método de selección de p k cambiará en iteraciones posteriores.
Ahora calculamos el escalar α 0 usando la relación
Ahora podemos calcular x 1 usando la fórmula
Este resultado completa la primera iteración, el resultado es una solución aproximada "mejorada" para el sistema, x 1 . Ahora podemos continuar y calcular el siguiente vector residual r 1 usando la fórmula
Nuestro siguiente paso en el proceso es calcular el β 0 escalar que eventualmente se utilizará para determinar la siguiente dirección de búsqueda p 1 .
Ahora, usando este escalar β 0 , podemos calcular la siguiente dirección de búsqueda p 1 usando la relación
Ahora calculamos el escalar α 1 usando nuestro p 1 recién adquirido usando el mismo método que el usado para α 0 .
Finalmente, encontramos x 2 usando el mismo método que se usó para encontrar x 1 .
El resultado, x 2 , es una aproximación "mejor" a la solución del sistema que x 1 y x 0 . Si en este ejemplo se usara 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 teóricamente 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 de gradiente conjugado es inestable con respecto a incluso pequeñas perturbaciones, por ejemplo, la mayoría de las direcciones no son en la práctica conjugadas, debido a la naturaleza degenerativa de generar los subespacios de Krylov.
Como método iterativo , el método de gradiente conjugado de forma monótona (en la norma energética) mejora las aproximacionesa 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 : el mas largo es decir, más lenta es la mejora. [7]
Si es grande, el preacondicionamiento se usa comúnmente para reemplazar el sistema original con tal que es más pequeña que , vea abajo.
Teorema de convergencia
Defina un subconjunto de polinomios como
dónde es el conjunto de polinomios de grado máximo.
Dejar ser las aproximaciones iterativas de la solución exacta , y defina los errores como . Ahora, la tasa de convergencia se puede aproximar como [4] [8]
dónde denota el espectro , ydenota el número de condición .
Tenga en cuenta, el límite importante cuando tiende a
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.
No round-off error se asume en el teorema de la convergencia, pero la convergencia unido es comúnmente válido en la práctica como se ha explicado teóricamente [5] por Anne Greenbaum .
Convergencia práctica
Si se inicializa al azar, 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 la convergencia teórica ligada 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 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 gran tamaño, 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.
El método de gradiente conjugado preacondicionado
En la mayoría de los casos, el preacondicionamiento es necesario para garantizar una rápida convergencia del método de gradiente conjugado. El método de gradiente conjugado preacondicionado adopta la siguiente forma: [9]
- repetir
- si r k +1 es suficientemente pequeño , salga del final del bucle si
- fin de repetir
- El resultado es x k +1
La formulación anterior es equivalente a aplicar el método de gradiente conjugado sin preacondicionamiento del sistema [10]
dónde
La matriz del preacondicionador M tiene que ser simétrica positiva definida y fija, es decir, no puede cambiar de una iteración a otra. Si se viola alguna de estas suposiciones sobre el preacondicionador, el comportamiento del método de gradiente conjugado preacondicionado puede volverse impredecible.
Un ejemplo de un preacondicionador de uso común es la factorización de Cholesky incompleta . [11]
El método de gradiente conjugado flexible preacondicionado
En aplicaciones numéricamente desafiantes, se utilizan preacondicionadores sofisticados, que pueden conducir a preacondicionamientos variables, cambiando 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 las pruebas prácticas, conduce a una desaceleración significativa de la convergencia del algoritmo presentado anteriormente. Usando la fórmula de Polak-Ribière
en lugar de la fórmula de Fletcher-Reeves
puede mejorar drásticamente la convergencia en este caso. [12] Esta versión del método de gradiente conjugado preacondicionado se puede llamar [13] flexible, ya que permite el preacondicionamiento variable. También se muestra que la versión flexible [14] 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,por 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 lento que el método localmente óptimo de descenso más empinado. [15]
Vs. el método de descenso más empinado localmente óptimo
Tanto en el método de gradiente conjugado original como en el preacondicionado, solo es necesario establecer para hacerlos localmente óptimos, utilizando la búsqueda de línea , los métodos de descenso más empinados . 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 empinado 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 (muy) variable y / o no SPD , ver más arriba.
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 . [16] En este enfoque, el método de gradiente conjugado cae como un controlador de retroalimentación óptimo ,
Gradiente conjugado en las ecuaciones normales
El método del gradiente conjugado se puede aplicar a un arbitrario n -by- m matriz aplicándolo a ecuaciones normales A T A y lado derecho vector A T b , ya que A T A es un simétrica positiva-semidefinida matriz para cualquier A . El resultado es un gradiente conjugado en las ecuaciones normales (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 las multiplicaciones matriz-vector. Por lo tanto, CGNR es particularmente útil cuando A es una matriz escasa, 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 κ 2 ( 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 acondicionado, 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 es extensible a resolver, dados la matriz A de valores complejos y el vector b, el sistema de ecuaciones lineales para el vector de valor complejo x, donde A es hermitiano (es decir, A '= A) y una matriz definida positiva , y el símbolo' denota la transposición conjugada usando el estilo de octava MATLAB / GNU . La modificación trivial es simplemente sustituir la transpuesta conjugada por la transpuesta real en todas partes. Esta sustitución es compatible con versiones anteriores, ya que la transpuesta conjugada se convierte en una transpuesta real en vectores y matrices de valor real. El código de ejemplo proporcionado anteriormente en MATLAB / GNU Octave, por lo tanto, ya funciona para matrices hermitianas complejas que no necesitan modificaciones.
Ver también
- Método de gradiente biconjugado (BiCG)
- Método residual conjugado
- Propagación de creencias gaussianas
- Método iterativo: sistemas lineales
- Subespacio de Krylov
- Método de gradiente conjugado no lineal
- Preacondicionamiento
- Multiplicación dispersa matriz-vector
Referencias
- ^ 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 .
- ^ Straeter, TA (1971). "Sobre la extensión de la clase de Davidon-Broyden de rango uno, métodos de minimización cuasi-Newton a un espacio de Hilbert de dimensión infinita con aplicaciones para problemas de control óptimos". Servidor de informes técnicos de la NASA . NASA. hdl : 2060/19710026200 .
- ^ 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.
- ^ a b c d Polyak, Boris (1987). Introducción a la optimización .
- ^ a b c Greenbaum, Anne (1997). Métodos iterativos para resolver sistemas lineales . doi : 10.1137 / 1.9781611970937 . ISBN 978-0898713961.
- ^ Shewchuk, Jonathan R (1994). Una introducción al método de gradiente conjugado sin dolor agonizante (PDF) .
- ^ Saad, Yousef (2003). Métodos iterativos para sistemas lineales dispersos (2ª ed.). Filadelfia, Pa .: Sociedad de Matemáticas Industriales y Aplicadas. págs. 195 . ISBN 978-0-89871-534-7.
- ^ Hackbusch, W. (21 de junio de 2016). Solución iterativa de grandes sistemas dispersos de ecuaciones (2ª ed.). Suiza: Springer. ISBN 9783319284835. OCLC 952572240 .
- ^ Barrett, Richard; Berry, Michael; Chan, Tony F .; Demmel, James; Donato, junio; 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. pag. 13 . Consultado el 31 de marzo de 2020 .
- ^ Golub, Gene H .; Van Loan, Charles F. (2013). Cálculos matriciales (4ª ed.). Prensa de la Universidad Johns Hopkins. segundo. 11.5.2. ISBN 978-1-4214-0794-4.
- ^ Concus, P .; Golub, GH; Meurant, G. (1985). "Preacondicionamiento 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 .
- ^ Golub, Gene H .; Ye, Qiang (1999). "Método de gradiente conjugado preacondicionado 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 .
- ^ 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 .
- ^ Henricus Bouwmeester, Andrew Dougherty, Andrew V Knyazev. Preacondicionamiento no simétrico para métodos de gradiente conjugado y descenso más pronunciado . Procedia Computer Science, Volumen 51, Páginas 276-285, Elsevier, 2015. doi : 10.1016 / j.procs.2015.05.241
- ^ Knyazev, Andrew V .; Lashuk, Ilya (2008). "Métodos de descenso más pronunciado y gradiente conjugado con preacondicionamiento variable". Revista SIAM sobre Análisis y Aplicaciones Matriciales . 29 (4): 1267. arXiv : math / 0605767 . doi : 10.1137 / 060675290 . S2CID 17614913 .
- ^ a b Ross, IM , "Una teoría de control óptimo para la optimización acelerada", arXiv : 1902.09004 , 2019.
Otras lecturas
- Atkinson, Kendell A. (1988). "Sección 8.9". Introducción al análisis numérico (2ª ed.). John Wiley e hijos. ISBN 978-0-471-50023-0.
- Avriel, Mordecai (2003). Programación no lineal: análisis y métodos . Publicación de Dover. ISBN 978-0-486-43227-4.
- Golub, Gene H .; Van Loan, 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.
enlaces externos
- "Gradientes conjugados, método de" , Enciclopedia de Matemáticas , EMS Press , 2001 [1994]