Método de gradiente conjugado - Conjugate gradient method

Una comparación de la convergencia del descenso del 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, asumiendo una 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-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 . Se atribuye comúnmente a Magnus Hestenes y Eduard Stiefel , quienes lo programaron en el Z4 y lo investigaron extensamente.

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.

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

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

Dado que es simétrico y positivo-definido, 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 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 , entonces es conjugado a . 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 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 por rendimientos ${\ Displaystyle \ mathbf {p} _ {k} ^ {\ mathsf {T}}}$

${\ Displaystyle \ 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} \ Rightarrow \ mathbf {p} _ {k} ^ {\ mathsf { T}} \ mathbf {b} = \ 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}} \ Rightarrow}$

${\ 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 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 . 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. ${\ Displaystyle \ mathbf {p} _ {k}}$${\ Displaystyle \ mathbf {x} _ {*}}$

Denotamos la estimación inicial de x _∗ por x ₀ (podemos suponer sin pérdida de generalidad que x ₀ = 0 ; de lo contrario, considere el sistema Az = b - Ax _{0 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 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

${\ 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 segunda derivada está dada por una matriz simétrica positiva definida.

${\ Displaystyle \ nabla ^ {2} f (\ mathbf {x}) = \ mathbf {A} \ ,,}$

y que el minimizador (use D f ( x ) = 0) resuelve el problema inicial es obvio a partir 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 suposición inicial x ₀ , esto significa que tomamos p ₀ = b - Ax ₀ . 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 ₀ es también el residuo 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 del gradiente requeriría moverse en la dirección r _k . Aquí, sin embargo, insistimos en que las direcciones se conjugan 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: ${\ 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}}$

(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

${\ 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 de se puede derivar si se sustituye la expresión de x _{k +1} en f y se minimiza wrt${\ Displaystyle \ mathbf {r} _ {k}}$${\ Displaystyle \ alpha _ {k}}$${\ Displaystyle \ alpha _ {k}}$

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

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 previas 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 el algoritmo avanza y 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 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}}$

El algoritmo se detalla a continuación para resolver Ax = b 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 {A}}$${\ Displaystyle \ mathbf {x} _ {0}}$

${\ displaystyle {\ begin {alineado} & \ mathbf {r} _ {0}: = \ mathbf {b} - \ mathbf {Ax} _ {0} \\ & {\ hbox {if}} \ mathbf {r } _ {0} {\ text {es lo suficientemente pequeño, luego devuelve}} \ mathbf {x} _ {0} {\ text {como resultado}} \\ & \ 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 {es lo suficientemente pequeño, luego sale del bucle}} \\ & \ 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 {como resultado}} \ end {alineado}}}$

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

Observamos que se calcula mediante el método de descenso de gradiente aplicado a . De manera similar, la configuración se calcularía mediante el método de descenso de gradiente de , es decir, se puede usar como una implementación simple de un reinicio de las iteraciones de gradiente conjugado. 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 residual explícito

Las fórmulas y , que tanto HOLD en aritmética exacta, hacen que las fórmulas y matemáticamente equivalentes. El primero se usa en el algoritmo para evitar una multiplicación adicional por ya que 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 errores de redondeo , por lo que se recomienda para una evaluación ocasional. ${\ 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}}$

Por lo general, se utiliza una norma del residual para detener los criterios. 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 sigue disminuyendo 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, α _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}}}}$

puesto . El β _k se elige de manera que se conjugue 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}}}}$

utilizando

${\ 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 la β en el algoritmo después de cancelar α _k .

Código de ejemplo en MATLAB / GNU Octave

function x = conjgrad(A, b, x)
    r = b - A * x;
    p = r;
    rsold = r' * r;

    for i = 1:length(b)
        Ap = A * p;
        alpha = rsold / (p' * Ap);
        x = x + alpha * p;
        r = r - alpha * Ap;
        rsnew = r' * r;
        if sqrt(rsnew) < 1e-10
              break
        end
        p = r + (rsnew / rsold) * p;
        rsold = rsnew;
    end
end

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

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}} = {\ begin {bmatrix} 0.2356 \\ 0.3384 \ end {bmatrix} }.}$

Este resultado completa la primera iteración, el resultado es una solución aproximada "mejorada" para el 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}} = {\ 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}}} = {\ frac {{\ begin {bmatrix} -0.2810 y 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} = {\ 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}}} = {\ 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} = {\ 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 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 del 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 la generación de los subespacios de Krylov.

Como método iterativo , el método de gradiente conjugado monótonamente (en la norma energética) mejora 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 más grande es, más lenta es la mejora. ${\ Displaystyle \ mathbf {x} _ {k}}$ ${\ Displaystyle \ kappa (A)}$${\ Displaystyle A}$${\ Displaystyle \ kappa (A)}$

Si es grande, acondicionamiento previo se utiliza comúnmente para reemplazar el sistema original con tal que es más pequeño que , ver más abajo. ${\ 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

Defina 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 defina los errores como . Ahora, la tasa de convergencia se puede aproximar como ${\ Displaystyle \ left (\ mathbf {x} _ {k} \ right) _ {k}}$${\ Displaystyle \ mathbf {x} _ {*}}$${\ Displaystyle \ mathbf {e} _ {k}: = \ mathbf {x} _ {k} - \ mathbf {x} _ {*}}$

${\ Displaystyle {\ begin {alineado} \ 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}} \ ,, \ end {alineado}}}$

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

Tenga en cuenta, 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 {para}} \ 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 asume 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 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 vinculada con , pero puede ser superlineal, según la distribución del espectro de la matriz y la distribución espectral del error. 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. ${\ Displaystyle {\ sqrt {\ kappa (\ mathbf {A})}}}$${\ Displaystyle A}$

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 toma la siguiente forma:

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

${\ Displaystyle \ mathbf {z} _ {0}: = \ mathbf {M} ^ {- 1} \ 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 , salga del final del bucle si

${\ Displaystyle \ mathbf {z} _ {k + 1}: = \ mathbf {M} ^ {- 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 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.

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

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

${\ 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. Esta versión del método de gradiente conjugado preacondicionado se puede llamar flexible, ya que permite el preacondicionamiento variable. La versión flexible también se muestra 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 localmente óptimo de descenso más empinado.

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 se necesita establecer para hacerlos localmente óptimos, utilizando la búsqueda de línea , métodos de descenso más empinado . 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. ${\ 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 . En este enfoque, el método de gradiente conjugado cae como un controlador de retroalimentación óptimo ,

$${\ 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 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 κ ² ( 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 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, 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 transposición conjugada usando el estilo MATLAB / GNU Octave . 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. ${\ Displaystyle \ mathbf {A} \ mathbf {x} = \ mathbf {b}}$

Ver también

Referencias

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]