MA&ALA3.10_LU分解 (The LU Factorization)

注：本文是对Matrix Analysis and Applied Linear Algebra一书3.10节The LU Factorization的学习笔记

$\mathbf{A}$的LU分解：$\mathbf{A}=\mathbf{L}\mathbf{U}$，其中$\mathbf{L}$是下三角矩阵，$\mathbf{U}$是上三角矩阵

LU分解其实是高斯消元的过程用矩阵描述

为了更好地理解这一点，我们首先引入一种特殊形式的矩阵(elementary lower-triangular matrix)：
$${\mathbf{T}}_k=\mathbf{I}-{\mathbf{c}}_k{\mathbf{e}}_k^T,$$
其中
$${\mathbf{c}}_k=\left(\begin{array}{c}0\\ 0\\ ⋮\\ {\mu }_{k+1}\\ ⋮\\ {\mu }_n\end{array}\right),$$
那么
$${\mathbf{T}}_k=\left(\begin{array}{ccccccc}1& 0& \cdots & 0& 0& \cdots & 0\\ 0& 1& \cdots & 0& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& & ⋮\\ 0& 0& \cdots & 1& 0& \cdots & 0\\ 0& 0& \cdots & -{\mu }_{k+1}& 1& \cdots & 0\\ ⋮& ⋮& & ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & -{\mu }_n& 0& \cdots & 1\end{array}\right).$$
观察到${\mathbf{e}}_k^T{\mathbf{c}}_k=0$，可以得到
$${\mathbf{T}}_k^{-1}=\mathbf{I}+{\mathbf{c}}_k{\mathbf{e}}_k^T=\left(\begin{array}{ccccccc}1& 0& \cdots & 0& 0& \cdots & 0\\ 0& 1& \cdots & 0& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& & ⋮\\ 0& 0& \cdots & 1& 0& \cdots & 0\\ 0& 0& \cdots & {\mu }_{k+1}& 1& \cdots & 0\\ ⋮& ⋮& & ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\mu }_n& 0& \cdots & 1\end{array}\right),$$
它也属于elementary lower-triangular matrix. 引入这个矩阵的目的是用它表示消除第$k$个pivot下的元素时需要的所有$\text{Type }\mathrm{I}\mathrm{I}\mathrm{I}$型行变换操作。具体来说，若经$k-1$步消元操作后的部分三角化矩阵${\mathbf{A}}_{k-1}$是如下形式：
$${\mathbf{A}}_{k-1}=\left(\begin{array}{ccccccc}\ast & \ast & \cdots & {\alpha }_1& \ast & \cdots & \ast \\ 0& \ast & \cdots & {\alpha }_2& \ast & \cdots & \ast \\ ⋮& ⋮& \ddots & ⋮& ⋮& & ⋮\\ 0& 0& \cdots & {\alpha }_k& \ast & \cdots & \ast \\ 0& 0& \cdots & {\alpha }_{k+1}& \ast & \cdots & \ast \\ ⋮& ⋮& & ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\alpha }_n& \ast & \cdots & \ast \end{array}\right),$$
那么
$$\begin{array}{rl}{\mathbf{T}}_k{\mathbf{A}}_{k-1}& =(\mathbf{I}-{\mathbf{c}}_k{\mathbf{e}}_k^T){\mathbf{A}}_{k-1}={\mathbf{A}}_{k-1}-{\mathbf{c}}_k{\mathbf{e}}_k^T{\mathbf{A}}_{k-1}\\ & =\left(\begin{array}{ccccccc}\ast & \ast & \cdots & {\alpha }_1& \ast & \cdots & \ast \\ 0& \ast & \cdots & {\alpha }_2& \ast & \cdots & \ast \\ ⋮& ⋮& \ddots & ⋮& ⋮& & ⋮\\ 0& 0& \cdots & {\alpha }_k& \ast & \cdots & \ast \\ 0& 0& \cdots & 0& \ast & \cdots & \ast \\ ⋮& ⋮& & ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 0& \ast & \cdots & \ast \end{array}\right),\text{where }{\mathbf{c}}_k=\left(\begin{array}{c}0\\ 0\\ ⋮\\ {\alpha }_{k+1}/{\alpha }_k\\ ⋮\\ {\alpha }_n/{\alpha }_k\end{array}\right)\end{array}$$
可以看到第$k$个pivot ${\alpha }_k$下的元素已经都清零了。

因此，如果不需要进行行交换的话，高斯消元的过程等价于左乘一系列${\mathbf{T}}_i$，也即${\mathbf{T}}_{n-1}\cdots {\mathbf{T}}_2{\mathbf{T}}_1\mathbf{A}=\mathbf{U}$，因此
$$\mathbf{A}={\mathbf{T}}_1^{-1}{\mathbf{T}}_2^{-1}\cdots {\mathbf{T}}_{n-1}^{-1}\mathbf{U}.$$
注意到$j\le k$时${\mathbf{e}}_j^T{\mathbf{c}}_k=0$，有
$$\begin{array}{rl}{\mathbf{T}}_1^{-1}{\mathbf{T}}_2^{-1}\cdots {\mathbf{T}}_{n-1}^{-1}& =(\mathbf{I}+{\mathbf{c}}_1{\mathbf{e}}_1^T)(\mathbf{I}+{\mathbf{c}}_2{\mathbf{e}}_2^T)\cdots (\mathbf{I}+{\mathbf{c}}_{n-1}{\mathbf{e}}_{n-1}^T)\\ & =\mathbf{I}+{\mathbf{c}}_1{\mathbf{e}}_1^T+{\mathbf{c}}_2{\mathbf{e}}_2^T+\cdots +{\mathbf{c}}_{n-1}{\mathbf{e}}_{n-1}^T.\end{array}$$
由于
$${\mathbf{c}}_k{\mathbf{e}}_k^T=\left(\begin{array}{ccccccc}0& 0& \cdots & 0& 0& \cdots & 0\\ 0& 0& \cdots & 0& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& & ⋮\\ 0& 0& \cdots & 0& 0& \cdots & 0\\ 0& 0& \cdots & l_{k+1,k}& 0& \cdots & 0\\ ⋮& ⋮& & ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & l_{n,k}& 0& \cdots & 0\end{array}\right),$$
其中$l_{ik}$是第$k$步消元时需要乘的系数。

那么有$\mathbf{A}=\mathbf{L}\mathbf{U}$，其中
$$\mathbf{L}=\mathbf{I}+{\mathbf{c}}_1{\mathbf{e}}_1^T+{\mathbf{c}}_2{\mathbf{e}}_2^T+\cdots +{\mathbf{c}}_{n-1}{\mathbf{e}}_{n-1}^T=\left(\begin{array}{ccccc}1& 0& 0& \cdots & 0\\ l_{21}& 1& 0& \cdots & 0\\ l_{31}& l_{32}& 1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ l_{n1}& l_{n2}& l_{n3}& \cdots & 1\end{array}\right).$$
注：也可以不引入elementary lower-triangular matrix直接一步步按高斯消元法左乘一系列$\text{Type }\mathrm{I}\mathrm{I}\mathrm{I}$型矩阵${\mathbf{G}}_{ji}=\mathbf{I}-{\mu }_i{\mathbf{e}}_j{\mathbf{e}}_i^T$(第$j$行=第$j$行-${\mu }_i$第$i$行)，可以发现${\mathbf{T}}_i={\mathbf{G}}_{ni}\cdots {\mathbf{G}}_{2i}{\mathbf{G}}_{1i}$，也即${\mathbf{T}}_i$表达了将第$i$个pivot下的元素清零需要的所有步骤。

例：对$\mathbf{A}=\left(\begin{array}{ccc}2& 2& 2\\ 4& 7& 7\\ 6& 18& 22\end{array}\right)$进行LU分解

原先$0$的位置可以用来记录系数$l_{ik}$，用加粗表示：
$$\begin{gathered} \left(\begin{array}{ccc}2& 2& 2\\ 4& 7& 7\\ 6& 18& 22\end{array}\right)\stackrel{R_2-2R_1 \\ {R}_3-3R_1}{⟶}\left(\begin{array}{ccc}2& 2& 2\\ 2& 3& 3\\ 3& 12& 16\end{array}\right)\stackrel{R_3-4R_2}{⟶}\left(\begin{array}{ccc}2& 2& 2\\ 2& 3& 3\\ 3& 4& 4\end{array}\right). \end{gathered}$$
有
$$\mathbf{L}=\left(\begin{array}{ccc}1& 0& 0\\ 2& 1& 0\\ 3& 4& 1\end{array}\right)\text{ and }\mathbf{U}=\left(\begin{array}{ccc}2& 2& 2\\ 0& 3& 3\\ 0& 0& 4\end{array}\right)$$
这也意味着LU分解不需要占用比原矩阵$\mathbf{A}$更多的内存。

LU Factorization

If $\mathbf{A}$ is an $n\times n$ matrix such that a zero pivot is never encountered when applying Gaussian elimination with Type III operations, then $\mathbf{A}$ can be factored as the product $\mathbf{A}=\mathbf{L}\mathbf{U}$, where the following hold.

• $\mathbf{L}$ is lower triangular and $\mathbf{U}$ is upper triangular.
• $l_{ii}=1$ and $u_{ii}\ne 0$ for each $i=1,2,...,n.$
• Below the diagonal of $\mathbf{L}$, the entry $l_{ij}$ is the multiple of row $j$ that is subtracted from row $i$ in order to annihilate the $(i,j)$ -position during Gaussian elimination.
• $\mathbf{U}$ is the final result of Gaussian elimination applied to $\mathbf{A}$.
• The matrices $\mathbf{L}$ and $\mathbf{U}$ are uniquely determined.

The decomposition of $\mathbf{A}$ into $\mathbf{A}=\mathbf{L}\mathbf{U}$ is called the LU factorization of $\mathbf{A}$, and the matrices $\mathbf{L}$ and $\mathbf{U}$ are called the LU factors of $\mathbf{A}$.

关于LU分解唯一性的证明

由于$\mathbf{L},\mathbf{U}$对角线元素都是非零的下/上三角矩阵，所以$\mathbf{L},\mathbf{U}$非奇异。假设存在${\mathbf{L}}_1{\mathbf{U}}_1=\mathbf{A}={\mathbf{L}}_2{\mathbf{U}}_2$，那么
$${\mathbf{L}}_2^{-1}{\mathbf{L}}_1={\mathbf{U}}_2{\mathbf{U}}_1^{-1}.$$
下/上三角矩阵的逆和乘积依然是下/上三角矩阵，因此${\mathbf{L}}_2^{-1}{\mathbf{L}}_1$是下三角矩阵，${\mathbf{U}}_2{\mathbf{U}}_1^{-1}$是上三角矩阵。要使两者相等，则必须有${\mathbf{L}}_2^{-1}{\mathbf{L}}_1={\mathbf{U}}_2{\mathbf{U}}_1^{-1}=\mathbf{I}$. 因此${\mathbf{L}}_1={\mathbf{L}}_2,{\mathbf{U}}_1={\mathbf{U}}_2$.

一旦矩阵$\mathbf{A}$的LU分解确定了，会极大简化线性方程组$\mathbf{A}\mathbf{x}=\mathbf{b}$的求解

• To solve a nonsingular system $\mathbf{A}\mathbf{x}=\mathbf{b}$ using the LU factorization $\mathbf{A}=\mathbf{L}\mathbf{U}$, first solve $\mathbf{L}\mathbf{y}=\mathbf{b}$ for $\mathbf{y}$ with the forward substitution algorithm , and then solve $\mathbf{U}\mathbf{x}=\mathbf{y}$ for $\mathbf{x}$ with the back substitution procedure.
• The advantage of this approach is that once the LU factors for $\mathbf{A}$ have been computed, any other linear system $\mathbf{A}\mathbf{x}=\tilde{\mathbf{b}}$ can be solved with only $n^2$ multiplications/divisions and $n^2-n$ additions/subtractions.

例：说明LU分解可以简化计算${\mathbf{A}}^{-1}$的过程

求${\mathbf{A}}^{-1}$也即求$\mathbf{A}\mathbf{X}=\mathbf{I}$. 注意到$\mathbf{A}\mathbf{X}=\mathbf{A}[{\mathbf{x}}_1{\mathbf{x}}_2\cdots {\mathbf{x}}_n]=[{\mathbf{e}}_1{\mathbf{e}}_2\cdots {\mathbf{e}}_n]$，则问题可以转化为分别求$\mathbf{A}{\mathbf{x}}_j={\mathbf{L}\mathbf{U}\mathbf{x}}_j={\mathbf{e}}_j$，这可以通过两步解决

1. 令${\mathbf{y}}_j={\mathbf{U}\mathbf{x}}_j$，然后通过forward substitution求${\mathbf{L}\mathbf{y}}_j={\mathbf{e}}_j$.
2. 接着通过back substitution求解${\mathbf{U}\mathbf{x}}_j={\mathbf{y}}_j$

---

不是所有矩阵都能进行LU分解

Existence of LU Factors

Each of the following statements is equivalent to saying that a nonsingular matrix ${\mathbf{A}}_{n\times n}$ possesses an LU factorization.

• Each leading principal submatrix ${\mathbf{A}}_k$ is nonsingular.
• A zero pivot does not emerge during row reduction to upper-triangular form with $\text{Type }\mathrm{I}\mathrm{I}\mathrm{I}$ operations.

关于第一点的证明：

a nonsingular matrix ${\mathbf{A}}_{n\times n}$ possesses an LU factorization $⟹$ each leading principal submatrix ${\mathbf{A}}_k$ is nonsingular

假设$\mathbf{A}$有LU分解：
$$\mathbf{A}=\mathbf{L}\mathbf{U}=\left(\begin{array}{cc}{\mathbf{L}}_{11}& 0\\ {\mathbf{L}}_{21}& {\mathbf{L}}_{22}\end{array}\right)\left(\begin{array}{cc}{\mathbf{U}}_{11}& {\mathbf{U}}_{12}\\ 0& {\mathbf{U}}_{22}\end{array}\right)=\left(\begin{array}{cc}{\mathbf{A}}_k& \ast \\ \ast & \ast \end{array}\right),$$
其中${\mathbf{L}}_{11}$和${\mathbf{U}}_{11}$都是$k\times k$的矩阵。那么有${\mathbf{A}}_k={\mathbf{L}}_{11}{\mathbf{U}}_{11}$，因为${\mathbf{L}}_{11}$和${\mathbf{U}}_{11}$都是非奇异的，${\mathbf{A}}_k$也是非奇异的。

a nonsingular matrix ${\mathbf{A}}_{n\times n}$ possesses an LU factorization $⟸$ each leading principal submatrix ${\mathbf{A}}_k$ is nonsingular

首先用数学归纳法进行证明每个leading principal submatrix都可以进行LU分解：

对$k=1$, 若${\mathbf{A}}_1=(a_{11})$是非奇异的，它存在LU分解${\mathbf{A}}_1=(1)(a_{11})$;

假设${\mathbf{A}}_k$存在LU分解${\mathbf{A}}_k={\mathbf{L}}_k{\mathbf{U}}_k$，那么
$${\mathbf{A}}_{k+1}=\left(\begin{array}{cc}{\mathbf{A}}_k& \mathbf{b}\\ {\mathbf{c}}^T& {\alpha }_{k+1}\end{array}\right)=\left(\begin{array}{cc}{\mathbf{L}}_k& 0\\ {\mathbf{c}}^T{\mathbf{U}}_k^{-1}& 1\end{array}\right)\left(\begin{array}{cc}{\mathbf{U}}_k& {\mathbf{L}}_k^{-1}\mathbf{b}\\ 0& {\alpha }_{k+1}-{\mathbf{c}}^T{\mathbf{A}}_k^{-1}\mathbf{b}\end{array}\right)\triangleq {\mathbf{L}}_{k+1}{\mathbf{U}}_{k+1},$$
因为${\mathbf{A}}_{k+1}$和${\mathbf{L}}_{k+1}$都是非奇异的，有${\mathbf{U}}_{k+1}$是非奇异的，即${\alpha }_{k+1}-{\mathbf{c}}^T{\mathbf{A}}_k^{-1}\mathbf{b}\ne 0.$ 因此${\mathbf{A}}_{k+1}$存在LU分解。

${\mathbf{A}}_n$本身也是一个leading principal submatrix，因此可以进行LU分解。

关于第二点的证明：

a nonsingular matrix ${\mathbf{A}}_{n\times n}$ possesses an LU factorization $⟹$ all pivots must be nonzero

首先根据
$${\mathbf{A}}_{k+1}=\left(\begin{array}{cc}{\mathbf{L}}_k& 0\\ {\mathbf{c}}^T{\mathbf{U}}_k^{-1}& 1\end{array}\right)\left(\begin{array}{cc}{\mathbf{U}}_k& {\mathbf{L}}_k^{-1}\mathbf{b}\\ 0& {\alpha }_{k+1}-{\mathbf{c}}^T{\mathbf{A}}_k^{-1}\mathbf{b}\end{array}\right)$$
得到第$(k+1)$个pivot：$p_{k+1}={\alpha }_{k+1}-{\mathbf{c}}^T{\mathbf{A}}_k^{-1}\mathbf{b}$.

借用第一点的结论：a nonsingular matrix ${\mathbf{A}}_{n\times n}$ possesses an LU factorization $⟹$ each leading principal submatrix ${\mathbf{A}}_k$ is nonsingular

已知${\mathbf{A}}_{k+1}$和${\mathbf{L}}_{k+1}$都非奇异，有${\mathbf{U}}_{k+1}$非奇异，$p_{k+1}\ne 0$.

a nonsingular matrix ${\mathbf{A}}_{n\times n}$ possesses an LU factorization $⟸$ all pivots must be nonzero

也即不需要进行行交换，满足从高斯消元过程中得到LU分解的前提条件

---

如果出现行交换应该怎么分解？

回忆之前的$\text{elementary lower-triangular matrix }{\mathbf{T}}_k=\mathbf{I}-{\mathbf{c}}_k{\mathbf{e}}_k^T$和$\text{Type }\mathrm{I}\text{ elementary interchange matrix }\mathbf{E}=\mathbf{I}-\mathbf{u}{\mathbf{u}}^{\mathrm{T}}\text{ with }\mathbf{u}={\mathbf{e}}_{\mathrm{k}+\mathrm{i}}-{\mathbf{e}}_{\mathrm{k}+\mathrm{j}}$，不难发现${\mathbf{E}}^2=\mathbf{I}$. 另外，有
$${\tilde{\mathbf{T}}}_k\triangleq \mathbf{E}{\mathbf{T}}_k\mathbf{E}={\mathbf{E}}^2-\mathbf{E}{\mathbf{c}}_k{\mathbf{e}}_k^T\mathbf{E}=\mathbf{I}-{\tilde{\mathbf{c}}}_k{\mathbf{e}}_k^T,$$
其中${\tilde{\mathbf{c}}}_k=\mathbf{E}{\mathbf{c}}_k$.

可以看出${\tilde{\mathbf{T}}}_k$也是一个elementary lower-triangular matrix，而且与${\mathbf{T}}_k$的差别仅仅在于${\mu }_{k+i}$和${\mu }_{k+j}$交换了一下位置。

现在我们假设在第$k$行高斯消元结束后必须进行一次行交换，也即用$\mathbf{E}{\mathbf{T}}_k\cdots {\mathbf{T}}_1$左乘$\mathbf{A}$，利用${\mathbf{E}}^2=\mathbf{I}$，有
$$\begin{array}{rl}\mathbf{E}{\mathbf{T}}_k{\mathbf{T}}_{k-1}\cdots {\mathbf{T}}_1& =\mathbf{E}{\mathbf{T}}_k{\mathbf{E}}^2{\mathbf{T}}_{k-1}{\mathbf{E}}^2\cdots {\mathbf{E}}^2{\mathbf{T}}_1{\mathbf{E}}^2\\ & =(\mathbf{E}{\mathbf{T}}_k\mathbf{E})(\mathbf{E}{\mathbf{T}}_{k-1}\mathbf{E})\cdots (\mathbf{E}{\mathbf{T}}_1\mathbf{E})\mathbf{E}\\ & ={\tilde{\mathbf{T}}}_k{\tilde{\mathbf{T}}}_{k-1}\cdots {\tilde{\mathbf{T}}}_1\mathbf{E}.\end{array}$$
可以发现$\mathbf{E}$被变换到了最右端，且${\tilde{\mathbf{T}}}_k{\tilde{\mathbf{T}}}_{k-1}\cdots {\tilde{\mathbf{T}}}_1$与${\mathbf{T}}_k{\mathbf{T}}_{k-1}\cdots {\mathbf{T}}_1$的区别仅在于第$k+i$和$k+j$行被交换了。

把高斯消元过程执行完，整体可以写成${\tilde{\mathbf{T}}}_{n-1}\cdots {\tilde{\mathbf{T}}}_2{\tilde{\mathbf{T}}}_1\mathbf{P}\mathbf{A}=\mathbf{U}$，其中$\mathbf{P}$是所有$\text{Type }\mathrm{I}$型行交换矩阵的乘积。令$\mathbf{L}={\tilde{\mathbf{T}}}_1^{-1}{\tilde{\mathbf{T}}}_2^{-2}\cdots {\tilde{\mathbf{T}}}_{n-1}^{-1}$，有$\mathbf{P}\mathbf{A}=\mathbf{L}\mathbf{U}$.

也即，对任意非奇异矩阵$\mathbf{A}$，存在置换矩阵$\mathbf{P}$使得$\mathbf{P}\mathbf{A}$存在LU分解。

（上面的推理过程也可以理解为，预先对$\mathbf{A}$做行交换，使交换后的矩阵所有pivot都不为零）

例：对矩阵$\mathbf{A}=\left(\begin{array}{cccc}1& 2& -3& 4\\ 4& 8& 12& -8\\ 2& 3& 2& 1\\ -3& -1& 1& -4\end{array}\right)$，确定$\mathbf{P}\mathbf{A}$的LU分解，其中$\mathbf{P}$是对应的置换阵

为了容易看出行交换过程，通常会加入permutation counter column $\mathbf{p}$，它一开始是$1,2,3,4$的自然排列：（下面的过程有一些交换不是必做的）
$$\begin{array}{rl}\left[\mathbf{A}|\mathbf{p}\right]=& \left(\begin{array}{cccc}1& 2& -3& 4\\ 4& 8& 12& -8\\ 2& 3& 2& 1\\ -3& -1& 1& -4\end{array}\right|\begin{array}{c}1\\ 2\\ 3\\ 4\end{array})⟶\left(\begin{array}{cccc}4& 8& 12& -8\\ 1& 2& -3& 4\\ 2& 3& 2& 1\\ -3& -1& 1& -4\end{array}\right|\begin{array}{c}2\\ 1\\ 3\\ 4\end{array})\\ ⟶& \left(\begin{array}{cccc}4& 8& 12& -8\\ 1/4& 0& -6& 6\\ 1/2& -1& -4& 5\\ -3/4& 5& 10& -10\end{array}\right|\begin{array}{c}2\\ 1\\ 3\\ 4\end{array})⟶\left(\begin{array}{cccc}4& 8& 12& -8\\ -3/4& 5& 10& 10\\ 1/2& -1& -4& 5\\ 1/4& 0& -6& 6\end{array}\right|\begin{array}{c}2\\ 4\\ 3\\ 1\end{array})\\ ⟶& \left(\begin{array}{cccc}4& 8& 12& -8\\ -3/4& 5& 10& 10\\ 1/2& -1/5& -2& 3\\ 1/4& 0& -6& 6\end{array}\right|\begin{array}{c}2\\ 4\\ 3\\ 1\end{array})⟶\left(\begin{array}{cccc}4& 8& 12& -8\\ -3/4& 5& 10& 10\\ 1/4& 0& -6& 6\\ 1/2& -1/5& -2& 3\end{array}\right|\begin{array}{c}2\\ 4\\ 1\\ 3\end{array})\\ ⟶& \left(\begin{array}{cccc}4& 8& 12& -8\\ -3/4& 5& 10& 10\\ 1/4& 0& -6& 6\\ 1/2& -1/5& 1/3& 1\end{array}\right|\begin{array}{c}2\\ 4\\ 1\\ 3\end{array}).\end{array}$$
因此有
$$\mathbf{L}=\left(\begin{array}{cccc}1& 0& 0& 0\\ -3/4& 1& 0& 0\\ 1/4& 0& 1& 0\\ 1/2& -1/5& 1/3& 1\end{array}\right),\mathbf{U}=\left(\begin{array}{cccc}4& 8& 12& -8\\ 0& 5& 10& 10\\ 0& 0& -6& 6\\ 0& 0& 0& 1\end{array}\right),\mathbf{P}=\left(\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 0& 1\\ 1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right).$$

LU Factorization with Row Interchanges

• For each nonsingular matrix $\mathbf{A}$, there exists a permutation matrix $\mathbf{P}$ such that $\mathbf{P}\mathbf{A}$ possesses an LU factorization $\mathbf{P}\mathbf{A}=\mathbf{L}\mathbf{U}$.
• To compute $\mathbf{L}$, $\mathbf{U}$, and $\mathbf{P}$, successively overwrite the array originally containing $\mathbf{A}$. Replace each entry being annihilated with the multiplier used to execute the annihilation. Whenever row interchanges such as those used in partial pivoting are implemented, the multipliers in the array will automatically be interchanged in the correct manner.
• Although the entire permutation matrix $\mathbf{P}$ is rarely called for, it can be constructed by permuting the rows of the identity matrix $\mathbf{I}$ according to the various interchanges used. These interchanges can be accumulated in a “permutation counter column” $\mathbf{p}$ that is initially in natural order $(1,2,...,n)$.
• To solve a nonsingular linear system $\mathbf{A}\mathbf{x}=\mathbf{b}$ using the LU decomposition with partial pivoting, permute the components in $\mathbf{b}$ to construct $\tilde{\mathbf{b}}$ according to the sequence of interchanges used—i.e.,
  according to $\mathbf{p}$ —and then solve $\mathbf{L}\mathbf{y}=\tilde{\mathbf{b}}$ by forward substitution followed by the solution of $\mathbf{U}\mathbf{x}=\mathbf{y}$ using back substitution.

---

The LDU factorization

在LU分解中，$\mathbf{L}$的对角线上是全$1$的，但$\mathbf{U}$并不是，这在结构上缺乏对称性。我们可以将原来的$\mathbf{U}$写为：
$$\left(\begin{array}{cccc}{u}_{11}& u_{12}& \cdots & u_{1n}\\ 0& u_{22}& \cdots & u_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & u_{nn}\end{array}\right)=\left(\begin{array}{cccc}{u}_{11}& 0& \cdots & 0\\ 0& u_{22}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & u_{nn}\end{array}\right)\left(\begin{array}{cccc}1& u_{12}/u_{11}& \cdots & u_{1n}/u_{11}\\ 0& 1& \cdots & u_{2n}/u_{22}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 1\end{array}\right),$$
令$\mathbf{D}=diag(u_{11},u_{22},...,u_{nn})$，得到LDU分解。由LU分解是唯一确定的，容易得到LDU分解也是唯一确定的。而且当$\mathbf{A}$是对称矩阵时，有$\mathbf{A}=\mathbf{L}\mathbf{D}{\mathbf{L}}^T$（根据唯一性容易得到）.

---

The Cholesky Factorization

A symmetric matrix $\mathbf{A}$ possessing an LU factorization in which each pivot is positive is said to be positive definite.

Prove that $\mathbf{A}$ is positive definite if and only if $\mathbf{A}$ can be uniquely factored as $\mathbf{A}={\mathbf{R}}^T\mathbf{R}$, where $\mathbf{R}$ is an upper-triangular matrix with positive diagonal entries.

This is known as the Cholesky factorization of $\mathbf{A}$, and $\mathbf{R}$ is called the Cholesky factor of $\mathbf{A}$.

$⟹$

若$\mathbf{A}$是正定的，那么存在LDU分解$\mathbf{A}=\mathbf{L}\mathbf{D}{\mathbf{L}}^T$，其中$\mathbf{D}=diag(p_1,p_2,\cdots ,p_n)$，$p_i>0$. 取$\mathbf{R}={\mathbf{D}}^{1/2}{\mathbf{L}}^T$即有$\mathbf{A}={\mathbf{R}}^T\mathbf{R}$，$\mathbf{R}$是上三角矩阵且对角元素均为正。

$⟸$

设$\mathbf{A}=\mathbf{R}{\mathbf{R}}^T$，$\mathbf{R}$是对角元素均为正的下三角矩阵，则它可以被分解为$\mathbf{R}=\mathbf{L}\mathbf{D}$，其中$\mathbf{L}$是对角线元素都为$1$的下三角矩阵，$\mathbf{D}$是对角线元素为$r_{ii}$的对角阵。因此$\mathbf{A}=\mathbf{L}{\mathbf{D}}^2{\mathbf{L}}^T$是$\mathbf{A}$的LDU分解，且所有pivot均为正值

唯一性可以从LDU分解的唯一性得到

---

Examples

If $\mathbf{A}$ is a matrix that contains only integer entries and all of its pivots are $1$, explain why ${\mathbf{A}}^{-1}$ must also be an integer matrix.

Note: This fact can be used to construct random integer matrices that possess integer inverses by randomly generating integer matrices $\mathbf{L}$ and $\mathbf{U}$ with unit diagonals and then constructing the product $\mathbf{A}=\mathbf{L}\mathbf{U}$.

由于所有的pivot都等于1，存在LU分解$\mathbf{A}=\mathbf{L}\mathbf{U}$，其中$\mathbf{L}$和$\mathbf{U}$分别是对角线全为1的下三角和上三角整数矩阵。对于这样的矩阵，它们的逆也是整数矩阵（因为在$[\mathbf{L}|\mathbf{I}]\stackrel{Gauss-Jordan}{⟶}[\mathbf{I}|{\mathbf{L}}^{-1}]$的过程中，每一次$\text{Type }\mathrm{I}\mathrm{I}$和$\text{Type }\mathrm{I}\mathrm{I}\mathrm{I}$型的操作都是乘以整数系数）。因此${\mathbf{A}}^{-1}={\mathbf{U}}^{-1}{\mathbf{L}}^{-1}$也是整数矩阵。

Consider the tridiagonal matrix $\mathbf{T}=\left(\begin{array}{cccc}{\beta }_1& {\gamma }_1& 0& 0\\ {\alpha }_1& {\beta }_2& {\gamma }_2& 0\\ 0& {\alpha }_2& {\beta }_3& {\gamma }_3\\ 0& 0& {\alpha }_3& {\beta }_4\end{array}\right)$. Assume that $\mathbf{T}$ possesses an LU fractorization, verify that it is given by
$$\mathbf{L}=\left(\begin{array}{cccc}1& 0& 0& 0\\ {\alpha }_1/{\pi }_1& 1& 0& 0\\ 0& {\alpha }_2/{\pi }_2& 1& 0\\ 0& 0& {\alpha }_3/{\pi }_3& 1\end{array}\right),\mathbf{U}=\left(\begin{array}{cccc}{\pi }_1& {\gamma }_1& 0& 0\\ 0& {\pi }_2& {\gamma }_2& 0\\ 0& 0& {\pi }_3& {\gamma }_3\\ 0& 0& 0& {\pi }_4\end{array}\right),\text{ where}$$

$${\pi }_1={\beta }_1\text{and}{\pi }_{i+1}={\beta }_{i+1}-\frac{{\alpha }_i{\gamma }_i}{{\pi }_i}.$$

代入验证即可，对任意大小的tridiagonal矩阵都成立。