Gyökkereső algoritmus

Legyenprel adott megbízhatósági szint és tekintsük az f(x_r) = p_rel

egyenletet, melynek pontos gyöke x_r – ez általánosan egy kvantilis meghatározási prob-léma. Feladatunk egy x_r közelítő gyök kiszámítása – ezt röviden gyökkeresési feladatnak nevezzük.

A leírásra kerülő heurisztikus algoritmus a t regressziós becslést használja közelítő gyök meghatározására, ahol t az előzőleg megadott t₁(x), t₂(x), t₃(x) becslések akárme-lyike lehet. Az eljárás működésének egy magyarázata a következő lehet. Ha a közelítő t függvény az f függvényhez közel van, akkor a t(x) = p_rel-ből meghatározható közelítő gyök is közel lesz a pontos gyökhöz. Az új x_i mintapontokat a közelítő gyök környékén vesszük fel, így a következő approximáló függvénynek remélhetőleg kisebb lesz a hibája az igazi gyök körül, tehát a következő közelítő gyök is közelebb kerül az igazi gyökhöz.

A továbbiakban egy iteratív eljárás segítségével keresünk egy, a pontos x_r gyökhöz tartó {x_r} közelítés–sorozatát; a közelítő gyök meghatározásához használt t(x) becslést fokozatosan egyre pontosabbá tesszük a közelítő gyök környékén. Ezt azáltal érjük el, hogy egy összehúzódó intervallum-sorozatot konstruálunk – ezeket a p értéket tartal-mazó, értékkészletekre vonatkozó [α_i, β_i] intervallumokat konfidencia-intervallumoknak nevezzük. A megfelelő,x értékeket tartalmazó[ai, bi]intervallumok a tⁱ(xr) =prel egyen-lőségből számítható x_r közelítő gyököket tartalmazzák, ahol tⁱ(a_i) = α_i, tⁱ(b_i) = β_i. Az intervallumok hosszát minden iterációs lépésben csökkentjük egy% <1 redukciós faktor-ral. Az algoritmus egyes lépéseire vonatkozó részletes megfontolásokat a [De 98a] cikkben közöltük, itt csak az algoritmusok végleges formáját közöljük.

A jelölések egyszerűsítése céljából az algoritmust csak at₁(x) =g₁(x) = a₁x²+b₁x+c₁ becslésre írjuk le (a szükséges apró módosítások elvégzése után az eljárás at₂ és at₃ bec-slésekre is alkalmazható), valamint a t₁, g₁, a₁, b₁, c₁ mennyiségek indexeit is elhagyjuk.

A redukciós faktor értékét % = 0.6-nek, az egyes iterációkban használt pontok számát K = 10-nek vettük. A kezdeti[a₀, b₀] intervallumot a tapasztalatok szerint szimmetriku-san kellene felvenni ax_r gyök körül, de mivel ezt nem ismerjük, egy nagyjából jó becslés elegendő itt. Az algoritmusban ezt a szimmetrizálást közelítőleg úgy érjük el, hogy a job-boldali és a baloldali deriváltak közelítő értékeineko hányadosát használjuk az algoritmus 2. lépésében.

Gyökkereső eljárás – az f(x) =p_rel egyenletre

0. [Előkészítés.] Tegyük fel, hogy a kezdeti [a₀, b₀] intervallum, az ebben lévő x_0j, j = 1,2, . . . , K pont, az ezekben felvett p_0j ∼f(x_0j) zajos függvényértékek adottak, valamint meghatároztuk már ezek alapján a g(x) becslést. Legyen kezdetben σ₀ =σ/%, az iterációs számláló i= 0, és a lépéshossz stl= (b0−a0)/10.

1. [A konfidencia-intervallum meghatározása.] Növeljük meg az

iterációs számlálót i=i+ 1, számítsuk ki a σ_i =%σi−1 értéket és az x_t=−b/(2a) csonkolási pontot.

2. [Az előzetes [a_i, b_i] intervallum kiszámítása.] Határozzuk meg a g(x_r) = p_rel egyenletből az x_r közelítő gyököt. Ha a <0, akkor elfogadjuk azt a gyököt, amely kisebb x_t-nél (ha nincs megoldás, akkor legyen xr=xt). Ha a >0, akkor elfogadjuk az xt-nél nagyobb gyököt (ha nincs megoldás, akkor legyen x_r =x_t). Számítsuk ki a δ= 2∗stl értéket és legyen x⁻=x_r−δ, x⁺ =x_r+δ. A jobb és baloldali deriváltak közelítő arányát az o= (p_rel−g(x⁻))/(g(x⁺)−p_rel) hányados mutatja. Legyen α_i = max(p_rel−o·σ_i,0.0001), β_i = min(p_rel+σ_i,0.9999).

3. [Az [a_i, b_i] intervallum meghatározása.] Számítsuk ki az a_i =g⁻¹(α_i), bi =g⁻¹(βi) gyököket. Ha a <0, akkor válasszuk ki azokat a gyököket, amelyek x_t-nél kisebbek, ha a >0, akkor válasszuk ki azokat a

gyököket, amelyek x_t-nél nagyobbak (ha nincs megoldás, akkor legyen a_i =x_t, b_i =a_i+ 0.01).

4. [A g(x) újraszámítása.] Válasszunk ki K = 10 darab új pontot az [a_i, b_i] intervallumban: x_ij = (b_i−a_i)(j−1)/(K −1) +a_i, j = 1, . . . , K és egy Monte Carlo eljárással határozzuk meg a pij ∼f(xij),

E(p_ij) = f(x_ij) zajos függvényértékeket. Az új pontok halmaza legyen S_i ={x_ij, p_ij}^K_j=1, és S=∪ⁱ_l=0 S_l. Határozzuk meg az új g(x)

approximációt az összes S-beli pont és függvényérték segítségével.

Határozzuk meg az új közelítő x^∗_r =g⁻¹(p_rel) gyököt (ha a <0, akkor vegyük a kisebbet, egyébként pedig a nagyobb értékű gyököt,

illetőleg legyen x^∗_r =xt, ha nincs megoldás).

5. [A konvergencia ellenőrzése.] Ha σ/p

(i+ 1)K < ε₀ fennáll, ahol ε₀ egy előírt hibatűrés, akkor STOP. Egyébként legyen x_r =x^∗_r, stl= 2(b_i−a_i)/K és menjünk vissza az 1. lépésre.

Néhány általános megjegyzést teszünk. Számítástechnikai szempontból az algoritmus-nak a következő két előnye látható:

(i) a számítógépes eredmények szerint [De 98a] az approximáció hibája azS={x_ij, p_ij} halmazban megadott pontok számának négyzetgyökével arányosan csökken, vagyis|f(x_r)−

f(x_r)| ≤σ/p

(N + 1)K, aholN az iterációk száma (ezt a tulajdonságot az algoritmus 5.

lépésében használtuk ki).

(ii) a g(x) becslésnek a 4. lépésben leírt újraszámolását nagyon egyszerűen lehet elvégezni, csak az eddig már meghatározott átlagokhoz kell hozzáadni azS_i-beli új pontok megfelelő hatványait (momentumait), vagyis nincs szükség arra, hogy minden iterációban teljesen előlről kezdjük a számításokat. Tehát a regressziós becslés meghatározásával kapcsolatos munka csak lineárisan függ azS-beli pontok (függvénykiszámítások) számától.

Ha a prel megbízhatósági szint 0.95-0.99 körül van, akkor a g(x) becslés lineáris csonkolását kell használni, mivel a becslés konstrukciója miatt a közelítő gyök nagyon közel lesz a csonkolási ponthoz, és a valódi gyök gyakran a csonkolási pontnál nagyobb.

Az 5. lépésben szereplő megállási kritérium helyett másmilyen megállási elvet is használhatunk; a legegyszerűbb az elvégzendő iterációk N számát előre lerögzíteni. Egy másik lehetőség lehet az, ha mintavétel hibáját számítjuk az egész algoritmus folyamán és a kívánt hibakorlát elérése esetén állítjuk le az algoritmust.

Eljárhatnánk úgy is a gyökkeresés algoritmusában, hogy a teljes m-dimenziós Φ füg-gvényre adunk meg egy m-dimenziós regressziós becslést, és ennek az egy egyenes menti (sokaságon felvett) egydimenziós függvényét használjuk a fenti eljárásban, de a szükséges munka mennyiségét nem indokolja meg a pontosságban elérhető növekedés, ezért ezzel nem foglalkozunk részletesen.

Néhány számítási eredményt közlünk annak megmutatására, hogyan működik a gyökkereső algoritmus rögzített N iterációszám esetén (további számítási eredmények találhatók a [De 98a] cikkben). A számítógépes futásokban a következőképpen állítottuk be az algo-ritmus állandóit: K = 10, N = 9, σ = 0.05, vagyis minden eljárásban 100 függvényérték kiszámítását végeztük el és a végső eredmény hibája a legtöbb esetben σ/10 = 0.005 körül volt. Az első példa esetében egy m = 2-dimenziós normális eloszlást használtunk, amelynek komponensei korreláltak voltak. A 3. példában pedig egy tíz-dimenziós nor-mális eloszlást használtunk, amely komponensei páronként korreláltak voltak (az első és a második komponensek egymással korreláltak voltak, de függetlenek az összes többi kom-ponenstől, a harmadik és a negyedik komponensek korreláltak voltak, de függetlenek a többitől, stb.) A 4. példában minden komponens független volt, a kapott numerikus eredményeket a 4.3 táblázatban közöljük.

Példa Dimenzió Megb. ”Pontos” gyök Becslés Gyök Hiba

száma m= p_rel = x_r x_r |f(x_r)−f(x_r)|

1. 2 0.80 1.1114 alap 1.099 0.0025

log 1.100 0.0035

f-log 1.131 0.0048

3. 10 0.95 5.895 alap 5.724 0.0024

log 5.827 0.0010

f-log 5.999 0.0012

4. 50 0.90 4.276 alap 4.201 0.0038

log 4.164 0.0050

f-log 4.216 0.0029

4.3. táblázat. A gyökkeresés eredményei.