Havaitaan oppositiossa olevaa asteroidia. Pyörimisestä huolimatta sen kirkkaus ei muutu. Projisoitu pinta-ala pysyy ilmeisesti vakiona. Onko asteroidi siis pyörähdyssymmetrinen?
Ei välttämättä, se voi olla esimerkiksi Reuleaux'n kolmio: \vskip7cm
Vastaavat ongelmat ovat tähtitieteessä tavallisia, koska havaintogeometriaa ei voi juuri muuttaa.
Suora menetelmä
Konstruoidaan malli havaittavalle ilmiölle esimerkiksi simuloimalla sitä tietokoneella. Lasketaan, millaisia havaintoja saadaan, kun mallia kuvaaville parametreille annetaan tietyt arvot. Parametreja muutellaan, kunnes mallista lasketut havainnot vastaavat todellisia.
Ongelma: onko ratkaisu yksikäsitteinen?
Toinen ongelma: ratkaisun stabiilisuus. Usein havaitaan integroitu suure, johon mallin parametrien suuretkin muutokset vaikuttavat vain vähän. Pienetkin havaintovirheet johtavat hyvin erilaisiin parametrien arvoihin.
Käänteinen menetelmä
Havainnoista lasketaan kohdetta kuvaavat parametrit. Usein matemaattisesti vaikeaa.
Havaintoprosessiin liittyy usein todellisen suureen integrointi, joten informaatiota häviää. (Tieto, että x+y=1, ei kerro mitään luvuista x ja y).
Integrointi on tasoittava operaatio, johon integroitavan pienet häiriöt eivät juuri vaikuta. Käänteinen operaatio vaatii aina jonkinlaista derivointia, joka puolestaan vahvistaa häiriöitä. Siksi pienistä havaintovirheistä voi aiheutua suuria virheitä havainnoista johdettuun malliin.
g(x) = \intab K(x,t) f(t) dt.
Tämä on Fredholmin ensimmäisen lajin integraaliyhtälö. Toisen lajin yhtälö on
f(x) = g(x) + \lambda \intab K(x,t) f(t) dt.
Integrointirajat ovat kiinteät. Jos integrointiväli riippuu x:stä, saadaan Volterran ensimmäisen ja toisen lajin integraaliyhtälöt.
Jos toisen lajin yhtälöissä \lambda on pieni, funktiot f ja g ovat likimain samoja, ja integraalitermi kuvaa pientä häiriötä. Yhtälö voidaan silloin ratkaista iteroimalla. Aluksi voidaan valita esim. f0=g.
fk+1(x) = g(x) + \lambda \int K(x,t) fk(t) dt.
Menetelmä ei toimi, jos \lambda on suuri. Silloin alkuperäisen funktion osuus jää pieneksi ja hallitsevaksi tulee integraalin tasoittava vaikutus. Yhtälö lähenee tällöin ensimmäisen lajin yhtälöä.
Hankalimman tapauksen muodostavat Fredholmin ensimmäisen lajin yhtälöt. Diskreetissä tapauksessa tällaiset yhtälöt voidaan kirjoittaa matriisimuotoon
g = H f .
Esimerkiksi kuvankäsittelyssä H on oleellisesti pistemäisen lähteen kuvan leviämistä esittävä funktio (PSF). Jos signaali pääsee laitteistosta lävitse muuttumattomana, H on yksikkömatriisi. Jos signaali leviää hiukan, H on nauhamatriisi, jossa nollasta poikkeavat alkiot sijaitsevat lävistäjällä ja sen lähistöllä. Mitä laajemmalle signaali leviää, sitä kauempaa H:n lävistäjästä löytyy nollasta eroavia alkioita ja samalla lävistäjän alkiot pienenevät.
Kun signaali on täysin levinnyt, ovat matriisin H kaikki alkiot suunnilleen yhtäsuuria. Matriisi on siten lähes singulaarinen ja sen kääntäminen on erittäin epästabiili operaatio. Tämä vastaa myös tilanteen fysikaalista tulkintaa: signaalit ovat levinneet ja sekoittuneet toisiinsa niin pahasti, että alkuperäinen muoto ei enää ole lainkaan näkyvissä.
Esimerkki:
\int01 (x+y) f(y) dy = x.
Jos tähän sijoitetaan yrite f(y)=ay+b, saadaan vasemman puolen arvoksi
(a/2 + b) x + a/3 + b/2.
Jotta tämä olisi identtisesti x, täytyy valita a=-6 ja b=4. Yhtälön ratkaisu on siten
f(y)=4-6y.
Yritetään ratkaista yhtälö numeerisesti. Lasketaan integraali trapetsisäännöllä, jolloin integraali on
(h/2) [ x f(0) + 2(x+h)f(h) + ... + 2(x+(n-1)h)f((n-1)h) + (x+1)f(1) ].
Lauseke sisältää n+1 tuntematonta, nimittäin funktion f arvot f(0), f(h), ..., f(1). Yhtälö voidaan kirjoittaa erikseen eri x:n arvoille, jolloin saadaan riittävä määrä yhtälöitä näiden tuntemattomien ratkaisemiseksi.
Valitaan h=0.2, jolloin tuntemattomia funktion arvoja on kuusi kappaletta. Käytetään x:n arvoja 0, 0.2, ..., 1.
(0.2/2) [0 + 0.4f(0.2)+0.8f(0.4)+1.2f(0.6)+1.6f(0.8)+f(1) ] = 0,
...
(0.2/2) [1 + 2.4f(0.2)+2.8f(0.4)+3.2f(0.6)+3.6f(0.8)+2f(1)] = 1
eli matriisimuodossa
| 0 0.4 0.8 1.2 1.6 1.0 | | f(0) | | 0 |
| 0.2 0.8 1.2 1.6 2.0 1.2 | | f(0.2) | | 0.2 |
0.2 | 0.4 1.2 1.6 2.0 2.4 1.4 | | f(0.4) | | 0.4 |
--- | 0.6 1.6 2.0 2.4 2.8 1.6 | | f(0.6) | = | 0.6 |
2 | 0.8 2.0 2.4 2.8 3.2 1.8 | | f(0.8) | | 0.8 |
| 1.0 2.4 2.8 3.2 3.6 2.0 | | f(1) | | 1 |
Kerroinmatriisi on melkein singulaarinen. Ratkaisulla ei välttämättä ole mitään tekemistä todellisuuden kanssa.
Regularisointi
Singulaarinen matriisi liittyy tilanteeseen, jossa informaatiota on menetetty. Kadonnutta informaatiota ei pystytä palauttamaan, mutta matriisin kääntämistä voidaan kuitenkin yrittää saada stabiilimmaksi.
Regularisoinnin ideana on saada käännettävän matriisin lävistäjälle muita suurempia alkioita.
Kirjoitetaan edellä olleen esimerkin yhtälö muotoon
\lambda f(x)+\int01 (x+y) f(y) dy = x,
missä regularisointiparametri \lambda on jokin positiivinen vakio. Kun tämä diskretoidaan samalla tavoin kuin edellä, saadaan yhtälöryhmä
\lambda f(0)+(0.2/2) [0 + 0.4f(0.2)+0.8f(0.4)+
1.2f(0.6)+1.6f(0.8)+f(1) ] = 0,
..
\lambda f(1)+(0.2/2)} [1 + 2.4f(0.2)+2.8f(0.4)+
3.2f(0.6)+3.6f(0.8)+2f(1)] = 1.
Kerroinmatriisin lävistäjälle tulevat luvut ovat nyt parametrin \lambda verran suurempia. Häiriöalttius on pienempi, ja yhtälö voidaan ratkaista ilman ongelmia.
Kun \lambda on suuri, uusi yhtälö poikkeaa huomattavasti alkuperäisestä; samoin sen ratkaisu on kaukana alkuperäisen yhtälön ratkaisusta.
On ilmeisesti olemassa jokin regularisointiparametrin \lambda optimiarvo, jolla ratkaisu on mahdollisimman lähellä alkuperäisen yhtälön ratkaisua, mutta matriisin singulaarisuus ei ole vielä alkanut aiheuttaa omia häiriöitään.
Esimerkkitapauksen ratkaisut eri $\lambda$:n arvoilla:
\lambda cond(H) f(0) f(1)
0.05 43 9.1 -5.6
0.01 123 4.3 -2.1
0.005 242 4.0 -1.9
0.001 1200 3.8 -1.8
Paras tulos saadaan, kun \lambda=0.005. Suuremmilla vakion arvoilla yhtälö poikkeaa liikaa alkuperäisestä ja pienemmillä arvoilla kerroinmatriisi alkaa lähestyä singulaarista.
Yksi mahdollinen keino on muuntaa alkuperäinen tehtävä optimointiongelmaksi. Yhtälöllä Hf=g ei välttämättä ole lainkaan ratkaisua, tai sen ratkaisu on hyvin herkkä havaintovirheille. Aina voidaan etsiä sellaisen ratkaisun f, joka antaa minimiarvon normille
|| H f - g ||.
Jos kerroinmatriisi on singulaarinen, ratkaisu ei ole yksikäsitteinen. Tämä ongelma voidaan korjata regularisoinnilla minimoimalla esimerkiksi normia
|| H f - g ||2 + \lambda2 || f ||2.
Tämän minimointitehtävän ratkaisuna saadaan sekä vektori f että parametri \lambda.
Haettava funktion f(x) minimi, x=(x1, ... , xN). Funktio f on kohdefunktio (objective function).
Jos tehtävään ei liity mitään lisäehtoja, kyseessä on rajoitteeton optimointitehtävä.
Rajoitetussa tehtävässä ratkaisu on lisäksi toteutettava annetut rajoitusehdot, jotka voivat olla
Usean muuttujan derivaattoja käyttävissä menetelmissä tarvitaan toisista derivaatoista muodostuvaa Hessen matriisia
| \pd2f(x)/\pd x12 ... \pd2f(x)/ \pd x1 \pd xn |
H(x) = | ... |
| \pd2f(x)/ \pd xn \pd x1 ... \pd2f(x)/\pd xn2 |
Hyvin monet tehtävät voidaan muuntaa optimointitehtäviksi:
1. Funktion derivaatat tunnetaan Esim. liittogradienttimenetelmät.
2. Derivaattoja ei tunneta tai ne ovat hankalia laskea. Esim. polytooppimenetelmä.
Globaali optimointi
Täysin luotettavia algoritmeja ei ole.
Tunnetaan funktion arvo kolmessa pisteessä a < b < c, s.e. f(b) < f(a) ja f(b) < f(c).
Valitaan piste x esimerkiksi väliltä (b, c). Jos f(x) < f(b), uusi väli on (b,c), muuten (a,x).
Toistetaan sama toimenpide, kunnes f(x) ei enää pienene.
Minimin kohdalla f'(x)=0, joten f:n muutos on korkeintaan toista kertalukua. Iterointi voidaan lopettaa, kun välin pituus on luokkaa \sqrt\epsilon, missä \epsilon on konevakio. Tätä pienemmät x:n muutokset eivät enää näy funktion arvoissa.
Kultaisen leikkauksen sääntö: välin keskimmäinen piste valitaan aina niin, että
(b-a)/(c-a) = (3-\sqrt 5) / 2 = 0.38197.
Seuraava piste valitaan pitemmältä osaväliltä s.e.
(x-b)/(c-b) = 0.38197.
Yksiulotteinen tapaus; derivaatta tunnetaan.
Oletetaan taas, että a < b < c ja f(b) < f(a) ja f(b) < f(c).
Lasketaan f'(b). Jos f'(b)>0, seuraava piste valitaan väliltä (a,b), muuten väliltä (b,c).
Kun on laskettu derivaatta kahdessa pisteessä, voidaan esim. sekanttimenetelmällä laskea seuraava piste, jossa derivaatan pitäisi olla nolla. Jos näin laskettu piste on välin ulkopuolella, otetaan uudeksi pisteeksi välin keskipiste.
N-ulotteinen tapaus
Edetään johonkin suuntaan niin kauan kuin kohdefunktion arvo pienenee, jonka jälkeen muutetaan suuntaa. Tätä viivahakua toistetaan, kunnes ei enää löydy pienempiä arvoja.
Kullakin askeleella ratkaistaan yksiulotteinen minimointitehtävä. Eteneminen tiettyyn suuntaan päättyy, kun kohdefunktion gradientilla ei ole etenemissuunnan suuntaista komponenttia.
Etenemissuunta voidaan valita eri tavoin.
1. Edetään vain koordinaattiakselien suuntaan; muutetaan vuorotellen koordinaatteja x1, ... , xN, x1, .... Usein tehotonta.
2. Edetään jyrkimmän muutoksen (eli kohdefunktion gradientin) suuntaan. Kunkin viivahaun jälkeen käännytään aina 90 astetta. Ei välttämättä johda kovin suoraan kohti minimiä. \vskip 8cm
3. Yritetään valita seuraava suunta siten, että otetaan huomioon aikaisemmat suunnat.
Liitto- eli konjugaattigradienttimenetelmä
Olkoon x pisteen paikka toistaiseksi parhaan minimikohdan p suhteen. Silloin on
f(x) \approx
f(p) + \sum_i \pd{f}{x_i}x_i
+{1\over2} \sum_{i,j} \pdd{f}{x_i}{x_j} x_ix_j
= c + {\bf b}.x+{1\over2}x.A\cdotx,
missä b on funktion f gradientti pisteessä p ja A Hessen matriisi pisteessä p. Tästä derivoimalla saadaan gradientille lauseke (pisteessä x)
\nabla f = A . x + b.
Kun pisteestä x lähdetään johonkin suuntaan v, gradientin muutos on A.v.
Oletetaan, että edellinen askel on kuljettu suuntaan u. Jotta seuraava askel ei heikentäisi tästä saatua hyötyä, gradientin täytyy olla kohtisuorassa vektoria u vastaan myös askelen v jälkeen:
u . A . v = 0.
Tällaisia suuntia u ja v sanotaan toistensa konjugaateiksi.
1. Valitaan jokin alkuarvo x1.
2. Toistetaan seuraavaa, kun k=1,2, ..., kunnes ||gk|| on riittävän pieni.
\betak = (gk . gk)/ (gk-1 . gk-1)
Polytooppimenetelmä
Nelderin ja Meadin polytooppimenetelmä; myös simpleksimenetelmä (eri asia kuin lineaaristen opitmointitehtävien simpleksimenetelmä). Kohdefunktion derivaattoja ei tarvita.
N-ulotteisen avaruuden simpleksi on monitahokas, jossa on N+1 kärkeä x1, ... , xN+1. Merkitään fi=f(xi), fl= min fi, fh= max fi.
Huonointa kärkeä vastaavan tahkon painopiste
x_c = (1/n) \sumi \ne hxi.
Peilauskerroin \alpha=1, pienennyskerroin \beta=0.5, laajennuskerroin \gamma=2.
Simpleksiä muunnellaan kolmenlaisilla operaatioilla:
xr = (1+\alpha)xc - \alpha xh.
xe = (1-\gamma)xc + \gamma xr.
xs = (1-\beta)xc + \gamma xh.
program amoeba
! Nelder-Mead polytooppi- eli simpleksimenetelma
! muunneltu Numerical Recipes -kirjasta
implicit none
integer, parameter :: NMAX = 1000, & ! askelten maksimimaara
ndim = 2, & ! avaruuden dimensio
mpts = 3 ! simpleksin karkien maara
real, parameter :: FTOL = 0.001, & ! lopetustarkkuus
ALPHA = 1.0, & ! peilauskerroin
BETA = 0.5, & ! pienennyskerroin
GAMMA = 2.0 ! venytyskerroin
integer i,j,ilo,ihi,inhi, nfunk
real ytry, ysave, sum, rtol, psum(ndim), &
p(ndim+1, ndim), &
y(ndim+1)
! ensimmaisen simpleksin karjet
p(1,:) = (/ 0.0, 0.0 /)
p(2,:) = (/ 0.0, 1.0 /)
p(3,:) = (/ 1.0, 0.0 /)
! funktion arvot simpleksin karjissa
do i=1,3
y(i) = funk(p(i,:))
end do
nfunk=0
! koordinaattien summat; karjen i vastakkaisen
! simpleksin painopiste on (psum-p(i))/ndim
do j=1,ndim
sum=0.0
do i=1,mpts
sum = sum+p(i,j)
end do
psum(j) = sum
end do
main: do
! haetaan paras (ilo), huonoin (ihi) ja
! toiseksi huonoin (inhi) karkipiste
ilo=1
if (y(1) > y(2)) then
inhi=2; ihi=1
else
inhi=1; ihi=2
end if
do i=1,mpts
if (y(i) < y(ilo)) ilo=i
if (y(i) > y(ihi)) then
inhi=ihi; ihi=i
else if (y(i) > y(inhi)) then
if (i /= ihi) inhi=i
end if
end do
! jos huonoin ja paras arvo lahes samoja, lopetetaan
rtol=2.0*abs(y(ihi)-y(ilo)) / (abs(y(ihi))+abs(y(ilo)))
if (rtol < FTOL) then
write(6,'("rtol=",f10.5)') rtol
write(6,'("nfunk=",i5)') nfunk
exit main
end if
! jos iteraatiokierroksia liikaa, lopetetaan
if (nfunk >= NMAX) then
write(6,'("nfunk=",i5)') nfunk
exit main
end if
! peilaus
ytry=amotry(-ALPHA)
! jos tulos parani, siirretaan uutta karkea viela kauemmas
if (ytry <= y(ilo)) then
ytry=amotry(GAMMA)
else if (ytry >= y(inhi)) then
! uusi piste huonompi kuin toiseksi huonoin;
! pienennetaan simpleksia
ysave=y(ihi)
ytry=amotry(BETA)
if (ytry >= ysave) then
! edelleen liian iso; pienennetaan simpleksia
! parhaan pisteen suhteen
do i=1,mpts
if (i /= ilo) then
do j=1,ndim
psum(j)=0.5*(p(i,j)+p(ilo,j))
p(i,j)=psum(j)
end do
y(i)=funk(psum)
end if
end do
nfunk = nfunk+ndim
! paivitetaan koordinaattien summat
do j=1,ndim
sum=0.0
do i=1,mpts
sum = sum+p(i,j)
end do
psum(j) = sum
end do
end if
end if
end do main
! tulostetaan karkipisteet ja funktion arvot; kaikkien pitaisi
! olla melkein samoja, silla ratkaisun loytyessa simpleksi on
! surkastunut optimipisteeksi
write(*,*) p(1,:), y(1)
write(*,*) p(2,:), y(2)
write(*,*) p(3,:), y(3)
contains
! optimoitava kohdefunktio
real function funk(p)
implicit none
real p(ndim)
funk=(p(1)-3.0)**2+(p(2)-2.0)**2+1.0
end function
! simpleksin muunnosaskel
! huonoin karki peilataan vastakkaisen sivun suhteen ja
! etaisyytta muutetaan kertoimella fac
real function amotry(fac)
implicit none
real fac, fac2
integer j
real fac1, ytry, ptry(ndim)
fac1=(1.0-fac)/ndim
fac2=fac-(1-fac)/ndim
do j=1,ndim
ptry(j)=psum(j)*fac1+p(ihi,j)*fac2
end do
ytry=funk(ptry)
nfunk = nfunk+1
if (ytry < y(ihi)) then
y(ihi)=ytry
do j=1,ndim
psum(j) = psum(j)+ptry(j)-p(ihi,j)
p(ihi,j)=ptry(j)
end do
end if
amotry=ytry
end function
end program
Rajoitteet
Ratkaisuvektoria päivitettävä niin, että rajoitteet pysyvät aina voimassa.
Rajoitteellista tehtävää voidaan käsitellä rajoittamattomien tehtävien ohjelmilla, jos rajoitteet korvataan sakko- ja estefunktioilla.
Lisätään optimoitavaan funktioon termi, joka on hyvin suuri sallitun alueen ulkopuolella. Jos esimerkiksi on optimoitava f rajoitteella g(x)<=0, voidaan kohdefunktioksi ottaa
f(x) + s max {0, g(x)},
missä vakio s on sakkoparametri.
Estefunktio on funktio, joka kasvaa rajatta lähestyttäessä sallitun alueen reunaa.
Mikään menetelmä ei takaa globaalin minimin löytymistä.
Valitaan riittävän suuri määrä alkuarvoja, joista lähtien suoritetaan lokaali optimointi. Toivotaan, että jokin löytyneistä lokaaleista minimeistä on myös globaali.
Varsinkin työläissä kombinatorisissa tehtävissä geneettiset algoritmit ja simuloitu jäähdytys saattavat olla käyttökelpoisia.
Analogia: kun nestettä jäähdytetään riittävän hitaasti, se jäätyy säännölliseksi kidehilaksi, jonka energia on minimissään.
Kohdefunktio vastaa systeemin energiaa. Konfiguraatiota muutellaan satunnaisesti. Uusi konfiguraatio valitaan lämpötilasta riippuvalla todennäköisyydellä.
Kun lämpötila on korkea, hyväksytään kohtalaisella todennäköisyydellä myös siirtymiset huonompiin (korkeamman energian) tiloihin. Tällä pyritään estämään juuttuminen lokaaleihin minimeihin.
Kun lämpötila laskee, todennäköisyys hyväksyä huonompi konfiguraatio pienenee.
Sopivia sovelluksia esim. kombinatoriset ongelmat, joille ei tunneta nopeita ratkaisumenetelmiä.
Löytää yleensä melko hyvän ratkaisun, mutta ei takaa sen optimaalisuutta.