Problém nejmenších čtverců

Máme zadány matice ${\bf A}$ a ${\bf B}$ a hledáme takovou matici ${\bf X}$, že pro každou matici ${\bf X'}$ platí

$$\vert\vert {\bf B} - {\bf A}{\bf X'} \vert\vert \geq \vert\vert {\bf B} - {\bf A}{\bf X} \vert\vert$$ (4.25)

Pokud neexistuje inverzní matice k matici $A$ použijeme k řešení jednu z následujících metod.

Moore-Penroseova pseudoinverse
Řešení spočteme

$${\bf X} = {\bf A}^{+}{\bf B},$$ (4.26)

kde ${\bf A^{+}}$ je pseudoinverze matice $A$. Pseudoinverzí matice ${\bf A}$ nazýváme takovou matici ${\bf A^{+}}$, pro kterou platí

$$\begin{array}{lcclc} 1.& {\bf A}{\bf A}^{+}{\bf A} = {\bf A... ....& ({\bf A}^{+}{\bf A})^T = {\bf A}^{+}{\bf A}. \\ \end{array}$$

Existuje-li k matici A inversní matice, je rovna její pseudoinversi. Pokud je matice $({\bf A}^T{\bf A})$ regulární, pak matice $({\bf A}^T{\bf A})^{-1}{\bf A}^T$ splňuje podmínky 1-4. Tato matice bývá označována jako Moore-Penroseova pseudoinverse.

SVD rozklad
Pro každou matici ${\bf A}\in\mathbb {R}^{m\times n}$ existují matice ${\bf U}$, ${\bf V}$ a ${\bf\Sigma}$ tak, že

$${\bf A} = {\bf U}{\bf\Sigma}{\bf V}^T$$ (4.27)

a ${\bf U}^T{\bf U}={\bf V}^T{\bf V}={\bf V}{\bf V}^T={\bf I}^n$ a ${\bf\Sigma}$ je diagonální matice. Prvky diagonály $(\sigma_{1}, \cdots ,\sigma_{n})$ se nazývají singulární čísla matice ${\bf A}$. Pokud je hodnost matice ${\bf A}$ rovna $r$, pak $\sigma_{r+1}, \cdots, \sigma_{n}$ jsou rovny nule.

Dosadíme-li za ${\bf A}$ její SVD rozklad dostáváme

$${\bf U}{\bf\Sigma}{\bf V}^T{\bf X} = {\bf B}$$ (4.28)

a tedy

$${\bf X} = {\bf V}{\bf\Sigma}^{+}{\bf U}^T{\bf B}$$ (4.29)

kde ${\bf\Sigma}^{+}=diag(\sigma_i^{+})$ a $\sigma_i^{+}=1/\sigma_i$, pro $\sigma_{i}\neq 0$ a $\sigma_i^{+}=0$, pro $\sigma_{i}=0$.

G.H.Golub a C.Reinsch [10, str. 134-151] navrhli algoritmus na výpočet SVD rozkladu. Nejprve upravíme matici na bidiagonální tvar a po té se provádí vlastní rozklad. Redukce na bidiagonální tvar se provádí pomocí Householderových transformací [8, str. 52-57]

$${\bf P}^{(k)} = {\bf I} - 2\vec{x}^{(k)}\vec{x}^{(k)T}$$ (4.30)

$${\bf Q}^{(k)} = {\bf I} - 2\vec{y}^{(k)}\vec{y}^{(k)T}$$ (4.31)

tak, aby ${\bf P}^{(n)}\cdots{\bf P}^{(1)}{\bf A}{\bf Q}^{(1)}\cdots{\bf Q}^{(n-2)} = {\bf J}^{0}$ byla v horním bidiagonálním tvaru. Matice ${\bf J}^0$ má stejná singulární čísla jako matice ${\bf A}$. Tedy je-li

$${\bf J}^0 = {\bf G}{\bf\Sigma}{\bf H}^T$$ (4.32)

pak

$${\bf A} = {\bf P}{\bf G}{\bf\Sigma}{\bf H}^T{\bf Q}^T$$ (4.33)

tak, že ${\bf U}={\bf P}{\bf Q},$ ${\bf V}={\bf Q}{\bf H}$ a ${\bf P}={\bf P}^1 \cdots {\bf P}^n$ a ${\bf Q}={\bf Q}^1 \cdots {\bf Q}^{n-2}$.

SVD rozklad bidiaogonální matice probíhá jako iterativní diagonalizace ${\bf J}^0 \rightarrow {\bf J}^1 \cdots \rightarrow {\bf\Sigma}$, kde ${\bf J}^{i+1} = {\bf S}^{iT}{\bf J}^i{\bf T}^i$ a ${\bf S}^i$ a ${\bf T}^i$ jsou transformační matice složené s Givensových rotací [8, str. 49-51].

QR rozklad
B.Businger a G.H.Golub [10, str. 111-118] navrhli řešení problému nejmenších čtverců pomocí Householderových transformací.

Máme matici ${\bf A}\in\mathbb {R}^{m\times n}$, kde $m\leq n$ a hodnost matice ${\bf A}$ je $n$, a matici ${\bf B}$. Hledáme ${\bf X}$ tak, aby pro každou matici ${\bf X'}$ platilo

$$\parallel{\bf B}-{\bf A}{\bf X'}\parallel \geq \parallel{\bf B}-{\bf A}{\bf X}\parallel$$ (4.34)

Platí, že $\parallel{\bf B}-{\bf A}{\bf X}\parallel = \parallel{\bf Q}{\bf B}-{\bf Q}{\bf A}{\bf X}\parallel$, pokud ${\bf Q}^T{\bf Q}=I.$

Zvolíme ${\bf Q}$ tak, aby

$${\bf Q}{\bf A} = {\bf R} = \left( \begin{array}{c} {\bf\tilde{R}} \\ {\bf0} \end{array}\right) ,$$ (4.35)

kde ${\bf\tilde{R}}$ je horní trojúhelníková matice. Pak ${\bf X}={\bf R}^{-1}({\bf Q}{\bf B})$.

Dekompozici lze realizovat pomocí Householderových transformací

$$\begin{array}{cr} {\bf P}^k={\bf I}-\beta_k\vec{u}^k\vec{u}^{kT} \qquad k=1, \cdots ,n \end{array}$$ (4.36)

Vektory $\vec{u}^k$ a parametry $\beta_k$ určíme tak, aby k-tá transformace vynulovala v k-tém sloupci všechny prvky pod diagonálou.

$$\sigma_k = \left( \sum_{i=k}^m ({a}_{ik}^k)^2 \right)$$ (4.37)

$$\beta_k = \left( \sigma_k(\sigma_k+\vert a_{kk}^k\vert)\right)$$ (4.38)

$$u_i^k = 0, \quad i<k$$ (4.39)

$$u_k^k = sgn(a_{kk}^k)\left( \sigma_k+\vert a_{kk}^k\vert \right)$$ (4.40)

$$u_i^k = a_{ik}^k, \quad i>k$$ (4.41)

Praktický výpočet probíhá iterativně. Po určení prvního řešení ${\bf\tilde{X}}$, položíme ${\bf X}={\bf\tilde{X}}+{\bf E}$. Potom $\parallel{\bf B}-{\bf A}{\bf X}\parallel = \parallel{\bf R}-{\bf A}{\bf E}\parallel$, kde ${\bf R}={\bf B}-{\bf A}{\bf\tilde{X}}$. Výpočet korekce ${\bf E}$ je snadný, neboť rozklad matice ${\bf A}$ už máme.