Pienimmän neliösumman menetelmä
Palataan hetkeksi ratkaisemaan yhtälöryhmää
Ax=b, missä A on matriisi ja
b tunnettu vektori. Insinöörien elämässä tulee usein
vastaan ongelmia, missä tällä yhtälöryhmällä ei yksinkertaisesti ole
ratkaisua. Yksi yleinen syy on se, että reaalimaailman mittaukset eivät
ole virheettömiä ja tämän vuoksi tuottavat ristiriitoja käytössä olevaan
malliin. Tällöin yhtälöryhmälle Ax=b voidaan
tarkan ratkaisun sijaan etsiä ”paras mahdollinen” ratkaisu, jota
kutsutaan pienimmän neliösumman ratkaisuksi.
Jos yhtälöryhmällä A\mathbf{x}=\mathbf{b} ei ole ratkaisua, niin
\mathbf{b}\not\in \mathcal{R}(A). Pienimmän neliösumman ratkaisu
on sellainen vektori \overline{\mathbf{x}}, että vektori
A\overline{\mathbf{x}}\in \mathcal{R}(A) on ikään kuin lähimpänä
vektoria \mathbf{b}. Nyt kun periaate on selvä, tarvitsee enää
tietää, miten \overline{\mathbf{x}} lasketaan.
Merkitään virhevektoria, eli residuaalia
\mathbf{e}= \mathbf{b}- A\mathbf{x}, jolloin komponenteittain
\begin{split}\mathbf{e}=
\begin{bmatrix}
e_1 \\ e_2 \\ \vdots \\ e_m
\end{bmatrix} =
\begin{bmatrix}
b_1 - (a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n) \\
b_2 - (a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n) \\
\vdots \\
b_m - (a_{m1}x_1 + a_{m2}x_2 + \cdots + a_{mn}x_n)
\end{bmatrix} = \mathbf{b}- A\mathbf{x}.\end{split}
Pienimmän neliösumman ratkaisun ajatuksena on minimoida normi
\|\mathbf{e}\|. Samanaikaisesti minimoidaan myös normin neliö
\|\mathbf{e}\|^2 = e_1^2 + e_2^2 + \cdots + e_m^2,
sillä normi on ei-negatiivinen reaaliluku. Tässä summassa jokainen termi
on ei-negatiivinen, ja siksi se saavuttaa minimiarvonsa täsmälleen
silloin, kun jokainen yhteenlaskettavista termeistä e_i^2,
i = 1, 2, \ldots, m minimoituu.
Pidetään summaa
\|\mathbf{e}\|^2 = e_1^2 + e_2^2 + \cdots + e_m^2 vuorotellen
muuttujavektorin \mathbf{x} komponenttien x_j,
j = 1, 2, \ldots, n funktiona. Nähdään, että
\frac{\mathrm{d}e_i}{\mathrm{d}x_j} = -a_{ij},
jolloin ketjusäännön nojalla
\frac{\mathrm{d}\|\mathbf{e}\|^2}{\mathrm{d}x_j} = 2e_1\frac{\mathrm{d}e_1}{\mathrm{d}x_j} + 2e_2\frac{\mathrm{d}e_2}{\mathrm{d}x_j} + \cdots + 2e_m\frac{\mathrm{d}e_m}{\mathrm{d}x_j} = -2(a_{1j}e_1 + a_{2j}e_2 + \cdots + a_{mj}e_m).
Matriisitulon avulla kirjoitettuna
\begin{split}\frac{\mathrm{d}\|\mathbf{e}\|^2}{\mathrm{d}x_j} = -2
\begin{bmatrix}
a_{1j} & a_{2j} & \cdots & a_{mj}
\end{bmatrix}
\begin{bmatrix}
e_1 \\ e_2 \\ \vdots \\ e_m
\end{bmatrix} = -2\mathbf{a}_j^T\mathbf{e},\end{split}
missä \mathbf{a}_j on matriisin A sarake j. Kun
kaikki tällaiset derivaatat kerätään pystyvektoriin, saadaan
\begin{split}\begin{bmatrix}
\frac{\mathrm{d}\|\mathbf{e}\|^2}{\mathrm{d}x_1} \\ \frac{\mathrm{d}\|\mathbf{e}\|^2}{\mathrm{d}x_2} \\ \vdots \\ \frac{\mathrm{d}\|\mathbf{e}\|^2}{\mathrm{d}x_n}
\end{bmatrix} =
\begin{bmatrix}
-2\mathbf{a}_1^T\mathbf{e}\\ -2\mathbf{a}_2^T\mathbf{e}\\ \vdots \\ -2\mathbf{a}_n^T\mathbf{e}
\end{bmatrix} = -2
\begin{bmatrix}
\mathbf{a}_1^T \\ \mathbf{a}_2^T \\ \vdots \\ \mathbf{a}_n^T
\end{bmatrix}\mathbf{e}= -2A^T\mathbf{e}= -2A^T(\mathbf{b}- A\mathbf{x}).\end{split}
Asiaa sen tarkemmin todistamatta väitetään, että virhevektorin normin
neliö \|\mathbf{e}\|^2 minimoituu silloin, kun kaikki derivaatat
ovat yhtä aikaa nollia, eli kun
-2A^T(\mathbf{b}- A\mathbf{x}) = \mathbf{0}. Tämä ehto on
yhtäpitävä yhtälön
A^TA\mathbf{x}= A^T\mathbf{b}
kanssa.
Pienimmän neliösumman menetelmää voidaan lähestyä myös teoreettisemmin
seuraavasta näkökulmasta. Jos matriisiyhtälössä
A\mathbf{x}= \mathbf{b} matriisi A ei ole kääntyvä,
voidaan silti yrittää löytää mahdollisimman hyvin käänteismatriisin
käyttäytymistä imitoiva matriisi. Tällaista matriisia kutsutaan
pseudoinverssiksi. Seuraava tulos otetaan käyttöön todistamatta.
Lause.
Jos A on m \times n-matriisi, niin on
olemassa yksikäsitteinen n \times m-matriisi A^+, joka
toteuttaa ehdot
- AA^+A = A,
- A^+AA^+ = A^+,
- (AA^+)^T = AA^+
- (A^+A)^T = A^+A.
Matriisia A^+ sanotaan matriisin A Moore-Penrosen
pseudoinverssiksi.
Lemma.
Jos A on m \times n-matriisi, niin
\mathcal{R}(A^TA) = \mathcal{R}(A^T).
Oletetaan, että \mathbf{y}\in \mathcal{R}(A^TA).
Tällöin \mathbf{y}= A^TA\mathbf{x} jollekin avaruuden
\mathbb R^n vektorille \mathbf{x} ja
A\mathbf{x}= \mathbf{z} on avaruuden \mathbb R^m
vektori. Täten \mathbf{y}= A^T\mathbf{z} jollekin avaruuden
\mathbb R^m vektorille \mathbf{z}, eli
\mathbf{y}\in \mathcal{R}(A^T).
Oletetaan sitten, että \mathbf{y}\in \mathcal{R}(A^T), eli
\mathbf{y}= A^T\mathbf{x} jollekin avaruuden \mathbb R^m
vektorille \mathbf{x}. Sijoitetaan matriisin A paikalle
Moore-Penrosen pseudoinverssin A^+ sisältävä matriisitulo
edellisen lauseen ensimmäisen ehdon mukaisesti ja sovelletaan muitakin
pseudoinverssin ominaisuuksia.
\mathbf{y}= (AA^+A)^T\mathbf{x}= A^T(AA^+)^T\mathbf{x}= A^T(AA^+)\mathbf{x}= (A^TA)A^+\mathbf{x}
Tässä A^+\mathbf{x}= \mathbf{z} on avaruuden \mathbb R^n
vektori, jolle \mathbf{y}= A^TA\mathbf{z}, ja täten
\mathbf{y}\in \mathcal{R}(A^TA).
Edellä olevista päättelyistä seuraa, että
\mathcal{R}(A^TA) \subseteq \mathcal{R}(A^T) ja
\mathcal{R}(A^T) \subseteq \mathcal{R}(A^TA), joten
\mathcal{R}(A^TA) = \mathcal{R}(A^T). \square
Nyt voidaan lopulta osoittaa, että yhtälöryhmän
A\mathbf{x}= \mathbf{b} normaaliryhmällä
A^TA\mathbf{x}= A^T\mathbf{b} on aina vähintään yksi ratkaisu.
Lause.
Olkoon A m\times n-matriisi ja
\mathbf{b}\in\mathbb R^m. Tällöin yhtälöryhmän
A\mathbf{x}=\mathbf{b} normaaliryhmällä on ratkaisu
\mathbf{x}= \overline{\mathbf{x}}, ja tämä on yhtälöryhmän
A\mathbf{x}= \mathbf{b} pienimmän neliösumman ratkaisu.
Vektori A^T\mathbf{b} kuuluu luonnollisesti
matriisin A^T sarakeavaruuteen \mathcal{R}(A^T).
Edellisen lemman nojalla siis
A^T\mathbf{b}\in \mathcal{R}(A^TA), ja täten normaaliryhmällä
A^TA\mathbf{x}= A^T\mathbf{b} on ratkaisu
\overline{\mathbf{x}}. Osion alussa johdetun mukaisesti
sijoituksella \mathbf{x}= \overline{\mathbf{x}} minimoidaan
virhevektorin \mathbf{b}- A\mathbf{x} normi, joten
\|\mathbf{b}- A\overline{\mathbf{x}}\| \leq \|\mathbf{b}- A\mathbf{x}\|
aina, kun \mathbf{x}\in \mathbb R^n. \square
Pienimmän neliösumman ratkaisu on siis aina olemassa, ja edellinen lause
tarjoaa keinon löytää sen.
Huomautus.
Pienimmän neliösumman ratkaisun olemassaolo ei tarkoita
sitä, että se olisi yksikäsitteinen. On hyvinkin mahdollista, että
matriisi normaaliryhmän kerroinmatriisi A^TA ei ole kääntyvä, ja
tällöin normaaliryhmällä on ääretön määrä ratkaisuja.
Varmistetaan vielä, että pienimmän neliösumman ratkaisu toimii kuten
yhtälön A\mathbf{x}= \mathbf{b} tarkkakin ratkaisu ainakin
seuraavan lauseen mielessä.
Lause.
Olkoon \overline{\mathbf{x}} yhtälöryhmän
A\mathbf{x}= \mathbf{b} pienimmän neliösumman ratkaisu. Jos
yhtälöryhmällä A\mathbf{x}=\mathbf{b} on ratkaisuja, niin
tällöin myös \overline{\mathbf{x}} on sen ratkaisu.
Koska yhtälöryhmällä on ratkaisuja, löydetään vektori
\mathbf{x}, jolle \|\mathbf{b}- A\mathbf{x}\| = 0. Koska
pienimmän neliösumman ratkaisu minimoi yhtälön vasemman puolen
lausekkeen, on oltava
\|\mathbf{b}- A\overline{\mathbf{x}}\| = 0. Ainoastaan
nollavektorin normi on nolla, joten
\mathbf{b}- A\overline{\mathbf{x}} = \mathbf{0} ja edelleen
A\overline{\mathbf{x}} = \mathbf{b}. \square
Normaaliryhmä A^TA\mathbf{x}=A^T\mathbf{b} voidaan
ratkaista Gaussin
eliminoinnilla, kuten mikä tahansa muukin yhtälöryhmä. Erityisen
kiintoisa tapaus syntyy silloin, kun ratkaisu on yksikäsitteinen, eli
kerroinmatriisi A^TA on kääntyvä. Tällöin pienimmän neliösumman
ratkaisuksi tulee
\overline{\mathbf{x}}=(A^TA)^{-1}A^T\mathbf{b}.
Moore-Penrosen pseudoinverssi sattuu tässä tapauksessa olemaan
täsmälleen se matriisi, joka ”kääntää” matriisin A: jos
A^TA on kääntyvä, niin A^+ = (A^TA)^{-1}A^T ja pienimmän
neliösumman ratkaisu on \overline{\mathbf{x}} = A^+\mathbf{b}.
Lause.
Olkoon A m \times n-matriisi. Tällöin
seuraavat väitteet ovat voimassa.
- Jos A^TA on kääntyvä, niin A^+A = I_n.
- Jos A on kääntyvä, niin A^+ = A^{-1}.
Jätetään harjoitustehtäväksi. \square
Käänteismatriisia (A^TA)^{-1} kutsutaan yhtälöryhmän
A\mathbf{x}= \mathbf{b} normalisoivaksi kertoimeksi, ja sen
olemassaololle löydetään riittävä ja välttämätön ehto.
Lause.
Normalisoiva kerroin (A^TA)^{-1} on olemassa
täsmälleen silloin, kun matriisin A sarakkeet ovat lineaarisesti
riippumattomat.
Tiedetään, että m \times n-matriisin A
sarakkeet ovat lineaarisesti riippumattomat täsmälleen silloin, kun
\operatorname{rank}(A) = n. Vastaavasti
n \times n-matriisi A^TA on kääntyvä jos ja vain jos
\operatorname{rank}(A^TA)=n. Mutta aikaisemman lemman nojalla
\mathcal{R}(A^TA) = \mathcal{R}(A^T), joten
\operatorname{rank}(A^TA) = \operatorname{rank}(A^T) = \operatorname{rank}(A) = n,
ja väite on todistettu. \square
Esimerkki.
Tutki yhtälöryhmän
\begin{split}\begin{cases}
x+5y = 3\\
2x-2y = 2\\
-x+y =5
\end{cases}\end{split}
ratkaisujen olemassaoloa. Jos ratkaisuja ei ole, etsi residuaalin normin
minimoiva approksimaatio ryhmän ratkaisulle. Jos tämä arvio on
yksikäsitteinen, etsi yhtälöryhmän kerroinmatriisille Moore-Penrosen
pseudoinverssi.
Viedään yhtälöryhmän kokonaismatriisi redusoituun
riviporrasmuotoon.
\begin{split}[A\mid\mathbf{b}]=
\begin{bmatrix}
1 & 5& 3\\
2 & -2 & 2\\
-1 & 1 &5
\end{bmatrix}\longrightarrow \operatorname{rref}[A\mid\mathbf{b}]=
\begin{bmatrix}
1 & 0& 0\\
0 & 1 & 0\\
0 & 0 &1
\end{bmatrix}\end{split}
Tästä nähdään, että yhtälöryhmällä ei ole tarkkaa ratkaisua, joten
etsitään pienimmän neliösumman ratkaisu. Normaaliryhmän kerroinmatriisi
ja vakiovektori
\begin{split}A^TA=
\begin{bmatrix}
6 & 0\\
0 & 30
\end{bmatrix} \qquad\text{ja}\qquad A^T\mathbf{b}=
\begin{bmatrix}
2\\16
\end{bmatrix}.\end{split}
Matriisi A^TA on selvästi kääntyvä, sillä se on kolmiomatriisi,
jolla ei ole nolla-alkiota diagonaalillaan. Täten normaaliryhmällä
A^TA\mathbf{x}=A^T\mathbf{b} on yksikäsitteinen pienimmän
neliösumman ratkaisu
\begin{split}\overline{\mathbf{x}}=(A^TA)^{-1}A^T\mathbf{b}=\frac{1}{180}
\begin{bmatrix}
30 & 0 \\ 0 & 6
\end{bmatrix}
\begin{bmatrix}
2 \\ 16
\end{bmatrix}=
\begin{bmatrix}
\frac{1}{3}\\\frac{8}{15}
\end{bmatrix}.\end{split}
Kerroinmatriisin A pseudoinverssiksi saadaan tällöin
\begin{split}A^+=(A^TA)^{-1}A^T=\frac{1}{180}
\begin{bmatrix}
30 & 0\\
0 & 6
\end{bmatrix}
\begin{bmatrix}
1 & 2 & -1\\
5 & -2 & 1
\end{bmatrix}=\begin{bmatrix}
\frac{1}{6} & \frac{1}{3} & -\frac{1}{6}\\
\frac{1}{6} & -\frac{1}{15} & \frac{1}{30}
\end{bmatrix}.\end{split}