Algoritmo di Gauss-Newton

Regressione di una curva con un modello a picco asimmetrico, utilizzando l'algoritmo di Gauss–Newton con un fattore di smorzamento $\alpha$ variabile.
Sopra: dati grezzi e curva del modello.
Sotto: evoluzione della somma normalizzata dei quadrati dei residui.

L'algoritmo di Gauss–Newton è un metodo iterativo per risolvere problemi di minimi quadrati e regressioni non lineari. È una versione modificata del metodo di Newton per trovare un minimo di una funzione. Diversamente da quest'ultimo, l'algoritmo di Gauss–Newton può essere utilizzato solo per minimizzare una somma di funzioni al quadrato, ma possiede il vantaggio che le derivate seconde, spesso faticose da calcolare, non sono richieste.

I problemi di minimi quadrati non lineari compaiono, ad esempio, nella regressione non lineare, in cui si cercano i parametri tali che il modello sia in buono accordo con le osservazioni disponibili.

Il nome del metodo deriva dai matematici Carl Friedrich Gauss e Isaac Newton.

Descrizione

Date $m$ funzioni ${\boldsymbol {r}}=(r_{1},\ldots ,r_{m})$ (spesso chiamati residui) di $n$ variabili ${\boldsymbol {\beta }}=(\beta _{1},\ldots ,\beta _{n})$, con $m\geq n$, l'algoritmo di Gauss–Newton trova iterativamente i valori delle variabili in modo da minimizzare la seguente somma di quadrati:

$S({\boldsymbol {\beta }})=\sum _{i=1}^{m}r_{i}^{2}({\boldsymbol {\beta }}).$

Iniziando con ${\boldsymbol {\beta }}^{(0)}$ come stima iniziale per il minimo, il metodo esegue iterativamente

${\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}-\left(\mathbf {J_{r}} ^{\mathsf {T}}\mathbf {J_{r}} \right)^{-1}\mathbf {J_{r}} ^{\mathsf {T}}\mathbf {r} ({\boldsymbol {\beta }}^{(s)}),$

dove, se ${\boldsymbol {r}}$ e ${\boldsymbol {\beta }}$ sono vettori colonna, gli elementi della matrice jacobiana sono

$(\mathbf {J_{r}} )_{ij}={\frac {\partial r_{i}({\boldsymbol {\beta }}^{(s)})}{\partial \beta _{j}}},$

e il simbolo $^{\mathsf {T}}$ indica la matrice trasposta.

Se $m=n$, l'iterazione si semplifica e diventa

${\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}-\left(\mathbf {J_{r}} \right)^{-1}\mathbf {r} ({\boldsymbol {\beta }}^{(s)}),$

che è una diretta generalizzazione in più dimensioni del metodo delle tangenti.

Nella regressione dei dati, in cui l'obiettivo è trovare i valori dei parametri ${\boldsymbol {\beta }}$ tali che una data funzione modello $y=f(x,{\boldsymbol {\beta }})$ sia il più possibile in accordo con la serie di punti $(x_{i},y_{i})$, le funzioni $r_{i}$ sono i residui:

$r_{i}({\boldsymbol {\beta }})=y_{i}-f(x_{i},{\boldsymbol {\beta }}).$

Allora, il metodo di Gauss–Newton può essere espresso in termini dello jacobiano ${\boldsymbol {J}}_{f}$ della funzione $f$ come

${\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}+\left(\mathbf {J_{f}} ^{\mathsf {T}}\mathbf {J_{f}} \right)^{-1}\mathbf {J_{f}} ^{\mathsf {T}}\mathbf {r} ({\boldsymbol {\beta }}^{(s)}).$

Da notare che $\left(\mathbf {J_{f}} ^{\mathsf {T}}\mathbf {J_{f}} \right)^{-1}\mathbf {J_{f}} ^{\mathsf {T}}$ è la pseudoinversa di $\mathbf {J_{f}}$. Nell'algoritmo, l'assunzione $m\geq n$ è necessaria, altrimenti la matrice $\mathbf {J_{r}} ^{\mathsf {T}}\mathbf {J_{r}}$ non è invertibile e le equazioni non possono essere risolte (almeno in modo unico).

L'algoritmo di Gauss–Newton si ottiene dall'approssimazione lineare del vettore delle funzioni $r_{i}$ utilizzando il teorema di Taylor. Infatti, a ogni iterazione si ottiene:

$\mathbf {r} ({\boldsymbol {\beta }})\approx \mathbf {r} ({\boldsymbol {\beta }}^{(s)})+\mathbf {J_{r}} ({\boldsymbol {\beta }}^{(s)})\Delta$

con $\Delta ={\boldsymbol {\beta }}-{\boldsymbol {\beta }}^{(s)}$. Trovare $\Delta$ che minimizza la somma dei quadrati nel membro destro, cioè

$\min \left\|\mathbf {r} ({\boldsymbol {\beta }}^{(s)})+\mathbf {J_{r}} ({\boldsymbol {\beta }}^{(s)})\Delta \right\|_{2}^{2},$

è un problema dei minimi quadrati lineare, che si risolve esplicitamente.

Le equazioni normali sono $n$ equazioni lineari simultanee nell'incremento $\Delta$ incognito. Si possono risolvere in un solo passaggio, usando la decomposizione di Cholesky, o, ancora meglio, la fattorizzazione QR di $\mathbf {J_{r}}$. Per sistemi grandi, può essere più efficiente un metodo iterativo, come quello del gradiente coniugato. Se esiste una dipendenza lineare tra le colonne di $\mathbf {J_{r}}$, le iterazioni falliranno a causa della singolarità di $\mathbf {J_{r}} ^{\mathsf {T}}\mathbf {J_{r}}$.

Esempio

Curva di best fit ottenuta (in blu), con ${\hat {\beta }}_{1}=0.362$ and ${\hat {\beta }}_{2}=0.556$, insieme ai dati osservati (in rosso).

In questo esempio, l'algoritmo di Gauss–Newton viene usato per la regressione della velocità $V$ di formazione del prodotto in una reazione catalizzata da un enzima in relazione alla concentrazione del substrato $[S]$, secondo il modello di Michaelis-Menten. I dati misurati sono riportati nella seguente tabella. Le incertezze di ogni misura sono state messe uguali a 1.

$i$	1	2	3	4	5	6	7
$[S]$	0.038	0.194	0.425	0.626	1.253	2.500	3.740
V	0.050	0.127	0.094	0.2122	0.2729	0.2665	0.3317

La funzione modello è della forma

$V={\frac {V_{\text{max}}[S]}{K_{M}+[S]}}$

con parametri $V_{\text{max}}$ e $K_{M}$ da determinare attraverso l'algoritmo.

Siano $x_{i}$ e $y_{i}$ i valori di $[S]$ e $V$ rispettivamente nella tabella, con $i=1,\dots ,7$. Siano $\beta _{1}=V_{\text{max}}$ e $\beta _{2}=K_{M}$. Si troveranno $\beta _{1}$ e $\beta _{2}$ tali che la somma dei quadrati dei residui

$r_{i}=y_{i}-{\frac {\beta _{1}x_{i}}{\beta _{2}+x_{i}}}\quad (i=1,\dots ,7)$

sia minima.

Lo jacobiano $\mathbf {J_{r}}$ del vettore dei residui $r_{i}$ rispetto alle incognite $\beta _{j}$ è una matrice $7\times 2$ in cui nell'$i$-esima riga si trova

${\frac {\partial r_{i}}{\partial \beta _{1}}}=-{\frac {x_{i}}{\beta _{2}+x_{i}}};{\frac {\partial r_{i}}{\partial \beta _{2}}}={\frac {\beta _{1}x_{i}}{(\beta _{2}+x_{i})^{2}}}.$

Iniziando con una stima iniziale $\beta _{1}^{(0)}=0.9$ e $\beta _{2}^{(0)}=0.2$, dopo cinque iterazioni dell'algoritmo di Gauss–Newton, si ottengono i valori ottimali ${\hat {\beta }}_{1}=0.362$ e ${\hat {\beta }}_{2}=0.556$. La somma dei quadrati dei residui descresce dal valore iniziali di $1.445$ a quello finale di $0.00784$. Il grafico in figura mostra i dati nella tabella insieme alla curva modello con i parametri ottimali ottenuti dall'algoritmo. Sotto è riportata una tabella dei valori intermedi di $\beta _{1}$ e $\beta _{2}$ durante l'algoritmo.

Iterazione $i$	$\beta _{1}^{(i)}$	$\beta _{2}^{(i)}$	$S(\mathbf {\beta ^{(i)}} )$
1	0.33266293	0.26017391	0.015072
2	0.34280925	0.42607918	0.008458
3	0.35777522	0.52950844	0.007864
4	0.36140546	0.5536581	0.007844
5	0.36180308	0.55607253	0.007844
6	0.36183442	0.55625246	0.007844

Convergenza del metodo

Si può mostrare che l'incremento $\Delta$ è una direzione di discesa per $S$, e, se l'algoritmo converge, che il limite è un punto stazionario di $S$. Tuttavia, la convergenza non è garantita, nemmeno quella locale come nel metodo delle tangenti, o sotto le comuni condizioni di Wolfe.

La velocità di convergenza di Gauss–Newton può diventare quadratica. L'algoritmo potrebbe anche convergere lentamente o affatto se la stima iniziale è lontana dal minimo oppure la matrice $\mathbf {J_{r}^{\mathsf {T}}J_{r}}$ è mal condizionata. Per esempio, si consideri il problema con $m=2$ equazioni e $n=1$ variabili, dato da

${\begin{aligned}r_{1}(\beta )&=\beta +1,\\r_{2}(\beta )&=\lambda \beta ^{2}+\beta -1.\end{aligned}}$

Il minimo è per $\beta =0$. (In realtà il minimo è per $\beta =-1$ se $\lambda =2$, poiché $S(0)=1^{2}+(-1)^{2}=2$, ma $S(-1)=0$.) Se $\lambda =0$, allora il problema diventa lineare e il metodo trova il minimo in una sola iterazione. Se $|\lambda |<1$, allora l'algoritmo converge linearmente e l'errore decresce asintoticamente di un fattore $|\lambda |$ a ogni iterazione. Tuttavia, se $|\lambda |>1$, non c'è convergenza nemmeno locale.

Derivazione dal metodo di Newton

In questo paragrafo, l'algoritmo di Gauss–Newton verrà derivato dal metodo di Newton per l'ottimizzazione di funzione. Come conseguenza, la velocità di convergenza dell'algoritmo di Gauss–Newton può essere quadratico sotto certe condizioni di regolarità. In generale (sotto condizioni più deboli), la convergenza è lineare.

La relazione di ricorrenza per il metodo di Newton per la minimizzazione della funzione $S$ di parametri ${\boldsymbol {\beta }}$ è

${\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}-\mathbf {H} ^{-1}\mathbf {g} ,$

dove $\mathbf {g}$ indica il vettore gradiente di $S$, e $\mathbf {H}$ la sua matrice hessiana. Poiché $S=\sum _{i=1}^{m}r_{i}^{2}$, il gradiente è dato da

$g_{j}=2\sum _{i=1}^{m}r_{i}{\frac {\partial r_{i}}{\partial \beta _{j}}}.$

Elementi dell'Hessiana si calcolano derivando le componenti del gradiente, $g_{j}$, rispetto a $\beta _{k}$:

$H_{jk}=2\sum _{i=1}^{m}\left({\frac {\partial r_{i}}{\partial \beta _{j}}}{\frac {\partial r_{i}}{\partial \beta _{k}}}+r_{i}{\frac {\partial ^{2}r_{i}}{\partial \beta _{j}\partial \beta _{k}}}\right).$

Il metodo di Gauss–Newton si ottiene trascurando i termini con le derivate seconde (il secondo nell'espressione). Cioè, la matrice Hessiana è approssimata come

$H_{jk}\approx 2\sum _{i=1}^{m}J_{ij}J_{ik},$

dove $J_{ij}={\frac {\partial r_{i}}{\partial \beta _{j}}}$ sono gli elementi del jacobiano $\mathbf {J_{r}}$. Si possono riscrivere il gradiente e l'Hessiana approssimata in notazione matriciale come

$\mathbf {g} =2\mathbf {J} _{\mathbf {r} }^{\mathsf {T}}\mathbf {r} ,\quad \mathbf {H} \approx 2\mathbf {J} _{\mathbf {r} }^{\mathsf {T}}\mathbf {J_{r}} .$

Si sostituiscono queste espressioni nella relazione di ricorrenza precedente, così da ottenere l'equazioni dell'algoritmo

${\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}+\Delta ;\quad \Delta =-\left(\mathbf {J_{r}} ^{\mathsf {T}}\mathbf {J_{r}} \right)^{-1}\mathbf {J_{r}} ^{\mathsf {T}}\mathbf {r} .$

La convergenza del metodo di Gauss–Newton non è garantita in tutte le situazioni. L'approssimazione

$\left|r_{i}{\frac {\partial ^{2}r_{i}}{\partial \beta _{j}\partial \beta _{k}}}\right|\ll \left|{\frac {\partial r_{i}}{\partial \beta _{j}}}{\frac {\partial r_{i}}{\partial \beta _{k}}}\right|,$

che serve per trascurare le derivate seconde può essere valida in due casi, così da aspettarsi la convergenza dell'algoritmo:

1. I valori della funzione $r_{i}$ sono piccoli, almeno intorno al minimo.
2. Le funzioni sono quasi-lineari, in modo che ${\frac {\partial ^{2}r_{i}}{\partial \beta _{j}\partial \beta _{k}}}$ sia relativamente piccolo.

Versioni migliorate dell'algoritmo

Con l'algoritmo di Gauss–Newton, la somma dei quadrati dei residui $S$ può non decrescere a ogni interazione. Tuttavia, poiché $\Delta$ è una direzione di discesa, a meno che $S({\boldsymbol {\beta }}^{s})$ sia un punto stazionario, vale che $S({\boldsymbol {\beta }}^{s}+\alpha \Delta )<S({\boldsymbol {\beta }}^{s})$ per ogni $\alpha >0$ sufficientemente piccolo. Quindi, se il metodo diverge, una soluzione è di impiegare una frazione $\alpha$ dell'incremento $\Delta$, utilizzando la seguente formula:

${\boldsymbol {\beta }}^{s+1}={\boldsymbol {\beta }}^{s}+\alpha \Delta .$.

In altre parole, il vettore incremento è troppo lungo, ma è diretto verso il basso, perciò avanzare soltanto di una frazione di cammino diminuirà il valore della funzione oggettiva $S$. Si può trovare il valore ottimale di $\alpha$ usando un algoritmo di line search, cioè il valore $\alpha$ viene determinato trovando quello che minimizza $S$, di solito con un metodo di ricerca diretta nell'intervallo $0<\alpha <1$.

Nei casi in cui la frazione ottimale $\alpha$ è vicina a zero, un metodo alternativo per il trattamento della divergenza è l'utilizzo dell'algoritmo di Levenberg-Marquardt, anche conosciuto come "metodo della regione di confidenza". Le equazioni normali sono modificate in modo che l'incremento sia rotato verso la direzione di massima decrescita,

$\left(\mathbf {J^{\mathrm {T} }J+\lambda D} \right)\Delta =-\mathbf {J} ^{\mathrm {T} }\mathbf {r} ,$

dove $\mathbf {D}$ è una matrice diagonale positiva. Da notare che quando $\mathbf {D}$ è la matrice identità $\mathbf {I}$ e $\lambda \to +\infty$, allora $\lambda \Delta =\lambda \left(\mathbf {J^{\mathrm {T} }J} +\lambda \mathbf {I} \right)^{-1}\left(-\mathbf {J} ^{\mathrm {T} }\mathbf {r} \right)=\left(\mathbf {I} -\mathbf {J^{\mathrm {T} }J} /\lambda +\cdots \right)\left(-\mathbf {J} ^{\mathrm {T} }\mathbf {r} \right)\to -\mathbf {J} ^{\mathrm {T} }\mathbf {r}$, perciò la direzione di $\Delta$ si avvicina alla direzione del gradiente negativo $-\mathbf {J} ^{\mathrm {T} }\mathbf {r}$.

Il parametro di Marquardt $\lambda$ può essere ottimizzato attraverso un line search, ma è molto inefficiente, dato che il vettore incremento deve essere ricalcolato a ogni modifica di $\lambda$. Una strategia più efficiente è questa: quando il metodo diverge, si incrementa il parametro di Marquardt fintanto che si ha una decrescita di $S$. Quindi si conserva il valore da una iterazione a quella successiva, ma si diminuisce fino a che non si raggiunge un valore limite, quando il parametro di Marquardt può essere posto uguale a 0; la minimizzazione di $S$ diventa pertanto una ottimizzazione standard di Gauss–Newton.

Ottimizzazione su larga scala

Per l'ottimizzazione su larga scala, l'algoritmo di Gauss–Newton è di particolare interesse perché in generale vale (sebbene non sempre) che la matrice $\mathbf {J} _{\mathbf {r} }$ è molto più sparsa dell'Hessiana approssimata $\mathbf {J} _{\mathbf {r} }^{\mathsf {T}}\mathbf {J_{r}}$. In questi casi, il passo dell'algoritmo viene fatto con un metodo iterativo approssimato adatto a problemi grandi e sparsi, come il metodo del gradiente coniugato.

Per far funzionare questo approccio, serve almeno un metodo efficiente per calcolare il calcolare il prodotto

$\mathbf {J} _{\mathbf {r} }^{\mathsf {T}}\mathbf {J_{r}} \mathbf {p}$

per un qualche vettore $\mathbf {p}$. Per l'immagazzinamento della matrice sparsa, in genere è pratico memorizzare le righe di $\mathbf {J} _{\mathbf {r} }$ in una forma compressa (cioè, senza gli elementi nulli), rendendo però il calcolo diretto del prodotto precedente alquanto complicato per via della trasposizione. Tuttavia, se si definisce $\mathbf {c_{i}}$ come la riga $i$-esima della matrice $\mathbf {J} _{\mathbf {r} }$, vale la seguente semplice relazione:

$\mathbf {J} _{\mathbf {r} }^{\mathsf {T}}\mathbf {J_{r}} \mathbf {p} =\sum _{i}\mathbf {c} _{i}(\mathbf {c} _{i}\cdot \mathbf {p} ),$

in modo che ogni riga contribuisce additivamente e indipendentemente al prodotto. In aggiunta alla memorizzazione molto pratica, questa espressione è adatta al calcolo parallelo. Da notare che ogni riga $\mathbf {c_{i}}$ è il gradiente del rispettivo residuo $r_{i}$; tenendo conto di questo, la forma precedente enfatizza il fatto che i residui contribuiscono al problema in modo indipendente uno dall'altro.

Algoritmi collegati

In un metodo quasi-Newton, come quello dovuto a Davidon, Fletcher e Powell oppure a Broyden–Fletcher–Goldfarb–Shanno (metodo BFGS), si calcola numericamente una stima della Hessiana ${\frac {\partial ^{2}S}{\partial \beta _{j}\partial \beta _{k}}}$ usando solo le derivate prime ${\frac {\partial r_{i}}{\partial \beta _{j}}}$, in modo che solo dopo $n$ cicli di perfezionamento il metodo si avvicina approssimativamente a quello di Newton in termini di prestazioni. Da notare che i metodi quasi-Newton possono minimizzare funzioni arbitrarie a valori reali, mentre Gauss–Newton, Levenberg–Marquardt, ecc. risolvono solo problemi di minimi quadrati non lineari.

Un altro metodo per risolvere problemi di minimo usando solo derivate prime è la discesa del gradiente. Tuttavia, quest'ultimo metodo non considera le derivate seconde nemmeno approssimativamente, perciò è altamente inefficiente per molte funzioni, specialmente se i parametri hanno una forte correlazione.

Bibliografia

• A. Björck, Numerical methods for least squares problems, SIAM, Philadelphia, 1996, ISBN 0-89871-360-9.
• Roger Fletcher, Practical methods of optimization, 2nd, New York, John Wiley & Sons, 1987, ISBN 978-0-471-91547-8.
• Jorge Nocedal e Wright, Stephen, Numerical optimization, New York: Springer, 1999, ISBN 0-387-98793-2.

Collegamenti esterni

Implementazioni

• Artelys Knitro è un risolutore non lineare con un'implementazione del metodo di Gauss–Newton. È scritto in linguaggio C e possiede interfacce per C++/C#/Java/Python/MATLAB/R.