数值代数 - LU 分解的充要条件，特殊矩阵、LDU 分解、列主元 Gauss 消元法

LU 分解的进一步讨论

可分解的条件

回顾 Gauss 消元法的过程，我们从消元过程当中得到了 $LU$ 分解，但是消元过程当中有一个限制条件是 $a_{ii}^{(i-1)}\ne 0$ ，这个条件并不是一开始就能从 $A$ 当中看出来的. 下面给出一个充要条件.

定理：主元非零的充要条件

主元 $a_{ii}^{(i-1)}\ne 0,\forall i=1,2,\cdots ,k$ 的充要条件是 $A$ 的 $i$ 阶顺序主子阵都是非奇异的.

证明： 考虑归纳法，当 $i=1$ 时结论是显然的，我们假设 $i⩽k-1$ 时 $A_1,\cdots ,A_{k-1}$ 都是非奇异的，现在仅需证明 $A_k$ 非奇异的等价条件为 $a_{kk}^{(k-1)}\ne 0$ .

现在我们考虑 $k-1$ 次消元之后的矩阵：

$$A^{(k-1)}=\left(\begin{array}{cc}{A}_{11}^{(k-1)}& A_{12}^{(k-1)}\\ O& A_{21}^{(k-1)}\end{array}\right)$$

这相当于将其前 $k-1$ 阶主子阵计算了 $LU$ 分解当中的 $L$，记 $(L_i{)}_k$ 为第 $i$ 个 Gauss 变换阵的前 $k$ 阶主子阵，于是

$$(L_{k-1}{)}_k(L_{k-2}{)}_k\cdots (L_2{)}_k(L_1{)}_k{A}_k=\left(\begin{array}{cc}{A}_{11}^{(k-1)}& \ast \\ O& a_{kk}^{(k-1)}\end{array}\right)$$

取行列式，由于 Gauss 变换阵对角线元素均为 $1$ ，从而

$$|A_k|=|A_{11}^{(k-1)}|a_{kk}^{(k-1)}$$

由于 $|A_{11}^{(k-1)}|$ 已经确定非零，所以剩下两者的等价性证明完毕. $\square$

从而根据这个定理有

定理：LU 分解的等价条件

矩阵 $A$ 存在唯一的 LU 分解当且仅当 $A$ 的所有顺序主子阵均非奇异.

证明： 仅需证明唯一性，假如存在两个 LU 分解：

$$A=L_1{U}_1=L_2{U}_2$$

那么有

$$L_2^{-1}{L}_1=U_2{U}_1^{-1}$$

左侧是单位下三角和单位下三角的乘积，右侧是两个上三角的乘积，说明两边应该均为单位阵，从而唯一性成立. $\square$

计算量

直接计算三角分解时，我们将加减乘除都放在一起考虑，先给出教材上的三角分解：

for k = 1:n-1
	A(k+1:n,k) = A(k+1:n,k)/A(k,k)
	A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k) * A(k,k+1:n)
end

那么我们首先看到第一步是 $n-k$ 次除法，然后发生了 $(n-k{)}^2$ 次乘法和 $(n-k{)}^2$ 次减法，因此总计算次数为

$$\sum _{k=1}^{n-1}[(n-k)+2(n-k{)}^2]=\frac{2}{3}{n}^3+O(n^2)$$

也就是说，计算复杂度是 $O(n^3)$ 的，相比较 Laplace 方法而言把阶乘复杂度变为多项式复杂度，是一个切实的飞跃.

但是可以注意到，LU 分解理应得到两个矩阵，为什么这里我们直接对 $A$ 作了就地操作？事实上，对 $A$ 作就地操作后，严格下三角部分是 $L$ 的严格下三角部分，上三角部分则是 $U$ 的上三角部分. 这样做的原因在于节省内存空间，当问题规模较大时，少存一个矩阵意味着空闲出相当多的内存空间.

特殊矩阵

特殊矩阵：对角占优矩阵

对角占优矩阵存在唯一的 LU 分解.

这个定理可以直接从下面的问题得到：出自数值分析第十四次作业.

练习 5.6

设 $\mathbf{A}$ 为 $n$ 阶按行严格对角占优矩阵，经 Gauss 消去法一步后 $\mathbf{A}$ 变为如下形式： ${\mathbf{A}}^{(2)}=\left(\begin{array}{cc}{a}_{11}& {\mathbf{A}}_{12}\\ 0& {\mathbf{A}}_{22}^{(2)}\end{array}\right)$ 试证明：${\mathbf{A}}_{22}^{(2)}$ 是 $n-1$ 阶按行严格对角占优矩阵.

即证明：

$$\sum _{j=2,j\ne i}^n\left|a_{ij}-\frac{a_{i1}{a}_{1j}}{a_{11}}\right|<\left|a_{ii}-\frac{a_{i1}{a}_{1i}}{a_{11}}\right|$$

从左到右有

$$\begin{array}{rl}\sum _{j=2,j\ne i}^n\left|a_{ij}-\frac{a_{i1}{a}_{1j}}{a_{11}}\right|& ⩽\sum _{j=2}^n|a_{ij}|+\left|\frac{a_{i1}}{a_{11}}\right|\sum _{j=2,j\ne i}^n|a_{1j}|\\ & =\sum _{j=1,j\ne i}^n|a_{ij}|-|a_{i1}|+\left|\frac{a_{i1}}{a_{11}}\right|\left(\sum _{j=2}^n|a_{1j}|-|a_{1i}|\right)\\ & <|a_{ii}|-|a_{i1}|+\left|\frac{a_{i1}}{a_{11}}\right|(|a_{11}|-|a_{1i}|)\\ & =|a_{ii}|-\left|\frac{a_{i1}{a}_{1i}}{a_{11}}\right|⩽\mathrm{RHS}\end{array}$$

于是结论成立. $\square$

LDU 分解

我们考虑 LU 分解进一步细化为如下形式：

$$A=LDU$$

这里的 $L$ 还是 LU 分解中的单位下三角矩阵，$U$ 为单位上三角矩阵，$D$ 为对角矩阵，可以看出，我们只不过就是把 $U$ 的对角线元素提出来了而已.

似乎我们可以到这里就结束对 LDU 分解的讨论了，但是一个有意思的问题在于：$D$ 当中的元素到底是什么呢？是 $A$ 的特征值，还是通过 $A$ 的某些量计算出来的呢？

首先 $D$ 的对角线元素并不是 $A$ 的特征值，考虑

$$A=\left(\begin{array}{cc}2& 1\\ 1& 2\end{array}\right),L=\left(\begin{array}{cc}1& 0\\ 0.5& 1\end{array}\right),U=\left(\begin{array}{cc}2& 1\\ 0& 1.5\end{array}\right)$$

我们可以知道 $A$ 的特征值为 $3,1$ ，但是可以看到 $U$ 的对角线元素不是特征值.

考虑从 $i$ 到 $i+1$ 步的 LDU 分解，取每个矩阵的前 $i+1$ 行，前 $i+1$ 列，设 $A_i$ 为 $A$ 的第 $i$ 个顺序主子式：

$$A_{i+1}=\left(\begin{array}{cc}{A}_i& A(1:i,i+1)\\ A(i+1,1:i)& A(i+1,i+1)\end{array}\right)$$

于是有

$$\begin{array}{rl}& L_{i+1}=\left(\begin{array}{cc}{L}_i& 0\\ L(i+1,1:i)& 1\end{array}\right),U_{i+1}=\left(\begin{array}{cc}{U}_i& U(1:i,i+1)\\ 0& 1\end{array}\right)\\ & D_{i+1}=\left(\begin{array}{cc}{D}_i& 0\\ 0& D(i+1,i+1)\end{array}\right)\end{array}$$

取行列式可得

$$|A_{i+1}|=|D_i|\cdot |D(i+1,i+1)|$$

而 $|D_i|=|A_i|$ ，因此

$$d_{i+1}=|D(i+1,i+1)|=\frac{|A_{i+1}|}{|A_i|},d_1=A_1$$

由 LU 分解求逆矩阵

从 LU 分解出发，考虑

$$L=(L_{n-1}{L}_{n-2}\cdots L_1{)}^{-1},U=A^{(n-1)}$$

这就让我们发现一个有意思的现象：$L$ 是很容易求逆的，而对 $U$ ，它也是一个三角阵，可以用类似的方法对其求逆，这就启发我们通过：

$$A^{-1}=U^{-1}{L}^{-1}$$

进行求逆. 对三角矩阵求逆是简单的（见第一次作业）.

当然，还有一种方案是，如果我们回顾高代当中的初等变换求逆方法：

$$\left(\begin{array}{cc}A& I\end{array}\right)\Rightarrow \left(\begin{array}{cc}I& A^{-1}\end{array}\right)$$

那么设 $A^{-1}$ 当中的第 $i$ 个列向量为 ${\mathbf{x}}_i$ ，对应单位列向量为 ${\mathbf{e}}_i$ ，则

$$A{\mathbf{x}}_i={\mathbf{e}}_i$$

解共 $n$ 个线性方程组即可.

列主元 Gauss 消元法

全主元 Gauss 消元法

现在来看 Gauss 消元法没能考虑到的一个方面，假定我们在消元过程中每个主元都非零，但是出现：

$$\left(\begin{array}{cccc}\ast & & & \\ & a_{ii}^{(i-1)}& \ast & \\ & \ast & \ast & \end{array}\right)$$

这其中的 $a_{ii}^{(i-1)}$ 绝对值较小，换言之，此时作为除数时，它将使得算法极为不稳定，当除数较小时，微小的扰动都将造成巨大的误差.

除此之外，LU 分解的条件也过于严格，我们研究线性方程组是针对全体可逆矩阵的，但是并不一定所有的可逆阵都有这么好的性质.

为解决上述问题，全主元 Gauss 消元法应运而生，它相较于 Gauss 消元法只有一个地方不同，就是在每次消元时，都查找右下角子矩阵当中绝对值最大的元替换到主元位置上.

有了置换过程，我们就能得到如下定理：

定理：推广 LU 分解

设 $A\in {\mathbb{R}}^{n\times n}$ ，则存在排列矩阵 $P,Q\in {\mathbb{R}}^{n\times n}$ ，以及单位下三角阵 $L\in {\mathbb{R}}^{n\times n}$ ，上三角阵 $U\in {\mathbb{R}}^{n\times n}$ 使得

$$PAQ=LU$$

我们只需理解全主元消元法即可知这个定理是显然的，其中

• $L$ 的所有元素均满足 $|l_{ij}|⩽1$ ；这个性质是简单的，因为我们选的主元是较大的，而且其充当除数.
• $U$ 的非零对角元个数正好等于 $\mathrm{rank}(A)$ .

相比于 Gauss 消元法，算法也只多出：

• 置换矩阵的记录；
• 主元的查找和置换.

但是，选主元在这里也消耗了大量时间，注意到每个情况下，每次查找主元时都要将右下角的矩阵搜索个遍，这就会多出：

$$\sum _{k=1}^{n-1}(n-k+1{)}^2=\frac{1}{3}{n}^3+O(n^2)$$

次比较操作.

列主元 Gauss 消元法

列主元是为了解决全主元代价过大的问题的，它的改进就是改为只搜索当列的绝对值最大元，这就使得代价变为：

$$\sum _{k=1}^{n-1}(n-k+1)=\frac{1}{2}{n}^2+O(n)$$

降低了整整一个数量级.