Algoritmul Gauss-Newton

Versiunea actuală a paginii nu a fost încă examinată de colaboratori experimentați și poate diferi semnificativ de versiunea revizuită la 25 ianuarie 2021; verificarea necesită 1 editare .

Algoritmul Gauss-Newton este folosit pentru a rezolva probleme prin metoda neliniară a celor mai mici pătrate . Algoritmul este o modificare a metodei lui Newton pentru găsirea minimului funcției . Spre deosebire de metoda Newton, algoritmul Gauss-Newton poate fi folosit doar pentru a minimiza suma pătratelor, dar avantajul său este că metoda nu necesită calculul derivatelor secunde, ceea ce poate reprezenta o dificultate semnificativă.

Problemele pentru care se aplică metoda celor mai mici pătrate neliniare apar, de exemplu, în regresia neliniară , în care sunt căutați parametrii modelului care sunt cel mai în concordanță cu valorile observate.

Metoda poartă numele matematicienilor Carl Friedrich Gauss și Isaac Newton .

Descriere

Având în vedere m funcții r = ( r 1 , …, r m ) (deseori numite reziduuri) a n variabile β = ( β 1 , …, β n ), pentru m ≥ n . Algoritmul Gauss-Newton găsește iterativ valorile variabilelor care minimizează suma pătratelor [1]

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

Pornind de la o aproximare inițială , metoda iterează ${\boldsymbol {\beta }}^{(0))$

{\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 }} ^{(e)})

Aici, dacă considerăm r și β ca vectori coloană, elementele matricei jacobiene sunt

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

iar simbolul înseamnă transpunerea matricei . $^{\mathsf {T}}$

Dacă m = n , iterațiile sunt simplificate la

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

care este o generalizare directă a metodei unidimensionale a lui Newton .

La ajustarea datelor, unde scopul este de a găsi parametrii β astfel încât un model dat de funcții y = f ( x , β ) să aproximeze cel mai bine punctele de date ( x i , y i ), funcțiile r i sunt erori reziduale

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

Atunci metoda Gauss-Newton poate fi exprimată în termenii jacobian J f al funcției f

{\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 }} ^{(e)}).

Rețineți că este o matrice pseudo -inversă pentru . $\left(\mathbf {J_{f}} ^{\mathsf {T}}\mathbf {J_{f}} \right)^{-1}\mathbf {J_{f}} ^{\mathsf {T}}$ $\mathbf {J_{f)}$

Note

Cerința m ≥ n în algoritm este necesară, deoarece altfel matricea J r T J r nu are inversă și ecuațiile normale nu pot fi rezolvate (cel puțin fără ambiguitate).

Algoritmul Gauss-Newton poate fi obținut folosind o aproximare liniară a vectorului funcție r i . Folosind teorema lui Taylor , putem scrie pentru fiecare iterație:

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

unde . Problema găsirii Δ minimizând suma pătratelor din partea dreaptă, i.e. $\Delta ={\boldsymbol {\beta }}-{\boldsymbol {\beta }}^{s}$

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

este o problemă liniară cu cele mai mici pătrate care poate fi rezolvată explicit, dând ecuații normale.

Ecuațiile normale sunt m ecuații liniare în incremente necunoscute Δ. Ecuațiile pot fi rezolvate într-o singură etapă folosind descompunerea Cholesky , sau mai bine, descompunerea QR a matricei J r . Pentru sistemele mari, metoda iterativă poate fi mai eficientă dacă sunt utilizate metode precum metoda gradientului conjugat . Dacă există o dependență liniară a coloanelor matricei J r , metoda iterației eșuează deoarece J r T J r devine degenerată.

Exemplu

Acest exemplu folosește algoritmul Gauss-Newton pentru a construi un model de date prin minimizarea sumei abaterilor pătrate ale datelor și ale modelului.

În biologia experimentală, studiul relației dintre concentrația substratului [ S ] și viteza de reacție în reacția de modulare a enzimei, s-au obținut următoarele date.

i	unu	2	3	patru	5	6	7
[ S ]	0,038	0,194	0,425	0,626	1.253	2.500	3.740
viteză	0,050	0,127	0,094	0,2122	0,2729	0,2665	0,3317

Este necesar să se găsească o curbă (funcție-model) a formei

viteza ,

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

care aproximează cel mai bine datele în sensul celor mai mici pătrate cu parametrii și de găsit. $V_{\text{max})$ $K_{M)$

Notați cu și valorile lui [ S ] și viteza din tabel, . Lasă și . Vom căuta și , astfel încât suma abaterilor pătrate $x_{i}$ $y_{i}$ $i=1,\dots ,7$ $\beta _{1}=V_{\text{max)}$ $\beta _{2}=K_{M)$ $\beta_{1}$ $\beta _{2}$

r_{i}=y_{i}-{\frac {\beta _{1}x_{i}}{\beta _{2}+x_{i}}}\;(i=1,\ puncte ,7)

minim.

Jacobianul vectorului de reziduuri peste necunoscute este o matrice cu rândul --lea având elementele $\mathbf {J_{r}}$ $r_{i}$ $\beta _{j)$ $7\times 2$ $i$

{\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}}{\left(\beta _{2} +x_{i}\dreapta)^{2}}}.

Pornind de la aproximarea inițială și după cinci iterații, algoritmul Gauss-Newton oferă valorile optime ale și . Suma reziduurilor pătrate scade de la valoarea inițială de 1,445 la 0,00784 cu a cincea iterație. Graficul din dreapta arată curba cu parametri optimi. $\beta _{1}=0,9$ $\beta _{2}=0,2$ ${\hat {\beta }}_{1}=0,362$ ${\hat {\beta }}_{2}=0,556$

Convergență

Se poate arăta [2] că direcția de creștere a Δ este direcția de descreștere pentru S , iar dacă algoritmul converge, limita va fi punctul staționar pentru S . Totuși, convergența nu este garantată chiar și atunci când punctul de plecare este aproape de soluția , ceea ce se întâmplă în metoda Newton sau metoda BFGS în condiții normale Volfe [3] .

Rata de convergență a algoritmului Gauss-Newton este apropiată de pătratică [4] . Algoritmul poate converge mai lent sau deloc dacă estimarea inițială este departe de minim sau dacă matricea este prost condiționată . De exemplu, imaginați-vă o problemă cu ecuații și o variabilă $\mathbf {J_{r}^{\mathsf {T}}J_{r}}$ $m=2$ $n=1$

{\begin{aligned}r_{1}(\beta )&=\beta +1\\r_{2}(\beta )&=\lambda \beta ^{2}+\beta -1.\ sfârşit{aliniat}}

Soluția optimă rezultată este . (Optimul real este pentru , deoarece , în timp ce .) Dacă , atunci problema este, de fapt, liniară și metoda găsește o soluție într-o iterație. Dacă |λ| < 1, atunci metoda converge liniar și eroarea scade cu o rată de |λ| la fiecare iterație. Totuși, dacă |λ| > 1, atunci metoda nu converge nici măcar local [5] . $\beta =0$ $\beta =-1$ $\lambda =2$ $S(0)=1^{2}+(-1)^{2}=2$ $S(-1)=0$ $\lambda=0$

Algoritm bazat pe metoda lui Newton

Următoarele presupune că algoritmul Gauss-Newton se bazează pe metoda lui Newton pentru minimizarea funcției prin aproximare. În consecință, rata de convergență a algoritmului Gauss-Newton poate fi pătratică dacă sunt îndeplinite anumite condiții. În cazul general (în condiții mai slabe), rata de convergență poate fi liniară [6] .

Relația de recurență a metodei lui Newton de minimizare a funcției S a parametrilor ${\boldsymbol {\beta })$

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

unde g reprezintă vectorul gradient al funcției S , iar H reprezintă Hessianul funcției S . Deoarece , gradientul este dat de egalitate $S=\sum _{i=1}^{m}r_{i}^{2)$

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

Elementele hessiene sunt calculate prin diferențierea elementelor de gradient în raport cu $g_{j)$ $\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}}}\dreapta).

Metoda Gauss-Newton se obține prin eliminarea derivatei a doua (al doilea termen din expresie). Adică Hessianul este aproximativ

H_{jk}\aprox 2\sum _{i=1}^{m}J_{ij}J_{ik)

unde sunt elemente ale jacobianului J r . Gradientul și Hessianul aproximativ pot fi scrise în notație matriceală $J_{ij}={\frac {\partial r_{i}}{\partial \beta _{j}}}$

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

Aceste expresii sunt substituite în relația de recursivitate de mai sus pentru a obține ecuațiile de operare

{\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} .

În general, convergența metodei Gauss-Newton nu este garantată. Apropiere

\left|r_{i}{\frac {\partial ^{2}r_{i}}{\partial \beta _{j}\partial \beta _{k}}}\right|\ll \ stânga|{\frac {\partial r_{i}}{\partial \beta _{j}}}{\frac {\partial r_{i}}{\partial \beta _{k}}}\right|

care trebuie să se țină pentru a putea renunța la termeni cu derivata a doua, se poate obține în două cazuri pentru care se așteaptă convergență [7]

Valorile funcției sunt mici ca mărime, cel puțin aproape de minim. $r_{i}$
Funcțiile sunt doar „puțin” neliniare, adică relativ mici ca mărime. ${\frac {\partial ^{2}r_{i}}{\partial \beta _{j}\partial \beta _{k)))$

Versiuni îmbunătățite

În metodele Gauss-Newton, suma reziduurilor pătrate S poate să nu scadă la fiecare iterație. Totuși, deoarece Δ este direcționat în direcția scăderii funcției, dacă nu este un punct staționar, inegalitatea este valabilă pentru suficient de mic . Astfel, dacă se găsește o divergență, se poate folosi fracția vectorului de creștere Δ în formula de actualizare: $S({\boldsymbol {\beta }}^{s})$ $S({\boldsymbol {\beta}}^{s}+\alpha \Delta )<S({\boldsymbol {\beta }}^{s})$ $\alpha >0$ $\alfa$

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

Cu alte cuvinte, vectorul de increment este prea lung, dar indică direcția „coborârii”, așa că dacă parcurgeți doar o parte din drum, puteți reduce valoarea funcției S . Valoarea optimă poate fi găsită folosind un algoritm de căutare unidimensional , adică valoarea este determinată prin găsirea valorii care minimizează S folosind o căutare unidimensională pe interval . $\alfa$ $\alfa$ $0<\alpha<1$

În cazurile în care fracția optimă este aproape de zero în direcția vectorului de increment, o metodă alternativă de calcul a divergenței este utilizarea algoritmului Levenberg-Marquardt , cunoscut și sub denumirea de „metoda regiunii de încredere” [1] . Ecuații normale modificate astfel încât vectorul de coborâre să se rotească în direcția celei mai abrupte coborâri , $\alfa$

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

unde D este o matrice diagonală pozitivă. Rețineți că dacă D este matricea de identitate a lui E și , atunci . Astfel direcția Δ aproximează direcția gradientului negativ . $\lambda \to +\infty$ $\lambda \Delta =\lambda \left(\mathbf {J^{E}J} +\lambda \mathbf {E} \right)^{-1}\left(-\mathbf {J} ^{ T}\mathbf {r} \right)=\left(\mathbf {E} -\mathbf {J^{T}J} /\lambda +\cdots \right)\left(-\mathbf {J} ^{ T}\mathbf {r} \right)\to -\mathbf {J} ^{T}\mathbf {r}$ $-\mathbf {J} ^{T}\mathbf {r}$

Așa-numitul parametru Marquardt poate fi optimizat și prin căutare liniară, dar nu are prea mult sens, deoarece vectorul de deplasare trebuie recalculat de fiecare dată când se modifică . O strategie mai eficientă este aceasta. Dacă se găsește o discrepanță, creșteți parametrul Marquardt pe măsură ce S scade. Apoi păstrăm valoarea între iterații, dar o reducem, dacă este posibil, până ajungem la o valoare în care parametrul Marquardt nu poate fi pus la zero. Minimizarea lui S devine apoi minimizarea standard Gauss-Newton. $\lambda$ $\lambda$

Optimizarea sarcinilor mari

Pentru optimizări de dimensiuni mari, metoda Gauss-Newton este deosebit de interesantă, deoarece adesea (deși cu siguranță nu întotdeauna) matricea este rară decât Hessianul aproximativ . În astfel de cazuri, etapa de calcul în sine necesită de obicei utilizarea unei metode de aproximare iterativă, cum ar fi metoda gradientului conjugat . $\mathbf {J} _{\mathbf {r} )$ $\mathbf {J} _{\mathbf {r} }^{\mathsf {T}}\mathbf {J_{r}}$

Pentru ca această abordare să funcționeze, aveți nevoie de cel puțin o metodă eficientă de calcul al produsului

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

pentru un vector p . Pentru a stoca o matrice rară, este practic să stocați rândurile matricei în formă comprimată (adică fără elemente zero), ceea ce face dificilă calcularea directă a produsului de mai sus (din cauza transpunerii). Totuși, dacă c i este definit ca rândul i al matricei , următoarea relație este valabilă: $\mathbf {J} _{\mathbf {r} )$ $\mathbf {J} _{\mathbf {r} )$

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

astfel încât orice rând contribuie aditiv și independent la produs. În plus, această expresie este bine studiată pentru aplicarea calculului paralel . Rețineți că orice rând c i este gradientul rezidualului corespunzător r i . Luând în considerare această împrejurare, formula de mai sus subliniază faptul că reziduurile contribuie la rezultat independent unele de altele.

Algoritmi înrudiți

În metodele cvasi-newtoniene , cum ar fi metodele lui Davidon, Fletcher și Powell sau Broyden-Fletcher-Goldfarb-Shanno ( metoda BFGSh ), aproximarea hessiană completă este construită folosind primele derivate, astfel încât după n rafinamente metoda să fie apropiată ca performanță de metoda Newton. Rețineți că metodele cvasi-newtoniene pot minimiza funcțiile reale de formă generală, în timp ce metodele lui Gauss-Newton, Levenberg-Marquardt etc. sunt aplicabile numai problemelor neliniare cu cele mai mici pătrate. ${\frac {\partial ^{2}S}{\partial \beta _{j}\partial \beta _{k))}$ ${\frac {\partial r_{i}}{\partial \beta _{j}}}$

O altă metodă de rezolvare a problemelor de minimizare folosind numai derivate primare este metoda coborârii gradientului . Cu toate acestea, această metodă nu ține cont de derivatele secunde, chiar și de cele aproximative. Ca urmare, metoda este extrem de ineficientă pentru multe funcții, mai ales în cazul unei influențe reciproce puternice a parametrilor.

Note

↑ 1 2 Björck, 1996 .
↑ Björck, 1996 , p. 260.
↑ Mascarenhas, 2013 , p. 253–276.
↑ Björck, 1996 , p. 341, 342.
↑ Fletcher, 1987 , p. 113.
↑ Gratton, Lawless, Nichols .
↑ Nocedal, Wright, 1999 , p. 259-262.

Literatură

A. Björck. Metode numerice pentru probleme cu cele mai mici pătrate. - Philadelphia: SIAM, 1996. - ISBN 0-89871-360-9 .
Roger Fletcher. Metode practice de optimizare . — al 2-lea. - New York: John Wiley & Sons , 1987. - ISBN 978-0-471-91547-8 .
Walter F. Mascarenhas. Divergența metodelor BFGS și Gauss Newton // Programare matematică. - 2013. - T. 147 , nr. 1 . - doi : 10.1007/s10107-013-0720-6 .
S. Gratton, AS Lawless, NK Nichols. Metode Gauss-Newton aproximative pentru probleme neliniare cu cele mai mici pătrate . RAPORT DE ANALIZĂ NUMERICALĂ 9/04 (engleză) (link nu este disponibil) . Universitatea din Reading (ianuarie 2007) . Preluat la 20 iulie 2017. Arhivat din original la 4 august 2016.
Jorge Nocedal, Stephen J. Wright. Optimizare numerică / Peter Glynn, Stephen M. Robinson. - New York: Springer, 1999. - (Seria Springer în cercetarea operațională). — ISBN 0-387-98793-2 .

Link -uri

Implementări

Artelys Knitro . Sistem de rezolvare a problemelor neliniare cu implementarea metodei Gauss-Newton. Sistemul este scris în C și are interfețe pentru C++/C#/Java/Python/MATLAB/R.

Metode de optimizare
Unidimensional	metoda secțiunii de aur Dihotomie Metoda parabolelor Căutare în grilă Metoda de căutare uniformă a blocurilor Metoda Fibonacci Căutare ternară metoda Piyavsky Metoda Strongin
Comanda zero	metoda Gauss Metoda Nelder-Mead Metoda Hook-Jeeves metoda Rosenbrock metoda Powell
Prima comanda	coborâre în gradient Metoda Zeutendijk Coordonarea coborârii Metoda gradientului conjugat Metode cvasi-newtoniene Algoritmul Levenberg-Marquardt
a doua comanda	metoda lui Newton Metoda Newton-Raphson Algoritmul Broyden-Fletcher-Goldfarb-Shanno (BFGS)
Stochastic	Metoda Monte Carlo Recoacere simulată Algoritmi evolutivi evolutie diferentiala Algoritmul furnicilor Metoda roiului de particule Algoritmul coloniilor de albine Metoda de mers aleatoriu
Metode de programare liniară	Metoda simplex algoritmul lui Gomori Metoda elipsoidă Metoda potențială
Metode de programare neliniară	Programare secvenţială pătratică