科列斯基分解

线性代数中，科列斯基分解（Cholesky decomposition 或 Cholesky factorization）是指将一个正定的埃尔米特矩阵分解成一个下三角矩阵与其共轭转置之乘积。这种分解方式在提高代数运算效率、蒙特卡罗方法等场合中十分有用。实数矩阵的科列斯基分解由安德烈-路易·科列斯基最先发明。实际应用中，科列斯基分解在求解线性方程组中的效率约两倍于LU分解。

## 描述

对正定埃尔米特矩阵$\mathbf{A}$进行科列斯基分解，即求矩阵$\mathbf{L}$使下式成立

$$\mathbf{A}=\mathbf{LL}^*$$

其中，$\mathbf{L}$是一个下三角矩阵且所有对角元素均为正实数，$\mathbf{L}^*$表示$\mathbf{L}$的共轭转置。每一个正定埃尔米特矩阵都有一个唯一的科列斯基分解。

当矩阵$\mathbf{A}$是一个半正定的埃尔米特矩阵，若允许$\mathbf{L}$的对角线元素为零，则$\mathbf{A}$也存在上述形式的分解。

当$\mathbf{A}$为实数矩阵，则$\mathbf{L}$也为实数矩阵且科列斯基分解可改写成$\mathbf{A}=\mathbf{LL}^\mathbf{T}$。

当$\mathbf{A}$是正定矩阵时，科列斯基分解是唯一的，即只存在一个对角元素均严格大于零的下三角矩阵，使$\mathbf{A}=\mathbf{LL}^*$成立。然而，当$\mathbf{A}$是半正定时，分解则不一定是唯一的。

定理的逆命题自然成立：对于某些可逆矩阵$\mathbf{L}$（下三角矩阵或其他矩阵），如果$\mathbf{A}$可被写成$\mathbf{LL}^*$，则$\mathbf{A}$是一个正定的埃尔米特矩阵。

## LDL分解

经典科列斯基分解的一个变形是LDL分解，即

$$\mathbf{A}=\mathbf{LDL}^*$$

其中，$\mathbf{L}$是一个单位下三角矩阵，$\mathbf{D}$是一个对角矩阵。

该分解与经典科列斯基分解犹有关联，如下：

$$\mathbf{A}=\mathbf{LDL}^*=\mathbf{LD}^{\frac{1}{2}}(\mathbf{D}^{\frac{1}{2}})^*\mathbf{L}^*=\mathbf{LD}^\frac{1}{2} (\mathbf{LD}^\frac{1}{2})^*$$

或者，由于$\mathbf{L}$的对角元素必须为$1$，且$\mathbf{L}^{Cholesky}$与$\mathbf{L}$都是下三角矩阵，故如果已知经典科列斯基分解$\mathbf{L}^{Cholesky}$，则$\mathbf{LDL} ^{T}$形式可依下式求出，设S是包含$\mathbf {L} ^{Cholesky}$的对角元素的对角矩阵。

$$\mathbf D = \mathbf S^2$$

$$\mathbf L = \mathbf L^{Cholesky} \mathbf S^{-1}$$

LDL变形如果得以有效执行，构造及使用时所需求的空间及计算的复杂性与经典科列斯基分解是相同的，但是可避免提取平方根。某些不存在科列斯基分解的不定矩阵，也可以执行LDL分解，此时矩阵$\mathbf{D}$中会出现负数元素。因此人们更倾向于使用LDL分解。对于实数矩阵，该种分解的形式可被改写成

$$\mathbf{A}=\mathbf{LDL}^\mathbf{T}$$

此形式通常称为LDLT分解（或LDLT分解）。它与实对称矩阵的特征分解密切相关，因为对于实对称矩阵，存在特征分解$\mathbf{A}=\mathbf{Q\Lambda Q}^\mathbf{T}$。

## 实例

以下乃一个实对称矩阵的科列斯基分解︰

$$\begin{align} \left( <pre> \begin{array}{*{3}{r}} 4 & 12 & -16 \\ 12 & 37 & -43 \\ -16 & -43 & 98 \\ \end{array} </pre> \right) = \left( <pre> \begin{array}{*{3}{r}} 2 & 0 & 0 \\ 6 & 1 & 0 \\ -8 & 5 & 3 \\ \end{array} </pre> \right) \left( <pre> \begin{array}{*{3}{r}} 2 & 6 & -8 \\ 0 & 1 & 5 \\ 0 & 0 & 3 \\ \end{array} </pre> \right) \end{align}$$

以下乃其LDLT分解︰

$$\begin{align} \left( <pre> \begin{array}{*{3}{r}} 4 & 12 & -16 \\ 12 & 37 & -43 \\ -16 & -43 & 98 \\ \end{array} </pre> \right) & = \left( <pre> \begin{array}{*{3}{r}} 1 & 0 & 0 \\ 3 & 1 & 0 \\ -4 & 5 & 1 \\ \end{array} </pre> \right) \left( <pre> \begin{array}{*{3}{r}} 4 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 9 \\ \end{array} </pre> \right) \left( <pre> \begin{array}{*{3}{r}} 1 & 3 & -4 \\ 0 & 1 & 5 \\ 0 & 0 & 1 \\ \end{array} </pre> \right) \end{align}$$

## 应用

科列斯基分解主要被用于线性方程组$\mathbf{Ax} = \mathbf{b}$的求解。如果A是对称正定的，我们可以先求出$\mathbf{A}=\mathbf{LL}^\mathbf{T}$，随后借向后替换法对y求解$\mathbf{Ly}=\mathbf{b}$，再以向前替换法对x求解$\mathbf{L}^\mathbf{T}\mathbf{x} =\mathbf{y}$即得最终解。

另一种可避免在计算$\mathbf{LL}^\mathbf{T}$时需要解平方根的方法就是计算$\mathbf{A} = \mathbf{LDL}^\mathrm{T}$，然后对y求解$\mathbf{Ly} = \mathbf{b}$，最后求解$\mathbf{DL}^\mathrm{T}\mathbf{x} = \mathbf{y}$

对于可以被改写成对称矩阵的线性方程组，科列斯基分解及其LDL变形是一个较高效率及较高数值稳定性的求解方法。相比之下，其效率几近为LU分解的两倍。

### 矩阵求逆

若欲对埃尔米特矩阵直接求逆，可以通过科列斯基分解，类似求解线性方程组一般求出逆矩阵，这需要$n^3$次运算（$\frac{1}{2}n^3$次乘法运算）。该方法即便要求逐步计算也非常有效率。

非埃尔米特矩阵$\mathbf{B}$也可以通过下例性质求逆，因为其中$\mathbf{BB}^*$总是埃尔米特矩阵︰

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

## 计算

有各种各样的方法用于计算科列斯基分解。常用的算法的计算复杂度在一般情况下为$O(n^3)$。

下面的算法何者较快取决于执行时的具体条件。总体而言，第一个算法会稍慢，因为其以一种不太寻常的方式读取数据。

### 科列斯基算法

用于计算矩阵$\mathbf{L}$的科列斯基算法，是以高斯消元法为基础而调整得来的。

递归算法由$i:=1$开始，令

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

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

$$\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},$$

其中，$\mathbf{I}_{i-1}$代表$i-1$维的单位矩阵。

此时定义矩阵$\mathbf{L}_i$为

$$\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},$$

随即$\mathbf{A}^{(i)}$可以被改写成

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

其中

$$\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}.$$

注意︰在此$\mathbf{b}_i\mathbf{b}^*_i$是一个外积。

重复此步骤直到$i$从$1$到$n$。$n$步之后，我们得到$A^{(n+1)} = \mathbf{I}$。因此，所求的下三角矩阵$L$为

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

### 科列斯基-巴那齐耶维茨及科列斯基-克劳特算法

<image>

写出等式$\mathbf{A}=\mathbf{LL}^*$。

$$\begin{align} {\mathbf{A=LL^T}} & = \begin{pmatrix} \mathbf{L}_{11} & 0 & 0 \\ <pre> \mathbf{L}_{21} & \mathbf{L}_{22} & 0 \\ \mathbf{L}_{31} & \mathbf{L}_{32} & \mathbf{L}_{33}\\ </pre> \end{pmatrix} \begin{pmatrix} \mathbf{L}_{11} & \mathbf{L}_{21} & \mathbf{L}_{31} \\ <pre> 0 & \mathbf{L}_{22} & \mathbf{L}_{32} \\ 0 & 0 & \mathbf{L}_{33} </pre> \end{pmatrix} \\ & = \begin{pmatrix} \mathbf{L}_{11}^2 & &(\text{symmetric}) \\ <pre> \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 </pre> \end{pmatrix} \end{align}$$

则得到矩阵$\mathbf{L}$的元素的计算公式如下︰

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

$$\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$$

只要$\mathbf{A}$是实数正定矩阵，则平方根下的表达式恒为正。

对于复埃尔米特矩阵，则适用如下公式：

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

$$\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.$$

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

- 科列斯基-巴那齐耶维茨算法从矩阵L的左上角开始，依行进行计算。
- 科列斯基-克劳特算法从矩阵L的左上角开始，依列进行计算。

若有需要，整个矩阵可以逐个元素计算得出，无论使用何种顺序读取。

### 计算的稳定性

假设我们要求解一个良置的线性方程组。采用了LU分解的算法，除非进行某种绕轴旋转，否则是不稳定的；如果算法进行了绕轴旋转，则其误差取决于矩阵的增长因子，这个数通常（但非总是）很小。

现在，假设算法适用科列斯基分解。如上所述，算法的效率将会是原来的两倍。此外，无必要进行绕轴旋转且误差总是很小。具体而言，若要求解$\mathbf{Ax}= \mathbf{b}$，$\mathbf{y}$表示已计算出的解，然后能解出干扰方程组$(\mathbf{A}+\mathbf{E})\mathbf{y}=\mathbf{b}$，其中

$$\|\mathbf{E}\|_2 \le c_n \varepsilon \|\mathbf{A}\|_2.$$

在这里，$\|\quad\|_2$是指矩阵2-范数，而$c_n$是一个取决于$n$的小常数，$\varepsilon$表示单位舍入。

值得一提的是，科列斯基分解与平方根的使用有关。如果被分解的矩阵是正定的，那么只要运算精确，矩阵中带有平方根的元素的平方根下的数字永远是正数。不幸的是，由于存在舍入误差，这些数字可能为负数，并使算法搁浅。然而，此种情况仅当矩阵为病态时才会出现。一种可解决此问题的方法，是增添一个对角修正矩阵至待分解矩阵，以增加矩阵的正定性。此法虽或将减少分解的准确性，但有在某些情况却颇有作为；譬如，执行应用于最优化的牛顿法时，若初期值距最优值较远，则此时引入对角矩阵可以提高算法的稳定性。

### LDL分解

科列斯基分解的另一种形式——LDL分解的计算方式如下所示。

$$\begin{align} {\mathbf{A=LDL}^\mathrm{T}} & = \begin{pmatrix} 1 & 0 & 0 \\ <pre> L_{21} & 1 & 0 \\ L_{31} & L_{32} & 1\\ </pre> \end{pmatrix} \begin{pmatrix} D_1 & 0 & 0 \\ <pre> 0 & D_2 & 0 \\ 0 & 0 & D_3\\ </pre> \end{pmatrix} \begin{pmatrix} 1 & L_{21} & L_{31} \\ <pre> 0 & 1 & L_{32} \\ 0 & 0 & 1\\ </pre> \end{pmatrix} \\ & = \begin{pmatrix} D_1 & &(\mathrm{symmetric}) \\ <pre> L_{21}D_1 & L_{21}^2D_1 + D_2 & \\ L_{31}D_1 & L_{31}L_{21}D_{1}+L_{32}D_2 & L_{31}^2D_1 + L_{32}^2D_2+D_3. </pre> \end{pmatrix} \end{align}$$

如果 \mathbf A 是实数矩阵，下述之递归计算式适用于矩阵 \mathbf D 及 \mathbf L 中的所有元素︰

$$D_j = A_{jj} - \sum_{k=1}^{j-1} L_{jk}^2 D_k$$

$$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.$$

如果 \mathbf A 是复埃尔米特矩阵，则矩阵 \mathbf D 及 \mathbf L 的计算适用以下公式：

$$D_{j} = A_{jj} - \sum_{k=1}^{j-1} L_{jk}L_{jk}^* D_k$$

$$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.$$

同样地，若有需求，该读取顺序可以逐一计算矩阵中的每一个元素。

### 分块矩阵变形

当$\mathbf{A}$是一个不定矩阵时，LDL分解在未经过谨慎的绕轴旋转的情况下是不稳定的；特别是，分解出的矩阵的元素可能是任意的。一个可取的改进手段是执行矩阵分块后再执行分解，通常将原矩阵分为包含$2\times2$的小矩阵的分块矩阵：

$$\begin{align} {\mathbf{A=LDL}^\mathrm{T}} & = \begin{pmatrix} <pre> \mathbf I & 0 & 0 \\ \mathbf L_{21} & \mathbf I & 0 \\ \mathbf L_{31} & \mathbf L_{32} & \mathbf I\\ </pre> \end{pmatrix} \begin{pmatrix} <pre> \mathbf D_1 & 0 & 0 \\ 0 & \mathbf D_2 & 0 \\ 0 & 0 & \mathbf D_3\\ </pre> \end{pmatrix} \begin{pmatrix} <pre> \mathbf I & \mathbf L_{21}^\mathrm T & \mathbf L_{31}^\mathrm T \\ 0 & \mathbf I & \mathbf L_{32}^\mathrm T \\ 0 & 0 & \mathbf I\\ </pre> \end{pmatrix} \\ & = \begin{pmatrix} <pre> \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 </pre> \end{pmatrix} \end{align}$$

在此，矩阵中每一个元素都是一个子矩阵。故此，可依照上述递归公式类比计算：

$$\mathbf D_j = \mathbf A_{jj} - \sum_{k=1}^{j-1} \mathbf L_{jk} \mathbf D_k \mathbf L_{jk}^\mathrm T$$

$$\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}$$

上述计算涉及矩阵相乘同精确的求逆，故而实践中不能使用过大的分块矩阵。

### 修正分解

在实践中经常有修正科列斯基分解的需求。即，经已计算一些矩阵$\mathbf{A}$的科列斯基分解$\mathbf{A}=\mathbf{LL}^*$，然后在$\mathbf{A}$上稍作修改而产生的矩阵$\tilde{\mathbf{A}}$，欲对其进行一个修正的科列斯基分解$\tilde{\mathbf{A}} = \tilde{\mathbf{L}} \tilde{\mathbf{L}}^*$。问题是，能否用已知的$\mathbf{A}$的分解去修正$\tilde{\mathbf{A}}$的分解。

#### 秩为1的修正（相加）

如果修正矩阵${\tilde{\mathbf{A}}}=\mathbf{A}+\mathbf{xx}^*$，则其修正的分解被称为秩为1的修正（相加）。

以下是一个 MATLAB 语言写的实现秩为1的修正的小程式：

```matlab
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
```

#### 秩为1的修正（相减）

如果$\tilde{\mathbf{A}} = \mathbf{A} - \mathbf{x} \mathbf{x}^*$，那么只有当$\tilde{\mathbf{A}}$仍然是正定的时候该方法才适用。

此条件下，上面的代码也可用于相减的情况，只需要将`r`和`L(k+1:n,k)`的赋值式的加号替换成减号。

## 证明半正定矩阵的科列斯基分解唯一

上面的算法表明，每一个正定矩阵$A$都有一个科列斯基分解。此结论可以加入一些限制条件后，推广到半正定的情况。但该条件并未被完全建立，例如，它未给出明确的数值算法以计算科列斯基因子。

如果$A$是一个$n\times n$的半正定矩阵，则序列$\{\mathbf{A}\}=\{\mathbf{A}+(1/k)\mathbf{I}_n\}$是一个由正定矩阵构成的序列。而且，在算子范数中

$$\mathbf{A}_k\rightarrow\mathbf{A}$$

在正定的情形下，每一个$\mathbf{A}_k$都有其科列斯基分解$\mathbf{A}_k = \mathbf{L}_k\mathbf{L}^*_k$。根据算子范数的性质。

$$\| \mathbf{L}_k \|^2 \le \| \mathbf{L}_k \mathbf{L}_k ^* \| = \| \mathbf{A}_k \|$$

因而$\{\mathbf{L}_k\}$是算子巴拿赫空间上的一个有界集，因此是相对紧空间（因为基础向量空间是有限维的）。因此，它有一个带有限制$\mathbf{L}$收敛子序列，也记作$\{\mathbf{L}_k\}$。容易验证，矩阵$\mathbf{L}$具有所需的特性，例如，满足$\mathbf{A}= \mathbf{LL}^*$以及$\mathbf{L}$是下三角矩阵且其对角元素非负。对于所有的$x$和$y$都有。

$$\langle \mathbf{A} x, y \rangle = \langle \lim \mathbf{A}_k x, y \rangle = \langle \lim \mathbf{L}_k \mathbf{L}_k^* x, y \rangle = \langle \mathbf{L} \mathbf{L}^*x, y \rangle$$

因此，$\mathbf{A}=\mathbf{LL}^*$。因为基础向量空间是有限维的，所以算子空间的所有拓扑结构都是等价的。故此，范数上$\mathbf{L}_k$收敛于$\mathbf{L}$，意味着，$\mathbf{L}_k$的每个元素都独立地收敛于$\mathbf{L}$。此同时暗示，由于每个$\mathbf{L}_k$都是对角元素非负的下三角矩阵，$\mathbf{L}$亦如此。

## 推广

科列斯基分解可以推广到元素中包含算子的矩阵（称为算子矩阵）。令$\{\mathcal{H}_n\}$是希尔伯特空间上的一个序列。考虑如下算子矩阵

$$\mathbf{A} = \begin{bmatrix} \mathbf{A}_{11} & \mathbf{A}_{12} & \mathbf{A}_{13} & \; \\ \mathbf{A}_{12}^* & \mathbf{A}_{22} & \mathbf{A}_{23} & \; \\ \mathbf{A} _{13}^* & \mathbf{A}_{23}^* & \mathbf{A}_{33} & \; \\ \; & \; & \; & \ddots \end{bmatrix}$$

满足直和

$$\mathcal{H} = \oplus _n \mathcal{H}_n,$$

其中每一个

$$\mathbf{A}_{ij} : \mathcal{H}_j \rightarrow \mathcal{H} _i$$

都是一个有界算子。如果$\mathbf{A}$是正定或半正定的，即对于有限的$k$及任意的

$$h \in \oplus _{n = 1 }^k \mathcal{H}_k ,$$

都有$\langle h,\mathbf{A}h\rangle\ge0$，则存在一个下三角的算子矩阵$\mathbf{L}$使得$\mathbf{A} = \mathbf{LL}^*$。此外也可以把$\mathbf{L}$的对角元素化为正数。

## 用程式语言实现

- C语言︰GNU科学库提供了几个科列斯基分解的实现。
- Maxima电脑代数系统︰函数Cholesky可用于计算科列斯基分解。
- GNU Octave数值计算系统提供了一些函数以计算，修正和应用科列斯基分解。
- LAPACK库提供了一个高性能的科列斯基分解的实现，可以以Fortran，C语言及其他大多数语言读取。
- 在Python中，numpy.linalg模块中的命令Cholesky可执行科列斯基分解。
- 在Matlab中，chol命令可以简单地对一个矩阵进行科列斯基分解。
- 在R语言中，chol函数可进行科列斯基分解。
- 在Wolfram Mathematica中，函数CholeskyDecomposition可以对一个矩阵执行科列斯基分解。
- 在C++中，armadillo库中的命令chol可执行科列斯基分解。特征库提供了稀疏矩阵和稠密矩阵的科列斯基分解。
- 在Analytica中，函数Decompose提供科列斯基分解。
- Apache Commons数学库中有科列斯基分解的实现，可应用于Java、Scala及任何其他JVM语言。

## 外部链接

### 科学史

- 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'

### 资讯

- 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，哥伦比亚大学

### 线上计算机

- 线上矩阵运算，可执行科列斯基分解。