Apri il menu principale

Algoritmo di Gauss-Newton

Algoritmo di regressione ai minimi quadrati
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, a 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.

DescrizioneModifica

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 l'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.

L'equazione 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  .

EsempioModifica

 
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 metodoModifica

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 NewtonModifica

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'algoritmoModifica

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 scalaModifica

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 collegatiModifica

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.

NoteModifica

  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.

BibliografiaModifica

Collegamenti esterniModifica

ImplementazioniModifica

  • 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