torstai 4. toukokuuta 2017

Funktion kaikkien paikalliset ääriarvojen haku Matlab-työkaluin

Lopullinen julkaisupäivä: 15.05.2017
Matlabin "Live editorilla" tuotettu, lyhennetty versio kirjoituksesta
(Lopussa mainittu "kaunis dokumentti").

Tiivistelmä:

Simo Kivelän blogissa huhtikuussa 2017 käsiteltiin tämän kevään lyhyen matematiikan ylioppilastehtävää, jossa "Kalle ja Leena tekevät biofysiikan perusteiden harjoitustyötä".  
Tässä tarina.

Yhdyn Simon mielipiteisiin. Ryhdyin kehittelemään diskreetin matematiikan ajattelutapoja, joiden kautta päästäisiin käsittelemään järkevällä tavalla tällaisia minimointi/maksimointi-tehtäviä. Kun nykyisessä signaalinkäsittelyssä jatkuvat funktiot digitoidaan, eli diskretoidaan riittävän hienojakoisesti, on ehkä perusteltua lähestyä koko problematiikkaa diskreeteistä funktioista lähtien. Tällä tavalla saadaan yhtenäisesti käsitellyksi diskreettien ja "sileiden" funktioiden tapaukset ja samantien numeeriset algoritmit ääriarvojen laskemiseen.

Matemaattisia pohjatietoja kirjoituksen ymmärtämiseksi tarvitaan varsin vähän. Lukion lyhyen linjan 1. luokan tiedot riittävät käsittääkseni vallan mainiosti.

Tarkoitus on kirjoittaa ensin tästä, ja jatkaa seuraavaksi käyrän sovittamisesta diskreettiin dataan, mikä on toinen osa yllä mainittua ylioppilastehtävää.

Erinomainen ohjelmisto näihin tarkoituksiin on Matlab, jonka ilmaisella "kloonilla" Octavella voidaan tehdä kaikki tässä esitettävät jutut yhtä hyvin. Editointiin ja julkaisuun tarvittavat työkalut ovat Matlab:ssa kehittyneempiä, mutta Octaven vaatimattomampi käyttöympäristö on loppujen lopuksi pieni haitta verrattuna siihen oivalliseen ajattelutapaan, johon korkean tason numeerinen vektorikieli antaa välineet. Octavella kehitetyt ohjelmat, m-tiedostot, toimivat sellaisenaan suoraan Matlab:ssa.

Kirjoitus on suunnattu myös lukijalle, jolla ei ole Matlab/Octave-kokemusta. Siksi annan yksityiskohtaisia opastustietoiskuja ohjelman käyttöön yksinkertaisin esimerkein, jotka Matlab-osaaja voi sivuuttaa. Ehkä antoisinta on lukea tarinaa avaamalla Matlab/Octave, ja suorittamalla samalla omia kokeiluja. Muistutan, että Octaven voi asentaa ohjelman webbisivulta eri käyttöjärjestelmissä itselleen, Ubuntu-Linux:ssa yksinkertaisimmin suoraan komennolla:
sudo apt-get install octave .

Matematiikan opetus ja ohjelmointi: Valitsemalla riittävän korkean tason kieli, jossa on matemaattisia käsitteitä ja sovelluksia tukevat rakenteet ja funktiot, voitaisiin matematiikan opetukseen tuoda ohjelmointiosuutta, joka tukisi itse matematiikan oppimista ja soveltamista ja samalla opettaisi ohjelmoinnin perusteita. Octave/Matlab voisi olla varteenotettava ehdokas.

Diskreetin ja sileän funktion paikalliset ääriarvot

 

Lähdetään liikkeelle diskreetissä pistejoukossa määritellystä funktiosta. Selvitetään juurta jaksain periaatteet, joilla päästään käsiksi ääriarvojen laskentaan, tarvitsematta mitään pohjatietoja asiasta. Kehitellään Matlab/Octave-työkalut ohjelmien perusominaisuuksia ja ajattelutapoja hyödyntäen, joilla periaatteet saadaan muutetuksi laskenta-algoritmiksi. Kun sileää (jatkuvasti derivoituvaa) funktiota f "digitoidaan", eli lasketaan arvoja rittävän tiheästi diskretoiden, saadaan samalla ohjelmalla lasketuksi halutulla tarkkuudella likiarvoja f:n ääriarvoille.

Kehitetään esimerkki, jossa on paljon sekä piikikkäitä että sileitä ääriarvoja ja sovelletaan kehittämäämme ohjelmaa.  Sen laskemia tuloksia alkuarvauksina käyttäen tarkennetaan tuloksia Matlab:n valmiin funktion fminsearch  tuottamiin, ja koetaan muutama yllätyskin.

Lopuksi valaistaan mahdollisuutta rinnakkaislaskentaan, josta tarkemmin tulevissa blogikirjoituksissa.

Matlab/Octave-oppia


Kirjoitus on suunnattu sellaiselle lukijalle, joka ei ole Matlab:ia aiemmin kohdannut. Siksi pyrin valaisemaan jokaista ohjelman käyttöä Matlabin perusteita opettavin selityksin ja esimerkein.
Matlab-istunnon muoto:
>> format compact   % Komento, jolla turhat välit tiivistetään.
>> 1+1              % Komento, ei sijoitusta muuttujaan.
ans =               % Ohjelman vastaus
     2              % Ohjelman vastaus
>> kaksi=1+1        % Komento, sijoitus muuttujaan
kaksi =             % Vastaus
     2              % Vastaus

Siis komentoa edeltää ohjelman kehote ( >> ) .
Jos seuraat juonta Matlab:n/Octaven äärellä, voit leikkaus/liimaus-tyyliin kopioida rivejä istuntoosi. Tietenkään ohjelman kehotetta ei kopioida. Suositeltavampaa on kopioda kokonaisuus tekstieditoriin ja poistaa kehotteet ja tulosteet. Tiedosto talletetaan vaikkapa nimelle opiaariarvot.m .
Tiedoston komennot voi suorittaa kokonaan tai osissa (kts. oppaista "skriptit", toki myös copy/paste-tyyliin). Työskentelysuosituksia on tuossa lyhyessä oppaassa ja muissa. Kts. myös tämän kirjoituksen viimeistä kappaletta.


Perustapaukset ääriarvoille


Lasketaan  ja piirretään Matlab:lla

>> x=[0 pi/2 4 5 2*pi]
x =
         0    1.5708    4.0000    5.0000    6.2832
>> y=sin(x)
y =
         0    1.0000   -0.7568   -0.9589   -0.0000
>> plot(x,y,'.-.',x,y,'or');grid on

Klikkaa kuva tarkemmaksi/suuremmaksi, jos siltä tuntuu. Samoin seuraavat kuvat.


Miten haetaan min ja max diskreetille funktiolle ?


Funktion kuvaajan muodostavat punaiset 'rinkulapisteet'.
>> [x;y]  % Vektorit x ja y allekkain:
ans =
         0    1.5708    4.0000    5.0000    6.2832
         0    1.0000   -0.7568   -0.9589   -0.0000
 Havainnot:
  •  Välillä $(x_1,x_2)$ funktio kasvaa, eli $y_2 - y_1 > 0$ 
  •  Välillä $(x_2,x_3)$ funktio pienenee, eli $y_3 - y_2 < 0$.
SIIS: $x_2$ on maksimikohta.
  •  Välillä $(x_2,x_3)$ funktio pienenee, samoin välillä $(x_3,x_4)$. SIIS:  $x_3$ ei ole min- eikä maxkohta.
  •  Välillä $(x_3,x_4)$ funktio pienenee ja välillä $(x_4,x_5)$ kasvaa. SIIS:  $x_4$ on minimkohta.
Lasketaan Matlab:lla (Octavella):

Avaimen onneen tarjoaa funktio diff, joka laskee argumenttivektorin peräkkäisten komponenttien erotusvektorin.
Ennenkuin jatketaan esimerkkiämme, katsotaan, miten diff toimii.
>> v=[-1 1 1 3 0]
v =
    -1     1     1     3     0
>> diff(v)
ans =
     2     0     2    -3
Jatketaan juonta, jossa siis
>> x=[0 pi/2 4 5 2*pi];  % Puolipiste estää tulostuksen
>> y=sin(x)      
   y =
         0    1.0000   -0.7568   -0.9589   -0.0000 
>> diff(y)              
ans =
    1.0000   -1.7568   -0.2021    0.9589
Erotukset mittaavat x-vektorin sisäpisteissä $x_2,x_3,x_4$ tapahtuvia muutoksia.
>> [x;y]
ans =
         0    1.5708    4.0000    5.0000    6.2832
         0    1.0000   -0.7568   -0.9589   -0.0000
diff(y) 
ans =
    1.0000   -1.7568   -0.2021    0.9589
Luemme diff(y)-vektorista sen, mikä näkyy y-vektorista (ja kuvasta)
  • Pisteessä $x_2$ erotus muuttuu $+ \rightarrow -$, siis max-piste.
  • Pisteessä $x_3$ $- \rightarrow -$ , siis ei ääriarvoa.
  • Pisteessä $x_4$ $- \rightarrow +$, siis min-piste.
Huomaa, että vektori diff(y)./diff(x) on erotusosamäärävektori, se edustaa siten diskreetin funktion derivaattaa. Kun selvitämme "derivaatan" merkinvaihteluita, ei tarvitse jakaa, koska diff(x) koostuu positiivisista luvuista, useissa tapauksissa vakioaskelista $\Delta x$.
(Miksi yllä jakolaskussa (./ ) on piste ? No tässähän jaetaan kaksi vektoria keskenään. Piste viittaa vastinalkioittain ("pisteittäin") tapahtuvaan jakoon, palataan kohta lähemmin.)
Olemme kiinostuneita erotusten merkeistä, joten kannattaa soveltaa vielä sign-funktiota:
>> sdy=sign(diff(y))
sdy =
     1    -1    -1     1
>> d2y=diff(sdy)  % Erotusvektorin (merkkivektorin) erotukset:
d2y =
    -2     0     2
Tästäpä näkyy suoraan, että $x_2$ on maksimipiste, $x_3$ ei ole ääriarvo, $x_4$ on minimipiste.

Huom: Jos erotusvektorissa diff(y) esiintyy 0, niin kaksi vierekkäistä y-arvoa yhtyy, jolloin "tangentti" on vaakasuora. Tällöin vastaavat y-arvot ovat "epäaitoja" ääriarvoja. Jätämme nämä tapaukset tässä ulkopuolelle, koska ne eivät käytännössä tule esiin sileiden funktioiden, eivätkä piikikkäisten ääriarvojen yhteydessä. Sinänsä niiden mukaan ottaminen ei ole vaikeaa, mutta lienee parempi suojella esitystä liiallisilta rönsyiltä. Jääköön lukijalle (Matlab-opettavaiseksi) harjoitustehtäväksi.
Jotta voisimme kehittää koodin, joka tekee valinnan puolestamme, tarvitsemme vielä yhden tehokkaan Matlab-työkalun, loogisen valinnan.

Indeksointi ja looginen valinta Matlab:ssa


 Kts. Lyhyen Matlab-oppaan opetusta 

Periaate yksinkertaisimmillaan valaistukoon tällä, oppaasta poimitulla esimerkillä: 

>> v=-3:3
v =
    -3    -2    -1     0     1     2     3
>> v>0  % v:n jokaiselle alkiolle tehdään vertailu 0:aan.
ans =
  1×7 logical array
   0   0   0   0   1   1   1   % 1 niille, joille ehto toteutuu.
>> v(v>0)                      % Valitaan ne.
ans =
     1     2     3             % Helppoa, mahtavaa!

Jatketaan Matlab-työtä esimerkkimme parissa


>> maxind=d2y < 0
maxind =
  1×3 logical array
  1   0   0
>> minind=d2y > 0
minind =
  1×3 logical array
   0   0   1  
Sisäpisteiden $x_2,x_3,x_4$ ja vastaavien y-pisteiden muodostamista vektoreista ....
>> xin=x(2:end-1)
xin =
    1.5708    4.0000    5.0000
>> yin=y(2:end-1)
yin =
    1.0000   -0.7568   -0.9589
... valitaan loogisten vektorien maxind ja minind avulla ne, joiden kohdalla on 1 (true).
>> xmax=xin(maxind)
xmax =
    1.5708
>> ymax=yin(maxind)
ymax =
     1
>> xmin=xin(minind)
xmin =
     5
>> ymin=yin(minind)
ymin =
   -0.9589
Kaksi viime vaihetta voidaan yhdistää helposti luettavaksi "idiomiksi" tyyliin:
>> xmin=xin(d2y>0)  % Valitaan ne, joissa toinen differenssi > 0
xmin =
     5
HUOM! Edellä suoritetut operaatiot ovat kaikki vektorioperaatioita. Niin ollen ne toimivat ilman muutoksia miten pitkille x- ja y-vektoreille tahansa, ja antavat kaikki (x,y)-vektoriparin määrittelemän diskreetin funktion max- ja min-pisteet.

Tehtävä: Määritä kaikki min-ja max-pisteet datalle:
>> x=[0 pi/2 4 5 10 15 19 23 8*pi]; % Puolipiste estää tulostuksen.
>> y=sin(x);                        % Hyödyllistä, kun dataa on paljon.
>> plot(x,y,'.-',x,y,'or');grid on;




Ohje: Tee yllä olevan mallin mukaan (kopioimalla koodirivit). Suorita sitten piirto tähän tapaan:
>> clf     % clear figure
>> plot(x,y,'.-.',xmin,ymin,'*r',xmax,ymax,'ob') % 'r'=red, 'b'=blue
>> legend('Datapisteet yhdistettyinä','minimit(red))','maksimit(blue)')

Kirjoitetaan yleispätevä Matlab-funktio (ohjelma)


Pitkän päälle koodirivien kopiointi käy turhauttavaksi. Siksipä kannattaa tässä vaiheessa kirjoittaa homman hoitava funktio, jolle annetaan argumentteina x ja y-vektorit. Funktion tulee palauttaa vektorit xmin, ymin, xmax, ymax

Jos tunnet Matlab:n perusteita oman funktiomäärittelyn verran, voit hypätä seuraavan kohdan yli.

Matlab-funktion määrittely ja käyttö


Tässä muutama johdatusviite aiheeseen:
  1. Syksyllä 2016 pitämäni kurssin luentokalvoja (englanniksi)
  2. Erään Matlab-oppaan ("YAGTOM") suomenkielinen mukaelma funktioaiheesta
  3. Ennen mainitun lyhyen oppaan tekstiä aiheesta
Otetaan tähän näkyville viitteen (1) esimerkkifunktio vielä hiukan yksinkertaistaen. Jospa haluamme laskea vektorin x alkioiden keskiarvon. Matlab-lauseke voidaan kirjoittaa seuraavan esimerkin tavoin:
>> x=1:10;
>> avg=sum(x)/length(x)  % Ei tarvita (./), koska jaetaan skalaareja,  
                            % ... eli lukuja, ei vektoreita.
Haluamme funktion "keskiarvo", jota voitaisiin kutsua millä tahansa vektorilla x tähän tapaan:
>> tulos=keskiarvo(x)
Tekstitiedostoon keskiarvo.m kirjoitetaan koodi, jonka 1. rivillä on avainsana function seuraavaan tapaan:
function y=keskiarvo(x)
% Lasketaan x-arvojen keskiarvo.
% Input: vektori x
% Output : x:n alkioiden keskiarvo
% Kutsuesimerkki: tulos=keskiarvo(1:10)
%
y=sum(x)/length(x);  % Ainoa koodirivi
Nyt voidaan komentaa Matlab-istunnossa vaikka:
>> help keskiarvo  
  Lasketaan x-arvojen keskiarvo.
  Input: vektori x
  Output : x:n alkioiden keskiarvo
  Kutsuesimerkki tulos=keskiarvo(1:10)
 
>> tulos=keskiarvo(1:10)
tulos =
    5.5000
>> keskiarvo(rand(1,100))
ans =
    0.5184
>> keskiarvo(rand(1,1000))
ans =
    0.5153
>> keskiarvo(rand(1,1000))
ans =
    0.5046
Toisinaan on tarvetta funktiolle, jolla on useita argumentteja lähtö- ja/tai tulospuolella. Voisimme haluta funktion, joka palauttaa keskiarvon lisäksi x-vektorin arvoalueen, eli vektorin [min(x),max(x)]. Edellistä muokattaisiin näin:

 function [ka,arvoalue]=ka_minmax(x)
 % Lasketaan x-arvojen keskiarvo ja arvoalue.
 % Input: vektori x
 % Out1 : x:n alkioiden keskiarvo
 % Out2 : [min(x) max(x)]
 % Kutsuesimerkki [ka,range]=keskiarvo(1:10)
 %
 ka=sum(x)/length(x);  
 arvoalue=[min(x),max(x)];
 end


HUOMAA: Funktiotiedosto on aivan tavallinen tekstitiedosto. Sen kirjoittamiseksi ei ole mitenkään pakko käyttää viitteessä (1) neuvottua klikkailua, ei edes Matlabin editoria, vaan mitä tahansa tekstieditoria. Tämä on tärkeää (ja helpottavaa) huomata erityisesti Octave-työskentelyä ajatellen.

Käyttöesimerkki:
>> x=linspace(-pi,4*pi); % Väli [-pi,4*pi] jaetaan 99 osaväliin.
>> y=x.*sin(x); % Pisteittäinen, vastinalkioittainen vektorikertolasku
>> [k,mima]=ka_minmax(y);
>> k
k =
   -0.5928
>> mima
mima =
  -11.0250    7.9160
>> plot(x,y)
>> title('x sin x')
>> hold on
>> plot(x([1 end]),[k k])
>> grid on
>> legend('x sin x','keskiarvo')
>> yrajat=mima+[-0.5 0.5]
yrajat =
  -11.5250    8.4160
>> ylim(yrajat)        % Väljennetään y-aluetta hiukan.

 

Ääriarvojuoni jatkuu



Kokoamme esitetyt palaset, funktiomäärittelysyntaksin ja kehittelemämme koodirivt yhteen. Tuloksena voimme kirjoittaa seuraavanlaisen funktion:

function [xmin,ymin,xmax,ymax] = lok_minmax(x,y)
% Etsitään diskreetillä x-akselilla annettujen y-arvojen paikalliset
% minimit ja maksimit
% Esim 1) 
% x=0:0.1:7; y=sin(5*x); [xmin,ymin,xmax,ymax]=lok_minmax(x,y)
% plot(x,y,xmin,ymin,'x');ylim([-1.1,1.1]);grid on;shg
% Esim 2)
% x=linspace(0,7);y=sin(5*x); [xmin,ymin,xmax,ymax]=lok_min(x,y)
%
% Voidaan siis käyttää myös sileän funktion ääriarvojen numeeriseen
% laskentaan.

dysign=sign(diff(y)); % Differenssivektorin (y(2)-y(1), y(3)-y(2), ...)
                      % ... "merkkivektori"
d2y=diff(dysign);     % 2. differenssi (+-2 tai 0)
xin=x(2:end-1);       % x-sisäpisteet
yin=y(2:end-1);       % vastaavat y-pisteet
xmin=xin(d2y>0);      % Poimitaan sisäpistevektorista xin ne, joissa 
                      % 2. differenssi > 0.;
ymin=yin(d2y>0);      % vastaavat y-pisteet
xmax=xin(d2y<0);      % Poimitaan sisäpistevektorista xin ne, joissa 
                      % 2. differenssi < 0.;
ymax=yin(d2y<0);      % vastaavat y-pisteet
(Jos haluat ottaa käyttöön, niin COPY/PASTE ja sijoita työhakemistoosi tiedostoksi lok_minmax.m ) tai kts. linkki viimeisessä kappaleessa.)

Nyt voidaan keskittyä rakentamiemme työkalujen tarjoamiin nautintoihin.

Lukuisa määrä sileitä ja piikikkäitä ääriarvoja



Kehitellään nopeasti värähtelevä funktio, paikoitellen sileä, toisin paikoin taas piikikäs, vaikkapa näin: $$f(x)=0.1\, x + (\cos x) |\sin(5\,x)|$$
Määritellään $f$ Matlab-funktiokahvaksi ("function handle"), kts. yllä olevat referenssit.

>> f=@(x)0.1*x+cos(x).*abs(sin(5*x))
f =
  function_handle with value:
    @(x)0.1*x+cos(x).*abs(sin(5*x))
Miksi (.*). Kts: Lyhyesti: cos(x) ja abs(sin(5*x)) ovat yhtä pitkiä vektoreita, jotka halutaan kertoa vastinalkioittain (eli "pisteittäin") keskenään. Tulossa 5*x, pistettä ei tarvita, koska 5 on skalaari.
Jatketaan Matlab-istuntoa:
>> vali=[0 10];
x=linspace(0,10,301);
plot(x,f(x));title('Piikkejä ja sileitä ääriarvoja');
grid on;ylim([-0.7 1.7]);shg  % shg:"show graphics"




Jaetaan väli [0,10] 300:aan osaväliin, joten osavälin pituus:
>> h=10/300
h =
    0.0333
(Kuva piirrettiin muutaman kokeilun jälkeen tällä "resoluutiolla", sen laadusta päätellen pisteitä on riittävästi, jotta päästään kaikkiin ääriarvoihin käsiksi.)
>> x=linspace(0,10,301);  % Kerrataan x:n määrittely.
>> y=f(x);                % Puolipiste estää tulostuksen.
x ja y ovat nyt 301:n alkion vektoreita, jotka yhdessä edustavat diskreettiä funktiota, jonka ääriarvoja haemme.
>> [xmin,ymin,xmax,ymax]=lok_minmax(x,y);
>> length(xmin)
ans =
    17
>> length(xmax)
ans =
    17
Sekä minimejä että maksimeja on 17 kpl.
Jatketaan piirtämistä edelliseen kuvaan: Piirretään minimipisteisiin punainen tähti ja maksimipisteisiin sininen rinkula.
>> hold on  % Jatka piirtoa edelliseen grafiikkaruutuun.
>> plot(xmin,ymin,'*r',xmax,ymax,'ob')
>> ylim([-0.7 1.7])  % Väljennetään y-aluetta.


Kuvan perusteella näyttää siltä, että olemme löytäneet approksimaatiot kaikille ääriarvoille. Approksimaatioiden absoluuttinen tarkkuuskin on tiedossamme: |virhe| on korkeintaan h= 0.0333
Katsotaan nyt numeerisia arvojakin, ehkä vaikka 7 ensimmäistä:
>> [xmin(1:7);ymin(1:7)]
ans =
    0.6355    1.2709    1.7391    2.2408    2.8428    3.4448    4.0134
    0.0922    0.1481    0.0623   -0.3835   -0.6686   -0.6085   -0.2023
>> format long    % Täysi näyttötarkkuus (=laskentatarkkuus)
>> [xmin(1:3);ymin(1:3)]  % Nyt vain 3 arvoa näytiksi.
ans =
   0.635451505016722   1.270903010033445   1.739130434782609
   0.092242289698521   0.148144547093171   0.062294731898389
>> format       % Takaisin oletusnäyttötarkkuuteen.


Onko kaikki löydetty, miten parannetaan tarkkuutta?


Helpoin tapa olisi hienontaa jakoa. Ensimmäiseen kysymyksen riittäisi kaksinkertaistaa kerran pari, tyyliin
>> N=300;
>> %% Toistetaan seuraavia rivejä:
>> N=2*N; x=linspace(0,10,N);y=f(x);   
>> [xmin,ymin,xmax,ymax]=lok_minmax(x,y);
>> minlkm=length(xmin);   % Muuttuuko edellisestä?
>> maxlkm=length(xmax);   %  ------- " ----------
>> minlkm                 % No katsotaan. Muuttujan arvon saa näkyviin
                          % kirjoittamalla sen nimen komentokehotteeseen.
minlkm =
    17                     
>> maxlkm                 
maxlkm =
    17
>> N                      % Tarkistetaan vielä N:n arvo
N =
   600

Pistemäärän kaksinkertaistus ei muuttanut lukumäärää. On siis hyvä syy uskoa, että kaikki ääriarvot on mainitulla tarkkuudella löydetty.
Tässä virheraja puolittui. Tarkkuutta saataisiin yksi dekadi lisää 10-kertaistamalla jakopisteiden lukumäärä (N=10*N).

Tässä voisi olla eräs luonteva lopetuskohta. Mutta luontevaa on myös jatkaa hiukan.

Tuloksien tarkentaminen kehittyneillä hakumenetelmillä


Numeerisia optimointimenetelmiä, samoin kuin nollakohdan hakumenetelmiä on kehitetty erityisesti yhden paikallisen ääriarvon tai nollakohdan laskentaan. Tällaisia ovat Fibonacci- ja kultaisen leikkauksen haku, sekä Newtonin menetelmä (ääriarvon tapauksessa derivaatan) nollakohdan etsinnässä. Yleensä nämä menetelmät keskittyvät yhden minimin hakuun, joka sijaitsee "riittävän lähellä" annettua alkuarvoa $x_0$
Matlab:ssa on valmis funktio fminsearch, joka toimii myös usean muuttujan funktioille. Se ei käytä derivaattaa, joten voidaan soveltaa myös piikikkäisiin tapauksiin. Yksinkertaisimmillaan funktion kutsu on
 >> xmin = fminsearch(Fun,x0)

Jos halutaan hakea kaikki ääriarvot jollain alueella, pitää ensin löytää sopivat alkuarvaukset. No nythän meillä on ohjelma juuri siihen tarkoitukseen. Voimme jatkaa seuraavanlaiselĺa koodilla, joka on luontevinta suorittaa skriptitiedostosta, siksi jätimme kehotteet pois.
 N=length(xmin);
 tarkennetutminimit=zeros(1,N);
 for k=1:N
   tarkennetutminimit(k)=fminsearch(f,xmin(k));
 end
 xminT=tarkennetutminimit;
 yminT=f(xminT);
Piirretään:
plot(xminT,yminT,'ok',xmax,ymax,'ob');ylim([-0.7 1.7]);
title('Rengastetut pisteet löytyvät, 3 jää saamatta')


Tutkitaanpa asiaa

Sijoitetaan allekkain kiinnostavat komponentit:
>> [xmin([7 8 9 13 14 15]);tarkennetutminimit([7 8 9 13 14 15])]
ans =
    4.0134    4.5151    5.0167    7.5251    8.0268    8.5284
    4.0296    3.4393    5.0265    9.1190    8.5254    8.5254
Ylärivillä on siis omalla ohjelmallamme saadut pisteet, joita käytetään alkuarvoina fminsearch-funktiolle. Kuten kuvasta näkyy, 3 minimiä jää saavuttamatta tässä "tarkennnuksessa". Kun katsotaan komponentteja 7 ja 8, huomataan, että alkuarvauksen kasvattaminen johtaa tuloksen pienenemiseen. Myös komponentti 13 (arvo 7.5251) suppenee "väärään" minimiin. Huomataan siis, että tarkentamalla saadaan oikeita arvoja, mutta osa suppenee samoihin arvoihin ja siten tässä tapauksessa 3 arvoa jää saavuttamatta.

Katsotaan vielä kokonaisuudessaan:
>>[xmin;tarkennetutminimit]
ans =
  Columns 1 through 7
    0.6355    1.2709    1.7391    2.2408    2.8428    3.4448    4.0134
    0.6283    1.2566    1.7324    2.2422    2.8359    3.4393    4.0296
  Columns 8 through 14
    4.5151    5.0167    5.6522    6.2876    6.9231    7.5251    8.0268
    3.4393    5.0265    5.6549    6.2832    6.9115    9.1190    8.5254
  Columns 15 through 17
    8.5284    9.1304    9.7324
    8.5254    9.1191    9.7225
>> diff(ans,1)    % Erotukset rivi-indeksin suunnassa
ans =
  Columns 1 through 7
   -0.0071   -0.0143   -0.0067    0.0014   -0.0069   -0.0055    0.0163
  Columns 8 through 14
   -1.0757    0.0098    0.0027   -0.0045   -0.0116    1.5939    0.4986
  Columns 15 through 17
   -0.0031   -0.0114   -0.0099
Virheraja näyttää pitävän paikkansa, paitsi komponenttien 8, 9, 13, 14, 15 kohdalla. MUTTA: Niistä lähtien fminsearch suppeneekin kohti "väärää" minimiä. Meidän menetelmämme on luotettavampi, joskin epätarkempi. Parannusta saadaan suorittamalla tarkennettu käsittely niille alkuarvoille, jotka fminsearch ajaa virherajan ulkopuolelle.
Jos tarkkuusvaatimus ei ole kovin suuri, voidaan tyytyä oman ohjelmamme käyttöön, jossa saadaan tarkkuutta lisää joko jakamalla koko väli tiheämmin (jokaista dekadia kohti dekadi lisää, aika tehotonta, tosin Matlab:ssa viuhuu lujaa). Parempi tapa: Jaetaan kukin osaväli (xmin(k)-h,xmin(k)+h) esim. 100:aan osaan.

Harjoitustehtävä


Tee vastaava selvitys maksimipisteille. Huomaa, että fminsearch sopii funktion $f$ maksimien etsintään tyyliin >> xmax=-fminsearch(@(x)-f(x),x0) , koska $\max f(x) = - min(-f(x))$

Alkuarvauksen valinnasta

Käsittelemme alkuarvon valintaproblematiikkaa tarkemmin seuraavissa blogikirjoituksissa funktion nollakohdan määritystehtävän ja Newtonin menetelmän yhteydessä. Senverran voidaan paljastaa, että pieni muutos alkuarvauksessa voi johtaa suppenemiseen eri ratkaisua kohti. Jos tutkitaan niiden pisteiden joukkoa, joista alkaen suppeneminen tapahtuu johonkin tiettyyn ratkaisuun, johdutaan usein fraktaaliseen rakenteeseen. Tämä näkyy erityisen havainnollisesti tutkittaessa esim. polynomifunktioiden nollakohtia kompleksitasossa.
Pääsääntönä toki on, että yleensä nopeasti suppenevat menetelmät toimivat, kun valitaan alkuarvo "riittävän läheltä" ratkaisua, joten jaon riittävä hienontaminen johtaa yleensä hyvään tulokseen.

Rinnakkaislaskentaa, parallell toobox


Viimeisen vuosikymmenen aikana on algoritmien suunnittelussa yleistynyt enenevässä määrin rinnakkaisuuden hyödyntäminen. Tämä tarkoittaa sitä, että ohjelma pyritään jakamaan toisistaan riippumattomiin osiin, jotka voidaan suorittaa samanaikaisesti. Fyysisesti sen tekee mahdolliseksi tietokonearkkitehtuuri, jossa ison muistin ja yhden tehokkaan prosessorin sijasta on jaettu muisti ja monta prosessoria. Tällaisia suuria "klustereita" on mm. Aalto-yliopiston Tritonus klusteri ja CSC-superlaskentakeskuksen Taito,
joissa kummassakin on Matlab ja seuraavassa puheena oleva parallel toolbox käytettävissä.
Matlabiin on siis saatavissa rinnakkaislaskentaan erikoistunut "parallell toolbox", jolla voidaan suorittaa tämänkaltainen rinnakkaistus hyvin helposti. Harjoittelu kotikoneella on vaivatonta, monessa tapauksessa koodia ei tarvitse muuttaa lainkaan siirryttäessä isompaan klusteriin. Kotiläppärin 2-ydinprosessori antaa vain 2 rinnakkaisprosessia, mutta se mahdollistaa alkutestaamisen, vaikkei vielä sanottavaa tehonlisäystä tuokaan.
Edellä esitetty algoritmi soveltuu erinomaisesti rinnakkaislaskentaan, kunkin alkuarvon laskemisen jälkeinen tarkentaminen jakautuu luonnollisesti riippumattomiin osiin.
Palaan aiheeseen seuraavissa blogikirjoituksissa. Viittaan tässä kurssiin, jonka pidin syksyllä 2016 yhdessä Juha Kuortin kanssa: Matlab course 2, "advanced", parallel

Optimointia ja nollakohtien hakua usemman muuttujan funktioille:

Tehokkuuden tarve lisääntyy huomattavasti, kun muuttujia tulee lisää, jolloin rinnakkaislaskennan edut moninkertaistuvat. Näistäkin aiheista lisää jatkossa.

Matlab/Octave-tiedostot, Matlab:n julkaisuvälineiden tuotokset


Viimeksi mainittu pdf-dokumentti on hiukan lyhennettynä sama, joka tässä on kirjoitettu Blogger-välineillä. Matlabin Live editor on kunnianhimoinen ja monin tavoin hieno työväline, jolla myös matemattiset kaavat voidaan editoida suoraan Latex-muodossa, kuvat tulevat automaattisesti ym. Omat pikku hankaluutensakin siinä toki on, ja suositeltavinta on monestakin syystä rakentaa Matlab-dokumentti yllä olevan esimerkin mukaisesti strukturoiduksi m-skriptiksi, ja ajaa (tai olla ajamatta) julkaisuksi. Käytön tietyistä epämukavuuksista kertovat omat korkeat versionumerot tiedostojen nimissä. Toisaalta tuloksena saatava dokumentti on kyllä vieläkin kauniimpi, kuin Bloggerilla tuotettu, vai mitä?.

Ei kommentteja:

Lähetä kommentti