对称正定矩阵的分解

上一节中，我们从高斯消元法出发，介绍了工程实际中更为常用的 LU 分解。 这是一种通用的方阵分解方法，基本上只要矩阵非奇异就可以将之分解为一个下三角和一个上三角矩阵的积。本节我们来看对称的正定矩阵，由于这类矩阵具有一些特殊的性质， 我们可以以更高效的方式将之分解，进而求解线性方程组和矩阵的逆。

1. 对称正定矩阵

一个 \(n \times n\) 的矩阵 \(\boldsymbol{A}\)，如果对于任意的非零向量 \(\boldsymbol{X}\)，都有 \(\boldsymbol{X}^T \boldsymbol{A} \boldsymbol{X} > \boldsymbol{0}\)， 那么该矩阵就是正定的。如果矩阵的转置等于它本身，即 \(\boldsymbol{A}^T = \boldsymbol{A}\)，那么该矩阵就是对称的。同时满足这两个条件的矩阵，就是对称正定矩阵。 这种矩阵具有很重要的物理意义，比如刚体动力学中惯性矩阵就是对称正定的。这类矩阵还有一些有用的数学特性。

一个对称正定矩阵 \(\boldsymbol{A}\) 与任意的可逆方阵 \(\boldsymbol{P}\) 构成的矩阵 \(\boldsymbol{B} = \boldsymbol{P}^T\boldsymbol{A}\boldsymbol{P}\) 也都是对称正定的。 因为 \(\boldsymbol{A}\) 是对称的，所以 \(\boldsymbol{B}^T = \boldsymbol{P}^T\boldsymbol{A}^T\boldsymbol{P} = \boldsymbol{B}\)，即 \(\boldsymbol{B}\) 也对称。 对于任意的非零向量 \(\boldsymbol{X}\)，由于 \(\boldsymbol{P}\) 可逆，所以有向量 \(\boldsymbol{Y} = \boldsymbol{P}\boldsymbol{X}\) 也非零。而方阵 \(\boldsymbol{A}\) 正定， 有 \(\boldsymbol{Y}^T \boldsymbol{A} \boldsymbol{Y} > 0\)，所以 \(\boldsymbol{X}^T\boldsymbol{B}\boldsymbol{X} > 0\)，即 \(\boldsymbol{B}\) 也正定。 利用这一特性，我们保证在例程中构建的矩阵是对称正定的。

如果对称矩阵 \(\boldsymbol{A}\) 是正定的，那么它可以唯一分解为 \(\boldsymbol{A} = \boldsymbol{L}\boldsymbol{L}^T\)，其中 \(\boldsymbol{L}\) 是一个对角元素都大于 0 的下三角矩阵。 这就是所谓的 Cholesky 分解定理。\(\boldsymbol{L}\) 的存在性可以通过数学归纳法来证明，唯一性可以通过反证法完成。 证明方法我们在后续介绍正定矩阵的二次型给出，这里我们主要关注矩阵 \(\boldsymbol{L}\) 的解法。

除此之外，对称的正定矩阵一定是可逆的，并且它的逆也是对称且正定的。对称正定矩阵的行列式和所有特征值都大于零。

2. Cholesky 分解

假设有一个 \(4 \times 4\) 的对称的正定矩阵 \(\boldsymbol{A}\) 如下式左侧所示，按照 Cholesky 定理，我们可以找到唯一的一个下三角矩阵 \(\boldsymbol{L}\) 将其分解为下式右侧的 \(\boldsymbol{L}\boldsymbol{L}^T\)形式。

$$ \underbrace{\begin{bmatrix} a_{11} & a_{21} & a_{31} & a_{41} \\ a_{21} & a_{22} & a_{32} & a_{42} \\ a_{31} & a_{32} & a_{33} & a_{43} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix}}_{\boldsymbol{A}} = \underbrace{\begin{bmatrix} l_{11} & 0 & 0 & 0 \\ l_{21} & l_{22} & 0 & 0 \\ l_{31} & l_{32} & l_{33} & 0 \\ l_{41} & l_{42} & l_{43} & l_{44} \end{bmatrix}}_{\boldsymbol{L}} \underbrace{\begin{bmatrix} l_{11} & l_{21} & l_{31} & l_{41} \\ 0 & l_{22} & l_{32} & l_{42} \\ 0 & 0 & l_{33} & l_{43} \\ 0 & 0 & 0 & l_{44} \end{bmatrix}}_{\boldsymbol{L^T}} $$

我们可以将上式，展开得到如下的形式。由于矩阵 \(\boldsymbol{A}\) 是对称的，所以我们省略了上三角的元素。

$$ \left[ \begin{array}{l} l_{11}l_{11} & & & \\ l_{21}l_{11} & l_{21}l_{21} + l_{22}l_{22} & & \\ l_{31}l_{11} & l_{31}l_{21} + l_{32}l_{22} & l_{31}l_{31} + l_{32}l_{32} + l_{33}l_{33} & \\ l_{41}l_{11} & l_{41}l_{21} + l_{42}l_{22} & l_{41}l_{31} + l_{42}l_{32} + l_{43}l_{33} & l_{41}l_{41} + l_{42}l_{42} + l_{43}l_{43} + l_{44}l_{44} \end{array} \right] = \begin{bmatrix} a_{11} & & & \\ a_{21} & a_{22} & & \\ a_{31} & a_{32} & a_{33} & \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix} $$

与 LU 分解的Crout's 算法类似，我们也可以按照特定的顺序依次求解出 \(l_{ij}\)。 对于第一列有 \(l_{11} = \sqrt{a_{11}}\)，因此 \(l_{i1} = a_{i1} / l_{11}\)。在此基础上对于第二列有，\(l_{22} = \sqrt{a_{22} - l_{21}l_{21}}\)， \(l_{i2} = \frac{1}{l_{22}}\left(a_{i2} - l_{i1}l_{21}\right)\)。可以从中归纳出如下的两式，用于求解 \(\boldsymbol{L}\)。

$$ l_{jj} = \left(a_{jj} - \sum_{k = 1}^{j - 1}l_{jk}^2 \right)^{\frac{1}{2}}, \qquad l_{ij} = \frac{1}{l_{jj}} \left(a_{ij} - \sum_{k = 1}^{j - 1}l_{ik}l_{jk} \right) $$

据说这种算法在数值上十分稳定，不需要像高斯消元法或者 LU 分解那样找一个主元。如果矩阵 \(\boldsymbol{A}\) 分解失败了，那么这个矩阵大概率就不是一个对称的正定矩阵。 Cholesky 分解算法有一个小缺陷，就是求解对角元素的时候有个开平方的操作，这个还是很耗时的。后来有了改进的 LDL' 分解。

3. LDL' 分解

LDL' 分解则是要把对称的正定矩阵 \(\boldsymbol{A}\) 分解成三个矩阵 \(\boldsymbol{L}\boldsymbol{D}\boldsymbol{L}^T\) 相乘的形式。 其中 \(\boldsymbol{D}\) 是一个对角矩阵，并且对角线的元素均大于零，\(\boldsymbol{L}\) 则是对角线全为 1 的下三角矩阵。一个 \(4 \times 4\) 的矩阵就可以写成如下的形式。

$$ \underbrace{\begin{bmatrix} 1 & 0 & 0 & 0 \\ l_{21} & 1 & 0 & 0 \\ l_{31} & l_{32} & 1 & 0 \\ l_{41} & l_{42} & l_{43} & 1 \end{bmatrix}}_{\boldsymbol{L}} \underbrace{\begin{bmatrix} d_1 & & & \\ & d_2 & & \\ & & d_3 & \\ & & & d_4 \end{bmatrix}}_{\boldsymbol{D}} \underbrace{\begin{bmatrix} 1 & l_{21} & l_{31} & l_{41} \\ 0 & 1 & l_{32} & l_{42} \\ 0 & 0 & 1 & l_{43} \\ 0 & 0 & 0 & 1 \end{bmatrix}}_{\boldsymbol{L^T}} = \underbrace{\begin{bmatrix} a_{11} & a_{21} & a_{31} & a_{41} \\ a_{21} & a_{22} & a_{32} & a_{42} \\ a_{31} & a_{32} & a_{33} & a_{43} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix}}_{\boldsymbol{A}} $$

整理一下，有:

$$ \underbrace{\begin{bmatrix} d_1 & 0 & 0 & 0 \\ l_{21}d_1 & d_2 & 0 & 0 \\ l_{31}d_1 & l_{32}d_2 & d_3 & 0 \\ l_{41}d_1 & l_{42}d_2 & l_{43}d_3 & d_4 \end{bmatrix}}_{\boldsymbol{LD}} \underbrace{\begin{bmatrix} 1 & l_{21} & l_{31} & l_{41} \\ 0 & 1 & l_{32} & l_{42} \\ 0 & 0 & 1 & l_{43} \\ 0 & 0 & 0 & 1 \end{bmatrix}}_{\boldsymbol{L^T}} = \underbrace{\begin{bmatrix} a_{11} & a_{21} & a_{31} & a_{41} \\ a_{21} & a_{22} & a_{32} & a_{42} \\ a_{31} & a_{32} & a_{33} & a_{43} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix}}_{\boldsymbol{A}} $$

与 LU 和 Cholesky 分解类似，我们只需要按照特定的顺序依次求解，就可以求得矩阵 \(\boldsymbol{L}\) 和 \(\boldsymbol{D}\)。它们的元素可以通过如下两式写出:

$$ d_{j} = a_{jj} - \sum_{k = 1}^{j - 1}l_{jk}^2 d_k, \qquad l_{ij} = \frac{1}{d_j} \left(a_{ij} - \sum_{k = 1}^{j - 1}l_{ik} d_k l_{jk} \right) $$

显然，这个分解就不存在开放运算了，效率上要比 Cholesky 分解高一些。

4. 完

在 XiaoTuMathBox 中，我们定义了类Cholesky 和类LDLT，在它们的构造函数中分别完成了 Cholesky 和 LDL' 分解。 并且还提供了Solve 和 Inverse 分别用于求解方程组和原矩阵的逆。