科列斯基分解

若欲對埃爾米特矩陣直接求逆，可以通過科列斯基分解，類似求解線性方程組一般求出逆矩陣，這需要${\displaystyle n^{3}}$次運算（${\displaystyle {\frac {1}{2}}n^{3}}$次乘法運算）。 該方法即便要求逐步計算也非常有效率。

非埃爾米特矩陣${\displaystyle \mathbf {B} }$也可以通過下例性質求逆，因為其中${\displaystyle \mathbf {BB} ^{*}}$總是埃爾米特矩陣︰

${\displaystyle \mathbf {B} ^{-1}=\mathbf {B} ^{*}(\mathbf {BB} ^{*})^{-1}}$

用於計算矩陣${\displaystyle \mathbf {L} }$的科列斯基算法，是以高斯消元法為基礎而調整得來的。

遞歸算法由${\displaystyle i:=1}$開始，令

${\displaystyle \mathbf {A} ^{(1)}:=\mathbf {A} }$

在步骤${\displaystyle i}$，矩陣${\displaystyle A^{(i)}}$的形式如下：

${\displaystyle \mathbf {A} ^{(i)}={\begin{pmatrix}\mathbf {I} _{i-1}&0&0\\0&a_{i,i}&\mathbf {b} _{i}^{*}\\0&\mathbf {b} _{i}&\mathbf {B} ^{(i)}\end{pmatrix}},}$

其中，${\displaystyle \mathbf {I} _{i-1}}$代表${\displaystyle i-1}$維的單位矩陣。

此時定義矩陣${\displaystyle \mathbf {L} _{i}}$為

${\displaystyle \mathbf {L} _{i}:={\begin{pmatrix}\mathbf {I} _{i-1}&0&0\\0&{\sqrt {a_{i,i}}}&0\\0&{\frac {1}{\sqrt {a_{i,i}}}}\mathbf {b} _{i}&\mathbf {I} _{n-i}\end{pmatrix}},}$

隨即${\displaystyle \mathbf {A} ^{(i)}}$可以被改寫成

${\displaystyle \mathbf {A} ^{(i)}=\mathbf {L} _{i}\mathbf {A} ^{(i+1)}\mathbf {L} _{i}^{*}}$

其中

${\displaystyle \mathbf {A} ^{(i+1)}={\begin{pmatrix}\mathbf {I} _{i-1}&0&0\\0&1&0\\0&0&\mathbf {B} ^{(i)}-{\frac {1}{a_{i,i}}}\mathbf {b} _{i}\mathbf {b} _{i}^{*}\end{pmatrix}}.}$

注意︰在此${\displaystyle \mathbf {b} _{i}\mathbf {b} _{i}^{*}}$是一個外積。

重複此步驟直到${\displaystyle i}$從${\displaystyle 1}$到${\displaystyle n}$。${\displaystyle n}$步之後，我們得到${\displaystyle A^{(n+1)}=\mathbf {I} }$。因此，所求的下三角矩陣${\displaystyle L}$為

${\displaystyle \mathbf {L}$ :=\mathbf {L} _{1}\mathbf {L} _{2}\dots \mathbf {L} _{n}.}

科列斯基-巴那齊耶維茨演算法以一個 5×5 矩陣執行的讀取順序（白色）及寫入順序（黄色）

寫出等式${\displaystyle \mathbf {A} =\mathbf {LL} ^{*}}$，

${\displaystyle {\begin{aligned}{\mathbf {A=LL^{T}} }&={\begin{pmatrix}\mathbf {L} _{11}&0&0\\\mathbf {L} _{21}&\mathbf {L} _{22}&0\\\mathbf {L} _{31}&\mathbf {L} _{32}&\mathbf {L} _{33}\\\end{pmatrix}}{\begin{pmatrix}\mathbf {L} _{11}&\mathbf {L} _{21}&\mathbf {L} _{31}\\0&\mathbf {L} _{22}&\mathbf {L} _{32}\\0&0&\mathbf {L} _{33}\end{pmatrix}}\\&={\begin{pmatrix}\mathbf {L} _{11}^{2}&&({\text{symmetric}})\\\mathbf {L} _{21}\mathbf {L} _{11}&\mathbf {L} _{21}^{2}+\mathbf {L} _{22}^{2}&\\\mathbf {L} _{31}\mathbf {L} _{11}&\mathbf {L} _{31}\mathbf {L} _{21}+\mathbf {L} _{32}\mathbf {L} _{22}&\mathbf {L} _{31}^{2}+\mathbf {L} _{32}^{2}+\mathbf {L} _{33}^{2}\end{pmatrix}}\end{aligned}}}$

則得到矩陣${\displaystyle \mathbf {L} }$的元素的計算公式如下︰

${\displaystyle \mathbf {L} _{j,j}={\sqrt {A_{j,j}-\sum _{k=1}^{j-1}\mathbf {L} _{j,k}^{2}}}}$

${\displaystyle \mathbf {L} _{i,j}={\frac {1}{\mathbf {L} _{j,j}}}\left(A_{i,j}-\sum _{k=1}^{j-1}\mathbf {L} _{i,k}\mathbf {L} _{j,k}\right),\qquad {\text{for }}i>j}$

只要${\displaystyle \mathbf {A} }$是實數正定矩陣，則平方根下的表達式恆為正。

對於複埃爾米特矩陣，則適用如下公式：

${\displaystyle \mathbf {L} _{j,j}={\sqrt {A_{j,j}-\sum _{k=1}^{j-1}\mathbf {L} _{j,k}\mathbf {L} _{j,k}^{*}}}.}$

${\displaystyle \mathbf {L} _{i,j}={\frac {1}{\mathbf {L} _{j,j}}}\left(\mathbf {A} _{i,j}-\sum _{k=1}^{j-1}\mathbf {L} _{i,k}\mathbf {L} _{j,k}^{*}\right),\qquad {\text{for }}i>j.}$

因此，要計算${\displaystyle \mathbf {L} _{i,j}(i\neq j)}$只需利用其的左、上方元素的值。計算通常是以以下其中一種順序進行的。

• 科列斯基-巴那齊耶維茨演算法從矩陣L的左上角開始，依行進行計算。
• 科列斯基-克勞特演算法從矩陣L的左上角開始，依列進行計算。

若有需要，整個矩陣可以逐個元素計算得出，無論使用何種順序讀取。

假設我們要求解一個良置的線性方程組。採用了LU分解的算法，除非進行某種繞軸旋轉，否則是不穩定的；如果算法進行了繞軸旋轉，則其誤差取決於矩陣的增長因子，這個數通常（但非總是）很小。

現在，假設算法適用科列斯基分解。如上所述，算法的效率將會是原來的兩倍。此外，無必要進行繞軸旋轉且誤差總是很小。具體而言，若要求解${\displaystyle \mathbf {Ax} =\mathbf {b} }$，${\displaystyle \mathbf {y} }$表示已計算出的解，然後能解出干擾方程組${\displaystyle (\mathbf {A} +\mathbf {E} )\mathbf {y} =\mathbf {b} }$，其中

${\displaystyle \|\mathbf {E} \|_{2}\leq c_{n}\varepsilon \|\mathbf {A} \|_{2}.}$

在這裡，${\displaystyle \|\quad \|_{2}}$是指矩陣2-範數，而${\displaystyle c_{n}}$是一個取決於${\displaystyle n}$的小常數，${\displaystyle \varepsilon }$表示單位捨入。

值得一提的是，科列斯基分解與平方根的使用有關。如果被分解的矩陣是正定的，那麽只要運算精確，矩陣中帶有平方根的元素的平方根下的數字永遠是正數。不幸的是，由於存在捨入誤差，這些數字可能為負數，並使算法擱淺。然而，此種情況僅當矩陣為病態時才會出現。一種可解決此問題的方法，是增添一個對角修正矩陣至待分解矩陣，以增加矩陣的正定性。[6] 此法雖或將減少分解的準確性，但有在某些情況卻頗有作為；譬如，執行應用於最優化的牛頓法時，若初期值距最優值較遠，則此時引入對角矩陣可以提高演算法的穩定性。

科列斯基分解的另一種形式——LDL分解的計算方式如下所示，

${\displaystyle {\begin{aligned}{\mathbf {A=LDL} ^{\mathrm {T} }}&={\begin{pmatrix}1&0&0\\L_{21}&1&0\\L_{31}&L_{32}&1\\\end{pmatrix}}{\begin{pmatrix}D_{1}&0&0\\0&D_{2}&0\\0&0&D_{3}\\\end{pmatrix}}{\begin{pmatrix}1&L_{21}&L_{31}\\0&1&L_{32}\\0&0&1\\\end{pmatrix}}\\&={\begin{pmatrix}D_{1}&&(\mathrm {symmetric} )\\L_{21}D_{1}&L_{21}^{2}D_{1}+D_{2}&\\L_{31}D_{1}&L_{31}L_{21}D_{1}+L_{32}D_{2}&L_{31}^{2}D_{1}+L_{32}^{2}D_{2}+D_{3}.\end{pmatrix}}\end{aligned}}}$

如果 ${\displaystyle \mathbf {A} }$ 是實數矩陣，下述之遞歸計算式適用於矩陣 ${\displaystyle \mathbf {D} }$ 及 ${\displaystyle \mathbf {L} }$ 中的所有元素︰

${\displaystyle D_{j}=A_{jj}-\sum _{k=1}^{j-1}L_{jk}^{2}D_{k}}$

${\displaystyle L_{ij}={\frac {1}{D_{j}}}\left(A_{ij}-\sum _{k=1}^{j-1}L_{ik}L_{jk}D_{k}\right),\qquad {\text{for }}i>j.}$

如果 ${\displaystyle \mathbf {A} }$ 是複埃爾米特矩陣，則矩陣 ${\displaystyle \mathbf {D} }$ 及 ${\displaystyle \mathbf {L} }$ 的計算適用以下公式：

${\displaystyle D_{j}=A_{jj}-\sum _{k=1}^{j-1}L_{jk}L_{jk}^{*}D_{k}}$

${\displaystyle L_{ij}={\frac {1}{D_{j}}}\left(A_{ij}-\sum _{k=1}^{j-1}L_{ik}L_{jk}^{*}D_{k}\right),\qquad {\text{for }}i>j.}$

同樣地，若有需求，該讀取順序可以逐一計算矩陣中的每一個元素。

當${\displaystyle \mathbf {A} }$是一個不定矩陣時，LDL分解在未經過謹慎的繞軸旋轉的情況下是不穩定的；[7] 特別是，分解出的矩陣的元素可能是任意的。一个可取的改進手段是執行矩陣分塊後再執行分解，通常將原矩陣分為包含${\displaystyle 2\times 2}$的小矩陣的分塊矩陣：[8]

${\displaystyle {\begin{aligned}{\mathbf {A=LDL} ^{\mathrm {T} }}&={\begin{pmatrix}\mathbf {I} &0&0\\\mathbf {L} _{21}&\mathbf {I} &0\\\mathbf {L} _{31}&\mathbf {L} _{32}&\mathbf {I} \\\end{pmatrix}}{\begin{pmatrix}\mathbf {D} _{1}&0&0\\0&\mathbf {D} _{2}&0\\0&0&\mathbf {D} _{3}\\\end{pmatrix}}{\begin{pmatrix}\mathbf {I} &\mathbf {L} _{21}^{\mathrm {T} }&\mathbf {L} _{31}^{\mathrm {T} }\\0&\mathbf {I} &\mathbf {L} _{32}^{\mathrm {T} }\\0&0&\mathbf {I} \\\end{pmatrix}}\\&={\begin{pmatrix}\mathbf {D} _{1}&&(\mathrm {symmetric} )\\\mathbf {L} _{21}\mathbf {D} _{1}&\mathbf {L} _{21}\mathbf {D} _{1}\mathbf {L} _{21}^{\mathrm {T} }+\mathbf {D} _{2}&\\\mathbf {L} _{31}\mathbf {D} _{1}&\mathbf {L} _{31}\mathbf {D} _{1}\mathbf {L} _{21}^{\mathrm {T} }+\mathbf {L} _{32}\mathbf {D} _{2}&\mathbf {L} _{31}\mathbf {D} _{1}\mathbf {L} _{31}^{\mathrm {T} }+\mathbf {L} _{32}\mathbf {D} _{2}\mathbf {L} _{32}^{\mathrm {T} }+\mathbf {D} _{3}\end{pmatrix}}\end{aligned}}}$

在此，矩陣中每一個元素都是一個子矩陣。故此，可依照上述遞歸公式類比計算：

${\displaystyle \mathbf {D} _{j}=\mathbf {A} _{jj}-\sum _{k=1}^{j-1}\mathbf {L} _{jk}\mathbf {D} _{k}\mathbf {L} _{jk}^{\mathrm {T} }}$

${\displaystyle \mathbf {L} _{ij}=\left(\mathbf {A} _{ij}-\sum _{k=1}^{j-1}\mathbf {L} _{ik}\mathbf {D} _{k}\mathbf {L} _{jk}^{\mathrm {T} }\right)\mathbf {D} _{j}^{-1}}$

上述計算涉及矩陣相乘同精確的求逆，故而實踐中不能使用過大的分塊矩陣。

在實踐中經常有修正科列斯基分解的需求。即，經已計算一些矩陣${\displaystyle \mathbf {A} }$的科列斯基分解${\displaystyle \mathbf {A} =\mathbf {LL} ^{*}}$，然後在${\displaystyle \mathbf {A} }$上稍作修改而產生的矩陣${\displaystyle {\tilde {\mathbf {A} }}}$，欲對其進行一個修正的科列斯基分解${\displaystyle {\tilde {\mathbf {A} }}={\tilde {\mathbf {L} }}{\tilde {\mathbf {L} }}^{*}}$。問題是，能否用已知的${\displaystyle \mathbf {A} }$的分解去修正${\displaystyle {\tilde {\mathbf {A} }}}$的分解。

如果修正矩陣${\displaystyle {\tilde {\mathbf {A} }}=\mathbf {A} +\mathbf {xx} ^{*}}$，則其修正的分解被稱為秩為1的修正（相加）。。

以下是一個 MATLAB 語言寫的實現秩為1的修正的小程式：

function [L] = cholupdate(L,x)
    n = length(x);
    for k=1:n
        r = sqrt(L(k,k)^2 + x(k)^2);
        c = r / L(k, k);
        s = x(k) / L(k, k);
        L(k, k) = r;
        L(k+1:n,k) = (L(k+1:n,k) + s*x(k+1:n)) / c;
        x(k+1:n) = c*x(k+1:n) - s*L(k+1:n,k);
    end
end

如果${\displaystyle {\tilde {\mathbf {A} }}=\mathbf {A} -\mathbf {x} \mathbf {x} ^{*}}$，那麽只有當${\displaystyle {\tilde {\mathbf {A} }}}$仍然是正定的時候該方法才適用。 此條件下，上面的代碼也可用於相減的情況，只需要將r和L(k+1:n,k)的賦值式的加號替換成减號。

• Sur la résolution numérique des systèmes d'équations linéaires, Cholesky's 1910 manuscript, online and analyzed on BibNum（页面存档备份，存于） （法文） （英文） [for English, click 'A télécharger']

• Hazewinkel, Michiel (编), , , Springer, 2001, ISBN 978-1-55608-010-4
• . PlanetMath.
• Cholesky Decomposition, The Data Analysis BriefBook
• Cholesky Decomposition（页面存档备份，存于） on www.math-linux.com
• Cholesky Decomposition Made Simple（页面存档备份，存于） on Science Meanderthal

• LAPACK（页面存档备份，存于），一個求解稠密線性代數問題的FORTRAN子程式的集合
• ALGLIB（页面存档备份，存于），包含部分從LAPACK到C++、C#、Delphi、Visual Basic等的埠
• libflame（页面存档备份，存于），帶有LAPACK功能的C函數庫
• 德克薩斯大學奧斯汀分校，有關科列斯基分解高性能實現的筆記及影片
• Google，Library "Ceres Solver"（页面存档备份，存于）.
• MATLAB程式，LDL分解
• Armadillo（页面存档备份，存于），一個C++線性函數包 (package)

• 生成相关的随机变量，Martin Haugh，哥倫比亞大學

• 線上矩陣運算，可執行科列斯基分解。