ADAPTIVNÍ ALGORITMUS PRO ODHAD PARAMETRŮ NELINEÁRNÍCH REGRESNÍCH MODELŮ

ADAPTIVNÍ ALGORITMUS PRO ODHAD PARAMETRŮ NELINEÁRNÍCH REGRESNÍCH MODELŮ Josef Tvrdík Katedra informatiky, Přírodovědecká fakulta, Ostravská univerzita, 30. dubna 22, 701 03 Ostrava e-mail: [email protected]

Abstrakt. Článek se zabývá stochastickými algoritmy pro globální optimalizaci a jejich využitím v odhadu parametrů nelineárních regresních modelů. Je popsán algoritmus řízeného náhodného prohledávání se soutěží heuristik lokálního vyhledávání a adaptivní podmínkou ukončení. Tento algoritmus se osvědčil v obtížných úlohách odhadu parametrů nelineárního regresního modelu a lze jej doporučit pro odhad parametrů v úlohách, kde standardní deterministické algoritmy pro lokální optimalizaci selhávají. Klíčová slova: Globální optimalizace, stochastické algoritmy, nelineární regrese, adaptivní algoritmy.

1. Nelineární regresní model Uvažujme data (y, X), kde y je náhodný vektor typu n × 1 a X je matice typu n × k deterministických hodnot vysvětlujících veličin (regresorů). Pak aditivní nelineární regresní model má tvar Yi = g(xi , β) + εi ,

i = 1, 2, . . . , n ,

(1)

kde xTi = (x1 , x2 , . . . , xk ) je i-tý řádek matice X, β = (β1 , β2 , . . . , βd )T je vektor parametrů, g je daná nelineární funkce parametrů a εi jsou stejně rozdělené nezávislé náhodné veličiny s nulovou střední hodnotou. Odhad parametrů nelineárního modelu metodou nejmenších čtverců znamená najít globální minimum reziduálního součtu čtverců n h i2 X ˆ = ˆ Q(β) Yi − g(xi , β) . (2) i=1

Je známo, že na těchto úlohách deterministické algoritmy užívané ve statistickém software někdy selhávají [8], [18]. Pak spolehlivější alternativou odhadu parametrů je využití stochastických algoritmů pro globální optimalizaci.

2. Globální optimalizace Problém globální optimalizace lze formulovat takto: Pro účelovou funkci f : D ⊆ Rd → R,

D 6= ∅

je hodnota f ∗ := f (x∗ ) > −∞ globální minimum právě tehdy, když ∀x ∈ D : f (x∗ ) ≤ f (x). Problém určení x∗ se nazývá problém globální optimalizace [2]. Tato formulace problému globální optimalizace jako nalezení minima není na újmu obecnosti, neboť globální maximum nalezneme jako globální minimum účelové funkce −f (x). Analýza problému globální optimalizace ukazuje, že neexistuje deterministický algoritmus řešící každou úlohu globální optimalizace v polynomiálním čase, problém globální optimalizace je NP-obtížný, viz Bäck [2], str. 35 – 57. Pro odhad parametrů nelineárních regresního modelů nám stačí omezit se na úlohy, kdy hledáme globální minimum v uzavřené souvislé oblasti Qd (3) D = j=1 haj , bj i, aj < bj , j = 1, 2, . . . , d , a účelovou funkci f (x) umíme vyhodnotit s požadovanou přesností v každém bodě x ∈ D. Omezení (3) se označuje jako box constraints, neboť oblast D je vymezena jako d-rozměrný kvádr. Jelikož neexistuje vhodný deterministický algoritmus pro řešení problému globální optimalizace, užívají se algoritmy stochastické, které hledají globální minimum heuristicky [14]. Většina takových algoritmů je 8.

NÁRODNÍ KONFERENCE

STATISTICKÉ DNY

V

BRNĚ, BRNO 27. - 28.

ČERVNA,

2006

založena na jednoduchých modelech Darwinovy vývojové teorie populací a poznatků genetiky. Proto se nazývají evoluční algoritmy [2], [7]. Jejich základním rysem je to, že pracují s populací jedinců (N bodů v D) a užívají tzv. evoluční operátory, což jsou modely selekce silnějších jedinců, křížení a mutace jejich genetické informace. V průběhu prohledávání prostoru D se populace vyvíjí tak, že postupně vznikají a přežívají jedinci s nižšími hodnotami účelové funkce. Lze ukázat, že za jistých podmínek se nejlepší jedinci populace dostatečně přibližují globálnímu minimu x∗ , viz např. [3], [13], [9]. Mezi evoluční algoritmy patří i řízené náhodné prohledávání, jehož adaptivní varianta užitá pro odhad parametrů nelineárních regresních modelů je popsána v tomto textu. Při aplikaci stochastických algoritmů v řešení praktických optimalizačních úloh narážíme na potíž nastavení vhodných hodnot řídících parametrů těchto algoritmů, což jsou např. velikost populace nebo výběr evolučních operátorů vhodných pro řešený problém. Pro rutinní využití stochastických algoritmů ve výpočetní statistice, tedy i pro odhad regresních parametrů, je nutné najít takové varianty stochastických algoritmů, které jsou adaptivní, tzn. v průběhu prohledávání prostoru D se adaptuje nejen populace bodů, ale také strategie prohledávání s ohledem na řešenou úlohu. Takové algoritmy jsou pak schopny řešit dosti širokou třídu problémů s implicitně nastavenými hodnotami řídících parametrů, tj. bez časově náročného dolaďování hodnot těchto parametrů pro řešenou úlohu, a přitom nacházet globální minimum s dostatečnou spolehlivostí a s přijatelnými časovými nároky.

3. Řízené náhodné prohledávání se soutěží lokálních heuristik Řízené náhodné prohledávání (Controlled Random Search, CRS) navrhl před třemi desetiletími Price [12] jako jednoduchý stochastický algoritmus pro hledání globálního minima. V řízeném náhodném prohledávání se v každém cyklu algoritmu nahrazuje nejhorší bod populace novým bodem, pokud tento nový bod populaci zlepšuje. Tak se populace postupně „smršťujeÿ do oblasti s nižšími hodnotami účelové funkce. Algoritmus CRS lze v pseudokódu zapsat takto: Algoritmus 1. Řízené náhodné prohledávání 1 generuj populaci P , tj. N bodů náhodně v D; 2 najdi xmax ∈ P , f (xmax ) ≥ f (x), x ∈ P ; 3 repeat 4 užij lokální heuristiku k vygenerování nového bodu y ∈ D; 5 if f (y) < f (xmax ) then 6 xmax := y; (nejhorší bod populace nahrazen novým) 7 najdi nové xmax ∈ P , f (xmax ) ≥ f (x), x ∈ P ; 8 endif ; 9 until podmínka ukončení; Lokální heuristikou na řádku 4 rozumíme libovolný nedeterministický předpis pro vygenerování nového bodu y. V původní variantě algoritmu CRS [12] se nový bod generuje reflexí v simplexu S, který je tvořen d + 1 body v D – viz Nelder a Mead [10], d + 1 bodů simplexu S vybírá z populace náhodně. Bod y lze však generovat jakoukoli jinou lokální heuristikou nebo dokonce tyto heuristiky střídat [16], [17]. Pokud v tomto střídání heuristik budou preferovány ty heuristiky, které byly v dosavadním průběhu hledání úspěšnější, je naděje, že algoritmus bude konvergovat rychleji. V takovém algoritmu lokální heuristiky soutěží o to, aby byly vybrány pro generování nového bodu. Mějme k dispozici h heuristik a v aktuálním kroku algoritmu vybíráme náhodně i-tou heuristiku s pravděpodobností qi , i = 1, 2, . . . , h. Pravděpodobnosti qi měníme v závislosti na úspěšnosti i-té heuristiky v dosavadním průběhu vyhledávacího procesu. Heuristiku považujeme za úspěšnou, když generuje takový nový bod y, že f (y) < f (xmax ). Pravděpodobnost qi můžeme vyhodnotit jako váženou relativní úspěšnost i-té heuristiky. Váha wi se určí jako míra zlepšení funkční hodnoty wi =

fmax − max(f (y), fmin ) , fmax − fmin

(4)

kde fmin , fmax jsou minimální a maximální hodnoty účelové funkce v populaci. Hodnoty wi jsou v intervalu (0, 1i. Pravděpodobnost qi se pak vyhodnotí jako qi = Ph

Wi + w 0

j=1 (Wj

8.

NÁRODNÍ KONFERENCE

STATISTICKÉ DNY

V

+ w0 )

,

BRNĚ, BRNO 27. - 28.

(5)

ČERVNA,

2006

kde Wi je součet vah wi v předcházejícím hledání a w0 > 0 je vstupní parametr algoritmu. Aby se zabránilo trvalému potlačení možnosti výběru některé z heuristik, lze zadat vstupní parametr δ > 0 a klesne-li kterákoli hodnota qi pod tuto hodnotu, jsou pravděpodobnosti výběru heuristik nastaveny na jejich počáteční hodnoty qi = 1/h. Pro různé varianty algoritmu můžeme volit různé množiny soutěžících lokálních heuristik, viz [16], [17], [19]. V odhadu parametrů nelineární regrese se osvědčil algoritmus CRS se soutěží následujících čtyř lokálních heuristik. Tři z těchto heuristik jsou založeny na znáhodněné reflexi simplexu, která byla navržena v [4]. Nový bod y se generuje ze simplexu S podle pravidla y = g + U (g − xH )

(6)

kde xH = arg maxx∈S f (x) a g je těžiště zbývajících d bodů simplexu S, U je náhodná veličina rovnoměrně rozdělená na [s, α − s), α > 0 a s jsou vstupní parametry, 0 < s < α/2. Ve dvou heuristikách se vybírá všech d + 1 bodů simplexu náhodně, a nastaví se různé hodnoty vstupních parametrů znáhodněné reflexe, α = 2, s = 0.5, resp. α = 5, s = 1.5. Ve třetí heuristice je jeden bod simplexu vždy bod s nejmenší funkční hodnotou v celé populaci P a zbývajících d se vybere náhodně, α = 2, s = 0.5. Čtvrtá heuristika vychází z diferenciální evoluce [15]. Nový bod y se vytváří křížením bodu (vektoru) u s náhodně vybraným bodem x. Bod u se generuje jako u = r1 + F (r2 − r3 ) (7) kde r1 , r2 a r3 jsou tři náhodně vybrané vzájemně různé body z populace P , F > 0 je vstupní parametr. Prvky yj , j = 1, 2. . . . , d vektoru y se vytvoří křížením x a u podle pravidla ( uj když Uj ≤ C nebo j = l (8) yj = xj když Uj > C a j 6= l , kde l je náhodně zvolená hodnota indexu z {1, 2, . . . , d}, U1 , U2 , . . . , Ud jsou nezávislé náhodné veličiny rovnoměrně rozdělené na [0, 1), C ∈ [0, 1] je vstupní parametr ovlivňující počet prvků vektoru y změněných křížením. Hodnota parametru F se určuje adaptivně postupem navrženým v [1] ( |<1 |) když | ffmax max(Fmin , 1 − | ffmax min min (9) F = fmin max(Fmin , 1 − | fmax |) jinak, kde fmin , fmax jsou minimální a maximální hodnoty účelové funkce v populaci a Fmin je vstupní parametr, který zajišťuje, aby F ∈ [Fmin , 1). Hodnoty vstupních parametrů této heuristiky jsou nastaveny Fmin = 0.4 a C = 0.9. V úlohách odhadu regresních parametrů je výhodné formulovat podmínku ukončení prohledávání pomocí rozdílu indexu determinace v aktuální populaci, neboť tím vyloučíme vliv jednotek, ve kterých jsou 2 2 měřeny hodnoty veličiny Y , na ukončení hledání. Podmínka je tedy Rmax − Rmin < ε, ε > 0 je vstupní 2 2 parametr algoritmu, Rmax a Rmin jsou maximální a minimální hodnoty indexu determinace R2 v aktuální populaci P . Index determinace je dán vztahem Q(β) . 2 i=1 (Yi − Y )

R2 = 1 − Pn

(10)

Pro další řídící parametry algoritmu CRS se soutěží těchto čtyř lokálních heuristik se osvědčily následující hodnoty: velikost populace N = 10 d, ε = 1 × 10−15 , w0 = 0.5 a δ = 0.05. V dalším textu je algoritmus s tímto nastavením označen CRS4HC.

4. Adaptace podmínky ukončení Další možnost, jak zvýšit adaptaci algoritmu pro odhad parametrů nelineární regrese, je v adaptivní formulaci podmínky ukončení. Ve většině úloh, kde data byla získána experimentem, jsou indexy determinace takové, že 1 − R2 > 1 × 10−4 . Pak doporučená hodnota ε = 1 × 10−15 je zbytečně nízká a jen zvyšuje časovou náročnost prohledávání. Bohužel správnou hodnotu indexu determinace neznáme a dopředu tedy nemůžeme vyloučit, že zvolená hodnota ε = 1 × 10−15 je naopak příliš vysoká. Hodnotu ε lze však nastavovat adaptivně podle aktuálně dosaženého indexu determinace. Hodnota ε na počátku je rovna zadané počáteční hodnotě ε0 a může být postupně snižována, pokud 1 − R2 < γ × ε, kde γ À 1 je další vstupní parametr, který ovlivňuje tuto adaptaci. Algoritmus takového adaptivního nastavení podmínky ukončení je zapsán v pseudokódu jako Algoritmus 2. Příkaz na řádku 14 zabraňuje nekonečnému 8.

NÁRODNÍ KONFERENCE

STATISTICKÉ DNY

V

BRNĚ, BRNO 27. - 28.

ČERVNA,

2006

smyčce vnějšího cyklu while. Ve variantě algoritmu, která je v dalším textu označena CRS4HCe, byly zvoleny hodnoty parametrů řídících ukončení ε0 = 1 × 10−9 a γ = 1 × 107 , ostatní řídící parametry jsou stejné jako v algoritmu CRS4HC. Algoritmus 2. Adaptace podmínky ukončení – CRS4HCe 1 generate population P ; 2 ε = ε0 ; 3 pass = false; 4 while 1 − R2 < γ × ε & pass is false (outer loop) 2 2 − Rmin 5 while Rmax > ε (inner loop) 6 generate a new trial point y ∈ D; 7 if f (y) < f (xmax ) then 8 xmax := y; 9 find new xmax ; 10 endif 11 pass = true; 12 endwhile; (end of inner loop) 13 if pass is false then γ = γ/10 endif ; 14 if 1 − R2 < γ × ε & pass is true then 15 ε = ε/10; pass = false; 16 endif ; 17 endwhile; (end of outer loop)

5. Porovnání algoritmů pro odhad parametrů nelineární regrese V tabulce 1 je uvedeno porovnání několika optimalizačních algoritmů užitých v odhadu parametrů nelineárních regresních modelů metodou nejmenších čtverců. Uvedeny jsou výsledky pro úlohy vyšší obtížnosti z databáze NIST [11]. Porovnávány byly deterministický Levenberg-Marquardtův algoritmus, který je jako funkce nlinfit součástí statistického toolboxu Matlab [5], a čtyři stochastické algoritmy. CRS1H je řízené náhodné prohledávání s reflexí simplexu jako jedinou lokální heuristikou, DER je standardní diferenciální evoluce a doporučeným nastavením vstupních řídících parametrů [15], CRS4HC a CRS4HCe jsou algoritmy se čtyřmi soutěžícími heuristikami popsané v předcházejícím textu. U stochastických algoritmů bylo pro každou úlohu provedeno sto opakování, ve kterých byla sledována časová náročnost potřebná k dosažení podmínky ukončení vyjádřená počtem vyhodnocení účelové funkce (2) a spolehlivost nalezení správného řešení vyjádřená počtem číslic, ve kterých se nalezené řešení reziduálního součtu čtverců shoduje se správným řešením [6]. V tabulce 1 jsou pak uvedeny průměrné počty vyhodnocení účelové funkce ve sloupci ne a a ve sloupcích R procento opakování, ve kterých bylo nalezeno řešení úlohy shodující se správným alespoň ve čtyřech číslicích. Pro snadnější porovnání časových nároků jsou Tabulka 1: Porovnání algoritmů v odhadu parametrů nelineární regrese Algoritmus Úloha bennett5 boxbod eckerle4 mgh09 mgh10 rat42 rat43 thurber

nlinfit výsledek L X L L L OK OK OK

CRS4HC

CRS4HCe

DER

CRSH1

R

ne

R

rne

R

rne

R

rne

100 100 100 100 100 100 100 100

41335 1308 2629 10422 20761 2942 4807 13915

100 100 100 100 100 100 100 100

-11 -37 -35 -15 1 -35 -39 -30

5 100 100 100 0 100 100 0

190 206 140 1427 478 474 904 1912

1 90 94 0 0 81 87 1

-69 -1 3 184 -50 42 83 824

ve sloupcích rne uvedeny průměrné hodnoty relativního rozdílu ne vzhledem k CRS4HC v procentech, takže záporné hodnoty znamenají menší časovou náročnost než má CRS4HC, kladné hodnoty naopak náročnost vyšší. Ve sloupci nlinfit výsledek L znamená ukončení iterace mimo globální minimum, X je numerický kolaps procedury a OK je správný výsledek. Vidíme, že Levenberg-Marquardtův algoritmus 8.

NÁRODNÍ KONFERENCE

STATISTICKÉ DNY

V

BRNĚ, BRNO 27. - 28.

ČERVNA,

2006

selhal v pěti úlohách z osmi a podobně neúspěšné byly i algoritmy DER a CRSH1, zatímco adaptivní algoritmy CRS4HC a CRS4HCe se soutěžícími lokálními heuristikami nacházely v těchto úlohách úspěšně globální minimum s přijatelnými časovými nároky (1000 vyhodnocení funkce na PC trvá okolo 1 vteřiny). Adaptace podmínky ukončení v CRS4HCe navíc ve většině úloh zmenšila časovou náročnost oproti CRS4HC zhruba o třetinu.

6. Závěr Abychom nevyvolávali přehnaně optimistická očekávání o efektivnosti a spolehlivosti adaptivních stochastických algoritmů, je nutno uvést, že pro úlohy nižší a střední obtížnosti nebyla převaha adaptivních algoritmů CRS4HC a CRS4HCe tak výrazná. Podrobnější výsledky stochastických algoritmů v odhadu parametrů nelineárních regresních modelů na všech úlohách NIST jsou uvedeny v [20] nebo v [21], kde jsou také definována ohraničení prohledávaného prostoru D pro jednotlivé úlohy. Úloha globální optimalizace nemá žádné definitivní řešení, a proto také doporučení týkající se stochastických algoritmů pro odhad parametrů nelineárních regresních modelů nemohou být jednoznačná. Je nutné si vždy uvědomovat, že stochastické algoritmy hledají řešení heuristicky, tzn. bez záruky nalezení správného řešení. To ovšem nezaručují ani dosud užívané deterministické algoritmy, které konvergují k nejbližšímu lokálnímu minimu. Důležité poznatky o vlastnostech stochastických algoritmů jsou založeny na experimentálních výsledcích, takže poznání postupuje pomalu a nikdy nebude uzavřeno. Úspěšnost adaptivního algoritmu CRS4HCe na testovaných úlohách je však dobrým důvodem, aby byl v aplikacích nelinární regrese využíván. Zdrojový text adaptivního algoritmu CRS4HCe pro odhad parametrů nelineární regresních modelů v Matlabu [5] je dostupný k volnému využití na webové stránce autora článku (http://albert.osu.cz/tvrdik/). Výzkum adaptivních stochastických algoritmů nepochybně bude pokračovat a zřejmě budou pak k dispozici algoritmy, které předčí CRS4HCe v úlohách odhadu parametrů nelineární regrese a budou vhodné i pro jiné aplikace ve výpočetní statistice. V tomto ohledu je velmi slibný adaptivní algoritmus diferenciální evoluce se soutěží různých nastavení řídících parametrů [22].

Literatura [1] Ali M. M., Törn A. (2004) Population set based global optimization algorithms: Some modifications and numerical studies. Computers and Operations Research 31, 1703 – 1725. [2] Bäck T. (1996) Evolutionary Algorithms in Theory and Practice. Oxford University Press, New York. [3] Holland J. H. (1975) Adaptation in Natural and Artificial Systems, The University of Michigan Press, Ann Arbor. [4] Křivý I., Tvrdík J. (1995) The Controlled Random Search Algorithm in Optimizing Regression Models. Comput. Statist. and Data Anal. 20, 229 – 234. [5] MATLAB, version 7, The MathWorks, Inc., 2005. [6] McCullough, B.D.,Wilson, B., (2005) On the accuracy of statistical procedures in Microsoft Excel 2003. Comput. Statist. and Data Anal. 49, 1244 – 1252. [7] Michalewicz Z. (1992) Genetic Algorithms + Data Structures = Evolution Programs. Berlin, Springer Verlag. [8] Militký, J. (1994) Nonlinear Regression on Personal Computers. In: R. Dutter and W. Grossmann (Eds.), COMPSTAT 1994. Proceedings in Computational Statistics, Physica-Verlag, Vienna, 395 – 400. [9] Mišík L., Tvrdík J., Křivý I. (2001) On Convergence of a Class of Stochastic Algorithms. In: Antoch, J. and Dohnal, G. (Eds.), ROBUST 2000, JČMF, Praha, 198 – 209. [10] Nelder J.A., Mead, R. (1964) A Simplex Method for Function Minimization. Computer J. 7, 308 – 313. [11] NIST Information Technology Laboratory. (2001) Statistical Reference Datasets. Nonlinear regression. http://www.itl.nist.gov/div898/strd/. [12] Price W.L. (1977) A Controlled Random Search Procedure for Global Optimization. Computer J. 20, 367 – 370.

8.

NÁRODNÍ KONFERENCE

STATISTICKÉ DNY

V

BRNĚ, BRNO 27. - 28.

ČERVNA,

2006

[13] Solis F. J. and Wets R. J-B., (1981) Minimization by random search techniques. Mathematics of Operations Research 6, 19 – 30. [14] Spall J. C. (2003) Introduction to Stochastic Search and Optimization, Wiley-Intersience. [15] Storn R., Price K. (1997) Differential Evolution – a Simple and Efficient Heuristic for Global Optimization. J. Global Optimization 11, 341 – 359. [16] Tvrdík J., Křivý I., Mišík L. (2001) Evolutionary Algorithm with Competing Heuristics. In: Matoušek R. and Ošmera P. (eds), MENDEL 2001, 7-th Int. Conference on Soft Computing, University of Technology, Brno, 58 – 64. [17] Tvrdík, J. Mišík L., Křivý I. (2002) Competing Heuristics in Evolutionary Algorithms. In: Sinčák P. et al. (Eds.), 2nd Euro-ISCI Intelligent Technologies - Theory and Applications, IOS Press, Amsterdam, 159 – 165. [18] Tvrdík J., Křivý I., (2004) Comparison of algorithms for nonlinear regression estimates. In: Antoch, J. (ed), COMPSTAT 2004. Physica-Verlag, Heidelberg New York, 1917 – 1924. [19] Tvrdík, J. (2004) Generalized controlled random search and competing heuristics. In: Matoušek R. and Ošmera P. (eds), MENDEL 2004, 10th International Conference on Soft Computing, University of Technology, Brno, 228 – 233. [20] Tvrdík J., Křivý I., Mišík L. (2006) Adaptive population-based algorithm for global optimization. COMPSTAT 2006 (v tisku). [21] Tvrdík J., Křivý I., Mišík L. (2006) Adaptive population-based search: Application to nonlinear regression. IRAFM Research Report No. 98, University of Ostrava. http://ac030.osu.cz/irafm/resrep.html [22] Tvrdík J. (2006) Competitive Differential Evolution. In: Matoušek R. and Ošmera P. (eds), MENDEL 2006, 12th International Conference on Soft Computing, University of Technology, Brno, 7 – 12.

Poděkování: Tato práce byla podporována grantem GAČR 201/05/0284.

8.

NÁRODNÍ KONFERENCE

STATISTICKÉ DNY

V

BRNĚ, BRNO 27. - 28.

ČERVNA,

2006