概述
对于 n n n
A x = b \mathbf A \mathbf x = \mathbf b Ax = b 其中 A ∈ R n × n \mathbf A\in\mathbb{R}^{n\times n} A ∈ R n × n x , b ∈ R n \mathbf x,\mathbf b\in\mathbb{R}^n x , b ∈ R n
直接法是指，通过有限次确定的算术运算得到舍入误差以内精确的数值解。
常见矩阵结构
对角矩阵：A = diag ( a 11 , … , a n n ) \mathbf A=\operatorname{diag}(a_{11},\dots,a_{nn}) A = diag ( a 11 , … , a nn )
三对角矩阵（tridiagonal）：仅在主对角线及其上下各一条对角线非零。
上/下三角矩阵：上三角 U \mathbf U U L \mathbf L L
单位三角矩阵：对角线元素全为1。
对称/正定矩阵：A ⊤ = A \mathbf A^\top=\mathbf A A ⊤ = A x ⊤ A x > 0 \mathbf x^\top \mathbf A \mathbf x>0 x ⊤ Ax > 0 ≥ 0 \ge0 ≥ 0
正交矩阵：Q ⊤ = Q − 1 \mathbf Q^\top = \mathbf Q^{-1} Q ⊤ = Q − 1
注 ：这些结构（稠密性）会影响解法选择与复杂度。
GAUSS 消去法
通过初等行变换将 A \mathbf A A U \mathbf U U x \mathbf x x
GAUSS 消去法思想
按列逐步消元。
在第 k k k k k k k k k k k k
对 k = 1 , 2 , … , n − 1 k=1,2,\dots,n-1 k = 1 , 2 , … , n − 1 U \mathbf U U b ~ \tilde {\mathbf b} b ~
数学基础：SWM 公式
Sherman–Morrison–Woodbury 公式：
( A + U V ⊤ ) − 1 = A − 1 − A − 1 U ( I + V ⊤ A − 1 U ) − 1 V ⊤ A − 1 . (\mathbf A+\mathbf U\mathbf V^\top)^{-1}
=
\mathbf A^{-1} - \mathbf A^{-1} \mathbf U
(\mathbf I + \mathbf V^\top \mathbf A^{-1} \mathbf U)^{-1} \mathbf V^\top \mathbf A^{-1}. ( A + U V ⊤ ) − 1 = A − 1 − A − 1 U ( I + V ⊤ A − 1 U ) − 1 V ⊤ A − 1 . 其中 A ∈ R n × n \mathbf A \in \mathbb R^{n \times n} A ∈ R n × n U , V ∈ R n × k \mathbf U, \; \mathbf V \in \mathbb R^{n \times k} U , V ∈ R n × k
证明 ：
构造块矩阵 M \mathbf M M
( A U V ⊤ − I ) \begin{pmatrix}
\mathbf A & \mathbf U \\ \mathbf V^\top & -\mathbf I
\end{pmatrix} ( A V ⊤ U − I ) 将其对角化，可得恒等式：
( A U V ⊤ − I ) = ( I 0 V ⊤ A − 1 I ) ( A 0 0 − I − V ⊤ A − 1 U ) ( I A − 1 U 0 I ) \begin{pmatrix}
\mathbf A & \mathbf U \\ \mathbf V^\top & -\mathbf I
\end{pmatrix}
=
\begin{pmatrix}
\mathbf I & \mathbf 0 \\ \mathbf V^\top \mathbf A^{-1} & \mathbf I
\end{pmatrix}
\begin{pmatrix}
\mathbf A & \mathbf 0 \\ \mathbf 0 & -\mathbf I - \mathbf V^\top \mathbf A^{-1} \mathbf U
\end{pmatrix}
\begin{pmatrix}
\mathbf I & \mathbf A^{-1}\mathbf U \\ \mathbf 0 & \mathbf I
\end{pmatrix} ( A V ⊤ U − I ) = ( I V ⊤ A − 1 0 I ) ( A 0 0 − I − V ⊤ A − 1 U ) ( I 0 A − 1 U I ) 对 M \mathbf M M
M 1 , 1 − 1 = ( A + U V ⊤ ) − 1 \mathbf M_{1, 1}^{-1} = (\mathbf A + \mathbf U\mathbf V^\top)^{-1} M 1 , 1 − 1 = ( A + U V ⊤ ) − 1 对恒等式右侧分别求逆并相乘得：
M − 1 = ( I − A − 1 U 0 I ) ( A − 1 0 0 − ( I + V ⊤ A − 1 U ) − 1 ) ( I 0 − V ⊤ A − 1 I ) \mathbf M^{-1}
=
\begin{pmatrix}
\mathbf I & -\mathbf A^{-1}\mathbf U \\ \mathbf 0 & \mathbf I
\end{pmatrix}
\begin{pmatrix}
\mathbf A^{-1} & \mathbf 0 \\ \mathbf 0 & -(\mathbf I + \mathbf V^\top \mathbf A^{-1} \mathbf U)^{-1}
\end{pmatrix}
\begin{pmatrix}
\mathbf I & \mathbf 0 \\ -\mathbf V^\top \mathbf A^{-1} & \mathbf I
\end{pmatrix} M − 1 = ( I 0 − A − 1 U I ) ( A − 1 0 0 − ( I + V ⊤ A − 1 U ) − 1 ) ( I − V ⊤ A − 1 0 I ) 故
M 1 , 1 − 1 = A − 1 − A − 1 U ( I + V ⊤ A − 1 U ) − 1 V ⊤ A − 1 \mathbf M_{1, 1}^{-1}
= \mathbf A^{-1} - \mathbf A^{-1} \mathbf U
(\mathbf I + \mathbf V^\top \mathbf A^{-1} \mathbf U)^{-1} \mathbf V^\top \mathbf A^{-1} M 1 , 1 − 1 = A − 1 − A − 1 U ( I + V ⊤ A − 1 U ) − 1 V ⊤ A − 1 矩阵视角下的 GAUSS 消去法
消元可以用矩阵乘法描述。定义第 k k k
M k = I − m k e k ⊤ , \mathbf M_k = \mathbf I - \mathbf m_k \mathbf e_k^\top, M k = I − m k e k ⊤ , 其中向量
m k = [ 0 , … , 0 , m k + 1 , k , … , m n , k ] ⊤ , m j , k = a j , k ( k ) a k , k ( k ) , j = k + 1 , … , n , \mathbf m_k = [0,\dots,0,\,m_{k+1,k},\dots,m_{n,k}]^\top,\quad m_{j,k}=\frac{a^{(k)}_{j,k}}{a^{(k)}_{k,k}},\; j=k+1,\dots,n, m k = [ 0 , … , 0 , m k + 1 , k , … , m n , k ] ⊤ , m j , k = a k , k ( k ) a j , k ( k ) , j = k + 1 , … , n , 且 e k \mathbf e_k e k k k k
第 k k k M k \mathbf M_k M k A ( k + 1 ) = M k A ( k ) \mathbf A^{(k+1)} = \mathbf M_k \mathbf A^{(k)} A ( k + 1 ) = M k A ( k ) n − 1 n-1 n − 1
A ( n ) = M n − 1 ⋯ M 2 M 1 A ( 1 ) ≡ U \mathbf A^{(n)} = \mathbf M_{n-1} \cdots \mathbf M_2 \mathbf M_1 \mathbf A^{(1)} \equiv \mathbf U A ( n ) = M n − 1 ⋯ M 2 M 1 A ( 1 ) ≡ U 且每个 M k \mathbf M_k M k
由 k = 1 k = 1 k = 1
( A + u v ⊤ ) − 1 = A − 1 − A − 1 u v ⊤ A − 1 1 + v ⊤ A − 1 u (\mathbf A+\mathbf u\mathbf v^\top)^{-1}
= \mathbf A^{-1}-\frac{\mathbf A^{-1} \mathbf u \mathbf v^\top \mathbf A^{-1}}{1+\mathbf v^\top \mathbf A^{-1} \mathbf u} ( A + u v ⊤ ) − 1 = A − 1 − 1 + v ⊤ A − 1 u A − 1 u v ⊤ A − 1 取 A = I \mathbf A=\mathbf I A = I u = − m k \mathbf u=-\mathbf m_k u = − m k v = e k \mathbf v=\mathbf e_k v = e k
M k − 1 = I + m k e k ⊤ \mathbf M_k^{-1}=\mathbf I + \mathbf m_k \mathbf e_k^\top M k − 1 = I + m k e k ⊤ 令
L = M 1 − 1 M 2 − 1 ⋯ M n − 1 − 1 \mathbf L = \mathbf M_1^{-1} \mathbf M_2^{-1}\cdots \mathbf M_{n-1}^{-1} L = M 1 − 1 M 2 − 1 ⋯ M n − 1 − 1 则得到标准的 LU 分解
A = L U \mathbf A = \mathbf L \mathbf U A = LU 其中 L \mathbf L L U \mathbf U U
注 ：
有趣的时，此处 L \mathbf L L ( j , k ) (j,k) ( j , k ) k k k a j , k a_{j,k} a j , k m j , k m_{j,k} m j , k L \mathbf L L
L = ( 1 m 21 1 m 31 m 32 1 ⋮ ⋮ ⋱ ⋱ m n 1 m n 2 ⋯ m n , n − 1 1 ) \mathbf L = \begin{pmatrix}
1 & & & & \\
m_{21} & 1 & & & \\
m_{31} & m_{32} & 1 & & \\
\vdots & \vdots & \ddots & \ddots & \\
m_{n1} & m_{n2} & \cdots & m_{n,n-1} & 1
\end{pmatrix} L = 1 m 21 m 31 ⋮ m n 1 1 m 32 ⋮ m n 2 1 ⋱ ⋯ ⋱ m n , n − 1 1 GAUSS 消去法代码
假设 A ∈ R n × n \mathbf A \in \mathbb{R}^{n \times n} A ∈ R n × n A[i, j]，b ∈ R n \mathbf b \in \mathbb{R}^{n} b ∈ R n b[i]
# 消去算法 for k in range (n - 1 ): # 枚举列 for i in range (k + 1 , n): # 枚举行 c = - A[i, k] / A[k, k] # 乘子 for j in range (k + 1 , n + 1 ): # 消元 A[i, j] = A[i, j] + c * A[k, j] b[i] = b[i] + c * b[k] # 回代算法 for i in range (n - 1 , - 1 , - 1 ): # 逆序回代 for j in range (i + 1 , n - 1 ): # 此时这些 b[j] 已经是准确的了 b[i] = b[i] - A[i, j] * b[j] b[i] = b[i] / A[i, i] # 因为没有规约化 注 ：
算法中没有必要对 A[k, k] 做归一化（节约运算），同时只需更新矩阵的上三角部分（下三角不再需要）。计算中没有考虑除 0 问题。
GAUSS 消去法复杂度分析
除法数：约 n ( n − 1 ) 2 \dfrac{n(n-1)}{2} 2 n ( n − 1 )
乘法数：约 n ( n − 1 ) ( n + 1 ) 3 \dfrac{n(n-1)(n+1)}{3} 3 n ( n − 1 ) ( n + 1 )
加法数：同样约 1 3 n 3 \dfrac{1}{3} n^3 3 1 n 3
时间复杂度：O ( n 3 ) O(n^3) O ( n 3 )
空间复杂度：原地修改 A \mathbf A A O ( n 2 ) O(n^2) O ( n 2 )
列主元消去法
列主元消去法思想
考虑下面问题：
[ ε 1 1 1 ] [ x 1 x 2 ] = [ 1 ε ] , \begin{bmatrix}\varepsilon & 1\\ 1 & 1\end{bmatrix}\begin{bmatrix}x_1\\x_2\end{bmatrix}=\begin{bmatrix}1\\\varepsilon\end{bmatrix}, [ ε 1 1 1 ] [ x 1 x 2 ] = [ 1 ε ] , 真实解：
x = [ − 1 1 + ε ] \mathbf x = \begin{bmatrix}-1\\ 1+\varepsilon\end{bmatrix} x = [ − 1 1 + ε ]
直接 Gauss 消去得:
[ ε 1 1 0 1 − 1 / ε ε − 1 / ε ] \begin{bmatrix} \varepsilon & 1 & 1 \\ 0 & 1-1/\varepsilon & \varepsilon-1/\varepsilon \end{bmatrix} [ ε 0 1 1 − 1/ ε 1 ε − 1/ ε ] 若 ε \varepsilon ε x = [ 0 1 ] \mathbf x = \begin{bmatrix}0\\ 1\end{bmatrix} x = [ 0 1 ]
也即当主元 a k , k ( k ) a^{(k)}_{k,k} a k , k ( k )
解决方案是在每一步选择该列中绝对值最大的元素作主元并交换行，称为部分选主元（partial pivoting）。由于 A \mathbf A A
注 ：还有完全选主元法（full pivoting），也即行列同时选，代价更高但更稳定。工业中部分选主元（partial pivoting）通常是均衡的选择。
数学基础：初等交换阵
交换第 i i i j j j j j j i i i
P i j = I − e i e i ⊤ − e j e j ⊤ + e i e j ⊤ + e j e i ⊤ \mathbf P_{ij}
= \mathbf I - \mathbf e_i \mathbf e_i^\top - \mathbf e_j \mathbf e_j^\top
+ \mathbf e_i \mathbf e_j^\top + \mathbf e_j \mathbf e_i^\top P ij = I − e i e i ⊤ − e j e j ⊤ + e i e j ⊤ + e j e i ⊤ 其中 e k \mathbf e_k e k k k k
性质 ：
P i j ⊤ = P i j \mathbf P_{ij}^\top = \mathbf P_{ij} P ij ⊤ = P ij P i j − 1 = P i j \mathbf P_{ij}^{-1} = \mathbf P_{ij} P ij − 1 = P ij
矩阵视角下的列主元消去法
这一部分对于算法设计用处不大，但还是看一眼喵~在列主元消去法中，每一步消元前需要先进行行交换。设第 k k k p p p p ≥ k p \geq k p ≥ k
A ( n ) = M n − 1 P n − 1 M n − 2 P n − 2 ⋯ M 2 P 2 M 1 P 1 A ( 1 ) \mathbf{\mathbf A}^{(n)}
= \mathbf{\mathbf M}_{n-1}\mathbf{P}_{n-1}\mathbf{\mathbf M}_{n-2}\mathbf{P}_{n-2}
\cdots\mathbf{\mathbf M}_2\mathbf{P}_2\mathbf{\mathbf M}_1\mathbf{P}_1\mathbf{\mathbf A}^{(1)} A ( n ) = M n − 1 P n − 1 M n − 2 P n − 2 ⋯ M 2 P 2 M 1 P 1 A ( 1 ) 其中 P k \mathbf{P}_k P k k k k p p p M k \mathbf{\mathbf M}_k M k A ( n ) ≡ U \mathbf{\mathbf A}^{(n)} \equiv \mathbf{\mathbf U} A ( n ) ≡ U
对于 M k = I − m k e k ⊤ \mathbf{\mathbf M}_k = \mathbf{\mathbf I} - \mathbf{m}_k\mathbf{e}_k^\top M k = I − m k e k ⊤
P k + 1 M k P k + 1 ≡ M ‾ k \mathbf{P}_{k+1}\mathbf{\mathbf M}_k\mathbf{P}_{k+1} \equiv \overline{\mathbf{\mathbf M}}_k P k + 1 M k P k + 1 ≡ M k 其中由于原消元矩阵的稀疏性，故 M ‾ k \overline{\mathbf{\mathbf M}}_k M k M k \mathbf{\mathbf M}_k M k ( k + 1 , k ) (k+1, k) ( k + 1 , k ) ( p , k ) (p, k) ( p , k )
因为
A ( n ) = M n − 1 P n − 1 ⋯ M 1 P 1 A \mathbf{\mathbf A}^{(n)}
= \mathbf{\mathbf M}_{n-1}\mathbf{P}_{n-1}\cdots\mathbf{\mathbf M}_1\mathbf{P}_1\mathbf{\mathbf A} A ( n ) = M n − 1 P n − 1 ⋯ M 1 P 1 A 可得：
M ‾ 1 − 1 M ‾ 2 − 1 ⋯ M ‾ n − 2 − 1 M n − 1 − 1 U = P n − 1 P n − 2 ⋯ P 2 P 1 ⏟ ≡ P A \overline{\mathbf{\mathbf M}}_1^{-1}\overline{\mathbf{\mathbf M}}_2^{-1}
\cdots\overline{\mathbf{\mathbf M}}_{n-2}^{-1}\mathbf{\mathbf M}_{n-1}^{-1}\mathbf{\mathbf U}
= \underbrace{\mathbf{P}_{n-1}\mathbf{P}_{n-2}
\cdots\mathbf{P}_2\mathbf{P}_1}_{\equiv \mathbf{P}}\mathbf{\mathbf A} M 1 − 1 M 2 − 1 ⋯ M n − 2 − 1 M n − 1 − 1 U = ≡ P P n − 1 P n − 2 ⋯ P 2 P 1 A 即 带列主元的 LU 分解 ：
P A = L U \mathbf{P\mathbf A} = \mathbf{LU} PA = LU 其中：
P \mathbf{P} P 排列矩阵 （permutation matrix），满足 P − 1 = P ⊤ \mathbf{P}^{-1} = \mathbf{P}^\top P − 1 = P ⊤ L = \mathbf{\mathbf L} = L = M ‾ 1 − 1 ⋯ M ‾ n − 2 − 1 M n − 1 − 1 \overline{\mathbf{\mathbf M}}_1^{-1}\cdots\overline{\mathbf{\mathbf M}}_{n-2}^{-1}\mathbf{\mathbf M}_{n-1}^{-1} M 1 − 1 ⋯ M n − 2 − 1 M n − 1 − 1 U \mathbf{\mathbf U} U
列主元消去法代码
假设 A ∈ R n × n \mathbf A \in \mathbb{R}^{n \times n} A ∈ R n × n A[i, j]，b ∈ R n \mathbf b \in \mathbb{R}^{n} b ∈ R n b[i]
# 消去算法 # 排列向量 p[i] = j 表示当前的第 i 行对应原始的第 j 行 for i in range (n): p[i] = i for k in range (n - 1 ): # 枚举列 s = k for i in range (k + 1 , n): if abs (A[p[i], k]) > abs (A[p[s], k]): s = i p[s], p[k] = p[k], p[s] for i in range (k + 1 , n): # 枚举行 c = - A[p[i], k] / A[p[k], k] # 乘子 for j in range (k + 1 , n + 1 ): # 消元 A[p[i], j] = A[p[i], j] + c * A[p[k], j] b[p[i]] = b[p[i]] + c * b[p[k]] # 回代算法 for i in range (n - 1 , - 1 , - 1 ): # 逆序回代 for j in range (i + 1 , n - 1 ): # 此时这些 b[p[j]] 已经是准确的了 b[p[i]] = b[p[i]] - A[p[i], j] * b[p[j]] b[p[i]] = b[p[i]] / A[p[i], i] # 因为没有规约化 注 ：操作上使用排列向量 p[i] 来避免昂贵的内存行交换。
列主元消去法复杂度分析
时间复杂度、空间复杂度与 Gauss 消元法量级相同，但是在时间复杂度上会多 O ( n 2 ) O(n^2) O ( n 2 ) O ( n ) O(n) O ( n )