Metoda nejmenších čtverců - připomenutí

V článku Metoda nejmenších čtverců: řešení rovnic, která nemají řešení jsme si představili metodu nejmenších čtverců, která se používá pro řešení následujícího problému: Máme matici $X \in \mathbb{R}^{n,p}$ a vektor $Y \in \mathbb{R}^{n}$, kde $n$ je větší nebo rovno $p$, a snažíme se nějak poprat se soustavou $$ X\beta = Y, $$ kde $\beta \in \mathbb{R}^{p}$ je vektor neznámých. Nesnažíme se jí vyřešit, neboť tušíme, že asi nemá řešení: soustava má více rovnic ($n$) než neznámých ($p$).

Namísto přesného řešení se snažíme najít $\beta$ takové, aby Euklidovská vzdálenost vektorů $X\beta$ a $Y$ byla co nejmenší. Jinými slovy, hledáme vlastně minimum následující funkce $p$ proměnných (složek vektoru $\beta$) $$ ||X\beta - Y||. $$ Takové minimum najít umíme za pomoci metod analogických metodám pro hledání extrému funkce jedné proměnné známých z prvácké analýzy. Výsledkem je vektor $$ \beta = (X^TX)^{-1}X^TY, $$ kterému se říká odhad metodou nejmenší čtverců.

Výpočet prostředky prvácké lineární algebry

Jestli jste v kurzu lineární algebra (BI-LIN) dávali pozor, tak si teď říkáte, co o tom jako bude psát, dyť je hotovo: Matici $X$ známe, transpozice a násobení je hračka, takže spočítáme $X^TX$ a pomocí Gaussovy eliminační metody (GEM) spočítám inverzi. Přesněji řečeno, upravíme pomocí GEM matici $$ \left( X^TX \mid E \right) $$ na matici $$ \left( E \mid B \right), $$ kde $E$ je jednotková matice z $\mathbb{R}^{p,p}$ a $B$ je, jak víme z BI-LIN, hledaná inverze k $X^TX$, tj. $B = (X^TX)^{-1}$.

Problém je, že takto to na počítači počítat nelze. Důvodem jsou velmi špatné vlastnosti GEM z pohledu numerické stability: jak si můžete přečíst v článku Nástrahy sčítání pomocí počítače, i triviální operace sčítání a odčítání může dělat počítači velké problémy. Velkým strašákem je tzv. *katastrofické krácení*, které nastává, když odečítáme dvě blízká čísla. Když si rozmyslíte operace prováděné v GEM, tak k takovému odčítání může snadno dojít. Navíc když v GEM uděláme takovou drobnou chybu, tak se tato chyba projeví i v počítání s dalšími řádky a sloupci a rozleze se nám tak po celé výsledné matici. Zkrátka GEM má z pohledu počítání na počítači velice nevhodné vlastnosti.

QR rozklad matice

Abychom si mohli ukázat, jak se tedy výpočet odhadu metodou nejmenších čtverců skutečně počítá, musíme si představit jeden pojem, který z BI-LIN neznáte a tím jsou ortogonální matice.

Ortogonální matice a její vlastnosti

Ze střední školy byste meli vědět, že pro (standardní) skalární součin dvou vektorů $u = (u_1\ u_2\ u_3)$ a $v = (v_1\ v_2\ v_3)$ platí $$ u \cdot v = u_1v_1 + u_2v_2 + u_3v_3 = ||u||\,||v||\,\cos \varphi, $$ kde $\varphi$ je úhel mezi těmito dvěma vektory a $||u|| = \sqrt{u_1^2 + u_2^2 + u_3^2}$ je Euklidovská norma vektoru. Díky této vlastnosti snadno poznáme, že dva vektory jsou na sebe kolmé, tj. $\varphi = \pi/2$. Protože $\cos(\pi/2) = 0$, musí být skalární součin takových dvou vektorů roven nule.

Tato vlastnost dovoluje rozšířit pojem kolmosti vektorů i do vyšších dimenzí, kde bychom jinak s představou toho, co je kolmé, dost bojovali. Abychom se právě od naší tří a dvojrozměrné představy odpoutali, říkáme místo "vektory jsou na sebe kolmé", že vektory jsou ortogonální.

Ortogonální matici lze definovat jako matici, jejíž sloupce mají jednotkovou velikost a jsou navzájem ortogonální. Když si uvědomíme, že skalární součit dvou (sloupcových) vektorů $u,v \in \mathbb{R}^{n}$ lze pomocí maticového násobení zapsat jako $$ u \cdot v = u_1v_1 + u_2v_2 + \cdots + u_nv_n = u^Tv, $$ můžeme objevit jednu zajímavou vlastnost ortogonálních matic.

Předpokládejme, že $Q \in \mathbb{R}^{n,n}$ je ortogonální matice. Zkusme spočítat součin $Q^TQ$. Maticové násobení je definováno tak, že v $i$tém řádku a $j$tém sloupci matice $Q^TQ$ bude vlastně skalární součin $i$tého a $j$tého sloupce matice $Q$. Jelikož jsou ale sloupce $Q$ navzájem ortogonální, bude pro $i \neq j$ tento skalární součit vždy nula. A v případě, že $i = j$, bude výsledkem jednička, neboť sloupce mají jednotkovou velikost a pro libovolný vektor $u \in \mathbb{R}^{n}$ a jeho velikost platí $$ ||u||^2 = u^Tu, $$ neboli vynásobíme-li skalární vektor sám se sebou, vyjde nám kvadrát jeho velikosti. Ta je ale pro sloupce ortogonální matice jednotková! Celkově můžeme říci, že výsledkem součinu je jednotková matice, tj. že platí \begin{equation} \label{eq:OG_vlastnost} Q^TQ = E. \end{equation}

Tuto vlastnost ortogonálních matic lze chápat i jako jejich definici. Platí totiž, že matice, které mají navzájem ortogonální sloupce velikosti jedna, jsou právě ty matice, které rovnost $\ref{eq:OG_vlastnost}$ splňují. Rovnost $\ref{eq:OG_vlastnost}$ má navíc zajímavé důsledky. První z nich je ten, že pro ortogonální matice je snadné spočítat inverzi, což je jinak numericky poměrně obtížný výpočet: Pro ortogonální matici platí, že její inverze je rovna její transpozici!

Další dvě důležité vlastnosti plynou z následujícího faktu platného pro součin dvou "vhodněrozměrných" matic $A$ a $B$ a jejich transpozici: $$ (AB)^T = B^TA^T, $$ tedy že *transpozice součinu je součin transpozic*. Zkusme si vzít dvě ortogonální matice $Q_1$ a $Q_2$ z $\mathbb{R}^{n,n}$. Potom platí toto: $$ (Q_1Q_2)^T Q_1Q_2 = (Q_2^TQ_1^T) Q_1Q_2 = Q_2^T( Q_1^TQ_1) Q_2 = Q_2^TEQ_2 = Q_2Q_2 = E, $$ tedy $Q_1Q_2$ je také ortogonální matice. Když to trošku zobecníme, můžeme říci, že *součin libovolného konečného počtu ortogonálních matic je opět ortogonální matice*.

A poslední vlastnost, kterou si ukážeme, je tato: pro libovolný vektor $x \in \mathbb{R}^{n}$ platí $$ ||Qx||^2 = (Qx)^TQx = (x^TQ^T)Qx = x^T(Q^TQ)x = x^TEx = x^Tx = ||x||^2. $$ Tedy velikosti vektorů $Qx$ a $x$ jsou stejné, neboli *násobení vektoru ortogonální maticí nemění velikost tohoto vektoru*!

QR rozklad matice

Rozkladem matice se obvykle myslí to, že zadanou matici napíšeme jako součin dvou a více matic, které mají nějaké speciální vlastnosti. V případě QR rozkladu se tím myslí to, že zadanou matici $X \in \mathbb{R}^{n,p}$ napíšeme jakou součin $$ X = QR, $$ kde $Q \in \mathbb{R}^{n,n}$ je ortogonální matice a $R \in \mathbb{R}^{n,p}$ je horní trojúhelníková matice, tj. má nenulové prvky pouze v šedé oblasti, jak je znázorněno na následujícím obrázku. Horní čtverec této matice si označme $R_1$ (na obrázku lehce modrá oblast).

Jelikož jsou (čtvercové) ortogonální matice vždy regulární, je hodnost matice $Q$ rovna $n$, a tedy musí být hodnosti matic $X$ a $R$ stejné. (Pokud jste absolventi BI-LIN a nechápete předchozí větu, napište prosím email autorovi těchto řádků, aby Vám mohl o stupeň zhoršit známku.) Metoda nejmenších čtverců funguje dobře, pouze pokud matice $X$ má hodnost $p$, neboli pokud jsou její sloupce lineárně nezávislé. Pouze tehdy totiž existuje inverze $(X^TX)^{-1}$. A jelikož hodnosti $X$ a $R$ jsou stejné, je hodnost matice $R$, a tedy i $R_1$, rovna $p$. To mimo jiné znamená, že $R_1$ má na diagonále samá nenulová čísla.

QR rozklad matice a výpočet odhadu metodou nejmenších čtverců

Předpokládejme, že jsme našli QR rozklad matice $X$ (jak se to dělá, si řekneme níže): $$ X = QR. $$ Z vlastnosti ortogonální matice plyne, že platí $$ Q^TX = Q^TQR = ER = R. $$ Odhad metodou nejmenších čtverců je definován jako bod minima funkce $$ ||Y - X\beta ||^2. $$ Je-li $Q$ ortogonální, je $Q^T$ také ortogonální, neb platí $QQ^T = (Q^T)^T Q^T = E$. Jelikož víme, že násobení ortogonální maticí $Q^T$ nemění velikost vektoru, platí $$ ||Y - X\beta ||^2 = ||Q^T(Y - X\beta) ||^2. $$ To můžeme dále upravit následovně: $$ ||Q^T(Y - X\beta) ||^2 = ||Q^TY - Q^TX\beta ||^2 = ||Q^TY - Q^TQR\beta ||^2 = ||Q^TY - R\beta ||^2. $$ Označíme-li si složky vektoru $Q^TY - R\beta$ z $\mathbb{R}^n$ jako $v_1,v_2, \ldots, v_n$, je $$ ||Q^TY - R\beta ||^2 = v_1^2 + v_2^2 + \cdots + v_n^2 = \sum_{i = 1}^p v_i^2 + \sum_{i = p+1}^n v_i^2. $$ Rozdělme si vektor $Q^TY \in \mathbb{R}^n$ takto: $$ \mathbb{R}^n = \begin{pmatrix} b_1 \\ b_2 \end{pmatrix}, $$ kde $b_1 \in \mathbb{R}^p$ je prvních $p$ složek a $b_2 \in \mathbb{R}^{n-p}$ jsou ty zbylé.

S tímto značením a s vědomím toho, že u matice $R$ jsou řádky $p+1$ až $n$ kompletně nulové, můžeme říci následující: složky $v_1, v_2, \ldots, v_p$ vektoru $Q^TY - R\beta$ tvoří vektor $b_ 1 - R_1\beta$ a zbylé složky vektor $b_2$. Můžeme tedy psát $$ \sum_{i = 1}^p v_i^2 + \sum_{i = p+1}^n v_i^2 = ||b_1 - R_1 \beta||^2 + || b_2 ||^2. $$ Hledat $\beta$ minimalizující výraz $$ ||Y - X\beta ||^2 $$ je tedy to samé, jako hledat $\beta$ minimalizující výraz $$ ||b_1 - R_1 \beta||^2 + || b_2 ||^2, $$ což je ale o poznání snazší úkol, neb minimum je nabyto pro $\beta$, pro které platí $$ R_1\beta = b_1. $$ A to je rovnice, která se vyřeší snadno, neboť $R_1$ je horní trojúhelníková matice. Lakonicky řečeno, potřebujeme vyřešit soustavu, kde už máme GEM hotovou!

Máme-li tedy $QR$ rozklad matice $X$, redukuje se hledání odhadu metodou nejmenších čtverců na vynásobení $Q^TY$, abychom spočítali $b_1$, a vyřešení jednoduché soustavy rovnic $R_1 \beta = b_1$. Zbývá tedy zjistit, jak QR rozklad najít.

Jedna metoda pro výpočet QR rozkladu

Existují dve numericky přívětivé metody pro výpočet QR rozkladu: Givensovy rotace, které si blíže představíme, a Householderovy reflexe. Obě metody fungují tak, že se konstruuje posloupnost ortogonálních matic $Q_1, Q_2, \ldots, Q_k$ taková, že násobení těmito maticemi z matice $X$ postupně vyrobí hledanou horní trojúhelníkovou matici $R$, tedy $$ Q_k \cdots Q_2 Q_1 X = R. $$ Jelikož součin ortogonálních matic, je opět ortogonální matice, můžeme položit $$ Q^T = Q_k \cdots Q_2 Q_1, $$ pak platí (pamatujte, že $Q^TQ = E$) $$ X = QR. $$

Jak tedy hledat matice $Q_i$? Myšlenka Givensových rotací je prostá a ukážeme si ji na případu, kdy matice $X$ je čtvercová matice $2 \times 2$: $$ X = \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix}. $$ Hledáme ortogonální matici $Q \in \mathbb{R}^{2,2}$ takovou, že $$ Q_1X = R = \begin{pmatrix} \alpha & ? \\ 0 & ? \end{pmatrix}, $$ kde $?$ značí, že na daném místě může být jakékoli číslo. To se redukuje na problém hledání matice $Q_1$ tak, že $$ Q_1\begin{pmatrix} a_{11} \\ a_{21} \end{pmatrix} = \begin{pmatrix} \alpha \\ 0 \end{pmatrix} $$ Jelikož $Q_1$ má být ortogonální, nemění tedy velikost vektorů a musí platit, že velikost $(a_{11}\ a_{21} )$ musí být stejná jako velikost vektoru $(\alpha\ 0)$. To znamená, že $$ \alpha = \sqrt{a_{11}^2 + a_{21}^2}. $$ Podíváme-li se na to geometricky, potřebujeme vlastně vektor $(a_{11}\ a_{21} )$ otočit tak, aby ležel na ose $x$ a $y$nová složka byla nulová.

Potřebujeme tedy matici rotace o úhel $\varphi$ (viz obrázek). Jednoduché geometrické cvičení využívající definici goniometrických funkcí pomocí pravoúhlého trojúhelníku (sinus je protilehlá ku přeponě, atp.) vede k výsledku $$ \cos \varphi = \frac{a_{11}}{\sqrt{a_{11}^2 + a_{21}^2}} \quad \text{a} \quad \sin \varphi = \frac{a_{21}}{\sqrt{a_{11}^2 + a_{21}^2}}. $$ Pokud si to zrovna nepamatujeme z BI-LIN, snadno si ověříme, že hledaná rotační matice je rovna $$ Q_1 = \begin{pmatrix} \cos \varphi & \sin \varphi \\ -\sin \varphi & \cos \varphi \end{pmatrix}, $$ skutečně: $$ \begin{pmatrix} \cos \varphi & \sin \varphi \\ -\sin \varphi & \cos \varphi \end{pmatrix} \begin{pmatrix} a_{11} \\ a_{21} \end{pmatrix} = \begin{pmatrix} \frac{a_{11}}{\sqrt{a_{11}^2 + a_{21}^2}} & \frac{a_{21}}{\sqrt{a_{11}^2 + a_{21}^2}} \\ -\frac{a_{21}}{\sqrt{a_{11}^2 + a_{21}^2}} & \frac{a_{11}}{\sqrt{a_{11}^2 + a_{21}^2}} \end{pmatrix} \begin{pmatrix} a_{11} \\ a_{21} \end{pmatrix} = \begin{pmatrix} \frac{a_{11}^2 + a_{21}^2}{\sqrt{a_{11}^2 + a_{21}^2}} \\ \frac{a_{21}a_{11} - a_{11}a_{21}}{\sqrt{a_{11}^2 + a_{21}^2}} \end{pmatrix} = \begin{pmatrix} \alpha \\ 0 \end{pmatrix}. $$

Příklad

Tak a máme vlastně hotovo. Podobnou myšlenku můžeme aplikovat i na větší matice: postupně pomocí jedoduchých matic vynulujeme všechny prvky $X$, které vynulovat chceme. Ukážeme si to na jednoduchém příkladě, kdy najdeme $QR$ rozklad matice $$ X = \begin{pmatrix} 1 & 0 \\ 2 & 1 \\ 1 & 3 \end{pmatrix}. $$ Nejprve najdeme matici $Q_1$, pomocí které vyrobíme nulu na místě, kde je nyní číslo $2$. Potřebujeme otočit vektor $(1\ 2)$ o úhel $\varphi_1$, který tento vektor svírá s osou $x$. Spočítáme $$ \cos \varphi_1 = \frac{1}{\sqrt{5}} \quad \text{a} \quad \sin \varphi_1 = \frac{2}{\sqrt{5}} $$ a získáme $$ Q_1X = \begin{pmatrix} \frac{1}{\sqrt{5}} & \frac{2}{\sqrt{5}} & 0 \\ -\frac{2}{\sqrt{5}} & \frac{1}{\sqrt{5}} & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 2 & 1 \\ 1 & 3 \end{pmatrix} = \begin{pmatrix} \sqrt{5} & \frac{2}{\sqrt{5}} \\ 0 & \frac{1}{\sqrt{5}} \\ 1 & 3 \end{pmatrix}. $$ Matici $Q_1$ jsme vyrobili z jednotkové matice tak, že jsme příslušné čtyři souřadnice nahradili danou rotační maticí. Výsledkem je, že matice rotuje první dva řádky a se zbytkem nedělá nic.

Nyní se potřebujeme zbavit $1$ dole v prvním sloupci, tedy otočit vektor $(\sqrt{5}\ 1)$ do směru osy $x$, tedy o úhel $\varphi_2$, pro který platí $$ \cos \varphi_2 = \frac{\sqrt{5}}{\sqrt{6}} \quad \text{a} \quad \sin \varphi_2 = \frac{1}{\sqrt{6}}. $$ Násobíme $$ Q_2Q_1X = \begin{pmatrix} \frac{\sqrt{5}}{\sqrt{6}} & 0 & \frac{1}{\sqrt{6}} \\ 0 & 1 & 0 \\ -\frac{1}{\sqrt{6}} & 0 & \frac{\sqrt{5}}{\sqrt{6}} \end{pmatrix} \begin{pmatrix} \sqrt{5} & \frac{2}{\sqrt{5}} \\ 0 & \frac{1}{\sqrt{5}} \\ 1 & 3 \end{pmatrix} = \begin{pmatrix} \sqrt{6} & \frac{5}{\sqrt{6}} \\ 0 & \frac{1}{\sqrt{5}} \\ 0 & \frac{13}{\sqrt{6}\sqrt{5}} \end{pmatrix}. $$

Teď zbývá vynulovat $\frac{13}{\sqrt{6}\sqrt{30}}$, resp. otočit vektor $\left(\frac{1}{\sqrt{5}}\ \frac{13}{\sqrt{5}\sqrt{6}}\right)$ o úhel $\varphi_3$, který tento vektor svírá s osou $x$. Výsledné výrazy už jsou poněkud rozšafné, tak je necháme na Vaší fantazii a napíšeme jen, že výsledný výpočet je $$ Q_3Q_2Q_1X = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos \varphi_3 & \sin \varphi_3 \\ 0 & -\sin\varphi_3 & \cos \varphi_3 \end{pmatrix} \begin{pmatrix} \sqrt{6} & \frac{5}{\sqrt{6}} \\ 0 & \frac{1}{\sqrt{5}} \\ 0 & \frac{13}{\sqrt{6}\sqrt{5}} \end{pmatrix} = \begin{pmatrix} \sqrt{6} & \frac{5}{\sqrt{6}} \\ 0 & \frac{\sqrt{35}}{\sqrt{6}} \\ 0 & 0 \end{pmatrix}, $$ kde $\frac{\sqrt{35}}{\sqrt{6}}$ je velikost vektoru $\left(\frac{1}{\sqrt{5}}\ \frac{13}{\sqrt{5}\sqrt{6}}\right)$. Hledaný $QR$ rozklad je tedy dán maticemi $$ Q = (Q_3Q_2Q_1)^T \quad \text{a} \quad R = \begin{pmatrix} \sqrt{6} & \frac{5}{\sqrt{6}} \\ 0 & \frac{\sqrt{35}}{\sqrt{6}} \\ 0 & 0 \end{pmatrix}. $$ Uf!

---

Tento článek byl vznikl v rámci řešení projektu Institucionální podpora Českého vysokého učení technického v Praze (CZ.02.2.69/0.0/0.0/16_015/0002382).

---