Algoritmo di Levenberg-Marquardt

L'algoritmo di Levenberg-Marquardt (LMA) è un algoritmo di ottimizzazione usato per la soluzione di problemi in forma di minimi quadrati non lineari, che trova comunemente applicazioni in problemi di curve fitting. LMA è un algoritmo iterativo, nel quale il vettore di aggiornamento della soluzione ad ogni iterazione è dato da un'interpolazione fra l'algoritmo di Gauss-Newton e il metodo di discesa del gradiente. LMA può essere considerato come una versione trust region dell'algoritmo di Gauss-Newton, rispetto al quale è più robusto ma, in generale, leggermente più lento. L'algoritmo è stato pubblicato nel 1944 da Kenneth Levenberg,[1] e fu riscoperto nel 1963 da Donald Marquardt[2] e, indipendentemente, da Girard,[3] Wynne[4] e Morrison.[5]

Formulazione

L'applicazione principale dell'algoritmo di Levenberg-Marquardt è il problema di curve fitting tramite minimi quadrati non lineari. Dato un insieme di m {\displaystyle m} osservazioni ( x i , y i ) {\displaystyle \left(x_{i},y_{i}\right)} , si vuole determinare il vettore di parametri β ^ {\displaystyle {\hat {\boldsymbol {\beta }}}} del modello f ( x , β ) {\displaystyle f\left(x,{\boldsymbol {\beta }}\right)} che minimizza la somma dei quadrati residui S ( β ) {\displaystyle S\left({\boldsymbol {\beta }}\right)}

β ^ ∈ argmin β ⁡ S ( β ) ≡ argmin β ⁡ ∑ i = 1 m [ y i − f ( x i , β ) ] 2 . {\displaystyle {\hat {\boldsymbol {\beta }}}\in \operatorname {argmin} \limits _{\boldsymbol {\beta }}S\left({\boldsymbol {\beta }}\right)\equiv \operatorname {argmin} \limits _{\boldsymbol {\beta }}\sum _{i=1}^{m}\left[y_{i}-f\left(x_{i},{\boldsymbol {\beta }}\right)\right]^{2}.}

L'algoritmo di Levenberg-Marquardt è un metodo iterativo che parte da una stima iniziale del vettore β {\displaystyle {\boldsymbol {\beta }}} . Nel caso di funzioni non-convesse con più minimi locali, la scelta di una stima iniziale sufficientemente vicina al punto di ottimo globale è importante per la convergenza. Ad ogni iterazione, la stima corrente della soluzione β {\displaystyle {\boldsymbol {\beta }}} viene aggiornata ad un nuovo valore β + δ {\displaystyle {\boldsymbol {\beta }}+{\boldsymbol {\delta }}} . Per determinare la scelta di δ {\displaystyle {\boldsymbol {\delta }}} , la funzione f {\displaystyle f} viene linearizzata con un polinomio di Taylor

f ( x i , β + δ ) ≈ f ( x i , β ) + J i δ , {\displaystyle f\left(x_{i},{\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)\approx f\left(x_{i},{\boldsymbol {\beta }}\right)+\mathbf {J} _{i}{\boldsymbol {\delta }},}

dove

J i = ∂ f ( x i , β ) ∂ β {\displaystyle \mathbf {J} _{i}={\frac {\partial f\left(x_{i},{\boldsymbol {\beta }}\right)}{\partial {\boldsymbol {\beta }}}}}

è il gradiente di f {\displaystyle f} rispetto a β {\displaystyle {\boldsymbol {\beta }}} .

Usando tale approssimazione, la somma dei quadrati residui S ( β ) {\displaystyle S\left({\boldsymbol {\beta }}\right)} diventa

S ( β + δ ) ≈ ∑ i = 1 m [ y i − f ( x i , β ) − J i δ ] 2 , {\displaystyle S\left({\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)\approx \sum _{i=1}^{m}\left[y_{i}-f\left(x_{i},{\boldsymbol {\beta }}\right)-\mathbf {J} _{i}{\boldsymbol {\delta }}\right]^{2},}

o, in notazione vettoriale

S ( β + δ ) ≈ ‖ y − f ( β ) − J δ ‖ 2 = [ y − f ( β ) − J δ ] T [ y − f ( β ) − J δ ] = [ y − f ( β ) ] T [ y − f ( β ) ] − [ y − f ( β ) ] T J δ − ( J δ ) T [ y − f ( β ) ] + δ T J T J δ = [ y − f ( β ) ] T [ y − f ( β ) ] − 2 [ y − f ( β ) ] T J δ + δ T J T J δ . {\displaystyle {\begin{aligned}S\left({\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)&\approx \left\|\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)-\mathbf {J} {\boldsymbol {\delta }}\right\|^{2}\\&=\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)-\mathbf {J} {\boldsymbol {\delta }}\right]^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)-\mathbf {J} {\boldsymbol {\delta }}\right]\\&=\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]-\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]^{\mathrm {T} }\mathbf {J} {\boldsymbol {\delta }}-\left(\mathbf {J} {\boldsymbol {\delta }}\right)^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]+{\boldsymbol {\delta }}^{\mathrm {T} }\mathbf {J} ^{\mathrm {T} }\mathbf {J} {\boldsymbol {\delta }}\\&=\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]-2\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]^{\mathrm {T} }\mathbf {J} {\boldsymbol {\delta }}+{\boldsymbol {\delta }}^{\mathrm {T} }\mathbf {J} ^{\mathrm {T} }\mathbf {J} {\boldsymbol {\delta }}.\end{aligned}}}

La somma dei quadrati residui S ( β ) {\displaystyle S\left({\boldsymbol {\beta }}\right)} ha un minimo in un punto dove il gradiente rispetto al vettore dei parametri si annulla. Derivando l'espressione precedente rispetto a δ {\displaystyle {\boldsymbol {\delta }}} ed imponendo l'uguaglianza a zero, si ottiene

( J T J ) δ = J T [ y − f ( β ) ] , {\displaystyle \left(\mathbf {J} ^{\mathrm {T} }\mathbf {J} \right){\boldsymbol {\delta }}=\mathbf {J} ^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right],}

dove J {\displaystyle \mathbf {J} } è la matrice jacobiana, la cui riga i {\displaystyle i} -esima è data da J i {\displaystyle \mathbf {J} _{i}} , e dove f ( β ) {\displaystyle \mathbf {f} \left({\boldsymbol {\beta }}\right)} e y {\displaystyle \mathbf {y} } sono vettori le cui righe i {\displaystyle i} -esime sono date rispettivamente da f ( x i , β ) {\displaystyle f\left(x_{i},{\boldsymbol {\beta }}\right)} e y i {\displaystyle y_{i}} . La matrice jacobiana ha dimensione m × n {\displaystyle m\times n} , dove n {\displaystyle n} è il numero di parametri, ovvero la dimensione del vettore β {\displaystyle {\boldsymbol {\beta }}} , e il prodotto ( J T J ) {\displaystyle \left(\mathbf {J} ^{\mathrm {T} }\mathbf {J} \right)} è una matrice quadrata di dimensione n × n {\displaystyle n\times n} .

Risolvendo tale sistema lineare rispetto a δ {\displaystyle {\boldsymbol {\delta }}} si ottiene il vettore di aggiornamento della soluzione secondo il metodo di Gauss-Newton. L'idea originale di Levenberg è di sostituire la precedente equazione con una versione smorzata

( J T J + λ I ) δ = J T [ y − f ( β ) ] , {\displaystyle \left(\mathbf {J} ^{\mathrm {T} }\mathbf {J} +\lambda \mathbf {I} \right){\boldsymbol {\delta }}=\mathbf {J} ^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right],}

dove I {\displaystyle \mathbf {I} } è la matrice identità. Il fattore λ {\displaystyle \lambda } determina il comportamento dell'algoritmo, e un valore ridotto corrisponde ad un comportamento prossimo al metodo di Gauss-Newton, mentre un valore elevato corrisponde a spostare la soluzione in direzione pressappoco opposta al gradiente, con un comportamento più simile al metodo di discesa del gradiente. Il valore viene adattato ad ogni iterazione, incrementandolo se la precedente iterazione ha prodotto una riduzione limitata della funzione obiettivo, o diminuendolo in caso di rapida diminuzione.

Uno degli svantaggi della formulazione di Levenberg è il fatto che il termine J T J + λ I {\displaystyle \mathbf {J} ^{\text{T}}\mathbf {J} +\lambda \mathbf {I} } è praticamente ignorato quando il parametro di smorzamento λ {\displaystyle \lambda } ha un valore elevato. Una variante proposta da Fletcher[6] sostituisce la matrice identità con la diagonale di J T J {\displaystyle \mathbf {J} ^{\text{T}}\mathbf {J} } , scalando ogni parametro rispetto alla curvatura e di conseguenza aumentando la velocità di convergenza lungo le direzioni nelle quali il gradiente è minore:

[ J T J + λ diag ⁡ ( J T J ) ] δ = J T [ y − f ( β ) ] . {\displaystyle \left[\mathbf {J} ^{\mathrm {T} }\mathbf {J} +\lambda \operatorname {diag} \left(\mathbf {J} ^{\mathrm {T} }\mathbf {J} \right)\right]{\boldsymbol {\delta }}=\mathbf {J} ^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right].}

Esistono diverse euristiche per la scelta del parametro di smorzamento λ {\displaystyle \lambda } . Marquardt suggerì di usare una scelta iniziale λ 0 {\displaystyle \lambda _{0}} e un fattore di aggiornamento ν {\displaystyle \nu } , e di calcolare la funzione obiettivo dopo un'iterazione dal valore iniziale ponendo λ = λ 0 {\displaystyle \lambda =\lambda _{0}} , e per un'iterazione dal valore iniziale con λ = λ 0 ν {\displaystyle \lambda ={\frac {\lambda _{0}}{\nu }}} . Se uno dei due valori produce un miglioramento maggiore della funzione costo rispetto all'altro, viene usato come nuovo valore di λ {\displaystyle \lambda } . Se in entrambi i casi la funzione costo ha un valore superiore a quello iniziale, λ {\displaystyle \lambda } è moltiplicato per ν {\displaystyle \nu } iterativamente k {\displaystyle k} volte, fino a quando non si ottiene un valore migliore, ponendo quindi λ = λ 0 ν k {\displaystyle \lambda =\lambda _{0}\nu ^{k}} .

Note

  1. ^ Kenneth Levenberg, A Method for the Solution of Certain Non-Linear Problems in Least Squares, in Quarterly of Applied Mathematics, vol. 2, n. 2, 1944, pp. 164–168, DOI:10.1090/qam/10666.
  2. ^ Donald Marquardt, An Algorithm for Least-Squares Estimation of Nonlinear Parameters, in SIAM Journal on Applied Mathematics, vol. 11, n. 2, 1963, pp. 431–441, DOI:10.1137/0111030.
  3. ^ André Girard, Excerpt from Revue d'optique théorique et instrumentale, in Rev. Opt., vol. 37, 1958, pp. 225–241, 397–424.
  4. ^ C. G. Wynne, Lens Designing by Electronic Digital Computer: I, in Proc. Phys. Soc. Lond., vol. 73, n. 5, 1959, pp. 777–787, Bibcode:1959PPS....73..777W, DOI:10.1088/0370-1328/73/5/310.
  5. ^ David D. Morrison, Methods for nonlinear least squares problems and convergence proofs, in Proceedings of the Jet Propulsion Laboratory Seminar on Tracking Programs and Orbit Determination, 1960, pp. 1–9.
  6. ^ Roger Fletcher, A modified Marquardt subroutine for non-linear least squares (technical report), Harwell, Atomic Energy Research Establishment, 1971.

Bibliografia

  • Jorge J. Moré e Daniel C. Sorensen, Computing a Trust-Region Step (PDF), in SIAM J. Sci. Stat. Comput., vol. 4, n. 3, 1983, pp. 553–572, DOI:10.1137/0904038.
  • Philip E. Gill e Walter Murray, Algorithms for the solution of the nonlinear least-squares problem, in SIAM Journal on Numerical Analysis, vol. 15, n. 5, 1978, pp. 977–992, Bibcode:1978SJNA...15..977G, DOI:10.1137/0715063.
  • Jorge Nocedal e Stephen J. Wright, Numerical Optimization, 2nd, Springer, 2006, ISBN 978-0-387-30303-1.
  • T. Strutz: Data Fitting and Uncertainty (A practical introduction to weighted least squares and beyond). 2nd edition, Springer Vieweg, 2016, ISBN 978-3-658-11455-8.
  • C. T. Kelley, Iterative Methods for Optimization, SIAM Frontiers in Applied Mathematics, no 18, 1999, ISBN 0-89871-433-8.

Collegamenti esterni

  • Numerical Recipes in C, Chapter 15.5: Nonlinear models (archiviato dall'url originale il 4 marzo 2012).
  • History of the algorithm in SIAM news (archiviato dall'url originale il 1º marzo 2014).
  • A tutorial by Ananth Ranganathan (PDF).
  • Methods for Non-Linear Least Squares Problems (PDF).
  • H. P. Gavin, The Levenberg-Marquardt method for nonlinear least squares curve-fitting problems (PDF). (con implementazione in MATLAB)
  Portale Matematica: accedi alle voci di Wikipedia che trattano di matematica