Levenberg-Marquardt

Gli algoritmi di risoluzione di sistemi non lineari normalmente si possono vedere divisi tra algoritmi di discesa del gradiente o algoritmi di Gauss-Newton. Una versione più recente di questi algoritmi tuttavia, quella proposta da Levenberg-Marquardt, cerca di lavorare nei punti di forza dei due approcci in modo da trarne vantaggio da entrambi.

L'algoritmo di Levenberg Marquardt (LM) è una tecnica di regressione iterativa ormai ritenuta standard per risolvere problemi non lineari multivariabili. Una ottima descrizione dell'algoritmo può essere trovata in [21]. L'algoritmo si può vedere come composto da una fase di discesa del gradiente, lenta ma che converge, seguita da un risolutore di tipo Gauss-Newton, più veloce.

Sia $f$ una funzione non lineare che trasforma un vettore di parametri $\boldsymbol\beta \in \mathbb{R}^m$ in un vettore misura $\hat{\mathbf{y}} = f(\boldsymbol\beta), \hat{\mathbf{y}} \in \mathbb{R}^n$.

Sia $\mathbf{y}$ un vettore di misure sul quale modellare i parametri della funzione, attraverso una minimizzazione di

$$s(\boldsymbol\beta) = \Vert \mathbf{y} - f(\boldsymbol\beta) \Vert^2 = \mathbf{r}^{\top} \mathbf{r}$$

con $\mathbf{r} = \mathbf{y} - f(\boldsymbol\beta)$ residuo dell'errore.

Per problemi non-convessi è necessario fornire una stima iniziale $\boldsymbol\beta_0$ della soluzione, abbastanza vicina al minimo da ricavare. Solitamente si esegue una regressione lineare che minimizza una qualche quantità algebrica, ottenendo un minimo assoluto rispetto a questa metrica, per ottenere il punto di partenza.

La funzione $f$ in un intorno $\delta_\beta$ di $\boldsymbol\beta$ può essere approssimata con una espansione in serie

$$f(\boldsymbol\beta + \delta_\beta) \approx f(\boldsymbol\beta) + \mathbf{J} \delta_\beta$$

con $\mathbf{J} = \frac{\partial f}{\partial \boldsymbol\beta}(\boldsymbol\beta)$ jacobiano della funzione $f$ calcolato in $\boldsymbol\beta$.

Come per ogni metodo iterativo, LM produce una serie di candidati $\boldsymbol\beta_i$ che convergono verso il minimo locale e di conseguenza per ogni iterazione è richiesto il calcolo di un $\boldsymbol\delta_\beta$ che minimizza la quantità

$$\Vert \mathbf{x} - f(\boldsymbol\beta) - \mathbf{J} \boldsymbol\delta_\beta \Vert = \Vert \mathbf{r} - \mathbf{J} \delta_\beta \Vert$$

Questo è un problema di minimo lineare che si può risolvere attraverso l'uso della normal equations:

$$\mathbf{H} \boldsymbol\delta_\beta = \mathbf{J}^{\top} \mathbf{r}$$ (3.11)

La matrice $\mathbf{H} = \mathbf{J}^{\top} \mathbf{J}$ è la Hessiana approssimata di $f$, matrice simmetrica e semidefinita positiva. L'algoritmo di Gauss-Newton indica come nuova stima per $\boldsymbol\beta$ la soluzione del sistema 3.11.

L'algoritmo di Levenberg-Marquardt risolve invece una versione leggermente differente dell'equazione 3.11 conosciuta come augmented normal equations:

$$\mathbf{N} \delta_\beta = \mathbf{J}^{\top} \mathbf{r}$$ (3.12)

dove $\mathbf{N} = \mathbf{H} + \mu \mathbf{I}$ con $\mu>0$ un fattore di attenuazione (damping factor). Quando il fattore $\mu$ è elevato, la matrice $\mathbf{N}$ è pressochè diagonale e l'algoritmo si avvicina a un metodo di discesa del gradiente (steepest descent gradient). Quando il termine $\mu$ è vicino a zero, l'algoritmo approssima il metodo di Gauss-Newton.

Come impostare e come modificare tra le iterazioni $\mu$ tuttavia è un problema lasciato al risolutore e diverse tecniche sono proposte in letteratura.

Siccome il massimo elemento della diagonale di $\mathbf{H}$ ha lo stesso ordine di grandezza del massimo autovalore si può prendere come $\mu_0$ un valore del tipo

$$\mu_0 = \tau \max \trace \mathbf{H}$$

con $\tau$ scelto liberamente dall'utente basandosi sulla propria fiducia rispetto al valore di $\boldsymbol\beta$.

La modifica di $\mu$ tra le iterazioni può essere controllata dal fattore di guadagno (gain ratio):

$$\rho = \frac{ s(\boldsymbol\beta) - s(\boldsymbol\beta + \delta_\... ...a^{\top}_\beta (\mu \boldsymbol\delta_\beta + \mathbf{J}^{\top} \mathbf{r} ) }$$

Un elevato valore di $\rho$ indica che la versione linearizzata di $f$ è molto buona e si può diminuire $\mu$. Viceversa se $\rho$ è elevato, allora il valore di $\mu$ è da aumentare. Caso limite, quando $\rho$ è negativo indica una soluzione peggiorativa da scartare e $\mu$ è da aumentare in modo da avvicinarsi a un metodo a discesa del gradiente.