Matriisitulkin laajennuksia

[viesti Survo-keskustelupalstalla (2001-2013)]

Kirjoittaja: Seppo Mustonen
Sähköposti:    -
Päiväys: 19.3.2003 8:49

Survon matriisitulkki on tarkoitettu tukemaan paitsi tilastotieteen
ja matriisilaskennan opetusta myös tutkimustyötä erityisesti
lineaaristen mallien ja monimuuttujamenetelmien alueella.
Koska useat Survon tilastolliset operaatiot tallettavat tuloksensa
(myös) matriisitiedostoihin, on aina mahdollista jatkaa laskentaa
matriisitulkin avulla ja saada vastauksia erityiskysymyksiin tai
jalostaa tuloksia käyttäjän itsensä haluamalla tavalla.
Useita tilastollisia menetelmiä voi myös "ohjelmoida" sukroina
Survon tilastollisten ja matriisioperaatioiden yhdistelminä.

Survon matriisitulkki on siinä mielessä ollut edelläkävijä, että
matriisit voidaan varustaa selväkielisillä rivi- ja sarakeotsikoilla,
jotka siirtyvät ja mukautuvat loogisella tavalla matriisioperaatioissa.
Samoin matriisin otsikkotiedoista näkyy automaattisesti sen
syntyhistoria "sisäisenä nimenä", joka on sisältöä vastaava
matriisilauseke.
Liikkuminen matriistiedostojen, Survon datatiedostojen ja toimitus-
kentän välillä on joustavaa eli "taulukkojen" esitystavan voi valita
aina tilanteen mukaan.

Periaatteellisella tasolla en näe matriisitulkin kohdalla mitään
uudistustarpeita, mutta aina on mahdollista harkita, mitä toimintoja
siihen tulisi lisätä.
Olen toteuttanut allamainitut laajennukset ja keskusteltuani
mm. Simo Puntasen kanssa.
Nämä tulevat mukaan versiosta 1.33 lähtien.
Olisi hyvä kuulla pikaisesti mahdollisista aihetta koskevista
lisätoiveista, ennenkuin lisäykset "pannaan pakettiin".

Uusia matriisioperaatioita:
(Muutamat ovat olleet saatavilla jo aikaisemmin toisessa muodossa.)
MAT C=NULL(A)   / C on A:n nolla-avaruuden ortonormaalinen kanta.
                  Sama tulos X=C saadaan myös ratkaisemalla yhtälö AX=0
                  komennolla MAT SOLVE X FROM A*X=0.

MAT C=BASIS(A)  / C on A:n sarakeavaruuden ortonormaalinen kanta.
                  Tästä aiheesta lisää vuoden 1999 artikkelistani
                  www.helsinki.fi/survo/matrix99.html .

MAT C=RANK(A)   / C(1,1)=A:n aste (C on a 1x1-matriisi)

Kaikki kolme edellä mainittua operaatiota perustavat ratkaisunsa
matriisin A singulaariarvohajotelmaan A=UDV', koska singulaari-
arvot D kertovat luotettavimmin matriisin asteesta ja U:n ensimmäiset
sarakkeet antavat sarakeavaruuden kannan sekä vastaavasti V:n viimeiset
sarakkeet nolla-avaruuden kannan.

MAT C=TRACE(A)  / C on tr(A).
MAT C=DET(A)    / C on A:n determinantti.
MAT C=LDET(A)   / C on A:n determinantin logaritmi.

MAT C=#INTSCAL(A) skaalaa matriisin A sarakkeet kokonaisluvuiksi
ketjumurtolukutekniikalla. Jos A on numeerinen vakio tai 1x1-matriisi,
muodostetaan astenevasti A:n ketjumurtolukuapproksimaatioita ja
niihin liittyviä muita suureita, jotka talletetaan matriisiksi C.

Tuo viimeinen (#INTSCAL) on hieman eksoottisempaa laatua ja
kaivannee selityksiä.

Yritetään ratkaista yhtälöryhmä
Y1: 1*X1+2*X2+0*X3+4*X4 = 5
Y2: 3*X1+1*X2+3*X3+7*X4 = 4   eli matriisimuodossa esim. A*X=B
Y3: 5*X1+5*X2+2*X3+0*X4 = 3
Y4: 7*X1+2*X2+1*X3+3*X4 = 0

ja saada numeerisen ratkaisun ohella myös "tarkka tulos".
Numeerinen ratkaisu syntyy matriisitulkilla mm. seuraavasti:
Muodostetaan yhtälölöiden kertoimista ja vakiotermeistä
matriisi YHT

MATRIX YHT

///  X1  X2  X3  X4   C
Y1    1   2   0   4   5
Y2    3   1   3   7   4
Y3    5   5   2   0   3
Y4    7   2   1   2   0

ja talletetaan se matriisitiedostoksi YHT.MAT:
MAT SAVE YHT

Erotetaan kerroinmatriisi A ja vakiovektori C:
MAT A=YHT(*,1:4)
MAT B=YHT(*,5)

Ratkaistaan yhtälö A*X=B:
MAT SOLVE X FROM A*X=B

ja tulostetaan ratkaisuvektori X=(X1,X2,X3,X4):
MAT LOAD X,12.123456789012345,CUR+1
MATRIX X
(X_from_A*X=B)
///                       C
X1       -0.529113924050633
X2        1.255696202531646
X3       -0.316455696202533
X4        0.754430379746836

Koska kaikki yhtälöryhmässä esiintyneet luvut ovat rationaalilukuja
(ja tässä tapauksessa jopa kokonaislukuja), ratkaisuvektori koostuu
rationaaliluvuista ja yllä ovat niiden 15-desimaaliset likiarvot.

Näitä likiarvoja vastaavat rationaaliluvut (eli "tarkat arvot")
voi "arvata" yksitellen esim. muunnoskomennolla
X2        1.255696202531646(15:ratio)=496/395 (4.4408920985006e-016)
eli X2:n tarkka arvo olisi 496/395, koska suluissa oleva virhe on
laskutarkkuuden rajoissa 0 (tai käyttämällä INTREL-komentoa).

#INTSCAL-komennolla tämä on tehtävissä kertaheitolla näin:
MAT C=#INTSCAL(X)

Matriisissa (vektorissa) C on nyt X skaalattuna kokonaisluvuiksi:
MAT LOAD C
MATRIX C
#INTSCAL((X_from_A*X=B))
///             C
X1           -209
X2            496
X3           -125
X4            298

Kun lasketaan vielä A:n determinantti
MAT D=DET(A)
MAT_D(1,1)=395

saadaan "tarkka" ratkaisu muodossa C/det(A) (tässä tapauksesssa).

Taustaa:
#INTSCAL etsii pystyriveittäin matriisin alkioille niitä lähinnä
vastaavat rationaaliluvut (kokonaislukujen suhteet) kehittämällä
likiarvot ketjumurtoluvuiksi (continued fraction).
Esim. luku  30/17=1.7647058823529411... ketjumurtolukuna on

              1
 30/17 = 1 + -----------
                  1
             1 + -------
                      1
                 3 + ---
                      4

On tapana merkitä ylläolevaa kehitelmää yksinkertaisesti näin:
30/17 = [1,1,3,4]
Luvuista 1,1,3,4 käytetään nimitystä "partial quotient", suomeksi
ehkä "osanimittäjä" (E.Lindelöf).

Rationaaliluvuilla (kuten 30/17) ketjumurtolukuesitys on päättyvä
(vaikka desimaaliesitys ei sitä aina ole). Irrationaalilukujen
kejumurtolukuesitykset eivät ole päättyviä mutta neliöjuurilausekkeilla
sikäli säännöllisiä, että ne ovat jaksollisia jostain vaiheesta lähtien.
Esim. luvun 3 neliöjuuri on ketjumurtolukuna
sqrt(3) = [1,1,2,1,2,1,2,1,2,...].
ja
(sqrt(5)-1)/2=[0,1,1,1,1,1,1,...]
(eli kultaisen leikkauksen suhde, joka on ketjumurtolukuna kaikkein
hitaimmin suppeneva!),
(5*sqrt(7)-10)/3 = [1,13,8,1,2,1,8,13,8,1,2,1,8,...] .

Kun matriisin sarakkeen alkioille on löydetty rationaaliesitykset,
#INTSCAL etsii niiden nimittäjien suurimman yhteisen jaettavan ja
tekee ne sitten "samannimisiksi", jolloin päästään haluttuun
kokonaislukuesitykseen erikseen kunkin sarakkeen osalta.

#INTSCAL:ia voi käyttää myös ketjumurtolukuesitysten muodostamiseen ja
tarkasteluun yksittäisillä luvuilla seuraavasti:

Halutaan löytää luvulle pi (likiarvo 3.141592653589793) hyviä
likimääräisiä rationaaliarvoja.

MAT C=#INTSCAL(3.141592653589793,10) / likiarvot termiä 10 myöten
LOADM C,(C11),CUR+1
3.1415926535897931_as_continued_fractions
            quotient           p           q         p/q       p/q-x
fract0             3           3           1 3.000000000 -0.14159265
fract1             7          22           7 3.142857143  0.00126449
fract2            15         333         106 3.141509434 -0.00008322
fract3             1         355         113 3.141592920  0.00000027
fract4           292      103993       33102 3.141592653 -0.00000000
fract5             1      104348       33215 3.141592654  0.00000000
fract6             1      208341       66317 3.141592653 -0.00000000
fract7             1      312689       99532 3.141592654  0.00000000
fract8             2      833719      265381 3.141592654 -0.00000000
fract9             1     1146408      364913 3.141592654  0.00000000
fract10            3     4272943     1360120 3.141592654 -0.00000000

Saadaan siis ketjumurtolukuesitykseksi
pi = [3,7,15,1,292,1,1,1,2,1,3,...] (ei jaksollinen)
ja askeleittain pi:n likiarvot (p/q) 3, 22/7, 333/106, 355/113, jne.,
jotka kaikki ovat tunnettuja matematiikan historiasta.

Kääntäen, matriisioperaatio #FRAC_TO_DEC muuntaa vektorin (matriisin
ensimmäisen sarakkeen), jonka oletetaan sisältävän ketjumurtoluku-
esityksen, desimaaliluvuksi. Siis äskeisen tilanteen jälkeen saadaan
                                              ACCURACY=16
MAT D=#FRAC_TO_DEC(C)
MAT_D(1,1)=3.141592653589389

- Seppo

Vastaukset:

Survo-keskustelupalstan (2001-2013) viestit arkistoitiin aika ajoin sukrolla, joka automaattisesti rakensi viesteistä (yli 1600 kpl) HTML-muotoisen sivukokonaisuuden. Vuoden 2013 alusta Survo-keskustelua on jatkettu entistäkin aktiivisemmin osoitteessa forum.survo.fi. Tervetuloa mukaan!

Etusivu  |  Keskustelu
Copyright © Survo Systems 2001-2013. All rights reserved.
Updated 2013-06-15.