Gauss-Newton algoritm

Den aktuella versionen av sidan har ännu inte granskats av erfarna bidragsgivare och kan skilja sig väsentligt från versionen som granskades den 25 januari 2021; verifiering kräver 1 redigering .

Gauss-Newton-algoritmen används för att lösa problem med den olinjära minsta kvadratmetoden . Algoritmen är en modifiering av Newtons metod för att hitta minimum av funktionen . Till skillnad från Newtonmetoden kan Gauss-Newton-algoritmen endast användas för att minimera summan av kvadrater, men dess fördel är att metoden inte kräver beräkning av andraderivator, vilket kan vara en betydande svårighet.

Problem för vilka den icke-linjära minsta kvadratmetoden tillämpas uppstår till exempel vid icke-linjär regression , där man söker efter modellparametrar som är mest överensstämmande med de observerade värdena.

Metoden är uppkallad efter matematikerna Carl Friedrich Gauss och Isaac Newton .

Beskrivning

Givet m funktioner r = ( r 1 , …, r m ) (ofta kallade residualer) av n variabler β = ( β 1 , …, β n ), för m ≥ n . Gauss-Newton-algoritmen hittar iterativt värdena för variabler som minimerar summan av kvadrater [1]

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

Med utgångspunkt från en initial uppskattning upprepas metoden ${\displaystyle {\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 }} ^{(s)})

Här, om vi betraktar r och β som kolumnvektorer, är elementen i den jakobiska matrisen

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

och symbolen betyder matristransponering . $^{\mathsf {T}}$

Om m = n förenklas iterationerna till

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

vilket är en direkt generalisering av Newtons endimensionella metod .

Vid anpassning av data, där målet är att hitta parametrar β så att en given modell av funktioner y = f ( x , β ) bäst approximerar datapunkter ( x i , y i ), är funktionerna r i restfel [

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

Då kan Gauss-Newton-metoden uttryckas i termer av Jacobian J f för funktionen 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 }} ^{(s)}).

Observera att det är en pseudo -invers matris till . $\left(\mathbf {J_{f}} ^{\mathsf {T}}\mathbf {J_{f}} \right)^{-1}\mathbf {J_{f}} ^{\mathsf {T}}$ $\mathbf {J_{f))$

Anteckningar

Kravet m ≥ n i algoritmen är nödvändigt, eftersom matrisen J r T J r annars inte har någon invers och normalekvationerna inte kan lösas (åtminstone entydigt).

Gauss - Newton-algoritmen kan erhållas genom att använda en linjär approximation av funktionsvektorn ri . Med hjälp av Taylors teorem kan vi skriva för varje iteration:

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

var . Problemet med att hitta Δ minimera summan av kvadrater på höger sida, dvs. $\Delta ={\boldsymbol {\beta }}-{\boldsymbol {\beta }}^{s}$

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

är ett linjärt minsta kvadratproblem som kan lösas explicit, vilket ger normala ekvationer.

Normalekvationer är m linjära ekvationer i okända steg Δ. Ekvationerna kan lösas i ett steg med hjälp av Cholesky-sönderdelningen , eller bättre, QR-sönderdelningen av matrisen Jr. För stora system kan den iterativa metoden vara effektivare om metoder som konjugatgradientmetoden används . Om det finns ett linjärt beroende av kolumnerna i matrisen Jr , misslyckas iterationsmetoden eftersom JrTJr blir degenererad .

Exempel

Det här exemplet använder Gauss-Newton-algoritmen för att bygga en datamodell genom att minimera summan av de kvadrerade avvikelserna för data och modellen.

Inom experimentell biologi, studiet av sambandet mellan koncentrationen av substratet [ S ] och reaktionshastigheten i enzymmoduleringsreaktionen, erhölls följande data.

i	ett	2	3	fyra	5	6	7
[ S ]	0,038	0,194	0,425	0,626	1,253	2 500	3,740
fart	0,050	0,127	0,094	0,2122	0,2729	0,2665	0,3317

Det är nödvändigt att hitta en kurva (funktionsmodell) av formen

hastighet ,

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

som bäst approximerar data i betydelsen minsta kvadrater med parametrarna och som ska hittas. $V_{\text{max))$ ${\displaystyle K_{M))$

Beteckna med och värdena för [ S ] och hastigheten från tabellen, . Låt och . Vi kommer att leta efter och , så att summan av de kvadratiska avvikelserna $x_{i}$ $y_{i}$ $i=1,\dots ,7$ $\beta _{1}=V_{\text{max))$ ${\displaystyle \beta _{2}=K_{M))$ $\beta_{1}$ $\beta _{2}$

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

minimal.

Jacobian för vektorn av residualer över okända är en matris där den -th raden har elementen $\mathbf {J_{r}}$ $r_{i}$ ${\displaystyle \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}\right)^{2}}}.

Från den initiala approximationen och efter fem iterationer ger Gauss-Newton-algoritmen de optimala värdena för och . Summan av kvadrerade residualer minskar från startvärdet 1,445 till 0,00784 med den femte iterationen. Grafen till höger visar kurvan med optimala parametrar. $\beta _{1}=0.9$ $\beta _{2}=0.2$ ${\hat {\beta }}_{1}=0,362$ ${\hat {\beta }}_{2}=0,556$

Konvergens

Det kan visas [2] att riktningen för ökande Δ är riktningen för fallande för S , och om algoritmen konvergerar kommer gränsen att vara den stationära punkten för S . Konvergens är dock inte garanterad även när utgångspunkten är nära lösningen , vilket sker i Newtonmetoden eller BFGS- metoden under normala Volfe-förhållanden [3] .

Konvergenshastigheten för Gauss-Newton-algoritmen är nära kvadratisk [4] . Algoritmen kan konvergera långsammare eller inte alls om den initiala gissningen är långt ifrån minimum, eller om matrisen är dåligt konditionerad . Föreställ dig till exempel ett problem med ekvationer och en variabel $\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.\ end{aligned}}

Den resulterande optimala lösningen är . (Det verkliga optimum är för , eftersom , medan .) Om , då är problemet i själva verket linjärt och metoden hittar en lösning i en iteration. Om |λ| < 1, då konvergerar metoden linjärt och felet minskar med en hastighet av |λ| vid varje iteration. Men om |λ| > 1, då konvergerar metoden inte ens lokalt [5] . $\beta =0$ $\beta =-1$ $\lambda =2$ $S(0)=1^{2}+(-1)^{2}=2$ $S(-1)=0$ $\lambda=0$

Algoritm baserad på Newtons metod

Följande antar att Gauss-Newton-algoritmen är baserad på Newtons metod för funktionsminimering genom approximation. Som en konsekvens kan konvergenshastigheten för Gauss-Newton-algoritmen vara kvadratisk om vissa villkor är uppfyllda. I det allmänna fallet (under svagare förhållanden) kan konvergenshastigheten vara linjär [6] .

Återkommande relation för Newtons metod för att minimera funktionen S för parametrar ${\boldsymbol {\beta ))$

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

där g betecknar gradientvektorn för funktionen S och H betecknar hessian för funktionen S . Eftersom , gradienten ges av jämlikheten ${\displaystyle S=\sum _{i=1}^{m}r_{i}^{2))$

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

De hessiska elementen beräknas genom att differentiera gradientelementen med avseende på ${\displaystyle g_{j))$ ${\displaystyle \beta _{k))$

H_{jk}=2\summa _{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}}}\höger).

Gauss-Newton-metoden erhålls genom att kassera den andra derivatan (den andra termen i uttrycket). Det vill säga, hessian är ungefärlig

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

var finns element av Jacobian J r . Gradienten och den ungefärliga Hessian kan skrivas i matrisnotation $J_{ij}={\frac {\partial r_{i}}{\partial \beta _{j}}}$

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

Dessa uttryck ersätts i rekursionsrelationen ovan för att erhålla driftsekvationerna

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

Konvergensen av Gauss-Newton-metoden är i allmänhet inte garanterad. Approximation

\left|r_{i}{\frac {\partial ^{2}r_{i}}{\partial \beta _{j}\partial \beta _{k}}}\right|\ll \ vänster|{\frac {\partial r_{i}}{\partial \beta _{j}}}{\frac {\partial r_{i}}{\partial \beta _{k}}}\right|

som måste hålla för att kunna kassera termer med den andra derivatan, kan erhållas i två fall för vilka konvergens förväntas [7]

Funktionsvärdena är små i storleksordningen, åtminstone nära minimum. $r_{i}$
Funktionerna är bara "något" icke-linjära, det vill säga relativt små till storleken. ${\displaystyle {\frac {\partial ^{2}r_{i}}{\partial \beta _{j}\partial \beta _{k))))$

Förbättrade versioner

I Gauss-Newton-metoder kanske summan av kvadratiska rester S inte minskar vid varje iteration. Men eftersom Δ är riktad i riktning mot att minska funktionen, om det inte är en stationär punkt, gäller olikheten för tillräckligt liten . Således, om en divergens hittas, kan man använda bråkdelen av inkrementvektorn Δ i uppdateringsformeln: $S({\boldsymbol {\beta }}^{s})$ $S({\boldsymbol {\beta }}^{s}+\alpha \Delta )<S({\boldsymbol {\beta }}^{s})$ $\alfa >0$ $\alfa$

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

Inkrementvektorn är med andra ord för lång, men den indikerar riktningen för "nedstigningen", så om du bara går en del av vägen kan du minska värdet på S -funktionen . Det optimala värdet kan hittas med en endimensionell sökalgoritm , det vill säga värdet bestäms genom att hitta värdet som minimerar S med en endimensionell sökning på intervallet . $\alfa$ $\alfa$ $0<\alfa<1$

I de fall där den optimala fraktionen är nära noll i inkrementvektorns riktning, är en alternativ metod för att räkna ut divergensen att använda Levenberg-Marquardt-algoritmen , även känd som "konfidensregionmetoden" [1] . Normalekvationer modifierade så att nedstigningsvektorn roterar i riktningen för den brantaste nedstigningen , $\alfa$

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

där D är en positiv diagonal matris. Observera att om D är identitetsmatrisen för E och , då . Således riktningen Δ approximerar riktningen för den negativa gradienten . $\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}$

Den så kallade Marquardt-parametern kan också optimeras genom linjär sökning, men det är inte så meningsfullt, eftersom skiftvektorn måste räknas om varje gång den ändras . En mer effektiv strategi är detta. Om en avvikelse hittas, öka Marquardt-parametern när S minskar. Sedan behåller vi värdet mellan iterationerna, men minskar det, om möjligt, tills vi når ett värde där Marquardt-parametern inte kan nollställas. Minimeringen av S blir då standard Gauss-Newton-minimeringen. $\lambda$ $\lambda$

Optimering av stora uppgifter

För storstorleksoptimeringar är Gauss-Newton-metoden särskilt intressant eftersom matrisen ofta (men absolut inte alltid) är gles än den ungefärliga Hessian . I sådana fall kräver själva beräkningssteget vanligtvis användning av en iterativ approximationsmetod, såsom den konjugerade gradientmetoden . ${\displaystyle \mathbf {J} _{\mathbf {r} ))$ $\mathbf {J} _{\mathbf {r} }^{\mathsf {T}}\mathbf {J_{r}}$

För att detta tillvägagångssätt ska fungera behöver du åtminstone en effektiv metod för att beräkna produkten

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

för någon vektor p . För att lagra en gles matris är det praktiskt att lagra matrisens rader i komprimerad form (dvs utan nollelement), vilket gör den direkta beräkningen av ovanstående produkt (på grund av transponering) svår. Men om c i definieras som rad i i matrisen gäller följande relation: ${\displaystyle \mathbf {J} _{\mathbf {r} ))$ ${\displaystyle \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} )

så vilken rad som helst bidrar additivt och oberoende till produkten. Dessutom är detta uttryck väl studerat för tillämpningen av parallell beräkning . Observera att varje rad c i är gradienten för motsvarande rest r i . Med hänsyn till denna omständighet understryker formeln ovan det faktum att restprodukter bidrar till resultatet oberoende av varandra.

Relaterade algoritmer

I kvasi-newtonska metoder , såsom metoderna av Davidon, Fletcher och Powell eller Broyden-Fletcher-Goldfarb-Shanno ( BFGSh-metoden ), konstrueras den fullständiga hessiska approximationen med hjälp av de första derivatorna så att metoden efter n förfinningar är i prestanda nära Newtonmetoden. Observera att kvasi-newtonska metoder kan minimera verkliga funktioner av en allmän form, medan metoderna för Gauss-Newton, Levenberg-Marquardt, etc. är endast tillämpliga på icke-linjära minsta kvadraters problem. ${\displaystyle {\frac {\partial ^{2}S}{\partial \beta _{j}\partial \beta _{k))))$ ${\frac {\partial r_{i}}{\partial \beta _{j}}}$

En annan metod för att lösa minimeringsproblem med endast förstaderivator är metoden för gradientnedstigning . Denna metod tar dock inte hänsyn till andraderivator, inte ens ungefärliga sådana. Som ett resultat är metoden extremt ineffektiv för många funktioner, speciellt vid stark ömsesidig påverkan av parametrar.

Anteckningar

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

Litteratur

A. Björck. Numeriska metoder för minsta kvadraters problem. - Philadelphia: SIAM, 1996. - ISBN 0-89871-360-9 .
Roger Fletcher. Praktiska metoder för optimering . — 2:a. - New York: John Wiley & Sons , 1987. - ISBN 978-0-471-91547-8 .
Walter F. Mascarenhas. Divergensen mellan BFGS och Gauss Newton-metoderna // Matematisk programmering. - 2013. - T. 147 , nr. 1 . - doi : 10.1007/s10107-013-0720-6 .
S. Gratton, AS Lawless, NK Nichols. Ungefärliga Gauss-Newton-metoder för olinjära minsta kvadraters problem . NUMERISK ANALYSRAPPORT 9/04 (engelska) (länk ej tillgänglig) . University of Reading (januari 2007) . Hämtad 20 juli 2017. Arkiverad från originalet 4 augusti 2016.
Jorge Nocedal, Stephen J. Wright. Numerisk optimering / Peter Glynn, Stephen M. Robinson. - New York: Springer, 1999. - (Springer Series in Operations Research). — ISBN 0-387-98793-2 .

Länkar

Implementeringar

Artelys Knitro . System för att lösa olinjära problem med implementering av Gauss-Newton-metoden. Systemet är skrivet i C och har gränssnitt för C++/C#/Java/Python/MATLAB/R.

Optimeringsmetoder _
En-dimensionell	gyllene snittmetoden Dikotomi Parabolmetoden Rutnätssökning Enhetlig blocksökningsmetod Fibonacci-metoden Ternär sökning Piyavsky-metoden Strongin metod
Noll ordning	Gauss metod Nelder-Mead metod Hook-Jeeves metod Rosenbrock-metoden Powell metod
Första beställning	lutning nedstigning Zeutendijk-metoden Koordinera nedstigning Konjugerad gradientmetod Kvasi-newtonska metoder Levenberg-Marquardts algoritm
andra beställning	Newtons metod Newton-Raphson-metoden Broyden-Fletcher-Goldfarb-Shanno-algoritm (BFGS)
Stokastisk	Monte Carlo metoden Simulerad glödgning Evolutionära algoritmer differentiell evolution Myralgoritm Partikelsvärmmetod Algoritm för bikoloni Random walk-metod
Linjära programmeringsmetoder _	Enkel metod Gomoris algoritm Ellipsoid metod Potentiell metod
Icke -linjära programmeringsmetoder	Sekventiell kvadratisk programmering