Metoda nejmenších čtverců

Graf zobrazuje proložení paraboly experimentálními daty zatíženými chybou. Koeficienty kvadratické funkce jsou nalezené metodou nejmenších čtverců.

Metoda nejmenších čtverců je matematicko-statistická metoda pro aproximaci řešení přeurčených soustav rovnic (tj. soustav, kde je více rovnic, než neznámých). +more „Nejmenší čtverce“ znamenají, že výsledné řešení má minimalizovat součet čtverců odchylek vůči každé rovnici. Metoda je v základní podobě určená pro řešení nekompatibilních soustav lineárních rovnic (v obecnější podobě hovoříme o nelineární metodě nejmenších čtverců), díky čemuž je fakticky ekvivalentní tzv. lineární regresi.

S nejjednodušší aplikací metody nejmenších čtverců se setkáváme například při prokládání (aproximaci) naměřených jednorozměrných dat přímkou. Nepatrně složitější aplikací je proložení dat parabolou, obecným polynomem předem daného stupně, nebo obecnou lineární kombinací předem daných bázových funkcí. +more Fakt, že proložení dat polynomem libovolného, ale předem daného stupně je stále lineární regresí, je častým zdrojem nedorozumění a terminologických nejasností. Další jednoduchou aplikací je nalezení nejpravděpodobnějšího průsečíku několika přímek (jejichž matematický popis je zatížen chybou) v rovině. Metoda nejmenších čtverců má velmi mnoho dalších aplikací v nejširším okruhu vědních oborů, ve kterých se setkáváme s nepřesnými daty, od statistiky a ekonomie, přes geodézii až po zpracování signálů a teorii řízení.

Obecně metoda nejmenších čtverců slouží k eliminaci chyb, kterou provádí optimálně vzhledem k pevně danému jednoznačnému kritériu (viz níže). Optimálně eliminovat chyby v datech lze i vzhledem k jiným kriteriím, takový postup může vést na metody převoditelné na metodu nejmenších čtverců (při použití různých typů vážení, např. +more když je známo, že chyba některých měření se výrazně liší od zbytku), nebo na metody obecně nepřevoditelné (nebo obtížně převoditelné) na metodu nejmenších čtverců (např. úplný problém nejmenších čtverců).

Historie

Metodu nejmenších čtverců poprvé použil Carl Friedrich Gauss v roce 1795 pro eliminaci chyb při geodetickém vyměřování.

Odvození metody

Uvažujme lineární aproximační problém

: \bold A\bold x\approx\bold b, \qquad \bold A \in\mathbb{R}^{n\times m}, \quad \bold x \in\mathbb{R}^{m}, \quad \bold b \in\mathbb{R}^{n}

takový, že \bold b\not\in\mathcal{R}(\bold A) (jinými slovy, neexistuje žádný vektor \bold x takový, aby nastala rovnost); symbol \mathcal{R}(\bold A) značí obor hodnot matice \bold A, tedy lineární obal jejích sloupců. Metoda nejmenších čtverců hledá vektor \bold x_{LS} (LS je zkratkou anglického least squares) splňující

: \|\bold A\bold x_{LS}-\bold b\|_2 = \min_{\bold x\in\mathbb{R}^m} \|\bold A\bold x-\bold b\|_2,

nebo ekvivalentně, opravu pravé strany \bold e_{LS} splňující

: \|\bold e_{LS}\|_2 = \min_{\bold e\in\mathbb{R}^n} \|\bold e\|_2, \quad \text{kde} \quad (\bold b+\bold e)\in\mathcal{R}(\bold A).

Zřejmě \bold A \bold x_{LS} = \bold b + \bold e_{LS} a oprava pravé strany \bold e = \bold A \bold x - \bold b je residuum odpovídající vektoru \bold x. Označíme-li e_1,\ldots,e_n jednotlivé komponenty vektoru \bold e, pak z definice Eukleidovské normy přímo plyne

: \|\bold A\bold x - \bold b\|_2=\|\bold e\|_2=(\sum_{j=1}^n e_j^2)^{1/2} \qquad \Longleftrightarrow \qquad \|\bold A\bold x - \bold b\|_2^2 = \|\bold e\|_2^2=\sum_{j=1}^n e_j^2.

Takže minimalizujeme součet čtverců jednotlivých komponent vektoru \bold e, tj. odchylek jednotlivých komponent pravé strany \bold b. +more Nalezení vztahu pro vektor \bold x_{LS} splňující podmínku minimality můžeme provést buď pomocí derivací (hledáním lokálního extrému) nebo přímo pomocí Pythagorovy věty.

Minimalizace kvadratického funkcionálu: Residuum \bold e = \bold A \bold x - \bold b je vektor, jehož normu minimalizujeme. Zřejmě

:\arg\min \|\bold e\|_2 = \arg\min \|\bold e\|_2^2,\quad \text{kde} \quad \|\bold e\|_2^2 = \bold e^T \bold e = \sum e_j^2,

Úlohu tedy můžeme zapsat jako minimalizaci kvadratického funkcionálu, kde volné parametry minimalizace jsou komponenty vektoru \bold x. Lokální extrém funkcionálu leží v bodě, kde se derivace funkcionálu podle všech komponent vektoru \bold x rovnají nule. +more Z vlastností normy je zřejmé, že se jedná o lokální minimum, které je zároveň minimem globálním. V maticovém zápisu lze formálně použít derivaci podle vektoru \bold x. Dostáváme.

:\left( \bold e^T \bold e \right)' = \left[ \left(\bold{Ax-b}\right)^T\left(\bold{Ax-b}\right) \right]' = \left[ \bold x^T \bold A^T \bold{Ax} - \bold x^T \bold A^T \bold b - \bold b^T \bold{Ax} + \bold b^T \bold b\right]' = 2 \bold A^T \bold{Ax} - 2 \bold A^T \bold b = 0.

Ortogonalita: Hledáme vektor \bold e s minimální normou, tj. eukleidovskou délkou. +more Z Pythagorovy věty vyplývá, že nejkratší vzdálenost dostaneme tehdy, pokud bude residuum \bold e ortogonální na obor hodnot \mathcal{R}(\bold A). Minimality tedy dosáhneme požadavkem.

:\bold e \perp \mathcal{R}(\bold A) \qquad \Longleftrightarrow \qquad \bold A^T\bold e = 0.

Dosazením za residuum dostáváme

:\bold A^T (\bold A \bold x- \bold b) = \bold A^T \bold A \bold x - \bold A^T \bold b = 0.

V obou případech tedy dostáváme, že hledaný vektor \bold x_{LS} je řešením tzv. normální soustavy rovnic

:\bold A^T \bold{Ax} = \bold A^T \bold b.

Za předpokladu, že \bold A má lineárně nezávislé sloupce (ekvivalentně, že \bold A^T \bold A je regulární), pak

:\bold x_{LS} = \left( \bold A^T \bold A \right)^{-1} \bold A^T \bold b,

Z předpokladu lineární nezávislosti sloupců přímo plyne, že řešení je jednoznačné.

Při praktickém výpočtu se ovšem normální soustava nikdy nesestavuje z důvodů numerické stability výpočtu. Použití zde odvozených vztahů pro praktický netriviální výpočet je jen obtížně ospravedlnitelné.

Pokud má matice \bold A lineárně závislé sloupce, je situace nepatrně komplikovanější. Řešení ve smyslu nejmenších čtverců tak jak je definované není jednoznačné. +more Zpravidla se pak uvažuje řešení ve smyslu nejmenších čtverců minimální v normě. Tedy mezi všemi přípustnými řešeními se hledá řešení s nejmenší normou. Opět užitím Pythagorovy věty můžeme minimalitu normy nahradit ortogonalitou. Řešení minimální v normě je pak ortogonální na nulový prostor matice \mathcal{N}(\bold A). Je tedy dáno vztahem.

:\bold x_{LS} = \bold A^{\dagger} \bold b,

kde \bold A^{\dagger} je Mooreova-Penroseova pseudoinverze. V tomto kontextu je třeba zmínit, že určení (numerické) hodnosti matice, tj. +more faktu, zda jsou sloupce matice \bold A (numericky) lineárně nezávisle, je samo o sobě velmi obtížná úloha.

V následujících odstavcích se seznámíme s nejdůležitějšími aplikacemi metody nejmenších čtverců, vesměs se však jedná pouze o návody jak sestavit matici \bold A a pravou stranu \bold b v konkrétních případech.

Řešení úlohy nejmenších čtverců

Metod k nalezení řešení \bold x_{LS} pro daný aproximační problém \bold A\bold x\approx\bold b je celá řada a dají se rozdělit podle dvou zásadních kritérií:

* velikost (a řídkost) problému (za velké považujeme zpravidla takové problémy, kdy s maticí \bold A nelze z paměťových a časových důvodů pracovat jako s polem obsahujícím všech mn čísel, u reálných problémů mohou být rozměry matic v řádech stovek milionů i vyšší), * hodnost matice \bold A, neboli počet lineárně nezávislých sloupců (pokud má matice lineárně závislé sloupce je řešení výrazně komplikovanější).

QR rozklad

Mezi základní metody řešení patří QR rozklad (počítaný buď ortogonálními transformacemi, nebo Gram-Schmidtovým ortogonalizačním procesem). Předpokládejme, že matice \bold A má lineárně nezávislé sloupce, pak existuje QR rozklad ve tvaru

:\bold Q\bold R=[\bold A,\bold b], \qquad \bold Q=[\bold Q_1,\bold q_{m+1}]\in\mathbb{R}^{n\times (m+1)}, \quad \bold R=\left[\begin{array}{cc}\bold{R}_{1,1}&\bold r_{1,m+1}\\&\rho_{m+1,m+1}\end{array}\right]\in\mathbb{R}^{(m+1)\times(m+1)}

kde \bold Q má ortonormální sloupce, \bold Q^T\bold Q=\bold I_m a \bold R je horní trojúhelníková. Zřejmě platí \bold A=\bold Q_1 \bold R_{1,1} a \bold b = \bold Q_1 \bold r_{1,m+1} + \bold q_{m+1}\rho_{m+1,m+1}. +more Formálním dosazením do normální soustavy rovnic dostáváme.

:\bold R_{1,1}^T \bold R_{1,1} \bold x = \bold R_{1,1}^T \bold r_{1,m+1}.

Po vynásobení maticí \bold R_{1,1}^{-T} zleva vidíme, že

:\bold x_{LS} = \bold R_{1,1}^{-1}\bold r_{1,m+1}.

Řešení tedy dostaneme snadno řešením soustavy s horní trojúhelníkovou maticí \bold R_{1,1} (ta je mimochodem Choleského faktorem matice \bold A^T\bold A ovšem spočteným numericky stabilní cestou).

Výpočet je použitelný pouze pokud má matice \bold A (numericky) lineárně nezávislé sloupce (postup lze rozšířit i pro matice s lineárně závislými sloupci pomocí tzv. RRQR rozkladu, obecně tak získáme pouze aproximaci řešení). +more Matice \bold Q při nalezení vlastního řešení potřeba není, můžeme tedy bez problémů použít levnější (ovšem méně přesný) Gram-Schmidtův proces. Pro rozsáhlé úlohy s velkým počtem sloupců může snadno dojít k rychlému zaplnění trojúhelníkového faktoru, proto je tento postup pro tyto typy úloh prakticky nepoužitelný.

Nebezpečí skryté v postupu hledání řešení problému nejmenších čtverců pomocí explicitního sestavení matice \bold A^T \bold A je dobře známe, viz např. Lauchliho matice.

Singulární rozklad

Má-li matice \bold A lineárně závislé sloupce, jedinou spolehlivou metodou je singulární rozklad a v podstatě explicitní sestavení Mooreovy-Penroseovy pseudoinverze. Pokud r=\mathrm{rank}(\bold A), a její ekonomický singulární rozklad je

:\bold A = \bold U_r \bold S_r \bold V_r^T, \qquad \text{pak} \qquad \bold A^{\dagger} = \bold V_r \bold S_r^{-1} \bold U_r^T

což umožňuje přímo nalézt řešení ve smyslu nejmenších čtverců minimální v normě.

Výpočet je použitelný pouze pro malé úlohy vzhledem k velké výpočetní náročnosti singulárního rozkladu.

Iterační metody

Pro řešení rozsáhlých úloh jsou standardem a jediným prakticky použitelným řešením iterační metody, zpravidla postavené na Krylovových podprostorech rostoucí dimenze. Klasickým představitelem je LSQR a další odvozené algoritmy CGLS, CGNE.

Aplikace: Regresní analýza

Metoda nejmenších čtverců bývá často používána při regresní analýze pro aproximaci zadaných (naměřených) hodnot nějakou funkcí z předepsaného prostoru. Nejjednodušším příkladem je proložení (aproximace) dat přímkou, tedy lineární funkcí. +more Protože klasická metoda nejmenších čtverců vede vždy na lineární aproximační model (de-facto soustavu rovnic) je často považována za ekvivalent pojmu lineární regrese. To souvisí s faktem, že libovolný aproximační problém \bold A \bold x \approx \bold b lze interpretovat jako lineární regresi v m-rozměrném prostoru, kde jsme v bodě o souřadnicích (a_{j,1},\ldots,a_{j,m}) naměřili nějakou veličinu b_j, j=1,\ldots,n a naměřené hodnoty prokládáme nadrovinou.

Pro první přiblížení uvažujme závislost nějaké nezávisle proměnné (např. fyzikální veličiny) y na nezávisle proměnné u, tj. +more y=f(u). přičemž závislost má nějaký předem daný tvar, např. f(u)=k\sin(u)+q\exp(-u), kde koeficienty k a q jsou známé. Pro identifikaci těchto koeficientů provedeme n měření veličiny y v přesně definovaných hodnotách nezávisle proměnné u. Získáme tak sadu n dvojic (u_j,y_j) přičemž hodnoty y_j jsou zatížené chybami měření. Dostaneme tak lineární aproximační problém.

:\bold A\bold x = \left[\begin{array}{cc}\sin(u_1)&\exp(-u_1)\\ \vdots & \vdots \\ \sin(u_n)&\exp(-u_n) \end{array}\right] \left[\begin{array}{c}k\\q\end{array}\right]\approx \left[\begin{array}{c}y_1\\\vdots\\y_n\end{array}\right]=\bold b.

Metoda nejmenších čtverců hledá takové hodnoty koeficientů k, q, aby součet čtverců „odchylek“ (reziduí) jejích funkčních hodnot od daných naměřených hodnot byl nejmenší možný.

Aproximace přímkou

Pro proložení daných hodnot přímkou (lineární funkcí, polynomem prvního řádu) s rovnicí y = f(u) = k_1u + k_0 dostaneme z normální soustavy rovnic vztahy

:k_1 = \frac{n \sum{u_i y_i} - \sum{u_i} \sum{y_i}}{n \sum{u_i^2} - \left( \sum{u_i} \right)^2}, \qquad k_0 = \frac{\sum{u_i^2}\sum{y_i} - \sum{u_i}\sum{u_i y_i}}{n \sum{u_i^2} - \left(\sum{u_i}\right)^2}.

Aproximace parabolou

Zadané hodnoty lze aproximovat parabolou (kvadratickou funkcí, polynomem druhého řádu) s rovnicí y = f(u) = k_2u^2 + k_1u + k_0. Optimální parametry k_i získáme opět řešením normální soustavy rovnic

: \bold A^T\bold A \bold x = \left[\begin{array}{ccc} \sum{u_i^4} & \sum{u_i^3} & \sum{u_i^2} \\ \sum{u_i^3} & \sum{u_i^2} & \sum{u_i} \\ \sum{u_i^2} & \sum{u_i} & n \end{array}\right] \left[\begin{array}{c}k_2\\k_1\\k_0\end{array}\right]= \left[\begin{array}{c}\sum{y_i u_i^2}\\\sum{y_i u_i} \\\sum{y_i}\end{array}\right]= \bold A^T \bold b.

Proložení bodů parabolou je lineární regresí kvadratické funkce (kvadratického modelu), zejména ve starší literatuře se můžeme setkat i s pojmem kvadratická regrese.

Aproximace polynomem

Ukázka aproximace zadaných bodů polynomem stupně s = 0,. +more,9. Obdobně jako v případě paraboly, i při aproximaci polynomem y = p_s(u) = k_s u^s + \ldots + k_1 u + k_0 přímým dosazením hodnot y_i, u_i získáme soustavu rovnic y_i=P_s(u_i). V maticovém zápisu.

:\bold A \bold x = \left[ \begin{matrix} u_1^s & \cdots & u_1 & 1 \\ \vdots & \ddots & \vdots & \vdots \\ u_n^s & \cdots & u_n & 1 \end{matrix} \right] \left[ \begin{matrix} k_s \\ \vdots \\ k_1 \\ k_0 \end{matrix} \right] \approx \left[ \begin{matrix} y_1 \\ \vdots \\ y_n \end{matrix} \right] = \bold b.

Přejdeme-li čistě formálně k normální soustavě rovnic (při řešení reálných problémů tuto soustavu nikdy nesestavujeme), dostáváme

: \bold A^T \bold A \bold x = \left[ \begin{matrix} \sum{u_i^{2s}} & \cdots & \sum{u_i^{s+1}} & \sum{u_i^s} \\ \vdots & \ddots & \vdots & \vdots \\ \sum{u_i^{s+1}} & \cdots & \sum{u_i^2} & \sum{u_i} \\ \sum{u_i^s} & \cdots & \sum{u_i} & n \end{matrix} \right] \left[ \begin{matrix} k_s \\ \vdots \\ k_1 \\ k_0 \end{matrix} \right] = \left[ \begin{matrix} \sum{y_i u_i^s} \\ \vdots \\ \sum{y_i u_i} \\ \sum{y_i} \end{matrix} \right] = \bold A^T \bold b,

která má vždy řešení. Pokud má matice \bold A lineárně nezávislé sloupce, pak (A^TA) je regulární a \bold x = \left( \bold A^T \bold A \right)^{-1} \bold A^T \bold b je jednoznačným řešením původního problému ve smyslu nejmenších čtverců.

Proložení bodů polynomem (předem) daného stupně s je lineární regresí polynomu (polynomiálního modelu), zejména ve starší literatuře se můžeme setkat i s pojmem polynomiální (polynomická) regrese.

Aplikace: Nalezení průsečíku několika nadrovin v prostoru

Nalezení průsečíku tří přímek (zatížených chybami) metodou nejmenších čtverců Opět se snažíme úlohu přeformulovat jako lineární aproximační problém \bold A\bold x\approx \bold b. +more Mechanismus lze nejsnáze ilustrovat na příkladu nalezení průsečíku několika přímek v rovině. Jsou dány tři přímky v rovině.

:u + v = 0, \quad 2u - v - 4 = 0, \quad u - y - 2 = 0,

které nemají společný průsečík. Zřejmě

: \bold A \bold x = \left[ \begin{matrix} 1 & 1 \\ 2 & -1 \\ 1 & -1 \end{matrix} \right] \left[ \begin{matrix} u \\ v \end{matrix} \right] \approx \left[ \begin{matrix} 0 \\ 4 \\ 2 \end{matrix} \right] = \bold b

a vektor \bold e = \bold A \bold x - \bold b je nenulový pro libovolné u a v. Za předpokladu, že přímky nejsou známy přesně (jejich popis je zatížen chybami) a chyby jsou obsaženy pouze v konstantních členech (tedy v pravé straně \bold b aproximačního problému), můžeme použít k nalezení nejpravděpodobnějšího místa průsečíku metodu nejmenších čtverců. +more Získáme tak bod.

:\bold x = \left( \bold A^T \bold A \right)^{-1} \bold A^T \bold b \doteq \left[ \begin{matrix} 1.286 \\ -1.143 \end{matrix} \right].

Všimněme si, že při přenásobení kterékoliv z rovnic definujících přímky nenulovou konstantou (různou od \pm1), pracujeme se stejnými přímkami ale hledaný bod vyjde jinde. Taková operace totiž není v problému nejmenších čtverců ekvivalentní operací, fakticky taková operace změní normu, ve které provádíme minimalizaci.

Odkazy

Reference

Související články

Lineární regrese * Kvadratická regrese * Polynomická regrese * Exponenciální regrese * Ortogonální regrese, error-in-variables modelling (EIV), úplný problém nejmenších čtverců (TLS) * Regresní analýza * Singular value decomposition (SVD), singulární rozklad

Externí odkazy

[url=http://herodes. feld. +morecvut. cz/mereni/mnc/mnc. php]http://herodes. feld. cvut. cz/…[/url] - aproximace přímkou, exponenciálou (linearizace) a polynomem * [url=http://amper. ped. muni. cz/jenik/nejistoty/html_tree/node10. html]http://amper. ped. muni. cz/…[/url] - aproximace obecnou funkcí, přímkou, polynomem, statistika * [url=http://cmp. felk. cvut. cz/cmp/courses/recognition/Lab_archive/RPZ_02-03l/crossval/crossval/node2. html]http://cmp. felk. cvut. cz/…[/url] - aproximace polynomem, soustava s obdélníkovou maticí * [url=http://www. vscht. cz/kot/era/]http://www. vscht. cz/kot/era/…[/url] - program ERA (autor Doc. Zámostný VŠCHT).

Kategorie:Lineární algebra Kategorie:Statistika Kategorie:Ekonometrie