[vastaus aiempaan viestiin]
| Kirjoittaja: | Seppo Mustonen |
|---|---|
| Sähköposti: | - |
| Päiväys: | 22.12.2003 14:04 |
Simo Puntanen kirjoitti 15.12.2003 12:58 :
> Se vaan on niin, että vaikka Survon suhteen "...sitä saa mitä
> tilaa...", niin silti siinä on jonkinmoinen kynnys ehdottaa niitä
> muutoksia, kun tietää että niissä helposti voi olla Sepolle kova homma
> toteuttaa, eikä yksinkertaisesti "kehtaa" (kai se paras sana on?)
> ottaa juttuja esille.
Onhan se jotenkin mukavaa, että joku sentään säälii näin meikäläistä,
mutta mielellään sitä aina yrittää tehdä jotain.
Simon ehdotukset Survon matriisitulkin lisätoiminnoiksi kiertyvät
etenkin yleistetyn käänteismatriisin (erityisesti ns. Moore-Penrose-
inverssin) ympärille.
Kuten hän toteaa, kyseinen käänteismatriisi on helppo muodostaa
matriisin singulaariarvohajotelman avulla, mutta uskaltaa nyt toivoa,
että se saataisiin Survossa suoraankin muotoa
MAT G=MPINV(A)
olevan komennon avulla.
Jotta nekin näitä keskusteluja seuraavat, jotka eivät riittävästi
tunne matriisilaskentaa, saisivat jonkinlaisen käsityksen siitä,
mistä on kysymys, yritän hiukan selittää taustaa esimerkein.
Tarkastellaan 3 yhtälön ja kolmen tuntemattoman (x1,x2,x3) lineaarista
yhtälöryhmää
2*x1+5*x2-3*x3=3
(1) 5*x1-8*x2+5*x3=4
3*x1+2*x2+7*x3=28
Käsin laskettaessa yhtälöryhmä ratkaistaan ns. eliminointimenetelmällä
ratkaisemalla ensin vaikka viimeinen yhtälö x3:n suhteen eli
x3=(28-3*x1-2*x2)/7 ja sijoittamalla tämä x3:n lauseke kahteen
ensimmäiseen yhtälöön. x3 on siis "eliminoitu" ja jäljellä on 2 yhtälön
ja kahden tuntemattoman (x1,x2) lineaarinen yhtälöryhmä. Soveltamalla
tähän vastaavaa menettelyä x2:n osalta päästään yhden tuntemattoman (x1)
lineaariseen yhtälöön, jonka ratkaisun jälkeen saadaan aikaisemmista
sijoituskaavoista lopullinen tulos eli tässä tapauksessa x1=1, x2=2 ja
x3=(28-3*x1-2*x2)/7=3. Ratkaisun voi tarkastaa sijoittamalla x-arvot
yhtälöiden vasemmanpuoleisiin osiin ja toteamalla, että ko. lausekkeet
ovat samanarvoisia oikealla puolella olevien vakioiden (3,4,28) kanssa.
"Matriisikielellä" sama yhtälöryhmä kirjoitetaan lyhyesti muodossa
AX=C
missä
| 2 5 -3 | |x1| | 3|
A= | 5 -8 5 | X= |x2| C= | 4|
| 3 2 7 | |x3| |28|
eli A on yhtälöryhmän kertoimien muodostama 3x3-matriisi, C on oikean
puolen vakioiden muodostama vektori (3x1-matriisi) ja X tuntemattomien
(x1,x2,x3) vektori.
Muodollisesti yhtälöryhmä ratkaistaan melkein kuten suureet A,X ja C
olisivat lukuja eli
-1
X = A *C (C kerrotaan siis A:n käänteisarvolla)
A:n käänteisarvo vain tässä yleisemmässä tilanteessa tarkoittaa
A:n käänteismatriisia, jota Survon matriisikielessä merkitään INV(A).
Survon matriisitulkilla lineaarisia yhtälöryhmiä on mahdollista
ratkaista tehokkaamminkin (esim. MAT SOLVE), mutta käänteismatriisin
avulla ratkaisu menee seuraavasti:
MATRIX A ///
2 5 -3
5 -8 5
3 2 7
MATRIX C ///
3
4
28
MAT SAVE A / Matriisin A talletus
MAT SAVE C / Matriisin C talletus
MAT X=INV(A)*C / Yhtälöryhmän ratkaisu
MAT LOAD X / Ratkaisuvektorin tulostus
MATRIX X
INV(A)*C
/// 1
1 1.000000
2 2.000000
3 3.000000
INV(A) vastaa A:n "käänteisarvoa" siinäkin mielessä, että esim.
MAT K=A*INV(A)
MAT LOAD K,12.123456789012345,CUR+1
MATRIX K
A*INV(A)
/// 1 2 3
1 1.000000000000000 -0.000000000000000 0.000000000000000
2 0.000000000000000 1.000000000000000 0.000000000000000
3 0.000000000000000 -0.000000000000000 1.000000000000000
eli K=I (yksikkömatriisi) vastaten ykköstä reaaliluvuilla.
Kaikilla matriiseilla ei ole kuitenkaan ole käänteismatriisia.
Jos tarkastellaan yhtälöryhmän (1) asemasta seuraavaa
2*x1+5*x2-3*x3=3
(2) 5*x1-8*x2+5*x3=4
5*x1-8*x2+5*x3=4
jossa kolmas yhtälö on korvattu toisella, joka siis esiintyy kahdesti,
on selvää, ettei tällä yhtälöryhmällä voi olla yksikäsitteistä
ratkaisua, koska kolmas yhtälö ei anna mitään lisäinformaatiota.
Tämä näkyy teknisesti tarkastelemalla kerroinmatriisin anatomiaa
ns. singulaariarvohajotelman (SVD) avulla
MATRIX A2 /// Tässä on uusi kerroinmatriisi
2 5 -3
5 -8 5
5 -8 5
MATRIX C2 /// ja tässä uusi vakiovektori.
3
4
4
MAT SAVE A2
MAT SAVE C2
MAT SVD OF A2 TO U,D,V / Singulaariarvohajotelman muodostaminen
MAT LOAD D / Singulaariarvojen tulostus
MATRIX D
Dsvd(A2)
/// sing.val
svd1 15.72724
svd2 4.31903
svd3 0.00000
Matriisille A2 on tehty SVD eli se esitetään kolmen matriisin tulona
A2=U*D*V'
missä U ja V ovat ns. ortogonaalisia matriiseja ja D on lävistäjä-
matriisi, jonka pelkkä lävistäjäosa on matriisitiedostossa D.
Edellä V' tarkoittaa, että V on transponoitu, eli rivit ja sarakkeet
vaihtavat paikkaa. Ortogonaalisuus merkitsee sitä, että U*U'=V*V'=I
eli ortogonaalisen matriisin käänteismatriisi saadaan yksinkertaisesti
transponoimalla.
Singulaariarvohajotelman muodostus (MAT SVD) on esim. matriisin kääntä-
miseen verrattuna aika raskas toimenpide (ei tosin tunnu missään
nykyisillä koneilla kohtuudimensioilla), mutta se lienee turvallisin
tapa tutkia matriisin rakennetta ja avain mitä moninaisimpiin
sovelluksiin.
Se, että ainakin yksi singulaariarvoista (D:n lävistäjäalkioista)
on 0, merkitsee säännöllisyyden menetystä eikä matriisilla voi silloin
olla käänteismatriisia; matriisi on vajaa-asteinen.
Muodollisesti matriisin A2 käänteismatriisi voitaisiin esittää SVD:n
kautta muodossa (Huom. tulon A*B*C inverssi on INV(C)*INV(B)*INV(A))
(3) INV(A2)=V*INV(D)*U'
ja tässä INV(D) saadaan yksinkertaisesti korvaamalla lävistäjäalkiot
käänteisarvoillaan. Näin ei kuitenkaan voi tehdä, jos yksikin on
nolla, kuten on asian laita äskeisellä A2:lla.
Nyt seuraa ratkaiseva askel eli yleistetyn käänteismatriisin
muodostaminen.
Koska nollia ei voi "kääntää", käännetään vain nollasta poikkeavat
alkiot, minkä Survossa voi tehdä näin: (arvot alle 10^(-14) katsotaan
nolliksi.)
MAT DI=DINV(DV(D),1e-14)
MAT LOAD DI
MATRIX DI
DINV(DV(Dsvd(A2)))
/// svd1 svd2 svd3
svd1 0.063584 0.000000 0.000000
svd2 0.000000 0.231533 0.000000
svd3 0.000000 0.000000 0.000000
ja sovelletaan kaavaa (3) käyttäen INV:in sijasta DINV:iä eli
A2:n yleistetty käänteismatriisi on silloin
MAT A2_INV=V*DI*U'
MAT LOAD A2_INV
MATRIX A2_INV
Vsvd(A2)*DINV(DV(Dsvd(A2)))*Usvd(A2)'
/// 1 2 3
1 0.19636 0.06068 0.06068
2 0.09103 -0.01712 -0.01712
3 -0.05072 0.01192 0.01192
Simon ensimmäistä toivomusta noudattaen olen nyt toteuttanut äskeisen
menettelyn suorana komentona
MAT G=MPINV(A2) / *G~MPINV(A2) 3*3
MAT LOAD G
MATRIX G
MPINV(A2)
/// 1 2 3
1 0.19636 0.06068 0.06068
2 0.09103 -0.01712 -0.01712
3 -0.05072 0.01192 0.01192
eli G täsmää äsken vaiheittain laskettuun A_INV-matriisiin. Laskentatapa
on käytännössä sama eli mennään SVD:n kautta. MAT MPINV on sikäli vielä
yleisempi, että sillä voi "kääntää mitä tahansa" eli mm. muitakin kuin
neliömatriiseja.
Katsotaan saatua yhtälöryhmän ratkaisua:
MAT X2=G*C2 / *X2~MPINV(B)*C2 3*1
MAT LOAD X2
MATRIX X2
MPINV(B)*C2
/// 1
1 1.07456
2 0.13611
3 -0.05678
Ratkaisu ei ole sama kuin yhtälöryhmälle (1) saatu (1,2,3), mutta
kuitenkin kelvollinen eli yhtälöt (2) toteutuvat mikä nähdään näin:
MAT E=A2*X2-C2
MAT LOAD E,12.12345678901234,CUR+1
MATRIX E
A2*MPINV(B)*C2-C2
/// 1
1 0.00000000000000
2 -0.00000000000000
3 -0.00000000000000
Miksi ei siis saatu sitä "alkuperäistä" ratkaisua?
Yhtälöryhmän (2) ratkaisu ei ole yksikäsitteinen, sillä eliminointi-
menettelyn ensimmäinen vaihe johtaa tässä tapauksessa yleiseen
ratkaisuun
(4) x1=(44-t)/41, x2=(7+25t)/41, x3=t,
missä t voi olla mikä tahansa reaaliluku. Jos valitaan t=3, saadaan
taas alkuperäinen ratkaisu (1,2,3), mutta mitä onkaan nyt saadun
ratkaisun (likim.) X1=1.07456, X2=0.13611, X3=-0.05678 salaisuus?
Se on kaikista mahdollisista ratkaisuista (4) "geometrisesti lyhin"
eli silloin ratkaisuvektorin alkioiden neliösumma on pienin mahdollinen.
Tämä on eräänlaista matematiikan ekonomiaa parhaimmillaan, mitä siis
SVD (ja yleistetty käänteismatriisi) tarjoaa.
Vastaavasti, kun tarkastellaan yhtälöryhmää
2*x1+5*x2-3*x3=3
5*x1-8*x2+5*x3=4
5*x1-8*x2+5*x3=28
jossa kerroinmatriisi on (1):stä ja vakiovektori (2):sta, syntyy risti-
riitainen tilanne: 2 viimeistä yhtälöä eivät voi olla voimassa saman-
aikaisesti.
Yritetään kuitenkin ratkaista se muodollisesti yleistetyn käänteismat-
riisin avulla:
MAT X3=MPINV(A2)*C
MAT LOAD X3
MATRIX X3
MPINV(A2)*C
/// 1
1 2.53099
2 -0.27482
3 0.22930
MAT E=A2*X3-C
MAT LOAD E
MATRIX E
A2*MPINV(A2)*C-C
/// 1
1 0.0000
2 12.0000
3 -12.0000
Nähdään, että vain ensimmäinen yhtälö toteutuu, mutta jälleen
metematiikan ekonomia toteutuu, sillä nyt "virhevektori" E=(0,12,-12)
on tälle ratkaisulle pienimmän neliösumman mielessä kaikkein lyhin.
Tämän voi päätellä mm. näin: Kahden viimeisen yhtälön vakiotermejä
koskeva ristiriita on kooltaan 28-4=24. Saatu ratkaisu yksinkertaisesti
puolittaa tämän riidan, jolloin E-vektorin alkioiden neliösumma on
pienin mahdollinen.
Puheena ollut matematiikan ekonomia on tässä yleisesti yhteydessä
pienimmän neliösumman keinoon lineaarisessa regressioanalyysissa.
Jos X on selittävien muuttujien matriisi ja y selitettävän muuttujan
vektori, niin yleisesti regressiokertoimet b saadaan suoraan lausekkeena
b=MPINV(X)*y
ja huomattakoon, että tässä X ei enää ole neliömatriisi, vaan siinä
on rivejä yleensä huomattavasti enemmän kuin sarakkeita ja lisäksi
X:n sarakkeet voivat olla lineaarisesti toisistaan riippuvia.
Viimeksi mainitusta vajavaisuudesta seuraa tiettyjä ongelmia mutta
ei laskuteknisiä.
Jos X on täysiasteinen, eli sen kaikki singulaariarvot ovat
positiivisia, MPINV(X) on sama kuin INV(X'X)*X' eli regressiokertoimille
saadaan oppikirjoista tuttu matriisilauseke.
Vajaa-asteisissa tilanteissa tuo tuttu kaava "räjähtää", mutta
MAT MPINV toimii.
Oma "pikku" juttunsa on se, mitä tarkoitetaan "nollalla" singulaari-
arvojen yhteydessä, koska "nollia" ei käännetä. Vaikka teoreettisesti
jollain matriisilla olisi singulaariarvoina tarkkoja nollia, koneen
laskiessa kaksoistarkkudella (noin 16 merkitsevän numeron) nollat
saattavat loksahtaa hieman paikoiltaan.
Näin on viisasta tulkita nolliksi aina tietyn rajan alittavat arvot.
MAT MPINV:issä noudatetaan seuraavaa sääntöä. Jos matriisin A suurin
singulaariarvo on s, nolliksi katsotaan sellaiset singulaariarvot,
jotka ovat pienempiä kuin 10^(-15)*s.
Tätä rajaa käyttäjä voi säädellä antamalla komennon muodossa
MAT G=MPINV(A,eps)
jolloin po. raja on eps*s.
Koska viestini on jo tässä vaiheessa ylipitkä, palaan vielä
Simon muihin toivomuksiin myöhemmin erikseen.
Joka tapauksessa kaikille (jotka ovat jaksaneet tunnollisesti kahlata
läpi tämän litanian) Hyvää Joulua!
-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!