Aceleración de Anderson

En matemáticas, la aceleración de Anderson , también llamada mezcla de Anderson , es un método para la aceleración de la tasa de convergencia de iteraciones de punto fijo . Introducida por Donald G. Anderson, ^[1] esta técnica se puede utilizar para encontrar la solución de ecuaciones de punto fijo que surgen a menudo en el campo de la ciencia computacional . ${\displaystyle f(x)=x}$

Dada una función , considere el problema de encontrar un punto fijo de , que es una solución a la ecuación . Un enfoque clásico para el problema es emplear un esquema de iteración de punto fijo ; ^[2] es decir, dada una estimación inicial para la solución, calcular la secuencia hasta que se cumpla algún criterio de convergencia. Sin embargo, la convergencia de dicho esquema no está garantizada en general; además, la tasa de convergencia suele ser lineal, lo que puede volverse demasiado lento si la evaluación de la función es computacionalmente costosa. ^[2] La aceleración de Anderson es un método para acelerar la convergencia de la secuencia de punto fijo. ^[2]${\displaystyle f:\mathbb {R} ^{n}\to \mathbb {R} ^{n}}$${\estilo de visualización f}$${\displaystyle f(x)=x}$$estilo de visualización x_{0}}$${\displaystyle x_{i+1}=f(x_{i})}$${\estilo de visualización f}$

Defina el residuo , y denote y (donde corresponde a la secuencia de iteraciones del párrafo anterior). Dado un valor inicial y un parámetro entero , el método se puede formular de la siguiente manera: ^{[3] [nota 1]}${\displaystyle g(x)=f(x)-x}$${\displaystyle f_{k}=f(x_{k})}$${\displaystyle g_{k}=g(x_{k})}$$Estilo de visualización x_{k}}$$estilo de visualización x_{0}}$${\displaystyle m\geq 1}$

$Estilo de visualización x_{1}=f(x_{0})}$

${\displaystyle \para todos k=1,2,\puntos}$

${\displaystyle m_{k}=\min\{m,k\}}$

${\displaystyle G_{k}={\begin{bmatrix}g_{k-m_{k}}&\dots &g_{k}\end{bmatrix}}}$

${\displaystyle \alpha _{k}=\operatorname {argmin} _{\alpha \in A_{k}}\|G_{k}\alpha \|_{2},\quad {\text{donde}}\;A_{k}=\{\alpha =(\alpha _{0},\dots ,\alpha _{m_{k}})\in \mathbb {R} ^{m_{k}+1}:\sum _{i=0}^{m_{k}}\alpha _{i}=1\}}$

${\displaystyle x_{k+1}=\sum _{i=0}^{m_{k}}(\alpha _{k})_{i}f_{k-m_{k}+i}}$

donde la multiplicación de la matriz por el vector , y es el elemento n de . Se pueden utilizar criterios de detención convencionales para finalizar las iteraciones del método. Por ejemplo, las iteraciones se pueden detener cuando cae por debajo de una tolerancia prescrita, o cuando el residuo cae por debajo de una tolerancia prescrita. ^[2]${\displaystyle G_{k}\alpha =\sum _{i=0}^{m_{k}}(\alpha )_{i}g_{k-m_{k}+i}}$${\displaystyle (\alpha )_{i}}$${\estilo de visualización i}$${\estilo de visualización \alpha}$${\displaystyle \|x_{k+1}-x_{k}\|}$$estilo de visualización g(x_{k})}$

Con respecto a la iteración de punto fijo estándar, se ha descubierto que el método converge más rápido y es más robusto, y en algunos casos evita la divergencia de la secuencia de punto fijo. ^{[3] [4]}

Para la solución , sabemos que , lo que equivale a decir que . Por lo tanto, podemos reformular el problema como un problema de optimización en el que queremos minimizar . ${\estilo de visualización x^{*}}$${\displaystyle f(x^{*})=x^{*}}$${\displaystyle g(x^{*})={\vec {0}}}$$estilo de visualización {\displaystyle g(x)\|_{2}}$

En lugar de ir directamente de a eligiendo como en la iteración de punto fijo , consideremos un punto intermedio que elegimos como la combinación lineal , donde el vector de coeficientes y es la matriz que contiene los últimos puntos, y elijamos de manera que minimice . Como los elementos de suman uno, podemos hacer la aproximación de primer orden y nuestro problema pasa a ser encontrar el que minimiza . Después de haber encontrado , podríamos en principio calcular . $Estilo de visualización x_{k}}$$estilo de visualización x_{k+1}}$${\displaystyle x_{k+1}=f(x_{k})}$$estilo de visualización x'_{k+1}}$${\displaystyle x'_{k+1}=X_{k}\alpha _{k}}$${\displaystyle \alpha_{k}\en A_{k}}$${\displaystyle X_{k}={\begin{bmatrix}x_{k-m_{k}}&\dots &x_{k}\end{bmatrix}}}$${\displaystyle m_{k}+1}$${\displaystyle x'_{k+1}}$${\displaystyle \|g(x'_{k+1})\|_{2}}$${\displaystyle \alpha _{k}}$${\displaystyle g(X_{k}\alpha _{k})=g\left(\sum _{i=0}^{m_{k}}(\alpha _{k})_{i}x_{k-m_{k}+i}\right)\approx \sum _{i=0}^{m_{k}}(\alpha _{k})_{i}g(x_{k-m_{k}+i})=G_{k}\alpha _{k}}$${\displaystyle \alpha }$${\displaystyle \|G_{k}\alpha \|_{2}}$${\displaystyle \alpha _{k}}$${\displaystyle x'_{k+1}}$

Sin embargo, dado que está diseñado para acercar un punto a , probablemente esté más cerca de que lo que es, por lo que tiene sentido elegir en lugar de . Además, dado que los elementos de suman uno, podemos hacer la aproximación de primer orden . Por lo tanto, elegimos ${\displaystyle f}$${\displaystyle x^{*}}$${\displaystyle f(x'_{k+1})}$${\displaystyle x^{*}}$${\displaystyle x'_{k+1}}$${\displaystyle x_{k+1}=f(x'_{k+1})}$${\displaystyle x_{k+1}=x'_{k+1}}$${\displaystyle \alpha _{k}}$${\displaystyle f(x'_{k+1})=f\left(\sum _{i=0}^{m_{k}}(\alpha _{k})_{i}x_{k-m_{k}+i}\right)\approx \sum _{i=0}^{m_{k}}(\alpha _{k})_{i}f(x_{k-m_{k}+i})=\sum _{i=0}^{m_{k}}(\alpha _{k})_{i}f_{k-m_{k}+i}}$

${\displaystyle x_{k+1}=\sum _{i=0}^{m_{k}}(\alpha _{k})_{i}f_{k-m_{k}+i}}$.

En cada iteración del algoritmo, se debe resolver el problema de optimización restringido , sujeto a . El problema se puede reformular en varias formulaciones equivalentes, ^[3] lo que produce diferentes métodos de solución que pueden dar como resultado una implementación más conveniente: ${\displaystyle \operatorname {argmin} \|G_{k}\alpha \|_{2}}$${\displaystyle \alpha \in A_{k}}$

• definir las matrices y , resolver y establecer ; ^{[3] [4]}${\displaystyle {\mathcal {G}}_{k}={\begin{bmatrix}g_{k-m_{k}+1}-g_{k-m_{k}}&\dots &g_{k}-g_{k-1}\end{bmatrix}}}$${\displaystyle {\mathcal {X}}_{k}={\begin{bmatrix}x_{k-m_{k}+1}-x_{k-m_{k}}&\dots &x_{k}-x_{k-1}\end{bmatrix}}}$${\displaystyle \gamma _{k}=\operatorname {argmin} _{\gamma \in \mathbb {R} ^{m_{k}}}\|g_{k}-{\mathcal {G}}_{k}\gamma \|_{2}}$${\displaystyle x_{k+1}=x_{k}+g_{k}-({\mathcal {X}}_{k}+{\mathcal {G}}_{k})\gamma _{k}}$
• Resuelve y luego establece . ^[1]${\displaystyle \theta _{k}=\{(\theta _{k})_{i}\}_{i=1}^{m_{k}}=\operatorname {argmin} _{\theta \in \mathbb {R} ^{m_{k}}}\left\|g_{k}+\sum _{i=1}^{m_{k}}\theta _{i}(g_{k-i}-g_{k})\right\|_{2}}$${\displaystyle x_{k+1}=x_{k}+g_{k}+\sum _{j=1}^{m_{k}}(\theta _{k})_{j}(x_{k-j}-x_{k}+g_{k-j}-g_{k})}$

Para ambas opciones, el problema de optimización se presenta en forma de un problema de mínimos cuadrados lineal sin restricciones , que se puede resolver mediante métodos estándar, incluida la descomposición QR ^[3] y la descomposición en valores singulares ^[4] , posiblemente incluyendo técnicas de regularización para abordar las deficiencias de rango y los problemas de condicionamiento en el problema de optimización. Resolver el problema de mínimos cuadrados mediante la resolución de ecuaciones normales generalmente no es aconsejable debido a las posibles inestabilidades numéricas y al alto costo computacional general. ^[4]

El estancamiento del método (es decir, iteraciones posteriores con el mismo valor, ) hace que el método deje de funcionar, debido a la singularidad del problema de mínimos cuadrados. De manera similar, el casi estancamiento ( ) da como resultado un mal condicionamiento del problema de mínimos cuadrados. Además, la elección del parámetro podría ser relevante para determinar el condicionamiento del problema de mínimos cuadrados, como se analiza a continuación. ^[3]${\displaystyle x_{k+1}=x_{k}}$${\displaystyle x_{k+1}\approx x_{k}}$${\displaystyle m}$

El algoritmo se puede modificar introduciendo un parámetro de relajación variable (o parámetro de mezcla) . ^{[1] [3] [4]} En cada paso, calcule la nueva iteración como La elección de es crucial para las propiedades de convergencia del método; en principio, podría variar en cada iteración, aunque a menudo se elige que sea constante. ^[4]${\displaystyle \beta _{k}>0}$$${\displaystyle x_{k+1}=(1-\beta _{k})\sum _{i=0}^{m_{k}}(\alpha _{k})_{i}x_{k-m_{k}+i}+\beta _{k}\sum _{i=0}^{m_{k}}(\alpha _{k})_{i}f(x_{k-m_{k}+i})\;.}$$${\displaystyle \beta _{k}}$${\displaystyle \beta _{k}}$

El parámetro determina cuánta información de iteraciones anteriores se utiliza para calcular la nueva iteración . Por un lado, si se elige demasiado pequeño, se utiliza muy poca información y la convergencia puede ser indeseablemente lenta. Por otro lado, si es demasiado grande, la información de iteraciones anteriores puede conservarse para demasiadas iteraciones posteriores, por lo que nuevamente la convergencia puede ser lenta. ^[3] Además, la elección de afecta el tamaño del problema de optimización. Un valor demasiado grande de puede empeorar el condicionamiento del problema de mínimos cuadrados y el costo de su solución. ^[3] En general, el problema particular a resolver determina la mejor elección del parámetro. ^[3]${\displaystyle m}$${\displaystyle x_{k+1}}$${\displaystyle m}$${\displaystyle m}$${\displaystyle m}$${\displaystyle m}$${\displaystyle m}$

Con respecto al algoritmo descrito anteriormente, la elección de en cada iteración puede modificarse. Una posibilidad es elegir para cada iteración (a veces denominada aceleración de Anderson sin truncamiento). ^[3] De esta manera, cada nueva iteración se calcula utilizando todas las iteraciones calculadas previamente. Una técnica más sofisticada se basa en elegir de manera que se mantenga un condicionamiento lo suficientemente pequeño para el problema de mínimos cuadrados. ^[3]${\displaystyle m_{k}}$${\displaystyle m_{k}=k}$${\displaystyle k}$${\displaystyle x_{k+1}}$${\displaystyle m_{k}}$

El método de Newton se puede aplicar a la solución de para calcular un punto fijo de con convergencia cuadrática. Sin embargo, dicho método requiere la evaluación de la derivada exacta de , lo que puede ser muy costoso. ^[4] Aproximar la derivada por medio de diferencias finitas es una alternativa posible, pero requiere múltiples evaluaciones de en cada iteración, lo que nuevamente puede volverse muy costoso. La aceleración de Anderson requiere solo una evaluación de la función por iteración y ninguna evaluación de su derivada. Por otro lado, la convergencia de una secuencia de puntos fijos acelerada por Anderson sigue siendo lineal en general. ^[5]${\displaystyle f(x)-x=0}$${\displaystyle f(x)}$${\displaystyle f(x)}$${\displaystyle f(x)}$${\displaystyle f(x)}$

Varios autores han señalado similitudes entre el esquema de aceleración de Anderson y otros métodos para la solución de ecuaciones no lineales. En particular:

• Eyert ^[6] y Fang y Saad ^[4] interpretaron el algoritmo dentro de la clase de métodos cuasi-Newton y multisecante, que generalizan el conocido método secante , para la solución de la ecuación no lineal ; también mostraron cómo el esquema puede verse como un método en la clase Broyden ; ^[7]${\displaystyle g(x)=0}$
• Walker y Ni ^{[3] [8]} demostraron que el esquema de aceleración de Anderson es equivalente al método GMRES en el caso de problemas lineales (es decir, el problema de encontrar una solución para alguna matriz cuadrada ), y por lo tanto puede verse como una generalización de GMRES al caso no lineal; Washio y Oosterlee encontraron un resultado similar. ^[9]${\displaystyle A\mathbf {x} =\mathbf {x} }$${\displaystyle A}$

Además, otros autores han desarrollado independientemente varios métodos equivalentes o casi equivalentes, ^{[9] [10] [11] [12] [13]} aunque con mayor frecuencia en el contexto de alguna aplicación específica de interés y no como un método general para ecuaciones de punto fijo.

El siguiente es un ejemplo de implementación en lenguaje MATLAB del esquema de aceleración de Anderson para hallar el punto fijo de la función . Observe que: ${\displaystyle f(x)=\sin(x)+\arctan(x)}$

• El problema de optimización se resolvió en la forma utilizando descomposición QR;${\displaystyle \gamma _{k}=\operatorname {argmin} _{\gamma \in \mathbb {R} ^{m_{k}}}\|g_{k}-{\mathcal {G}}_{k}\gamma \|_{2}}$
• El cálculo de la descomposición QR es subóptimo: de hecho, en cada iteración se agrega una sola columna a la matriz , y posiblemente se elimina una sola columna; este hecho se puede aprovechar para actualizar eficientemente la descomposición QR con menos esfuerzo computacional; ^[14]${\displaystyle {\mathcal {G}}_{k}}$
• El algoritmo se puede hacer más eficiente en el uso de la memoria almacenando sólo las últimas iteraciones y residuos, si no se necesita todo el vector de iteraciones;${\displaystyle x_{k}}$
• El código se generaliza directamente al caso de un valor vectorial .${\displaystyle f(x)}$

f = @( x ) sin ( x ) + atan ( x ); % Función cuyo punto fijo se va a calcular. x0 = 1 ; % Estimación inicial.      
   

k_max = 100 ; % Número máximo de iteraciones. tol_res = 1e-6 ; % Tolerancia en el residuo. m = 3 ; % Parámetro m.   
   
   

x = [ x0 , f ( x0 )]; % Vector de iteraciones x. g = f ( x ) - x ; % Vector de residuos.    
     

G_k = g ( 2 ) - g ( 1 ); % Matriz de incrementos en residuos. X_k = x ( 2 ) - x ( 1 ); % Matriz de incrementos en x.     
     

k = 2 ; mientras k < k_max && abs ( g ( k )) > tol_res m_k = min ( k , m ); % Resolver el problema de optimización por descomposición QR. [ Q , R ] = qr ( G_k ); gamma_k = R \ ( Q ' * g ( k )); % Calcular nueva iteración y nuevo residuo. x ( k + 1 ) = x ( k ) + g ( k ) - ( X_k + G_k ) * gamma_k ; g ( k + 1 ) = f ( x ( k + 1 )) - x ( k + 1 ); % Actualizar matrices de incremento con nuevos elementos. X_k = [ X_k , x ( k + 1 ) - x ( k )]; G_k = [ G_k , g ( k + 1 ) - g ( k )]; n = tamaño ( X_k , 2 ); si n > m_k X_k = X_k (:, n - m_k + 1 : fin ); G_k = G_k (:, n - m_k + 1 : fin ); fin k = k + 1 ; fin  
       
       
 
    
       
          
 
    
                
              
 
    
           
           
 
       
       
               
               
    
 
        


% Imprime el resultado: Punto fijo calculado 2.013444 después de 9 iteraciones 
fprintf ( "Punto fijo calculado %f después de %d iteraciones\n" , x ( fin ), k );  

• Punto fijo
• Iteración de punto fijo
• Algoritmos de búsqueda de raíces