http://math.aalto.fi/~apiola/matlab/opas/lyhyt/LA.html
Päivitetty 5.4.2010 HA
matriisijako
(\),det,inv,cond,rref |
Takakeno-operaattori, eli vasemmanpuoleinen matriisijako (mldivide, "matrix left divide")
suorittaa yhtälöryhmän ratkaisun tapauksessa, jossa tehtävällä on yksikäsitteinen ratkaisu. Tällöin matriisin
A on käytännössä oltava neliömatriisi, jonka determinantti ei ole 0.
Takakeno toimii neliömatriisin lisäksi "pystymatriisille" eli
ylimääräytyvän yhtälösysteemin tapauksessa (enemmän yhtälöitä kuin tuntemattomia).
Tällöin x = A \ b määrittää vektorin x, joka ratkaisee yhtälön
A x = b "pienimmän neliösumman (PNS, LSQ) mielessä. Aiheesta lisää
tässä ja tässä ,
(Mutta nämä vielä puuttuvat
[Tue May 31 16:31:29 EEST 2011])
Yhtälöryhmä
voidaan esittää matriisimuodossa Ax=b. Muodostetaan kerroinmatriisi A ja oikean puolen vektori b:
10x1-7x2 = 7 -3x1+2x2+6x3 = 4 5x1-x2+5x3 = 6
» A=[10,-7,0;-3,2,6;5,-1,5]
A =
10 -7 0
-3 2 6
5 -1 5
» b=[7,4,6]'
b =
7
4
6
Yhtälöryhmä ratkeaa takakenolla
» x=A\b
x =
0
-1
1
ja ratkaisun voi tarkistaa matriisikertolaskulla:
» A*x
ans =
7
4
6
Ennen ratkaisua voi laskea determinantin ja "häiriöalttiuden".
>> det(A)
ans =
-155 % det ei 0, ==> yksikäs, ratk.
>> cond(A)
ans =
9.4070 % cond on luokkaa 10 ==> hyvänlaatuinen ongelma.
Determinanttiehto on teoriassa oikea, mutta käytännön laskut antavat
yleensä jotain, joka on lähellä nollaa. Determinantti riippuu tehtävän skaalauksesta; determinantti voi olla lähellä nollaa ja tehtävä silti
"hyvänlaatuinen" ("well-conditioned").
Toisaalta det voi olla kaukana nollasta ja tehtävä silti "häiriöaltis"
("ill-conditioned"). Näistä lisää jatkossa, ainakin viitteissä.
>> B=[1 2 3;4 5 6;7 8 9] % Matriisi, jonka ihminen aina
B = % ensimmäiseksi Matlabille tarjoaa
1 2 3 % sattuu olemaan singulaarinen.
4 5 6
7 8 9
>> x=B\b % Mitä sanoo "takakeno"?
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.541976e-018.
x = % rcond on estimaatti alla selitetyn cond-luvun
1.0e+016 * % käänteisluvulle. rcond on laskentatyön ja tarkkuuden
-2.2518 % kannalta edullisempi. Kts. help rcond
4.5036
-2.2518
>> cond(B) % cond - "condition number", "häiriöluku". Nyrkkisääntö:
ans = % Luvun 10-kantainen logaritmi (tässä luokkaa 16)
3.8131e+016 % ilmaisee, miten monta numeroa suhteellisesta tarkkuudesta
% häviää, tässä käytännöllisesti katsoen kaikki.
>> det(B)
ans =
0
x =
1.0e+016 *
-2.2518
4.5036
-2.2518
>> B*x
ans =
8
0
0
>> b
b = % Yhtälö ei ole lähelläkään toteutumista.
7
4
6
Yllä matriisi B on singulaarinen. Lineaaristen yhtälösysteemien teoria sanoo,
että tällaisella systeemillä on joko ääretön määrä ratkaisuja tai ei ollenkaan.
Gaussin eliminaatio johtaa muotoon, josta näkyy, kumpi on kyseessä. Gaussin
algoritmi on helppo ohjelmoida Matlabiin, se on valmiiksi toteutettu funktiossa rref ("reduced row echelon form", redusoitu riviporrasmuoto).
Jatketaan esimerkkiä:
>> [B b] % Liitetään b-vektori matriisin B perään.
ans =
1 2 3 7
4 5 6 4
7 8 9 6
>> rref(ans)
ans =
1 0 -1 0
0 1 2 0
0 0 0 1 % Ristiriitainen yhtälö (0 = 1) ==> ei ratkaisua.
Ratkaisujen olemassaolo riippuu siitä, kuuluuko b-vektori matriisin B
sarakeavaruuteen. Otetaan (tässä vaikka tuurilla) vektori, joka koostuu
ykkösistä:
>> e=ones(3,1)
e =
1
1
1
>> rref([B e])
ans =
1 0 -1 -1
0 1 2 1
0 0 0 0 % yhtälö 0 = 0 toteutuu, yksi vapaasti valittava
% muuttuja (x3). Siis ääretön määrä ratkaisuja.
>> B\e
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.541976e-018.
ans =
-2.5000
4.0000
-1.5000
>> B*ans
ans =
1
1
1
Matlab päättelee aivan oikein lähes singulaarisuuden kerroinmatriisista.
Sen antama ratkaisu on kuitenkin oikea, tosin äärettömästä ratkaisuparvesta
se poimii yhden.
Miten tuo onnekas valinta muuttuu hallituksi?
Koska matriisi B on singulaarinen, sen rangi (aste) on korkeintaan 2.
Matlabilla:
>> rank(B)
ans =
2
Matriisin sarakeavaruus on siis 2-ulotteinen. Jotta yhtälösysteemillä olisi
ratkaisuja, on oikean puolen vektori (e) voitava esittää sarakevektorien
lineaarikombinaationa. Tämä tarkoittaa, että jos matriisiin B liitetään
sarakevektori e, niin rangi ei kasva, eli tässä tapauksessa on = 2.
No, katsotaan:
>> rank([B e])
ans =
2
No niinpä on. rref([B e]) . Lähinnä nyt palautimme
mieleen/esittelimme
lineaarialgebran perusasioita Matlab-avusteisesti. )
Matlabin takakenon laatu on alansa huippua. RCOND-arvio lähes singulaarisessa tapauksessa antaa hyvän pohjan ratkaisun luotettavuuden ja tarkkuuden arvioinnille.
>> inv(A)
ans =
-0.1032 -0.2258 0.2710
-0.2903 -0.3226 0.3871
0.0452 0.1613 0.0065
>> A*ans
ans =
1.0000 0 0.0000
-0.0000 1.0000 0.0000
-0.0000 -0.0000 1.0000
>> format long
>> ans
ans =
1.000000000000000 0 0.000000000000000
-0.000000000000000 1.000000000000000 0.000000000000000
-0.000000000000000 -0.000000000000000 1.000000000000000
% Tulos on yksikkömatriisi, miinus-merkit johtuvat pyöristysvirheistä,
% jotka tässä ovat minimaalisia.
>> inv(B)
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.541976e-018.
ans =
1.0e+016 *
-0.450359962737050 0.900719925474099 -0.450359962737049
0.900719925474100 -1.801439850948198 0.900719925474099
-0.450359962737050 0.900719925474099 -0.450359962737050
>> B*ans
ans =
2 0 2 % Käänteismatriisia ei ole, joten
8 0 0 % tulos on mitä sattuu.
16 0 8
Huom! Käänteismatriisia ei yleensä lasketa. Tämä johtuu siitä, että käänteismatriisin muodostaminen ja sillä kertominen on laskentatyön kannalta raskaampaa ja myös numeerisesti epätarkempaa kuin systeemin ratkaisu suoraan Gaussin algoritmilla, siinäkin tapauksessa, että ratkaistavana on useita yhtälöitä, joilla on sama kerroinmatriisi. (Käänteismatriisikin toki lasketaan Gaussin algoritmilla.)
Esim: Käänteismatriisilla kertomalla ja suoralla "matriisijakolaskulla" operointien eroa niin suoritusaikojen kuin laskentatarkkuuden kannalta valaistaan tässä Mathworksin inv -viitteessä esimerkillä.
Enemmän yhtälöitä kuin tuntemattomia ==> PNS-ratkaisu lu, chol, qr, svd, ...
[ov, oa]=eig(A);
help sparse, doc sparse