Sinter MATLAB opas, lineaarialgebraa
http://math.tkk.fi/~apiola/matlab/opas/lyhyt/LA.html
Päivitetty 5.4.2010 HA


[
Edellinen] [Seuraava] [Sisällys]

Lineaarialgebraa/matriisilaskentaa

Lineaarinen yhtälöryhmä

Viitteitä
  1. Mathworks-doku: "Systems of Linear Equations"
  2. Mathworks-doku: "Overdetermined Systems, Least squares"

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])

Perusesimerkki, yksikäsitteinen ratkaisu

Yhtälöryhmä

10x1-7x2 = 7
-3x1+2x2+6x3= 4
5x1-x2+5x3 = 6
voidaan esittää matriisimuodossa Ax=b. Muodostetaan kerroinmatriisi A ja oikean puolen vektori b:
» 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ä.

Singulaarinen/häiriöaltis tapaus

Katsotaanpa, miten Matlab reagoi singulaariseen neliömatriisiin.
>> 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.
( Tämän kirjoittaminen Matlabille ei ollut vähemmän työlästä kuin 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.

Käänteismatriisi inv


>> 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

Kyseessä on siis "kapea ja korkea" kerroinmatriisi, eli enemmän rivejä kuin sarakkeita. Yksi Matlabin hienous on, että tällöin muodostetaan ns. pienimmän neliösumman ratkaisu, jota käsitellään ohimennen kohdassa interpolaatio ja käyrän sovitus sekä perusteellisemmin kohdassa ...

Matriisihajotelmia

lu, chol, qr, svd, ...

Ominaisarvot

[ov, oa]=eig(A);

Harvat ("sparse") matriisit

help sparse, doc sparse
[Edellinen] [Seuraava] [Alkusivu]