Algoritmo di Gauss-Newton

Da Wikipedia, l'enciclopedia libera.
Vai alla navigazione Vai alla ricerca
Regressione di una curva con un modello a picco asimmetrico, utilizzando l'algoritmo di Gauss–Newton con un fattore di smorzamento 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.

Date funzioni (spesso chiamati residui) di variabili , con , l'algoritmo di Gauss–Newton trova iterativamente i valori delle variabili in modo da minimizzare la seguente somma di quadrati:[1]

Iniziando con come stima iniziale per il minimo, il metodo esegue iterativamente

dove, se e sono vettori colonna, gli elementi della matrice jacobiana sono

e il simbolo indica la matrice trasposta.

Se , l'iterazione si semplifica e diventa

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

Nella regressione dei dati, in cui l'obiettivo è trovare i valori dei parametri tali che una data funzione modello sia il più possibile in accordo con la serie di punti , le funzioni sono i residui:

Allora, il metodo di Gauss–Newton può essere espresso in termini dello jacobiano della funzione come

Da notare che è la pseudoinversa di . Nell'algoritmo, l'assunzione è necessaria, altrimenti la matrice 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 utilizzando il teorema di Taylor. Infatti, a ogni iterazione si ottiene:

con . Trovare che minimizza la somma dei quadrati nel membro destro, cioè

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

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

Curva di best fit ottenuta (in blu), con and , insieme ai dati osservati (in rosso).

In questo esempio, l'algoritmo di Gauss–Newton viene usato per la regressione della velocità di formazione del prodotto in una reazione catalizzata da un enzima in relazione alla concentrazione del substrato , 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.

1 2 3 4 5 6 7
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

con parametri e da determinare attraverso l'algoritmo.

Siano e i valori di e rispettivamente nella tabella, con . Siano e . Si troveranno e tali che la somma dei quadrati dei residui

sia minima.

Lo jacobiano del vettore dei residui rispetto alle incognite è una matrice in cui nell'-esima riga si trova

Iniziando con una stima iniziale e , dopo cinque iterazioni dell'algoritmo di Gauss–Newton, si ottengono i valori ottimali e . La somma dei quadrati dei residui descresce dal valore iniziali di a quello finale di . 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 e durante l'algoritmo.

Iterazione
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

[modifica | modifica wikitesto]

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

La velocità di convergenza di Gauss–Newton può diventare quadratica.[4] L'algoritmo potrebbe anche convergere lentamente o affatto se la stima iniziale è lontana dal minimo oppure la matrice è mal condizionata. Per esempio, si consideri il problema con equazioni e variabili, dato da

Il minimo è per . (In realtà il minimo è per se , poiché , ma .) Se , allora il problema diventa lineare e il metodo trova il minimo in una sola iterazione. Se , allora l'algoritmo converge linearmente e l'errore decresce asintoticamente di un fattore a ogni iterazione. Tuttavia, se , non c'è convergenza nemmeno locale.[5]

Derivazione dal metodo di Newton

[modifica | modifica wikitesto]

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.[6]

La relazione di ricorrenza per il metodo di Newton per la minimizzazione della funzione di parametri è

dove indica il vettore gradiente di , e la sua matrice hessiana. Poiché , il gradiente è dato da

Elementi dell'Hessiana si calcolano derivando le componenti del gradiente, , rispetto a :

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

dove sono gli elementi del jacobiano . Si possono riscrivere il gradiente e l'Hessiana approssimata in notazione matriciale come

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

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

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

  1. I valori della funzione sono piccoli, almeno intorno al minimo.
  2. Le funzioni sono quasi-lineari, in modo che sia relativamente piccolo.

Versioni migliorate dell'algoritmo

[modifica | modifica wikitesto]

Con l'algoritmo di Gauss–Newton, la somma dei quadrati dei residui può non decrescere a ogni interazione. Tuttavia, poiché è una direzione di discesa, a meno che sia un punto stazionario, vale che per ogni sufficientemente piccolo. Quindi, se il metodo diverge, una soluzione è di impiegare una frazione dell'incremento , utilizzando la seguente formula:

.

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 . Si può trovare il valore ottimale di usando un algoritmo di line search, cioè il valore viene determinato trovando quello che minimizza , di solito con un metodo di ricerca diretta nell'intervallo .

Nei casi in cui la frazione ottimale è 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".[1] Le equazioni normali sono modificate in modo che l'incremento sia rotato verso la direzione di massima decrescita,

dove è una matrice diagonale positiva. Da notare che quando è la matrice identità e , allora , perciò la direzione di si avvicina alla direzione del gradiente negativo .

Il parametro di Marquardt può essere ottimizzato attraverso un line search, ma è molto inefficiente, dato che il vettore incremento deve essere ricalcolato a ogni modifica di . Una strategia più efficiente è questa: quando il metodo diverge, si incrementa il parametro di Marquardt fintanto che si ha una decrescita di . 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 diventa pertanto una ottimizzazione standard di Gauss–Newton.

Ottimizzazione su larga scala

[modifica | modifica wikitesto]

Per l'ottimizzazione su larga scala, l'algoritmo di Gauss–Newton è di particolare interesse perché in generale vale (sebbene non sempre) che la matrice è molto più sparsa dell'Hessiana approssimata . 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

per un qualche vettore . Per l'immagazzinamento della matrice sparsa, in genere è pratico memorizzare le righe di 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 come la riga -esima della matrice , vale la seguente semplice relazione:

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 è il gradiente del rispettivo residuo ; tenendo conto di questo, la forma precedente enfatizza il fatto che i residui contribuiscono al problema in modo indipendente uno dall'altro.

Algoritmi collegati

[modifica | modifica wikitesto]

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 usando solo le derivate prime , in modo che solo dopo 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.

  1. ^ a b Björck (1996)
  2. ^ Björck (1996), p. 260.
  3. ^ Mascarenhas, The divergence of the BFGS and Gauss Newton Methods, in Mathematical Programming, vol. 147, n. 1, 2013, pp. 253–276, DOI:10.1007/s10107-013-0720-6, arXiv:1309.7922.
  4. ^ Björck (1996), p. 341, 342.
  5. ^ Fletcher (1987), p. 113.
  6. ^ Copia archiviata, su henley.ac.uk. URL consultato il 2 novembre 2018 (archiviato dall'url originale il 4 agosto 2016).
  7. ^ Nocedal (1999), p. 259.

Collegamenti esterni

[modifica | modifica wikitesto]

Implementazioni

[modifica | modifica wikitesto]
  • 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.
  Portale Matematica: accedi alle voci di Wikipedia che trattano di matematica