biscuitの博客

高斯消元法#

分为两步

1. 前向消元：逐渐将第i行到第n行的 $x_j,j<i$消去，得到上三角矩阵 $A'$ 和变换后的向量 $b'$。
2. 回代求解：从最后一行开始，逐行向上求解未知数。

把系数矩阵变成上三角形的过程，等价于第 $j$ 步用初等变换把第 $j$ 列的主对角线下方元素消为零。

设 $A^{(0)} = A$ 是原始矩阵，$A^{(1)}, A^{(2)}, \dots, A^{(n-1)}$ 是每一步消元后的矩阵。每一步的初等变换都可以表示为一个初等矩阵 $E^{(j)}$，使得：

$E^{(j)}$ 是单位阵初等变换来的，唯一的区别就是他的j列j行的元素( $E^{(j)}_{jj}$ )被修改成 $-\dfrac{A^{(j-1)}_{ij}}{A^{(j-1)}_{jj}}$

定义一个向量 $\mathbf{l}_j$，其第 $i$ 个元素为 $-\dfrac{A^{(j-1)}_{ij}}{A^{(j-1)}_{jj}}$，其他元素为 $0$，则 $E^{(j)}$ 可以表示为：

其中 $\mathbf{e}_j$ 是第 $j$ 个标准基向量。

可能的问题：某一行的主元过小，不适合做分母,某一行的应该等于零的值在数值计算中很可能因为舍入误差变成非零，导致出错。

解决：主元法

在第k步前向消除时，选择第k列中绝对值最大的元素所在的行作为主元行，与第k行交换位置。这样可以避免数值不稳定的问题，提高算法的鲁棒性。

这个不是换一次就可以，而是每一步都需要选一次列主元，直到消元完成。

置换的操作可以通过一个置换矩阵 $P$ 来表示，$P$ 是一个单位矩阵经过行交换得到的矩阵。对于每一步的主元选择和行交换，我们都有一个对应的置换矩阵 $P^{(j)}$，使得：

所以总结下来完整的带主元法的高斯消元过程可以表示为：

矩阵的LU分解#

TIP

高斯消元的前向消元过程等价于矩阵分解

其中 $L$ 是一个单位下三角矩阵（对角线元素全为 1，上三角部分全为 0），$U$ 是一个上三角矩阵。通过这种分解，我们可以将求解 $A\mathbf{x} = \mathbf{b}$ 转化为两个简单的三角矩阵方程：

首先求解方程 (1) 来得到 $\mathbf{y}$，然后将 $\mathbf{y}$ 代入方程 (2) 来求解 $\mathbf{x}$。

LU分解的手动计算#

LU 分解的核心思想是将一个矩阵 $A$ 分解为一个单位下三角矩阵 $L$（Lower triangular，主对角线全为 1）和一个上三角矩阵 $U$（Upper triangular）。它的本质其实就是高斯消元法的矩阵表达。我们用一个简单的 $3 \times 3$ 矩阵来演示手算过程。

待分解矩阵示例假设我们要分解矩阵 $A$：

手算步骤我们在对 $A$ 进行 高斯消元 （按列顺序处理，从左到右。）得到 $U$ 的过程中，将消元所用的系数填入 $L$ 的对应位置。

系数 $l_{ij}$ 是在消元过程中用来消去 $A(U)$ 第 $i$ 行第 $j$ 列元素的系数，也是存入 $L(i,j)$ 的元素。

对于第 $j$ 列的消元，系数 $l_{ij}$ 的计算公式为：

• 选取第 $j$ 列的主元元素 $A_{jj}$。
• 对于第 $i$ 行（$i > j$），计算消元系数：$l_{ij} = \dfrac{A_{ij}}{A_{jj}}$。

然后用这个系数 $l_{ij}$ 来更新第 $i$ 行：

第一步：处理第一列，我们要把第 2 行和第 3 行的第一列元素化为 0。

$R_2 \leftarrow R_2 - (2)R_1$：消掉 4。系数 $l_{21} = 2$。

$R_3 \leftarrow R_3 - (-1)R_1$：消掉 -2。系数 $l_{31} = -1$。

变换后的矩阵：

第二步：处理第二列我们要把第 3 行的第二列元素化为 0。

$R_3 \leftarrow R_3 - (5)R_2$：消掉 5。系数 $l_{32} = 5$。

变换后的矩阵（这就是 $U$）：

构建矩阵 $L$矩阵 $L$ 的主对角线全是 1，其余位置填入我们刚才记录的系数 $l_{ij}$：

验证结果我们将 $L$ 和 $U$ 相乘，看是否回到 $A$：

💡 手算小技巧总结

保持 $L$ 的对角线为 1：这是 Doolittle 分解法的标准。

符号不要搞反：如果在消元时是 $R_i - k R_j$，那么 $L$ 矩阵对应位置存的就是这个 $k$；如果是 $R_i + k R_j$，存的就是 $-k$。

什么时候不能直接分解？：如果在消元过程中发现主元（对角线元素）为 0，普通的 LU 分解就会失效，这时候需要进行行交换，也就是引入置换矩阵 $P$，变为 $PA = LU$ 分解。

条件数的由来：误差e与残差r#

设 $\mathbf{x}$ 是 $A\mathbf{x} = \mathbf{b}$ 的一个精确解（难以得到），$\tilde{\mathbf{x}}$ 是数值计算得到的近似解，误差 $\mathbf{e}$ 和残差 $\mathbf{r}$ 定义如下，消去其中精确解的部分得到关系：

在实际计算中，单纯用误差和残差不能反映矩阵本身的稳定性，因此对绝对误差的上限的范数进行表示：

所以相对误差的上限为：

想办法把右边化成和精确解$\mathbf{x}$无关：

这说明相对误差的上限是条件数 $\kappa(A)$ 和相对残差 $\frac{||r||}{||\mathbf{b}||}$ 的乘积。因此，矩阵的条件数越大，即使残差很小，误差也可能很大，这就是矩阵不稳定的表现。

优化误差和残差：迭代#

求解迭代的基本步骤：

从一个初始猜测 $\mathbf{x}^{(0)}$ 开始，

对于第k步计算残差 $\mathbf{r}^{(k)} = \mathbf{b} - A\mathbf{x}^{(k)}$。

根据残差，更新修正量 $\mathbf{\delta x}^{(k)}$，使得 $A\mathbf{\delta x}^{(k)} \approx \mathbf{r}^{(k)}$。

用修正量更新解 $\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \mathbf{\delta x}^{(k)}$。

重复，直到残差足够小或者达到预设的迭代次数。

Hager算法：不求逆来估计逆的范数#

应该不会考，放着

优化条件数1：建模时进行无量纲化#

因为量纲问题导致大条件数的矩阵，一般可以改成