高斯消去法

置换矩阵（permutaion matrix)： 一个单位矩阵的行重新排列后形成的矩阵。

$$\begin{pmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{pmatrix}$$

Lemma 2.2: 转换矩阵的性质

1. $PX$ 和 $XP$ 的意义
2. $P^{-1} = P^T$
3. $det(P) = \pm 1$
4. $P_1\cdot P_2$ 仍然是转换矩阵。

Alg 2.1: 高斯法求解 $Ax = b$ 的流程如下：

1. 分解 $A = PLU$
   • $P$ 为转换矩阵。
   • $L$ 为单位下三角矩阵。
   • $U$ 非奇异的上三角矩阵。
2. $LUx = P^{-1}b = P^Tb$.
3. $Ux = L^{-1}(P^{-1}b)$
4. $x = U^{-1}(L^{-1}P^{-1}b)$

定义 2.2：顺序主子矩阵 第 $j$ 个顺序主子矩阵为 $A(1:j，1:j)$。

定理 2.4 下面的两个声明等价：

1. $A = LU$ 分解是唯一的。
2. $A$ 的所有顺序主子矩阵都是非奇异的。

证明分析：

1. 从分解唯一出发证明，对矩阵 $A$ 进行分块，分出任意一个顺序主子矩阵。
2. 从顺序主子矩阵非奇异出发，归纳法证明，假设 $n-1\times n-1$ 的矩阵可分解

定理 2.5 对于非奇异的矩阵 $A$， 一定存在如下分解：

$$P_1AP_2 = LU$$

其中 $P_1$ 和 $P_2$ 只有一个是必须的。

证明分析： 用归纳法证明

1) 对于 $1\times 1$ 矩阵 $A$ ， $P_1=P_2=L=1$ 和 $U = A$ 。 2) 假设对于 $n-1\times n-1$ 矩阵， 上述结论成立。 3) 下面证明对 $n\times n$ 矩阵成立，首先利用 $A$ 非奇异 ，把问题转化为 $n-1\times n-1$ 的问题。

$$\begin{aligned} P_{1}^{\prime} A P_{2}^{\prime} &=\left[ \begin{array}{cc}{a_{11}} & {A_{12}} \\ {A_{21}} & {A_{22}}\end{array}\right] &= \left[ \begin{array}{cc}{1} & {0} \\ {L_{21}} & {I}\end{array}\right] \cdot \left[ \begin{array} {cc}{u_{11}} & {U_{12}} \\ {0} & {\tilde{A}_{22}} \end{array} \right] \\ &=\left[ \begin{array}{cc}{u_{11}} & {U_{12}} \\ {L_{21} u_{11}} & {L_{21} U_{12}+\tilde{A}_{22}}\end{array}\right] \end{aligned}$$

比较等式两边的系数， 可得：

• $u_{11}=a_{11}$
• $U_{12}=A_{12}$
• $L_{21} = \frac{A_{21}}{u_{11}} = \frac{A_{21}}{a_{11}}$
• $L_{21}U_{12}+\tilde A_{22}=A_{22}$ 推出 $\tilde A_{22} = A_{22} - L_{21}U_{12} = A_{22} - \frac{A_{21}}{a_{11}}A_{12}$
• $\tilde A_{22}$ 是非奇异的, 验证如下： $$\operatorname{det} P_{1}^{\prime} A P_{2}^{\prime}=\operatorname{det} \left[ \begin{array}{cc}{1} & {0} \\ {L_{21}} & {I}\end{array}\right] \cdot \operatorname{det} \left[ \begin{array}{cc}{u_{11}} & {U_{12}} \\ {0} & {\tilde{A}_{22}}\end{array}\right]=1 \cdot\left(u_{11} \cdot \operatorname{det} \tilde{A}_{22}\right)$$
• 上述转化过程实际上给出的高斯消去法的主要算法步骤。

4) $n-1 \times n-1$ 的分解假设存在

$$\tilde P_1 \tilde{A}_{22} \tilde P_2 = \tilde L\tilde U$$

最后推出 $A$ 的分解

$$\begin{align} P_1'AP_2' &= \begin{bmatrix} 1 & 0 \\ L_{21} & I \end{bmatrix} \begin{bmatrix} u_{11} & U_{12} \\ 0 & \tilde P_1^T \tilde L \tilde U \tilde P_2^T \end{bmatrix}\\ &= \begin{bmatrix} 1 & 0 \\ L_{21} & I \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0 & \tilde P_1^T\tilde L \end{bmatrix} \begin{bmatrix} u_{11} & U_{12} \\ 0 & \tilde U\tilde P_2^T \end{bmatrix}\\ &= \begin{bmatrix} 1 & 0 \\ L_{21} & \tilde P_1^T\tilde L \end{bmatrix} \begin{bmatrix} u_{11} & U_{12}\tilde P_2 \\ 0 & \tilde U \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0 & \tilde P_2^T \end{bmatrix}\\ &= \begin{bmatrix} 1 & 0 \\ 0 & \tilde P_1^T \end{bmatrix} \begin{bmatrix} 1 & 0 \\ \tilde P_1 L_{21} & \tilde L \end{bmatrix} \begin{bmatrix} u_{11} & U_{12}\tilde P_2 \\ 0 & \tilde U \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0 & \tilde P_2^T \end{bmatrix}, \end{align}$$

最终得到 $A$ 的分解：

$$\begin{align} P_1AP_2 &= \left( \begin{bmatrix} 1 & 0 \\ 0 & \tilde P_1 \end{bmatrix} P_1' \right) A \left( P_2' \begin{bmatrix} 1 & 0 \\ 0 & \tilde P_2 \end{bmatrix} \right)\\ &= \begin{bmatrix} 1 & 0 \\ \tilde P_1 L_{21} & \tilde L \end{bmatrix} \begin{bmatrix} u_{11} & U_{12}\tilde P_2 \\ 0 & \tilde U \end{bmatrix} \end{align}$$ $\square$

推论 2.1（GEPP） 列主元高斯消去法

推论 2.2 （GECP） 完全主元高斯消去法

算法 2.2 主元 $LU$ 分解 $$\begin{align} &\text{for } i=1 \text{ to } n-1 \\ &\quad \text{apply permutation so } a_{ii} \not= 0 \\ &\quad \text{for } j = i + 1 \text{ to } n \\ &\quad\quad l_{ji} = a_{ji}/a_{ii}\\ &\quad \text{end for} \\ &\quad \text{for } j = i \text{ to } n \\ &\quad\quad u_{ij} = a_{ij} \\ &\quad \text{end for}\\ &\quad \text{for } j=i+1 \text{ to } n \\ &\quad\quad \text{for } k=i+1 \text{ to } n\\ &\quad\quad\quad a_{jk} = a_{jk} - l_{ji}*u_{ik} \\ &\quad\quad \text{end for} \\ &\quad \text{end for}\\ &\text{end for} \end{align}$$

算法 2.3 主元 $LU$ 分解，利用 $A$ 的存储空间存放 $L$ 和 $U$ $$\begin{align} &\text{for } i=1 \text{ to } n-1 \\ &\quad \text{apply permutations} \\ &\quad \text{for } j=i+1 \text{ to } n \\ &\quad\quad a_{ji} = a_{ji}/a_{ii} \\ &\quad \text{end for} \\ &\quad \text{for } j=i+1 \text{ to } n \\ &\quad\quad \text{for } k=i+1 \text{ to } n\\ &\quad\quad\quad a_{jk} = a_{jk} - a_{ji}*a_{ik} \\ &\quad\quad \text{end for} \\ &\quad \text{end for}\\ &\text{end for}\\ \end{align}$$

讲课笔记

• 用一个 4 阶的矩阵演示 $LU$ 分解的过程

Python 实现 下面给出上面算法的 Python 实现

def algorithm_2_2(A):
    n = len(A)
    P, L, U = np.eye(n), np.zeros((n, n)), np.zeros((n, n))
    for i in range(n - 1):
        m = np.argmax(np.abs(A[i:, i])) + i
        if A[m, i] == 0:
            raise ValueError("matrix is singular.")
        else:
            if m != i:
                A[[i, m], :] = A[[m, i], :]
                L[[i, m], :] = L[[m, i], :]
                P[[i, m], :] = P[[m, i], :]
            for j in range(i + 1, n):
                L[j, i] = A[j, i] / A[i, i]

            for j in range(i, n):
                U[i, j] = A[i, j]

            for j in range(i + 1, n):
                for k in range(i + 1, n):
                    A[j, k] -= L[j, i] * U[i, k]

    P = P.T
    L += np.eye(n)
    U[-1, -1] = A[-1, -1]
    return A, P, L, U

def algorithm_2_3(A):
    n = len(A)
    P = np.eye(n)
    for i in range(n - 1):
        m = np.argmax(np.abs(A[i:, i])) + i
        if A[m, i] == 0:
            raise ValueError("matrix is singular.")
        else:
            if m != i:
                A[[i, m], :] = A[[m, i], :]
                P[[i, m], :] = P[[m, i], :]

            for j in range(i + 1, n):
                A[j, i] = A[j, i] / A[i, i]

            for j in range(i + 1, n):
                for k in range(i + 1, n):
                    A[j, k] -= A[j, i] * A[i, k]

    P = P.T
    L = np.tril(A, -1) + np.eye(n)
    U = np.triu(A, 0)
    return A, P, L, U

def algorithm_2_4(A):
    n = len(A)
    P = np.eye(n)

    for i in range(n - 1):
        m = np.argmax(np.abs(A[i:, i])) + i

        if A[m, i] == 0:
            raise ValueError("matrix is singular.")
        else:
            if m != i:
                A[[i, m], :] = A[[m, i], :]
                P[[i, m], :] = P[[m, i], :]

            A[(i + 1):, i] /= A[i, i]
            A[(i + 1):, (i + 1):] -= A[(i + 1):, i].reshape(-1, 1)*A[i, (i + 1):]

    P = P.T
    L = np.tril(A, -1) + np.eye(n)
    U = np.triu(A, 0)
    return A, P, L, U

$LU$ 分解的算法复杂度

$$\begin{align} &\text{for } i=1 \text{ to } n-1 \\ &\quad \text{apply permutations} \\ &\quad \text{for } j=i+1 \text{ to } n \\ &\quad\quad l_{ji} = a_{ji}/a_{ii} \\ &\quad \text{end for} \\ &\quad \text{for } j=i+1 \text{ to } n \\ &\quad\quad \text{for } k=i+1 \text{ to } n\\ &\quad\quad\quad a_{jk} = a_{jk} - l_{ji}*a_{ik} \\ &\quad\quad \text{end for} \\ &\quad \text{end for}\\ &\text{end for}\\ \end{align}$$

把上面的 for 循环用 $\sum_{}^{}$ 替换，可得

$$\begin{align} &\sum_{i=1}^{n-1} \left( \sum_{j=i+1}^{n} 1 + \sum_{j=i+1}^{n} \sum_{k=i+1}^{n}2 \right) \\ =& \sum_{i=1}^{n-1} \left( ( n-i ) + 2(n-i)^2\right) \\ =& \sum_{i=1}^{n-1} \left( n-i + 2n^2 + 2i^2 - 4ni \right) \\ =& \sum_{i=1}^{n-1} ( n + 2n^2 ) + 2\sum_{i=1}^{n-1} i^2 - (4n+1)\sum_{i=1}^{n-1} i\\ =& ( n - 1 ) ( n + 2n^2 ) + \frac{(n-1)n(2n - 1)}{3} - ( 4n+1 )\frac{ n(n+1) }{ 2 }\\ =& \frac{2}{3}n^3 + O(n^2) \end{align}$$

$$1^{2}+2^{2}+\ldots +k^{2}=\frac{k(k+1)(2k+1)}{6}$$

$$1 + 2 + 3 + \ldots + k = \frac{ ( k+1 ) ( k+2 ) }{ 2 }$$

实验报告

给定一个 100 阶的非奇异矩阵，比较上述几种算法的效率。