C++ Bc. 33

Choleskyho rozklad

Napište funkci, která počítá Choleskyho rozklad pozitivně semidefinitní matice.

Matice $\mathbf {A}$ je pozitivně definitní, pokud pro každý nenulový vektor $\mathbf {x}$ platí $\mathbf {x} ^{t}\mathbf {A} \mathbf {x} >0.$ Pro každou pozitivně definitní matici existuje jednoznačný symetrický rozklad (Choleskyho rozklad) $\mathbf {A} =\mathbf {L} \mathbf {L} ^{t},$ kde $\mathbf {L}$ je dolní trojúhelníková matice. Například

${\begin{pmatrix}225&210&180&135&75\\210&365&311&230&122\\180&311&365&266&134\\135&230&266&230&110\\75&122&134&110&55\\\end{pmatrix}}={\begin{pmatrix}15&0&0&0&0\\14&13&0&0&0\\12&11&10&0&0\\9&8&7&6&0\\5&4&3&2&1\end{pmatrix}}{\begin{pmatrix}15&14&12&9&5\\0&13&11&8&4\\0&0&10&7&3\\0&0&0&6&2\\0&0&0&0&1\end{pmatrix}}$

První sloupec Choleskyho rozkladu můžeme vypočítat jako

$\mathbf {A} =\mathbf {A} ^{(1)}={\begin{pmatrix}a_{11}&\mathbf {a} _{*1}^{t}\\\mathbf {a} _{*1}&\mathbf {A} ^{(2)}\end{pmatrix}}\longrightarrow {\begin{pmatrix}{\sqrt {a_{11}}}&\mathbf {0} \\\mathbf {b} _{*1}&\mathbf {B} ^{(2)}\end{pmatrix}},$

kde $\mathbf {b} _{*1}={1 \over {\sqrt {a_{11}}}}\mathbf {a} _{*1}$ (tj. prvky pod diagonálou vydělíme odmocninou z diagonálnho prvku) a $\mathbf {B} ^{(2)}=\mathbf {A} ^{(2)}-\mathbf {b} _{*1}\mathbf {b} _{*1}^{t}$. V našem příkladu tedy

$\mathbf {B} ^{(2)}={\begin{pmatrix}169&143&104&52\\143&221&158&74\\104&158&149&65\\52&74&65&30\end{pmatrix}}={\begin{pmatrix}365&311&230&122\\311&365&266&134\\230&266&230&110\\122&134&110&55\end{pmatrix}}-{\begin{pmatrix}14\times 14&14\times 12&14\times 9&14\times 5\\12\times 14&12\times 12&12\times 9&12\times 5\\9\times 14&9\times 12&9\times 9&9\times 5\\5\times 14&5\times 12&5\times 9&5\times 5\end{pmatrix}}.$

Submatice $\mathbf {B} ^{(2)}$ je také pozitivně definitní a stejným způsobem můžeme vypočítat druhý sloupec choleskyho rozkladu a obdobně i zbývající sloupce.

Explicitní vzorce pro výpočet koeficentů matice $\mathbf {L}$ Choleskyho rozkladu jsou

$l_{i,j}={\frac {1}{l_{j,j}}}\left(a_{i,j}-\sum _{k=1}^{j-1}l_{i,k}l_{j,k}\right),\qquad {\mbox{pro }}i>j.$

$l_{i,i}={\sqrt {a_{i,i}-\sum _{k=1}^{i-1}l_{i,k}^{2}}}.$

[ Zpět | C++ | Další ]