Accelerazione di Anderson

L'accelerazione di Anderson (nota anche in inglese come Anderson mixing) è un metodo per l'accelerazione della convergenza delle iterazioni di punto fisso. Introdotto da Donald G. Anderson, tale metodo può essere utilizzato per trovare la soluzione di equazioni di punto fisso $f(x)=x$ che emergono spesso nel campo della scienza computazionale.

Definizione

Data una funzione $f:\mathbb {R} ^{n}\to \mathbb {R} ^{n}$, si consideri il problema di trovare un punto fisso di $f$, ovvero una soluzione all'equazione $f(x)=x$ . Un approccio classico al problema è quello di impiegare uno schema di iterazione di punto fisso, ovvero data una prima ipotesi $x_{0}$ per la soluzione, calcolare la sequenza $x_{i+1}=f(x_{i})$ fino a quando non viene soddisfatto un criterio di convergenza. Tuttavia, la convergenza dello schema non è garantita in generale. Inoltre, la convergenza è in generale lineare, e può quindi risultare troppo lenta, soprattutto se la valutazione della funzione $f$ è costosa. L'accelerazione di Anderson è un metodo per accelerare la convergenza delle iterazioni di punto fisso.

Si definisca $g(x)=f(x)-x$, e $g_{k}=g(x_{k})$. Data un'iterata iniziale $x_{0}$ e un parametro intero $m\geq 1$, il metodo può essere formulato come segue:

$x_{1}=f(x_{0})$

$\forall k=1,2,\dots$

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

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

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

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

Per terminare le iterazioni del metodo possono essere usati i criteri di arresto tipici dei metodi iterativi. Ad esempio, le iterazioni possono essere interrotte quando $||x_{k+1}-x_{k}||$ è al di sotto di una tolleranza prescritta o quando il residuo $g(x_{k})$ è al di sotto di una tolleranza prescritta.

Rispetto alla semplice iterazione di punto fisso, è stato riscontrato che il metodo di accelerazione di Anderson converge più velocemente, più robusto e, in alcuni casi, evita la divergenza della sequenza di punto fisso.

Soluzione del problema di minimizzazione

Ad ogni iterazione dell'algoritmo, il problema di ottimizzazione vincolata $\operatorname {argmin} ||G_{k}\alpha ||_{2}$, soggetto a $\alpha \in A_{k}$ deve essere risolto. Il problema può essere espresso in alcune formulazioni equivalenti, cui corrispondono diversi metodi di soluzione che possono risultare in un'implementazione più efficiente:

• definendo le matrici ${\mathcal {G}}_{k}={\begin{bmatrix}g_{k-m_{k}+1}-g_{k-m_{k}}&\dots &g_{k}-g_{k-1}\end{bmatrix}}$ e ${\mathcal {X}}_{k}={\begin{bmatrix}x_{k-m_{k}+1}-x_{k-m_{k}}&\dots &x_{k}-x_{k-1}\end{bmatrix}}$, risolvere $\gamma _{k}=\operatorname {argmin} _{\gamma \in \mathbb {R} ^{m_{k}}}||g_{k}-{\mathcal {G}}_{k}\gamma ||_{2}$ e porre $x_{k+1}=x_{k}+g_{k}-({\mathcal {X}}_{k}+{\mathcal {G}}_{k})\gamma _{k}$ ;
• risolvere $\theta _{k}=\{(\theta _{k})_{i}\}_{i=1}^{m_{k}}=\operatorname {argmin} _{\theta \in \mathbb {R} ^{m_{k}}}||g_{k}+\sum _{i=1}^{m_{k}}\theta _{i}(g_{k-i}-g_{k})||_{2}$, quindi porre $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})$.

Per entrambe le scelte precedenti, il problema di ottimizzazione si presenta sotto forma di un problema ai minimi quadrati lineare non vincolato, che può essere risolto con metodi standard tra cui decomposizione QR e decomposizione ai valori singolari, possibilmente includendo tecniche di regolarizzazione per affrontare problemi legati al rango e al condizionamento del problema di ottimizzazione. Risolvere il problema dei minimi quadrati risolvendo le equazioni normali non è generalmente consigliabile a causa di potenziali instabilità numeriche e costi di calcolo generalmente elevati.

La stagnazione del metodo (ovvero il calcolo di iterazioni successive con lo stesso valore, $x_{k+1}=x_{k}$) ne causa il breakdown per via della singolarità del problema dei minimi quadrati. Allo stesso modo, la quasi-stagnazione ($x_{k+1}\approx x_{k}$) provoca il mal condizionamento del problema ai minimi quadrati. Inoltre, la scelta del parametro $m$ potrebbe essere rilevante nel determinare il condizionamento del problema dei minimi quadrati, come illustrato più avanti.

Rilassamento

L'algoritmo può essere modificato introducendo un parametro di rilassamento variabile (o parametro di mixing) $\beta _{k}>0$. Ad ogni passaggio, la nuova iterata va calcolata come $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})\;.$ La scelta del valore di $\beta _{k}$ è cruciale per determinare le proprietà di convergenza del metodo; in linea di principio, $\beta _{k}$ potrebbe variare ad ogni iterazione, sebbene sia spesso scelto costante.

Scelta di m

Il parametro $m$ determina quante informazioni dalle iterazioni precedenti vengono utilizzate per calcolare la nuova iterazione $x_{k+1}$. Da un lato, se $m$ viene scelto troppo piccolo, vengono utilizzate poche informazioni e la convergenza può essere lenta. D'altra parte, se $m$ è troppo grande, le informazioni delle vecchie iterazioni possono essere conservate per troppe iterazioni successive, anche in questo caso rallentando la convergenza. Inoltre, la scelta di $m$ influisce sulla dimensione del problema di ottimizzazione. Un valore troppo grande di $m$ può peggiorare il condizionamento del problema dei minimi quadrati e il costo della sua soluzione. In generale, quale sia una buona scelta di $m$ dipende dal particolare problema da risolvere.

Scelta di $m_{k}$

Rispetto all'algoritmo sopra descritto, la scelta di $m_{k}$ ad ogni iterazione può essere modificata. Una possibilità è scegliere $m_{k}=k$ per ogni iterazione $k$ (questa scelta è a volte indicata come accelerazione di Anderson senza troncamento). In questo modo, ogni nuova iterazione $x_{k+1}$ viene calcolata utilizzando tutte le iterazioni precedenti. Una tecnica più sofisticata si basa sullo scegliere $m_{k}$ in modo da mantenere un condizionamento abbastanza piccolo per il problema ai minimi quadrati.

Confronto con altre classi di metodi

Il metodo di Newton può essere applicato alla soluzione di $f(x)-x=0$ per calcolare un punto fisso di $f(x)$ con convergenza quadratica. Tuttavia, tale metodo richiede la valutazione della derivata esatta di $f(x)$, cosa può essere molto costosa. L'approssimazione della derivata mediante differenze finite è una possibile alternativa, ma richiede valutazioni multiple di $f(x)$ ad ogni iterazione, che di nuovo può essere molto costoso. L'accelerazione di Anderson richiede solo una valutazione della funzione $f(x)$ per iterazione e nessuna valutazione della sua derivata. D'altra parte, la convergenza di una sequenza di punto fisso con accelerazione di Anderson è ancora lineare in generale.

Diversi autori hanno messo in evidenza somiglianze tra lo schema di accelerazione di Anderson e altri metodi per la soluzione di equazioni non lineari. In particolare:

• Eyert e Fang e Saad hanno interpretato l'algoritmo all'interno della classe dei metodi quasi-Newton e multisecante, che generalizzano il noto metodo delle secanti, per la soluzione dell'equazione non lineare $g(x)=0$; hanno anche mostrato come lo schema può essere visto come un metodo nella di Broyden;
• Walker e Ni hanno dimostrato che lo schema di accelerazione di Anderson è equivalente al metodo GMRES in caso di problemi lineari (ovvero il problema di trovare una soluzione a $A\mathbf {x} =\mathbf {x}$ per una matrice quadrata $A$ e può quindi essere visto come una generalizzazione di GMRES nel caso non lineare; un risultato simile è stato trovato da Washio e Oosterlee.

Inoltre, altri metodi equivalenti o quasi equivalenti sono stati sviluppati in modo indipendente da altri autori, sebbene il più delle volte nel contesto di un'applicazione specifica di interesse piuttosto che come metodo generale per equazioni di punto fisso.

Esempio di implementazione in MATLAB

Di seguito è riportato un esempio di implementazione in linguaggio MATLAB dello schema di accelerazione di Anderson per trovare il punto fisso della funzione $f(x)=\sin(x)+\arctan(x)$. Si noti che:

• il problema di ottimizzazione è stato risolto nella forma $\gamma _{k}=\operatorname {argmin} _{\gamma \in \mathbb {R} ^{m_{k}}}||g_{k}-{\mathcal {G}}_{k}\gamma ||_{2}$ usando la decomposizione QR;
• il calcolo della decomposizione QR non è ottimale: infatti, ad ogni iterazione viene aggiunta una singola colonna alla matrice ${\mathcal {G}}_{k}$ ed eventualmente una singola colonna viene rimossa: questo può essere sfruttato per aggiornare in modo efficiente la decomposizione QR con minori costi computazionali;
• l'algoritmo può essere reso più efficiente in termini di memoria salvando solo le ultime iterazioni e residui, se l'intero vettore di iterazioni $x_{k}$ non è necessario;
• il codice è facilmente generalizzabile al caso di una $f(x)$ a valori vettoriali.

f = @(x) sin(x) + atan(x); % Funzione di cui calcolare il punto fisso.
x0 = 1; % Iterata iniziale.

k_max = 100; % Numero massimo di iterazioni.
tol_res = 1e-6; % Tolleranza sul residuo.
m = 3; % Parametro m.

x = [x0, f(x0)]; % Vettore di tutte le iterazioni.
g = f(x) - x; % Vettore di tutti i residui.

G_k = g(2) - g(1); % Matrice degli incrementi dei residui.
X_k = x(2) - x(1); % Matrice degli incrementi nelle iterazioni.

k = 2;
while k < k_max && abs(g(k)) > tol_res
  m_k = min(k, m);
  
  % Soluzione del problema ai minimi quadrati con fattorizzazione QR.
  [Q, R] = qr(G_k);
  gamma_k = R \ (Q' * g(k));
  
  % Calcolo della nuova iterata e del residuo corrispondente.
  x(k+1) = x(k) + g(k) - (X_k + G_k) * gamma_k;
  g(k+1) = f(x(k+1)) - x(k+1);
  
  % Aggiornamento delle matrici degli incrementi.
  X_k = [X_k, x(k+1) - x(k)];
  G_k = [G_k, g(k+1) - g(k)];
  
  n = size(X_k, 2);
  if n > m_k
    X_k = X_k(:,n-m_k+1:end);
    G_k = G_k(:,n-m_k+1:end);
  end
  
  k = k + 1;
end

% Stampa il risultato: Trovato punto fisso 2.013444 dopo 9 iterazioni
fprintf ("Trovato punto fisso %f dopo %d iterazioni\n", x(end), k);

Voci correlate

• Iterazioni di punto fisso
• Calcolo di uno zero di una funzione