6. SRA a sztochasztikus programozásban 99
6.6. A kétlépcsős feladat hatékony numerikus megoldása
6.6.1. Numerikus stabilitás
A számítások numerikus stabilitásának elérése lényegében két, egymásnak ellentmondó cél megvalósítását jelenti. Az egyik cél az, hogy azSk halmaz „elég széles” legyen azért, hogy az előállítottxi, i= 0,1, . . . , k pontok (illetőleg az általuk felvett függvényértékek) segít-ségével meghatározandó függvény egy pozitív definit mátrix által megadott kvadratikus függvény alakját vegye fel (kifeszítsen egy konvex függvényt). A másik célunk az, hogy a pontok nagy része lehetőleg az optimális megoldás (megoldások) egy szűk környezetében legyenek, hiszen ezáltal lesz a közelítésünk pontos az optimum környékén. A két cél együttes megvalósításának kulcsa abban rejlik, hogy néhány „messze lévő” pont elég a konvex függvény kifeszítéséhez, míg az összes többi pont lehet az optimum környékén.
A kezdeti ponthalmaz jó megadásának, a pontrendszer adaptív kiszélesítésének és a degeneráció elleni védekezésnek a célja az, hogy a ponthalmaz elég széles legyen (a közelítés stabil legyen), míg a legkisebb négyzetek módszerében alkalmazott súlyok használata azt eredményezi, hogy elég szűk lesz a ponthalmazban szereplő pontok nagy része által alkotott halmaz. Ebben a részben azxelső lépcsős döntési változó dimenzióját egyszerűen n-nel jelöljük.
Kezdeti halmaz
Az algoritmus elindításához egy kezdeti ponthalmazra van szükségünk, ezt a halmazt például a következőképpen lehet előállítani. Tekintjük a kétlépcsős feladat definiálásánál használt determinisztikus alapfeladatot, ahol aξ valószínűségi változót helyettesítettük a várható értékével. Tehát megoldjuk a
min c0x +qy
f.h. Ax = b,
Tx +Wy = Eξ, x, y ≥ 0,
feladatot. Jelölje ennek az optimális megoldásátx0, ezt nevezzük gyökérpontnak. A többi pontot a gyökérpont körül generáljuk véletlenszerűen: xi = x0 +δdi, i = 1, . . . , kin−1, aholδ tetszőleges (de a választás után rögzített), 0.1−0.5értékű volt, a di pontok pedig véletlenszerűen +1, vagy -1 értéket tartalmaztak minden komponensükben. Ezt az eljárást addig folytattuk, amíg a megfelelő kin számú pontot nem állítottuk elő: az eljárás az esetek többségében megfelelő kezdeti halmazt állított elő. Nagyméretű feladataink esetén (6.7 szakasz) a kezdeti pontoknak mintegy 10-20%-nyi részét megengedett megoldásoknak vettük, a többit véletlenszerűen és egymástól függetlenül generáltuk ezek környékén.
A kvadratikus regresszióban szereplő Dk mátrix pozitív definitségét a 6.2.4 alatt leírt módon biztosítottuk – ez még rosszul előállított kezdeti ponthalmazok esetén is működött, csak akkor valamivel több „messze lévő” pontot kellett előállítani (ezek száma a 10-et nem haladta meg). Tehát a kezdeti halmaz jó megadása nem volt annyira fontos, mint kezdet-ben gondoltuk, mert a ponthalmaz alább leírt adaptív kiszélesítése pótolja a hiányossá-gokat.
Pontok adaptív megadása
Mint láttuk a 6.2.4 pontban, a kvadratikus regressziós függvény meghatározása lényegében egynq×nq méretű mátrix invertálását jelenti, aholnq= 1 +n+n(n+ 1)/2. Az invertálás során néha stabilitási problémák léptek fel, amelyek kezelésére a következő adaptív mech-anizmust használtuk. Az előállított M mátrixot és annak inverzét összeszoroztuk, ha ez a szorzat az egységmátrixtól egy 10−3 hibahatárnál jobban tért el, akkor néhány újabb, a pillanatnyi optimális megoldástól egyre távolabb lévő pontot adtunk hozzá a ponthal-mazhoz, addig, amíg az invertálás nem lett sikeres.
Ezeket a pontokat a következő módon állítottuk elő. Generáltunk két, az egységgömb felületén egyenletes eloszlású, egymásra ortogonális u1,u2 vektort, és az Sk halmazhoz négy pontot csatoltunk: xk+1 =xk+Csu1,xk+2 =xk−Csu1,xk+3 = xk+Csu2,xk+4 = xk−Csu2, aholCs egy skálázási konstans volt (kezdetben a Cs értékének 0.05-öt válasz-tottunk). Ha az így kibővített ponthalmazon meghatározott mátrix, illetőleg az inverze megfelelő volt, akkor folytattuk az SRA algoritmust. Ha az MM−1 szorzat még mindig a hibahatárnál jobban különbözött az egységmátrixtól, akkor aCs értékét megdupláztuk, és újabb négy pontot adtunk hozzá az Sk halmazhoz. Ezt addig ismételtük, amíg az in-vertálás nem lett sikeres, aCskonstans értékét, amely bizonyos értelemben megadta a sta-bilitáshoz szükséges, „elég széles” ponthalmaz sugarát, későbbi felhasználásra megőriztük.
Az SRA algoritmus használatában, az általunk vizsgált példák esetében a ponthalmaz ilyen adaptív szélesítésére általában 2-8 alkalommal volt szükség, az algoritmus első 1000 lépésében – a későbbiekben már nem kellett ezt használni.
Degeneráció elleni védekezés
Az algoritmus végrehajtása során közelítő kvadratikus optimalizálási feladatokat kell megoldani, egy optimális megoldás esetén némely feltételek aktívak, mások pedig nem aktívak lesznek – ezek az aktív feltételek vagy nemnegatívitási, vagy pedig egyenlőtlenséggel megadott feltételek, vagy ezek valamilyen keveréke lesz. Emiatt (hacsak az aktív halmaz az al-goritmus előrehaladása során nem változik nagyon gyakran) a kapott közelítő optimális megoldások egy, az aktív feltételek által meghatározott L lineáris sokaságon fekszenek.
Mivel ezen lineáris sokaság nL dimenziós (legfeljebb (n − m1) dimenziós), így a qk(x) kvadratikus közelítésünk határértékben degenerálttá válik. Ugyan lehetséges lenne egy, aznL-dimenziósLsokaságban előállítani egy kvadratikus függvényt, de előre nem tudjuk
eldönteni, hogy adott esetben mi lesz az aktív halmaz az optimum közelében, illetőleg az hogyan változhat meg, ezért a megoldáshoz más utat követtünk.
Az ilyen, degenerációhoz közeli helyzetek elkerülése érdekében nem csak egyetlen xk pontot adunk hozzá az Sk halmazhoz az algoritmus 3. lépésében, hanem több pontot.
Ezeket a pontokat az utolsó közelítő optimális megoldás kis környezetében vesszük fel.
Jelölje a hozzáadott pontok számátnadd, ennek értékét a számítógépes futtatások alatt 2 és2n+ 2 közötti értéknek választottuk; az alábbi leírásban nadd = 2n+ 2pontot használ-tunk, bár a gyakorlatban minden nadd ≥ 5 értékre jó eredményeket kaptunk. Ennek a módosításnak egy további számítástechnikai előnye is volt: csökkent az elvégzendő munka mennyisége. Ugyanis nem minden egyes pont esetén határoztuk meg újra a közelítő kvadratikus függvényt, csak mindennadd-edik pont hozzáadása után.
Az nadd számú {xl}k+2n+1l=k+1 pótpontot a következő módon konstruáltuk meg. Ahogy szokásos, ak-adik iterációs lépésben megoldott közelítő feladat optimális megoldását xk -val jelöljük. Legyen adva egy véletlen ortonormalizált U ponthalmaz, ahol U={ui}ni=1, U egyenletes az ortogonális rendszerek között, uiuj = δij, ahol δii = 1, δij = 0, i 6=
j, i, j = 1,2, . . . , n. Miután az xk pontot hozzáadtuk Sk-hoz, csatoljuk még az xk+j = xk +γkuj, xk+n+j = xk −γkuj pótpontokat is, ahol j = 1, . . . , n, és γk ∼ Cs/√
k. Itt γk egy olyan változó, amelyet csak lassan engedünk nullához tartani, ennek az értéke kezdetben az előzőleg leírtCs konstans, vagyis ez az alacsonyabb dimenzióba való eséstől (degenerációtól) óvó állandó. Másrészt az1/√
k faktor azt eredményezi, hogy a pótpontok halmaza egyre közelebb kerül a pillanatnyi optimális megoldáshoz – így nem akadályozza a pontsorozat konvergenciáját.
A tapasztalatok szerint a stabilitást viszonylag könnyen el lehetett érni azzal, hogy a generált ortonormált vektorokat két csoportban állítottuk elő. Az első csoport tartal-mazta af(x) =c0x+qk(x)célfüggvény∇f(xk)gradiensét és néhány további vektort (itt n−1−nL darab, az L lineáris sokaságra merőleges vektort választottunk), a második csoportban pedig nL számú, L-beli vektort választottunk (természetesen a két csoport-ban lévő vektorok egymásra is merőlegesek voltak). Ennek a csoportosításnak az volt az értelme, hogy az első csoportban lévő vektorok a gradiens pontosabb meghatározását segítették elő, míg a második csoportban lévő vektorok segítségével az egymásutáni iterá-ciókban meghatározott közelítő optimum jobban tudott azL sokaságon belül elmozdulni.
Aznadd = 5számú pótpont esetén a gradienst és negatívját választottuk az első csoportba, a második csoportba pedig a két utolsó közelítő optimális megoldást összekötő egyenes irányvektorát + és - előjellel.
Hozzáadhatunk még egy pótpontot a következők szerint. Az algoritmusban állandóan számoltuk a közelítő minimális függvényértéket, az fappr.min. értékét, amelyet az xa.min pontban értünk el, ezért ez a pont a rendelkezésünkre áll, és nehézség nélkül csatolható azSk halmazhoz.
Tehát az nadd = 2n + 2 számú pótpont a következőképpen áll össze: az xk közelítő megoldás, az ortonormált rendszer vektorjai (mind +, mind - előjellel) és a közelítő min-imum pontja. Minél nagyobb az nadd értéke, annál kevesebb első lépcsős feladatot kell
megoldani, viszont ez lassíthatja a konvergenciát. Ezért érdemes az algoritmus elején csak néhány pontot hozzáadni a ponthalmazhoz, míg a vége felé akár10n−100n számú pótpont is használható.
A közelítő minimum
Az algoritmus működése során könnyen számon tarthatjuk a közelítő minimum értékét – ezt az algoritmus 3. lépésében, közvetlenül az új függvényérték kiszámítása után tesszük meg. Az algoritmusban előállított közelítő optimális megoldások indexét jelöljük K-val, ugyanis az Sk halmazból csak az első lépcsős feladat determinisztikus feltételeinek eleget tevő xk, k ∈ K pontok függvényértékeit vesszük itt figyelembe. Jelölje az xk, k ∈ K pontban kiszámított zajos qk függvényérték szórását σk, amelyet a tényleges mintából számított empírikus szórás értékével azonosítunk. Az fappr.min. közelítő minimum értéke kezdetben legyen egy nagy szám, és utána ennek értékét minden egyesk ∈K index esetén felfrissítjük a
fappr.min. =qk+ 3σk, ha qk+ 3σk≤fappr.min.
egyenlet alapján; ha az egyenletben megadott feltétel nem áll fenn, akkor megnöveljük azfappr.min. értékét (lásd lentebb). Jelöljexa.minazt a pontot, amely esetén ezt a közelítő minimum értékét találtuk; azfappr.min.a tényleges minimum egy felső korlátjának tekinthető.
Egy kis (0.05) valószínűséggel a kiszámított qk függvényérték hibája nagyobb lehet 3σk-nál is, ha ezt nem vesszük figyelembe, akkor ugyanazt a hibás minimum pontot újra és újra csatolnánk az Sk halmazhoz. Ezért rugalmasabbá tesszük a felfrissítést: ha a jelenlegi közelítő minimumot nem csökkentettük egy k ∈ K értékre, akkor megnöveltük egy kis értékkel az
fappr.min. =fappr.min.+ 0.1σk
egyenlet alapján, aholk az az index, amely esetén az utolsó alkalommal csökkentettük a közelítő minimum értékét. Ezzel a módosítással előbb, vagy utóbb kimozdulunk az aktuális xa.min pontból (persze még akkor is, ha ez tényleg a minimum pont volt).
Az SRA viselkedése és a súlypont
Néhány heurisztikus gondolattal világítjuk meg az algoritmus viselkedését. Ez a SUMT eljárásokra emlékeztet – egy korlátot alkot a xk, k ∈ K pontok számára. Az algoritmus egyxkközelítő megoldást számít ki, és aztán újra ellenőrzi, aqk∼ Q(xk,ξk)függvényérték kiszámításával a pont optimalitásának helyességét. Ha qk kisebb, mint a regressziós füg-gvény értéke ezen a helyen, akkor lefelé módosítja ebben a pontban aqk(x)függvényt (ha nagyobb, akkor felfelé). Vagyis egy idő után egy többé-kevésbé stabil tálat alakít ki, és a xk, k ∈K pontok csak ebben tudnak mozogni.
Véletlent használó algoritmusok esetében általában jó ötlet a súlypont vizsgálata.
Jelöléseinkkel ez az M1 = (1/P λi)P
λixi pontot jelenti, ahol az összegezés az i =
0, . . . , k indexekre történik. A közelítés stabilizálása alatt létrehozott „messze” lévő pon-tok hatása miatt ez a súlypont lényegesen eltérhet az optimális megoldástól, és elég sok időt vehet igénybe ennek a konvergenciája az optimális megoldáshoz. Ezek szerint érdemes egy módosított súlypontot alkotni, amelyben nem szerepeltetjük az algoritmus elején meghatározott, főleg a stabilizálást szolgáló „messze” lévő pontokat, vagyis legyen
s= (1/
A legtöbb esetben ez az s pont jobb eredményeket adott (kisebb célfüggvényértéket), mint az utolsó xk közelítő megoldás, hiszen az előbbi a tál közepén lassan változtatta a helyét, míg az utóbbi jóval nagyobb ingadozásokkal próbált meg a tál szélének nekimenni.