Pienimmän neliösumman menetelmä¶

Palataan hetkeksi ratkaisemaan yhtälöryhmää \(A\mathbf{x}= \mathbf{b}\), missä \(A\) on matriisi ja \(\mathbf{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 \(A\mathbf{x}=\mathbf{b}\) voidaan tarkan ratkaisun sijaan etsiä ”paras mahdollinen” ratkaisu, jota kutsutaan pienimmän neliösumman ratkaisuksi.

Määritelmä.

Olkoon \(A\) jokin \(m\times n\)-matriisi ja \(\mathbf{b}\) avaruuden \(\mathbb R^m\) vektori. Tällöin yhtälöryhmän \(A\mathbf{x}=\mathbf{b}\) pienimmän neliösumman ratkaisu (least squares solution) on avaruuden \(\mathbb R^n\) vektori \(\overline{\mathbf{x}}\), joka toteuttaa epäyhtälön

\[\|\mathbf{b}-A\overline{\mathbf{x}}\|\le \|\mathbf{b}-A\mathbf{x}\|\]

aina, kun \(\mathbf{x}\in\mathbb R^n\).

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.

Määritelmä.

Yhtälöryhmän \(A\mathbf{x}= \mathbf{b}\) normaaliryhmä on \(A^TA\mathbf{x}= A^T\mathbf{b}\).

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)\).

Todistus.

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.

Todistus.

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.

Todistus.

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}\).

Todistus.

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.

Todistus.

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.

Ratkaisu.

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}\]