04_teorie_chyb:0419_metody_reseni_normalnich_rovnic | IngGeo

« 18. Metody robustního odhadu
» 20. Kritéria pro testování opakovaných měření

19. Metody řešení normálních rovnicb

Úvod

Metodu volíme podle počtu rovnic a možnosti použití výpočetní techniky. K dispozici je velké množství metod, ať už vhodných pro případné ruční nebo počítačové zpracování. Zde budou zmíněny pouze některé vybrané metody, které se tradičně využívají a dále některé metody poskytující vyšší spolehlivost výsledků zejména z hlediska numerické přesnosti. Metody detailně popsané lze nalézt např. v [1] a [2], algoritmy jsou podrobně uvedeny (včetně zdrojových kódů) také v [3] (případně v [4]). Základní třídění prezentovaných metod řešení vyrovnávacích úloh:

a) Výpočet inverzní matice (jak odpovídá odvozením):
1. Jordanův algoritmus (Gaussova eliminace)
2. LU rozklad
3. QR rozklad
4. Pseudoinverze
b) Přímé řešení normálních rovnic:
1. Gaussova eliminační metoda
2. Choleskyho metoda
3. Metoda postupné iterace
  - Newtonovo řešení
  - Gaussovo řešení
c) Přímé řešení rovnic oprav pseudoinverzí
d) Řešení inverze singulární matice normálních rovnic

Jednotlivé metody nyní budou postupně popsány tak, aby byl zřejmý jejich princip a možnosti. K využití pseudoinverze je vhodné poznamenat, že jednotlivá využití vyplývají z jejích vlastností a proto jí bude věnována samostatný odstavec, kde budou popsány vlastnosti, postupy výpočtu a posléze jednotlivé možnosti využití. Vzhledem k tomu, že v současné době je běžné řešení velkých a špatně podmíněných soustav, lze tuto kapitolu doporučit zejména vzhledem k možnostem, které poskytuje singulární rozklad.

Při použití mnohých metod je z numerického hlediska pro zajištění proveditelnosti a kvality výpočtu vhodné provést přeskládání řešené matice nebo celé soustavy částečným nebo úplným výběrem hlavního prvku, s využitím permutační matice. Tomuto tématu bude věnován samostatný odstavec na konci kapitoly.

Výpočet inverzní matice

Ve všech předchozích odvozeních se řešení normálních rovnic provádí inverzí čtvercové matice. Numerické řešení inverze není triviální záležitostí, v mnoha případech je také numericky nestabilní a představuje největší problém při řešení vyrovnání. K dispozici je nepřeberné množství metod, např. v [5] lze nalézt metodu založenou na výpočtu determinantů, inverze matice rozdělením na bloky a další, které se pro inverzi matice o větším počtu řádků/sloupců nehodí. Dále proto budou uvedeny vybrané metody vhodné pro počítačové řešení o reálně využitelném počtu neznámých.

Pro inverzní matici ${{\mathbf N}}^{-1}$ k čtvercové matici ${\mathbf N}$ platí:

${\rm \ }{\mathbf N}{\rm \ }\cdot {{\mathbf N}}^{-1}={{\mathbf N}}^{-1}\cdot {\mathbf N}={\mathbf E} ,$

kde ${\mathbf E}$ je jednotková matice (tj. matice s jedničkami na diagonále a nulami všude jinde). Inverzní matice existuje vždy, pokud je matice ${\mathbf N}$ regulární, tj. pokud je čtvercová, její řádky jsou lineárně nezávislé, její sloupce jsou lineárně nezávislé, její hodnost je rovna jejímu řádu. Ekvivalentně pokud platí, že její determinant je různý od nuly.

Jordanův algoritmus

Podle [6] se invertovaná čtvercová matice ${\mathbf N}$ rozšíří o stejný počet sloupců jednotkové matice ${\mathbf E}$ . Matice ${\mathbf N}$ se odečítáním násobků lineárních kombinací ostatních řádků upravuje na matici jednotkovou (praktický postup Gaussovy eliminace, popsáno dále), stejné operace se provádějí na přidružené jednotkové matici. Po úpravě ${\mathbf N}$ na jednotkovou matici na místě ${\mathbf E}$ vznikla inverzní matice k ${\mathbf N}$ , tj. ${{\mathbf N}}^{-1}$ . Symbolicky lze zapsat ve tvaru:

${\mathbf (}{\mathbf N}{\mathbf |}{\mathbf E}{\mathbf )}\to {\mathbf (}{\mathbf E}{\mathbf |}{{\mathbf N}}^{-1}{\mathbf )}.$

Stejně jako u všech metod založených na Gaussově eliminaci je vhodné aplikovat výběr hlavního prvku před začátkem ekvivalentních úprav.

LU rozklad

Základem algoritmů inverze matic je vždy nějaká forma rozkladu (faktorizace) invertované matice na součin jednoduše invertovatelných matic. V případě LU rozkladu se jedná o trojúhelníkový rozklad čtvercové invertované matice tak, aby např. podle [7] platilo:

${\mathbf N}={\mathbf L}\cdot {\mathbf U} ,$

kde ${\mathbf L}$ je dolní trojúhelníková matice a ${\mathbf U}$ horní trojúhelníková matice. Rozklad existuje, pokud je matice ${\mathbf N}$ regulární. Matici ${\mathbf L}$ lze získat například Gaussovou eliminací obdobně jako u předchozí metody, matici ${\mathbf U}$ pak operací:

${\mathbf U}\ ={\mathbf N}\cdot {{\mathbf L}}^{-1} .$

Inverzní matice se pak vypočítá podle vzorce

${{\mathbf N}}^{-1}={{\mathbf U}}^{-1}\cdot {{\mathbf L}}^{-1} .$

Výpočet inverze dolní trojúhelníkové matice je jednoduchý, vždy je to také dolní trojúhelníková matice, na diagonále jsou reciproké hodnoty odpovídajících si prvků původní matice a ostatní prvky lze postupně dopočítat. Zde uvedeme pouze pro dolní trojúhelníkovou matici:

${\left({{\mathbf L}}^{-1}\right)}_{r,r}=\frac{1}{{{\mathbf L}}_{i,i}} ,$

${\left({{\mathbf L}}^{-1}\right)}_{r,s}=-\frac{1}{{{\mathbf L}}_{s,s}}\cdot \sum^r_{i=s+1}{\left({{{\mathbf L}}^{-1}}_{r,i}\cdot {{\mathbf L}}_{i,s}\right)} .$

V rovnici $r$ značí pořadí řádku, $s$ pořadí sloupce. Odvození vychází ze vztahů jednoduše ilustrovaných v následujících rovnicích:

$\begin{array}{ccccc} {\mathbf L} & \cdot & {{\mathbf L}}^{{\mathbf -}{\mathbf 1}} & = & {\mathbf E} \\ \left( \begin{array}{ccc} a & 0 & 0 \\ b & c & 0 \\ d & e & f \end{array} \right) & \cdot & \left( \begin{array}{ccc} 1/a & 0 & 0 \\ x & 1/c & 0 \\ y & z & 1/f \end{array} \right) & = & \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right) \end{array},$

$\left( \begin{array}{ccc} b & c & 0 \end{array} \right)\cdot \left( \begin{array}{c} 1/a \\ x \\ y \end{array} \right)=0 \to x=-\frac{1}{c}\cdot \left(b\cdot \frac{1}{a}\right),$

$\left( \begin{array}{ccc} d & e & f \end{array} \right)\cdot \left( \begin{array}{c} 1/a \\ x \\ y \end{array} \right)=0 \to y=-\frac{1}{f}\cdot \left(d\cdot \frac{1}{a}+e\cdot x\right).$

$\left( \begin{array}{ccc} d & e & f \end{array} \right)\cdot \left( \begin{array}{c} 0 \\ 1/c \\ z \end{array} \right)=0 \to z=-\frac{1}{f}\cdot \left(e\cdot \frac{1}{c}\right).$

LU rozklad je založen na stejném principu jako Gaussova eliminace, neřeší se však pouze soustava rovnic, ale rozkládá se matice na součin dvou snadno invertovatelných matic.

QR rozklad

Rozklad matice na ortonormální matici ${\mathbf Q}$ a horní trojúhelníkovou matici ${\mathbf R}$ , využívá se často právě při řešení úloh metodou nejmenších čtverců a také při výpočtu vlastních čísel.

Každá čtvercová matice ${\mathbf N}$ , jejíž prvky jsou reálná čísla, může být podle [8] rozložena do tvaru

${\mathbf N}={\mathbf Q}\cdot {\mathbf R} .$

Jestliže ${\mathbf N}$ není singulární, rozklad je jedinečný za předpokladu, že diagonální prvky matice ${\mathbf R}$ jsou kladné.

Pro výpočet inverzní matice pak platí:

${{\mathbf N}}^{-1}={\left({\mathbf Q}\cdot {\mathbf R}\right)}^{-1}={{\mathbf R}}^{-1}\cdot {{\mathbf Q}}^{-1}={{\mathbf R}}^{-1}\cdot {{\mathbf Q}}^T .$

Protože ${\mathbf Q}$ je ortonormální, platí:

${{\mathbf Q}}^{-1}={{\mathbf Q}}^T , {{\mathbf Q}}^T\cdot {\mathbf Q}={\mathbf E} ,$

kde ${\mathbf E}$ je jednotková matice.

Výpočet ortonormální matice ${\mathbf Q}$ lze provést více postupy, např. Householderovou transformací, pomocí Givensových rotací, zde bude podle [9] popsán postup Gram-Schmidtovy ortogonalizace. Jestliže se ortonormalizovaná matice ${\mathbf D}$ skládá z $n$ vektorů ${\mathbf D}=[ \begin{array}{ccc} {{\mathbf d}}_1 & \dots & {{\mathbf d}}_n \end{array} ]$ , ortonormální matice ${\mathbf O}{\mathbf = }\left[ \begin{array}{ccc} {{\mathbf o}}_1 & \dots & {{\mathbf o}}_n \end{array} \right]$ , která se skládá z $n$ jednotkových vzájemně ortogonálních vektorů se získá:

${{\mathbf o}}_1=\frac{{{\mathbf d}}_1}{\left|{{\mathbf d}}_1\right|} ,$

${{\mathbf p}}_2={{\mathbf d}}_2-({{\mathbf d}}_2,{{\mathbf o}}_1)\cdot {{\mathbf o}}_1 , {{\mathbf o}}_2=\frac{{{\mathbf p}}_2}{\left|{{\mathbf p}}_2\right|} ,$

pro $k=1\dots n$ :

${{\mathbf p}}_k{\rm =}{{\mathbf d}}_k{\rm -}\left({{\mathbf d}}_k,{{\mathbf o}}_{{\rm 1}}\right)\cdot {{\mathbf o}}_{{\rm 1}}{\rm -\dots -}\left({{\mathbf d}}_k,{{\mathbf o}}_{k{\rm -1}}\right)\cdot {{\mathbf o}}_{k{\rm -1}} , {{\mathbf o}}_k{\rm =}\frac{{{\mathbf p}}_k}{\left|{{\mathbf p}}_k\right|} ,$

kde symbol $\left(\ ,\ \right)$ značí skalární součin a $\left|\right|$ velikost vektoru.

Rozklad lze provést také pro obdélníkovou matici ${{\mathbf A}}_{(m\times n)}$ , kde $m\ge n$ , ve tvaru:

${{\mathbf A}}_{\left(m\times n\right)}={{\mathbf Q}}_{\left(m\times m\right)}\cdot {{\mathbf R}}_{\left(m\times n\right)} ,$

${\mathbf A}{\mathbf =}{\mathbf Q}\cdot \left[ \begin{array}{c} {{\mathbf R}}_1 \\ {\mathbf 0} \end{array} \right]{\mathbf =}\left[ \begin{array}{cc} {{\mathbf Q}}_{1_{(m\times n)}} & {{{\mathbf Q}}_2}_{(m\times (m-n))} \end{array} \right]\cdot \left[ \begin{array}{c} {{{\mathbf R}}_1}_{(n\times n)} \\ {\mathbf 0} \end{array} \right]{\mathbf =}{{\mathbf Q}}_{1_{(m\times n)}}\cdot {{{\mathbf R}}_1}_{(n\times n)} .$

Jestliže ${\mathbf A}$ má maximální možný počet lineárně nezávislých řádků a diagonální prvky ${\mathbf R}$ jsou kladné, pak ${{\mathbf R}}_1$ a ${{\mathbf Q}}_1$ jsou unikátní, ${{\mathbf Q}}_2$ nikoli. ${\mathbf \ }{{\mathbf R}}_1$ je rovno hornímu trojúhelníkovému tvaru Choleskyho dekompozice matice ${\mathbf N}={{\mathbf A}}^T\cdot {\mathbf A}$ .

Numerickou stabilitu výpočtu lze zlepšit pomocí permutační matice ${\mathbf P}$ , kdy se přeskládají sloupce tak, aby diagonální prvky tvořily neklesající řadu. Rozklad pak má tvar:

${\mathbf A}\cdot {\mathbf P}={\mathbf Q}\cdot {\mathbf R}\Rightarrow {\mathbf A}={\mathbf Q}\cdot {\mathbf R}\cdot {{\mathbf P}}^T .$

Z toho inverzní matice:

${{\mathbf A}}^{-1}={\mathbf P}\cdot {{\mathbf R}}^{{\mathbf -}1}\cdot {{\mathbf Q}}^T .$

Poznámka: Permutační matice slouží k přeskládání původní matice, je tvořena ortonormálními vektory (vždy v každém řádku i sloupci je právě jedna nenulová hodnota, která je zároveň rovna jedné). Protože je to matice ortonormální, její inverze je zároveň rovna transpozici.

Přímé řešení normálních rovnic

Gaussova eliminační metoda

Podstata Gaussovy eliminační metody je řečena již v názvu. Jednoduchými matematickými operacemi (násobení, sčítání) - redukcemi, dosáhneme postupné eliminace jednotlivých neznámých. Výpočet je mechanicky upraven do algoritmu.

Postup redukcí bude ukázán na jednoduchém příkladu. Zvolme soustavu tří normálních rovnic:

$\begin{array}{c} \left[{pa}_{{\rm 1}}a_{{\rm 1}}\right]{dx}_{{\rm 1}}{\rm +}\left[pa_{{\rm 1}}a_{{\rm 2}}\right]dx_{{\rm 2}}{\rm +}\left[pa_{{\rm 1}}a_{{\rm 3}}\right]dx_{{\rm 3}}{\rm +[}pa_{{\rm 1}}l{\rm ']=0,} \\ \left[{pa}_{{\rm 1}}a_{{\rm 2}}\right]{dx}_{{\rm 1}}{\rm +}\left[{pa}_{{\rm 2}}a_{{\rm 2}}\right]{dx}_{{\rm 2}}{\rm +}\left[{pa}_{{\rm 2}}a_{{\rm 3}}\right]{dx}_{{\rm 3}}{\rm +[}{pa}_{{\rm 2}}l{\rm ']=0,} \\ \left[{pa}_{{\rm 1}}a_{{\rm 3}}\right]{dx}_{{\rm 1}}{\rm +}\left[{pa}_{{\rm 2}}a_{{\rm 3}}\right]{dx}_{{\rm 2}}{\rm +}\left[{pa}_{{\rm 3}}a_{{\rm 3}}\right]{dx}_{{\rm 3}}{\rm +[}{pa}_{{\rm 3}}l{\rm ']=0.} \end{array}$

První neznámou eliminujeme tak, že násobíme první rovnici postupně zlomky:

nejdříve $-[pa_1a_2]/[pa_1a_1]$ a přičteme ji k rovnici druhé,

potom $-[pa_1a_3]/[pa_1a_1]$ a přičteme ji k rovnici třetí.

Je to tzv. první redukce normálních rovnic, kde dostaneme počet rovnic zmenšený o jednu:

$\left(\left[pa_2a_2\right]-\left(\frac{\left[pa_1a_2\right]}{\left[pa_1a_1\right]}\right)\cdot \left[pa_1a_2\right]\right)\cdot dx_2+\left(\left[pa_2a_3\right]-\left(\frac{\left[pa_1a_2\right]}{\left[pa_1a_1\right]}\right)\cdot \left[pa_1a_3\right]\right)\cdot dx_3+([pa_2l']-\{[pa_1a_2]/[pa_1a_1])\cdot [pa_1l'])=0,$

$\left(\left[pa_2a_3\right]-\left(\frac{\left[pa_1a_3\right]}{\left[pa_1a_1\right]}\right)\cdot \left[pa_1a_2\right]\right)\cdot dx_2+\left(\left[pa_3a_3\right]-\left(\frac{\left[pa_1a_3\right]}{\left[pa_1a_1\right]}\right)\cdot \left[pa_1a_3\right]\right)dx_3+([pa_3l']-\ -([pa_1a_3]/[pa_1a_1])\cdot [pa_1l'])=0.$

Pro redukované členy zavedeme Gaussovo označení, např. $[pa_2a_2.1]$ , což znamená „ $[pa_2a_2]$ po první redukci“. Dostaneme tak dvě jednou redukované normální rovnice o dvou neznámých:

$\begin{array}{c} \left[{pa}_2a_2.1\right]\cdot {dx}_2+\left[{pa}_2a_3.1\right]\cdot {dx}_3+[{pa}_2l'.1]=0,\ \\ \left[{pa}_2a_3.1\right]\cdot {dx}_2+\left[{pa}_3a_3.1\right]\cdot {dx}_3+[{pa}_3l'.1]=0.\ \end{array}$

Podílem $-[pa_2a_3.1]/[pa_2a_2.1]$ násobíme nyní první redukovanou rovnici, přičteme k druhé a dostaneme dvakrát redukovanou normální rovnici:

$\left(\left[pa_3a_{3.1}\right]-\left(\frac{\left[pa_2a_{3.1}\right]}{\left[pa_2a_{2.1}\right]}\right)\cdot \left[pa_2a_{3.1}\right]\right)\cdot dx_3+\left(\left[pa_3l'.1\right]-\left(\frac{\left[pa_2a_{3.1}\right]}{\left[pa_2a_{2.1}\right]}\right)\cdot \left[pa_2l'.1\right]\right)=0$

kterou napíšeme s Gaussovými symboly:

$\left[pa_3a_3.2\right]\cdot dx_3+[pa_3l'.2]=0.$

Tato druhá redukce dává již jen jednu rovnici o jedné neznámé. Při čtyřech neznámých by přibyla čtvrtá normální rovnice a byly by nutné tři redukce. K výpočtu neznámých použijeme vždy první rovnice z původních a postupně redukovaných rovnic:

$\left[{pa}_1a_1\right]\cdot {dx}_1+\ \left[{pa}_1a_2\right]\cdot {dx}_2+\left[{pa}_1a_3\right]\cdot {dx}_3+[{pa}_1l']=0 ,$

$\left[pa_2a_2.1\right]\cdot dx_2+\left[pa_2a_3.1\right]\cdot dx_3+[pa_2l'.1]=0$ (po první redukci),

$\left[pa_3a_3.2\right]\cdot dx_3+[pa_3l'.2]=0$ (po druhé redukci),

z nichž určíme postupně neznámé:

$dx_3=-[pa_3l'.2]/[pa_3a_3.2],$

$dx_2=-([pa_2l'.1]+\left[pa_2a_3.1\right]\cdot dx_3)/[pa_2a_2.1],$

$dx_1=-([pa_1l']+\left[pa_1a_3\right]\cdot dx_3+\ \left[pa_1a_2\right]\cdot dx_2)/[pa_1a_1].$

Postupným dosazením za $dx_3$ , $dx_2$ na pravé straně bychom mohli případně každou neznámou vyjádřit obecně redukovanými měřeními $l_i'$ .

Algoritmus řešení lze tedy jednoduše vyjádřit pro soustavu normálních rovnic v maticovém vyjádření ${\mathbf A}\cdot {\mathbf x}-{\mathbf b}{\mathbf =}{\mathbf 0}{\mathbf \ }$ následujícím způsobem. Prvky matice ${\mathbf b}$ zahrneme jakoby do posledního sloupce matice ${\mathbf A}$ , tedy $b_i=a_{i,n+1}$ :

$pro\ k=1,...,n-1:$

$pro\ i=k+1,...,n:$

$pro\ j=k+1,...,n+1:$

$a_{ij}:=a_{ij}-a_{ik}\cdot a_{kj}/a_{kk}$

$x_n=a_{n,n+1}/a_{nn}, pro\ i=n-1,...,1:x_i=\frac{\left(a_{i,n+1}-\sum^n_{j=i+1}{a_{ij}}\cdot x_j\right)}{a_{ji}} .$

Pro dosažení správného a přesného řešení je vhodné před řešením soustavy provést přeskládání (pivoting) soustavy.

Choleskyho metoda

Je vhodná pro řešení tzv. špatně podmíněných soustav rovnic. Také je významně rychlejší, vyžaduje $1/2$ operací oproti Gaussově eliminaci a pro získání stabilního a přesného řešení není třeba přeskládání prvků před řešením (pivoting). Řešení odpovídá rozkladu symetrické pozitivně definitní matice ${\mathbf N}$ na dvě trojúhelníkové, vzájemně transponované matice ${\mathbf C}$ :

${\mathbf N} {\mathbf = }{{\mathbf C}}^T\cdot {\mathbf C} .$

Rozklad je unikátní, existuje pouze jedna dolní trojúhelníková matice s pozitivními hodnotami na diagonále, pokud jsou hodnoty matice ${\mathbf N}$ reálná čísla, pro prvky matice ${\mathbf C}$ platí totéž. Řešený systém normálních rovnic pak má tvar

${{\mathbf C}}^T\cdot {\mathbf C}\cdot {\mathbf x}{\rm \ }-{{\mathbf A}}^T\cdot {\mathbf P}\cdot {\mathbf l}{\rm =}{\mathbf 0} ,$

který se rozdělí na dvě na ekvivalentní rovnice:

${\mathbf C}\cdot {\mathbf x}{\mathbf =}{\mathbf y}, {{\mathbf C}}^T\cdot {\mathbf y}{\mathbf =}{{\mathbf A}}^T\cdot {\mathbf P}\cdot {\mathbf l} .$

Nejprve se určí neznámé ${\mathbf y}$ , poté ${\mathbf x}$ . Získat trojúhelníkovou matici lze např. algoritmem založeným na Gaussově eliminaci ([10]), zde bude uveden princip Cholesky-Banachiewicz a Cholesky-Crout algoritmu (rozdíl je pouze ve směru výpočtu v matici, princip je stejný). Rozklad matice po roznásobení pro matici o rozměru $3\times 3$ :

${\mathbf N}={{\mathbf C}}^T\cdot {\mathbf C}{\mathbf =}\left( \begin{array}{ccc} C_{11} & 0 & 0 \\ C_{21} & C_{22} & 0 \\ C_{31} & C_{32} & C_{33} \end{array} \right)\cdot \left( \begin{array}{ccc} C_{11} & C_{21} & C_{31} \\ 0 & C_{22} & C_{32} \\ 0 & 0 & C_{33} \end{array} \right) ,$

${\mathbf N}=\left( \begin{array}{ccc} C^2_{11} & C_{21}\cdot C_{11} & C_{31}\cdot C_{11} \\ C_{21}\cdot C_{11} & C^2_{21}+C^2_{22} & C_{31}\cdot C_{21}+C_{32}\cdot C_{22} \\ C_{31}\cdot C_{11} & C_{31}\cdot C_{21}+C_{32}\cdot C_{22} & C^2_{31}+C^2_{32}+C^2_{33} \end{array} \right) .$

Označí-li se prvky matice ${\mathbf N}$ obecně ${{\mathbf N}}_{ij}$ , pak obecný algoritmus výpočtu pro matici $r\times r$ vypadá takto:

$C_{11}=\sqrt{N_{11}} , C_{1j}=\frac{N_{1j}}{C_{11}} (j=2,\dots ,r) .$

pro $i\ =\ 2,\dots ,r$ :

$C_{ii}={\left(N_{ii}-\sum^{i-1}_{k=1}{C^2_{ki}}\right)}^{\frac{1}{2}},$

$C_{ij}=\frac{\left(N_{ij}-\sum^{i-1}_{k=1}{C_{ki}\cdot C_{kj}}\right)}{C_{ii}}\$ , pro $j=i+1,\ \dots ,\ r$ .

Dále výpočet ekvivalentních rovnic:

Pro $i=1,...,r$ :

$y_i=\frac{{\left({{\mathbf A}}^T\cdot {\mathbf P}\cdot {\mathbf l}\right)}_i-\sum^{i-1}_{k=1}{C_{ki}\cdot y_k}}{C_{ii}} ,$

pro $i=r,...,1$ :

$x_i=\frac{y_i-\sum^r_{k=i+1}{C_{ik}\cdot x_k}}{C_{ii}} .$

Zápis rozkladu matice ${\mathbf C}$ v symbolice Gaussovy eliminace:

${\mathbf C}=\left( \begin{array}{cccc} \sqrt{\left[pa_1a_1\right]} & \frac{\left[pa_1a_2\right]}{\sqrt{\left[pa_1a_1\right]}} & \frac{\left[pa_1a_3\right]}{\sqrt{\left[pa_1a_1\right]}} & \dots \\ 0 & \sqrt{\left[pa_2a_2.1\right]} & \frac{\left[pa_2a_3.1\right]}{\sqrt{\left[pa_2a_2.1\right]}} & \dots \\ 0 & 0 & \sqrt{\left[pa_3a_3.2\right]} & \dots \\ \vdots & \vdots & \vdots & \ddots \end{array} \right).$

Pro odstranění nutnosti výpočtu odmocniny lze použít podobný rozklad:

${\mathbf N}={{\mathbf C}}^T\cdot {\mathbf D}\cdot {\mathbf C}=\left( \begin{array}{ccc} 1 & 0 & 0 \\ C_{21} & 1 & 0 \\ C_{31} & C_{32} & 1 \end{array} \right)\cdot \left( \begin{array}{ccc} D_1 & 0 & 0 \\ 0 & D_2 & 0 \\ 0 & 0 & D_3 \end{array} \right)\cdot \left( \begin{array}{ccc} 1 & C_{21} & C_{31} \\ 0 & 1 & C_{32} \\ 0 & 0 & 1 \end{array} \right) ,$

${\mathbf N}=\left( \begin{array}{ccc} D_1 & C_{21}\cdot D_1 & C_{31}\cdot D_1 \\ C_{21}\cdot D_1 & {C^2_{21}\cdot D_1+D}_2 & C_{31}\cdot C_{21}\cdot D_1+C_{32}\cdot D_2 \\ C_{31}\cdot D_1 & C_{31}\cdot C_{21}\cdot D_1+C_{32}\cdot D_2 & {C^2_{31}\cdot D_1+C^2_{32}\cdot D_2+D}_3 \end{array} \right) ,$

$D_j=N_{jj}-\sum^{j-1}_{k=1}{C^2_{jk}\cdot D_k}$ , pro $j=1,\dots ,3$ ,

$C_{ij}=\frac{1}{D_j}\cdot \left(N_{ij}-\sum^{j-1}_{k=1}{C_{ik}\cdot C_{jk}\cdot D_k}\right)$ , pro $i>j$ .

Pro rozklad pak platí:

${\mathbf N}={{\mathbf C}}^T\cdot {\mathbf D}\cdot {\mathbf C}={{\mathbf C}}^T\cdot {{\mathbf D}}^{\frac{1}{2}}\cdot {{\mathbf D}}^{\frac{1}{2}}\cdot {\mathbf C}={{\mathbf C}}^T\cdot {{\mathbf D}}^{\frac{1}{2}}\cdot {\left({\mathbf C}\cdot {{\mathbf D}}^{\frac{1}{2}}\right)}^T .$

Další výpočet probíhá obdobně, jako v předchozím případě, případně lze místo dvou soustav postupně řešit tři (více např. v [11]).

Matici ${\mathbf C}$ lze získat také QR rozkladem.

Metoda postupné iterace

Tato metoda je vhodná pro řešení velikých systémů rovnic, kde přesné metody již nelze použít buď pro rozsah rovnic, nebo pro potřebnou přesnost. Princip metody je takový, že řešený systém ${\mathbf N}\cdot {\mathbf x}{\mathbf =}{{\mathbf A}}^T\cdot {\mathbf P}\cdot {\mathbf l}$ převedeme na ekvivalentní:

${\mathbf x}{\mathbf =}{\mathbf B}\cdot {\mathbf x}{\mathbf +}{\mathbf g} ,$

kde

${\mathbf B}{\rm =-}{\left(diag\left({\mathbf N}\right)\right)}^{{\rm -1}}\cdot \left({\mathbf N}{\rm -}{\mathbf E}\right){\rm =}\left( \begin{array}{cccc} 0 & b_{{\rm 12}} & \cdots & b_{{\rm 1}k} \\ b_{{\rm 21}} & 0 & \cdots & b_{{\rm 2}k} \\ \vdots & \vdots & \vdots & \vdots \\ b_{k{\rm 1}} & b_{k{\rm 2}} & \cdots & 0 \end{array} \right),$

$g={\left(diag\left({\mathbf N}\right)\right)}^{{\rm -1}}\cdot {{\mathbf A}}^T\cdot {\mathbf P}\cdot {\mathbf l} .$

Takový systém nelze řešit přímo, ale zpravidla lze použít iterace. Při prvém přiblížení dosadíme libovolný vektor ${{\mathbf x}}_0$ do pravé strany rovnice a obdržíme „zlepšené“ hodnoty ${{\mathbf x}}_1$ . Při druhém přiblížení dosadíme ${{\mathbf x}}_1$ a dostaneme ${{\mathbf x}}_2$ atd. Tím obdržíme posloupnost vektorů:

${{\mathbf x}}_1={\mathbf B}\cdot {{\mathbf x}}_0+{\mathbf g} ,$

${{\mathbf x}}_2={\mathbf B}\cdot {{\mathbf x}}_1+{\mathbf g} ,$

$\vdots$

${{\mathbf x}}_m={\mathbf B}\cdot {{\mathbf x}}_m+{\mathbf g} .$

Po vzájemném postupném dosazení dostaneme:

${{\mathbf x}}_m{\mathbf =}{{\mathbf B}}^m\cdot {{\mathbf x}}_0{\mathbf +}\left({\mathbf E}{\mathbf +}{{\mathbf B}}^1+{{\mathbf B}}^{{\mathbf 2}}{\mathbf +\dots +}{{\mathbf B}}^{m-1}\right)\cdot {\mathbf g}$

Aby ${{\mathbf x}}_m$ konvergovalo ke správné hodnotě, je nutné, aby konvergovala posloupnost:

${\mathbf E}{\mathbf +}{{\mathbf B}}^1+{{\mathbf B}}^{{\mathbf 2}}{\mathbf +\dots +}{{\mathbf B}}^{m-1}.$

To nenastává vždy. Podmínka nutná a postačující požaduje, aby všechna vlastní čísla matice ${\mathbf B}$ byla menší než 1. Vlastní čísla ale většinou nebudeme počítat. Podmínkami postačujícími pro konvergenci dané posloupnosti je, aby kterákoliv norma matice ${\mathbf B}$ byla menší než 1, což prakticky znamená, aby členy na diagonále původní matice ${\mathbf N}$ svojí velikostí převyšovaly součet absolutních hodnot ostatních prvků v témže řádku matice.

Nedostatkem této metody je, že konverguje jen při splnění určitých podmínek a není proto použitelná pro všechny systémy. Výpočetní postup iterace může být dvojí:

a) Důsledné dosazování hodnot z předchozího kola (Jacobi),
b) dosazování vždy nejnověji vypočtených hodnot (Gauss).

U druhého způsobu je zpravidla rychlejší konvergence k limitám, avšak záleží také na pořadí počítaných neznámých.

Pseudoinverze

Pseudoinverzní matice a její výpočet má při řešení rozsáhlejších úloh v geodézii stále větší uplatnění, tento odstavec je věnován popisu vlastností a možnostem výpočtu. Vyplynou z toho praktické aplikace při řešení vyrovnání metodou nejmenších čtverců, které budou podrobněji popsány v následujících odstavcích.

Pojem pseudoinverze je rozšířením inverze na matice, které invertovat nelze, tj. i na matice singulární a obdélníkové. Pseudoinverze jsou různé, tj. splňují různé podmínky, pro uváděné výpočty bude dále pojmem pseudoinverze označována pouze Moore-Penroseova pseudoinverze, která byla popsána a propracována v literatuře [12], [13] a [14]. Pro další studium i souvisejících problémů lze doporučit také [15]. Zde budou uvedeny jen některé postupy a algoritmy, které ilustrují výpočet a jeho vlastnosti.

Definice a vlastnosti pseudoinverze

Pokud pseudoinverze ${{\mathbf A}}^{{\mathbf +}}$ matice ${\mathbf A}$ splňuje následující relace, jedná se o Moore-Penrose pseudoinverzi:

${\mathbf A}\cdot {{\mathbf A}}^{{\mathbf +}}\cdot {\mathbf A}{\mathbf =}{\mathbf A} ,$

${{\mathbf A}}^{{\mathbf +}}\cdot {\mathbf A}\cdot {{\mathbf A}}^{{\mathbf +}}{\mathbf =}{{\mathbf A}}^{{\mathbf +}} ,$

${\left({\mathbf A}\cdot {{\mathbf A}}^{{\mathbf +}}\right)}^T{\mathbf =}{\mathbf A}\cdot {{\mathbf A}}^{{\mathbf +}} ,$

${\left({{\mathbf A}}^{{\mathbf +}}\cdot {\mathbf A}\right)}^T{\mathbf =}{{\mathbf A}}^{{\mathbf +}}\cdot {\mathbf A} .$

Tato pseudoinverzní matice existuje vždy a je jedinečná. Pokud nesplňuje některé uvedené vlastnosti, může se jednat o některý typ generalizované inverze, kterou však nelze obecně použít pro všechny dále popsané účely.

Vybrané základní vlastnosti Moore-Penroseovi pseudoinverze:

Jestliže jsou prvky matice ${\mathbf A}$ reálná čísla, totéž platí pro prvky matice ${{\mathbf A}}^{{\mathbf +}}$ .
Pseudoinverze je komutativní s transpozicí.
Jestliže ${\mathbf A}$ je invertovatelná, platí: ${{\mathbf A}}^{{\mathbf -}1}{\mathbf =}{{\mathbf A}}^+$ .
Pseudoinverze pseudoinverzní matice je původní matice: ${\left({{\mathbf A}}^{{\mathbf +}}\right)}^{{\mathbf +}}{\mathbf =}{\mathbf A}{\mathbf \ }$ .
Pro transpozici inverzní matice platí: ${\left({{\mathbf A}}^T\right)}^+{\mathbf =}{\left({{\mathbf A}}^+\right)}^T$ .
Pro pseudoinverzi skalárního násobení nenulovým reálným číslem $\alpha$ platí: ${\left(\alpha \cdot {\mathbf A}\right)}^+={{\alpha }^{-1}\cdot {\mathbf A}}^+$ .
Pseudoinverze nulové matice je transponovaná nulová matice.
Pokud má ${\mathbf A}$ ortonormální sloupce nebo řádky, platí: ${{\mathbf A}}^{{\mathbf +}}{\mathbf =}{{\mathbf A}}^T$ .
Dále platí: ${{\mathbf A}}^{{\mathbf +}}{\mathbf =}{\left({{\mathbf A}}^T\cdot {\mathbf A}\right)}^{{\mathbf +}}\cdot {{\mathbf A}}^T$ , ${{\mathbf A}}^{{\mathbf +}}{\mathbf =}{{{\mathbf A}}^T\cdot \left({\mathbf A}{\cdot {\mathbf A}}^T\right)}^{{\mathbf +}}$ .

Metody výpočtu pseudoinverzní matice

Výpočet pseudoinverzní matice lze provést mnoha různými metodami, zde bude uvedena metoda rozkladu na submatice, dále aplikačně velmi jednoduchá iterační metoda Ben Israela a Cohena a jako poslední nejvyužívanější metoda založená na singulárním rozkladu.

Rozklad na součin matic

Pseudoinverzní matici k matici ${\mathbf N}$ , která je součinem dvou matic daných libovolných rozkladem

${\mathbf N}={\mathbf B}\cdot {\mathbf C} ,$

lze vypočítat podle vztahu:

${{\mathbf N}}^+={{\mathbf C}}^T\cdot {\left({\mathbf C}\cdot {{\mathbf C}}^T\right)}^{-1}\cdot {\left({{\mathbf B}}^T\cdot {\mathbf B}\right)}^{-1}{{\mathbf B}}^T,$

Obdélníkové matice ${\mathbf C}$ , ${\mathbf B}$ dostaneme libovolným rozkladem matice ${\mathbf N}{\mathbf =}{\mathbf B}\cdot {\mathbf C}$ tak, aby platilo $h_B=h_N$ ( $h$ je hodnost příslušných matic).

Iterační metoda Ben-Israela a Cohena

Metoda je velmi jednoduchá, má zajištěnou konvergenci při dodržení triviálního pravidla. Nevýhodou je vyšší počet iterací. Jestliže ${\mathbf A}$ je obdélníková matice, určí se první přiblížení ${{\mathbf Y}}_0$ :

${{\mathbf Y}}_0=\alpha \cdot {{\mathbf A}}^T ,$

kde $\alpha$ je reálné číslo. Dále se opakovaně určuje:

${{\mathbf Y}}_{k+1}={{\mathbf Y}}_k\cdot \left(2\cdot {\mathbf E}-{\mathbf A}\cdot {{\mathbf Y}}_k\right)=2\cdot {{\mathbf Y}}_k-{{\mathbf Y}}_k\cdot {\mathbf A}\cdot {{\mathbf Y}}_k .$

${{\mathbf A}}^+$ je limitou posloupnosti $\left\{{{\mathbf Y}}_k\right\}$ , jestliže reálné $\alpha$ splňuje podmínku:

$0<\alpha <\frac{2}{{\lambda }^2_{max}\ }=\frac{2}{{\lambda }_1\left({{\mathbf A}}^T\cdot {\mathbf A}\right)} ,$

kde ${\lambda }^2_{max}$ je kvadrát maximálního vlastního čísla matice ${\mathbf A}$ , ${\lambda }_1\left({{\mathbf A}}^T\cdot {\mathbf A}\right)$ je maximální vlastní číslo symetrické matice ${{\mathbf A}}^T\cdot {\mathbf A}$ . Tato podmínka je obtížněji vyčíslitelná nežli následující „nižší“ horní hranice intervalu:

$0<\alpha <\frac{2}{k\ } ,$

kde

$k={\mathop{\max }_{i=1\dots r} \sum \left|{{\mathbf b}}_i\right|\ }={\mathop{{\rm ma}{\rm x}}_{i=1\dots r} \sum^s_{j=1}{\left|B_{ij}\right|}\ } ,$

kde

${\mathbf B}={{\mathbf A}}^T\cdot {\mathbf A}{\mathbf =}\left( \begin{array}{cccc} {{\mathbf b}}_1 & {{\mathbf b}}_2 & {\mathbf \dots } & {{\mathbf b}}_n \end{array} \right) ,$

kde $r,s$ je počet řádků a sloupců matice ${\mathbf B}$ (platí $r=s$ ).

Pokud podmínka není splněna, posloupnost nemusí konvergovat. Metoda konverguje kvadraticky v „blízkosti“ výsledné pseudoinverzní matice, než se však dostane do této blízkosti, je to obvykle velmi výpočetně náročné. Metodu lze použít pro zpřesnění pseudoinverzní matice určené jinou metodou. Iterace se ukončí při dosažení požadované přesnosti, např. pokud je splněno kritérium:

${\max \left(\left|{{\mathbf Y}}_i\cdot {\mathbf A}\cdot {{\mathbf Y}}_i-{{\mathbf Y}}_i\right|\right)\ }<\ \varepsilon ,$

kde $\varepsilon$ je zvolené (malé kladné) číslo. $\alpha$ je vhodné pro zajištění konvergence volit ve tvaru $2/k-\delta$ , kde $\delta$ je malé kladné číslo. Větší hodnota $\alpha$ zajišťuje mírně rychlejší konvergenci. Blíže k metodě v [16].

Metoda využívající singulární rozklad

Základem tohoto výpočtu pseudoinverzní matice je výpočet vlastních vektorů a vlastních čísel matice (jak bude uvedeno dále), nejpoužívanějším postupem je vzhledem k jeho stabilitě je Golub-Reinschův postup, který je založen na transformaci do bidiagonálního tvaru pomocí Householderovy transformace, dále vytvoření tridiagonální symetrické matice a následně aplikováním QR rozkladu, blíže např. viz [17]. Zde bude uveden jednodušší a méně kvalitní postup.

Singulární rozklad je rozklad matice ${\mathbf A}$ na součin tří matic se speciálním tvarem a vlastnostmi:

${\mathbf A}{\mathbf =}{\mathbf U}\cdot {\mathbf S}\cdot {{\mathbf V}}^T .$

${\mathbf U}$ a ${\mathbf V}$ jsou ortonormální matice (matice levých, resp. pravých vlastních vektorů matice ${\mathbf A}$ ) a ${\mathbf S}$ je diagonální matice tvořená vlastními čísly matice ${\mathbf A}$ . Pro pseudoinverzní matici pak platí:

${{\mathbf A}}^{{\mathbf +}}{\mathbf =}{\left({\mathbf U}\cdot {\mathbf S}\cdot {{\mathbf V}}^T\right)}^{{\mathbf +}}{\mathbf =}{\left({{\mathbf V}}^T\right)}^{{\mathbf +}}\cdot {{\mathbf S}}^{{\mathbf +}}\cdot {{\mathbf U}}^{{\mathbf +}} .$

Protože pseudoinverze ortogonální matice je rovná její transpozici, platí:

${{\mathbf A}}^+{\mathbf =}{\mathbf V}\cdot {{\mathbf S}}^{{\mathbf +}}\cdot {{\mathbf U}}^T .$

Pseudoinverzní matice diagonální matice se vypočítá také velmi jednoduše, výsledkem je diagonální matice, která se z nenulových prvků původní matice určí takto $S^+_{ii}={1/S}_{ii}$ . Nulové prvky zůstanou nulové.

Samotný výpočet pseudoinverzní matice je tedy triviální, nyní k principu výpočtu singulárního rozkladu. Zde se určují vlastní čísla matice ${\mathbf A}$ , využívá se QR rozkladu aplikovaného iteračně na symetrickou matici ${\mathbf N}={{\mathbf A}}^T\cdot {\mathbf A}$ a vlastnosti, že pro vlastní čísla matice ${\mathbf N}$ a matice ${\mathbf A}$ platí $\lambda {\left({\mathbf A}\right)}^2_i=\lambda {\left({\mathbf N}\right)}_i$ . Výpočtem QR rozkladu ${\mathbf N}={\mathbf Q}\cdot {\mathbf R}$ se získají matice ${\mathbf Q}$ (ortonormální) a ${\mathbf R}$ (horní trojúhelníková). Dále se aplikuje opačná ortonormální transformace ve tvaru

${{\mathbf T}}_1={{\mathbf R}}_1\cdot {{\mathbf Q}}_1 ,$

jejímž výsledkem je matice ${\mathbf T}$ , která je podobná ${\mathbf N}$ se stejnými vlastními čísly (ortonormální transformací se vlastní čísla nemění). Iteračně se opakuje:

${{\mathbf Q}}_{k+1},{{\mathbf R}}_{k+1}\leftarrow {{\mathbf T}}_k ,$

${{\mathbf T}}_{k+1}={{\mathbf R}}_{k+1}\cdot {{\mathbf Q}}_{k+1} .$

Posloupnost $\left\{{{\mathbf T}}_k\right\}$ konverguje k trojúhelníkové matici s vlastními čísly na diagonále, jsou to vlastní čísla symetrické matice ${\mathbf N}$ , nikoli matice ${\mathbf A}$ a pro sestavení diagonální matice ${\mathbf S}$ je třeba je odmocnit. Iterace se zastaví po $n$ cyklech, pokud jsou ideálně nulové prvky trojúhelníkové matice ${\mathbf R}$ dostatečně malé.

Matice ${\mathbf V}$ se získá součinem ortonormálních transformačních matic ${\mathbf Q}$ :

${\mathbf V}=\prod^n_{k=1}{{{\mathbf Q}}_k}={{\mathbf Q}}_1\cdot {{\mathbf Q}}_2\cdot \dots \cdot {{\mathbf Q}}_n .$

Matice ${\mathbf U}$ se pak získá:

${\mathbf U}{\mathbf =}{\mathbf A}\cdot {\mathbf V}\cdot {{\mathbf S}}^{{\mathbf +}}.$

V případě, že vstupní matice pro výpočet pseudoinverze již je symetrická matice normálních rovnic, lze s výhodou vynechat první krok řešení, tj. tvorbu symetrické matice ${\mathbf N}={{\mathbf A}}^T\cdot {\mathbf A}$ , toto proběhlo již při tvorbě normálních rovnic. Při přítomnosti matice vah ${\mathbf P}$ , která je vždy symetrická, se na symetričnosti pseudoinvertované matice nic nemění.

Využití pseudoinverze při řešení vyrovnání MNČ

Využití pseudoinverze při řešení úloh metodou nejmenších čtverců je několikeré, některé možnosti budou dále popsány.

Výpočet inverzní matice pseudoinverzí

Jak již bylo uvedeno, pokud existuje k dané matici matice inverzní, je totožná s maticí pseudoinverzní a lze tedy využít algoritmy pro výpočet pseudoinverzní matice, která na rozdíl od inverzní matice existuje vždy a výpočetní algoritmy zvládnou vyřešit i špatně podmíněné matice, u kterých výpočet inverze selhává. Zároveň lze s výhodou využít vlastností singulárního rozkladu, kdy se vlastní čísla matice v matici ${\mathbf S}$ posoudí z hlediska jejich velikosti a velikosti odhadu zaokrouhlovacích chyb těchto vlastních čísel, vlastní čísla menší než takto určená hranice se považují za nedůvěryhodná a nahradí se nulami. Tímto postupem lze potlačit nestabilitu výpočtu a v mnohých hraničních případech získat výsledky i tam, kde to pomocí inverze nelze. Podle [18] lze hranici $tol$ stanovit následujícím způsobem:

$tol=\varepsilon \cdot {\rm max}(m,n)\cdot {\max \left(\lambda \right)\ },$

kde $\varepsilon$ je počítačové (strojové) epsilon, tj relativní chyba vznikající zaokrouhlováním při aritmetických operacích, pro ilustraci pro typ double v jazyce C++ bývá uváděna hodnota $\varepsilon =1\cdot {10}^{-16}$ , $m$ a $n$ jsou rozměry matice a ${\max \left(\lambda \right)\ }$ značí maximální vlastní číslo matice.

Přímé řešení rovnic oprav pseudoinverzí

Jednou z vlastností Moore-Penroseovy pseudoinverze velmi dobře využitelnou pro řešení vyrovnání metodou nejmenších čtverců je minimalizace normy přeurčených systémů lineárních rovnic ve tvaru ${\mathbf A}\cdot {\mathbf x}={\mathbf b}$ ve smyslu euklidovské normy tak, že pro všechna řešení ${\mathbf x}$ platí:

${\left\|{\mathbf A}\cdot {\mathbf x}-{\mathbf b}\right\|}_{{\rm 2}}\ge {\left\|{\mathbf A}\cdot {\mathbf z}-{\mathbf b}\right\|}_{{\rm 2}}{\rm \ } ,$

kde

${\mathbf z}={{\mathbf A}}^+\cdot {\mathbf b} ,$

tj. splňuje podmínku metody nejmenších čtverců

${\left\|{\mathbf A}\cdot {\mathbf x}-{\mathbf b}\right\|}_{{\rm 2}}=min.$

Při využití vah pro nezávislá měření při vyrovnání se provede rozklad matice normálních rovnic ve smyslu

${\mathbf N}={{\mathbf A}}^T\cdot {\mathbf P}\cdot {\mathbf A}={{\mathbf A}}^T\cdot {{\mathbf P}}^{\frac{{\mathbf 1}}{{\mathbf 2}}}\cdot {{\mathbf P}}^{\frac{{\mathbf 1}}{{\mathbf 2}}}\cdot {\mathbf A}$

a řeší se rovnice ve tvaru

${\left\|{{\mathbf P}}^{\frac{{\mathbf 1}}{{\mathbf 2}}}\cdot {\mathbf (}{\mathbf A}\cdot {\mathbf x}-{\mathbf b}{\mathbf )}\right\|}_{{\rm 2}}=min.$

Řešení inverze singulární matice normálních rovnic

V případě, že matice ${\mathbf N}$ je singulární, nemá systém normálních rovnic jednoznačné řešení. V podstatě to znamená, že některé řádky matice ${\mathbf N}$ jsou lineární kombinací dalších řádků. Tyto závislé řádky by měly být vyškrtnuty a dostane se tím méně rovnic, než je hledaných neznámých. O řešení takovýchto systémů je známo, že je jich nekonečně mnoho a lze je vyjádřit parametricky. Tato situace nastává např. u tzv. volných výškových a polohových sítí, kde jsou za neznámé voleny souřadnice všech bodů.

Řešení se provádí pomocí pseudoinverzní matice ${{\mathbf N}}^+$ . Řešení systému normálních rovnic splní podmínku ${{\mathbf x}}^T\cdot {\mathbf x}=min$ a zapíšeme jej ve tvaru:

${\mathbf x}{\mathbf =}{{\mathbf N}}^+{\cdot {\mathbf A}}^T\cdot {\mathbf P}\cdot {\mathbf l}.$

Pro výpočet matice ${{\mathbf Q}}_x$ použijeme zákon hromadění vah s tím, že ${\mathbf N}$ je symetrická matice:

${{\mathbf Q}}_x{\mathbf =}{{\mathbf N}}^+{\cdot {\mathbf A}}^T\cdot {\mathbf P}\cdot {{\mathbf P}}^{-1}\cdot {\mathbf P}\cdot {\mathbf A}\cdot {\left({{\mathbf N}}^+\right)}^T,$

${{\mathbf Q}}_x{\mathbf =}{{\mathbf N}}^+\cdot {\mathbf N}\cdot {\left({{\mathbf N}}^+\right)}^T{\mathbf =}{\mathbf N}\cdot {\left({\mathbf N}\cdot {\mathbf N}\right)}^{-1}\cdot {\mathbf N}\cdot {\left({\mathbf N}\cdot {\mathbf N}\right)}^{-1}\cdot {\mathbf N}.$

Je to vzorec odlišný, než bychom čekali z analogie s regulární maticí. Pouze při ${\mathbf N}$ regulární platí: ${{\mathbf Q}}_x{\mathbf =}{{\mathbf N}}^{-1}$ . V opačném případě je ${{\mathbf Q}}_x$ maticí sice čtvercovou a symetrickou, ale singulární, vyhovující podmínce $Sp({{\mathbf Q}}_x)\ =\ min.$ , kde $Sp({{\mathbf Q}}_x)$ je symbol pro stopu matice ${{\mathbf Q}}_x$ .

^{[1],
[15],
[17]} Björck, A.: Numerical Methods for Least Squares Problems. Society for Industrial and Applien Mathematics, Philadelphia, 1996. ISBN 0-89871-360-9.

^{[2],
[10],
[11]} Golub, G. H. - Van Loan, Ch. F.: Matrix computations. John Hopkins University Press, London, 1996. (3rd ed.) ISBN 0-8018-5414-8.

^[3] Press, W. H. - Flannery, B. P. - Teukolsky, S. A. - Vetterling, W. T.: Numerical Recipes. The Art of Scientific Computing. Cambridge University Press, 2007. ISBN 0-5218-8068-8.

^[4],
[18] Press, W. H. - Flannery, B. P. - Teukolsky, S. A. - Vetterling, W. T.: Numerical Recipes in Pascal. Cambridge University Press, 1989. ISBN 0-5213-7516-9.

^{[5],
[6],
[7],
[9]} Bubeník, F. - Pultar, M.: Matematické vzorce a metody. Vydavatelství ČVUT, Praha, 1994.

^[8] http://en.wikipedia.org/wiki/QR_decomposition; 6.1.2011.

^[12] Moore, E. H.: On the reciprocal of the general algebraic matrix. Bulletin of the American Mathematical Society 26: 394-395, 1920.

^[13] Bjerhammar, A.: Application of calculus of matrices to method of least squares; with special references to geodetic calculations. Trans. Roy. Inst. Tech. Stockholm 49, 1951.

^[14] Penrose, R.: A generalized inverse for matrices. In: Proceedings of the Cambridge Philosophical Society 51: 406-413, 1955.

^[16] Ben Israel, A. - Cohen, D.: On Iterative Computation of Generalized Inverses and Associated Projections. SIAM Numer. Anal., vol. 3, no. 3, USA, 1966.

« 18. Metody robustního odhadu
» 20. Kritéria pro testování opakovaných měření