1. 引言
几乎所有的物理系统或科学问题都涉及一个线性代数方程组。例如，拟合数据点曲线、优化成本函数、分析电路网络等问题，都需要求解线性方程。
计算线性代数中一些基本任务包括：
- 求解线性方程组 $ A\mathbf{x} = \mathbf{b} $
- 求一个方阵 $ A $ 的逆
- 计算矩阵的行列式
接下来将介绍一些用于执行这些基本任务的实用算法。
我们将学习如何将一个矩阵分解成更简单的“积木块”。对这些较小的“积木块”进行操作通常比直接处理原始矩阵更容易。
2. 部分选主元的高斯消去法
考虑一个包含三个未知数的线性方程组：
$$ \begin{align*} x + 2y + z &= 2\ 2x + 6y + z &= 7\ x + y + 4z &= 3 \end{align*} $$
我们可以将其写成矩阵形式 $ A\mathbf{x} = \mathbf{b} $：
$$ \begin{bmatrix} 1 & 2 & 1 \ 2 & 6 & 1 \ 1 & 1 & 4 \end{bmatrix} \begin{bmatrix} x \ y \ z \end{bmatrix} = \begin{bmatrix} 2 \ 7 \ 3 \end{bmatrix} $$
高中阶段学过的解法是通过初等行变换对方程进行处理。由于行变换会影响方程右边的值，因此我们通常使用增广矩阵 $ M=[A|\mathbf{b}] $ 来跟踪所有变化。
我们将 $ A $ 中第 $ (1,1) $ 位置的元素称为第一个主元（pivot），第 $ (2,2) $ 位置的元素称为第二个主元，依此类推。在高斯消去法中，核心思想是：在第 $ j $ 列中，将主元 $ a_{jj} $ 下方的所有元素变为 0。
因此，一个关键要求是：主元必须非零。 有时我们需要交换第 $ j $ 行和第 $ k $ 行（其中 $ k > j $），以确保主元非零。最常见的做法是选择 $ k $，使得 $ |a_{kj}| $ 最大化。这被称为部分选主元（Partial Pivoting）。
高斯消去法的伪代码如下：
algorithm GaussianElim(A, b):
// INPUT
// A = 非奇异矩阵，大小 n x n
// b = 右侧向量
// OUTPUT
// U = 上三角形式的 A
// c = 对应 U 的变换后的 b 向量
if A[j][j] == 0:
big = 0.0
kRow = j
for k from j+1 to n-1:
if abs(A[k][j]) > big:
big = abs(A[k][j])
kRow = k
// 交换行
for l from j to n-1:
swap A[j][l] and A[kRow][l]
swap b[j] and b[kRow]
pivot = A[j][j]
if pivot == 0:
throw exception (矩阵奇异)
for i from j+1 to n-1:
mult = A[i][j] / pivot
for l from j to n-1:
A[i][l] -= mult * A[j][l]
b[i] -= mult * b[j]
return A, b
消去完成后，我们得到一个上三角系统 $ U\mathbf{x} = \mathbf{c} $，可以使用回代法（Back Substitution）求解：
algorithm BackSubstitution(U, c):
// INPUT
// U = 上三角矩阵
// c = 右侧向量
// OUTPUT
// 解向量 x
for i from n-1 downto 0:
sum = 0
for j from i+1 to n-1:
sum += x[j] * U[i][j]
x[i] = (c[i] - sum) / U[i][i]
return x
例如，求解上述电路问题得到的解向量为 $ \mathbf{x} = [-3, 2, 1] $。
3. LU 分解
假设我们可以将矩阵 $ A $ 分解为两个矩阵的乘积：
$$ A = L \cdot U $$
其中 $ L $ 是下三角矩阵，$ U $ 是上三角矩阵。我们可以用这个分解来解线性方程组：
$$ A\mathbf{x} = (L \cdot U)\mathbf{x} = L \cdot (U\mathbf{x}) = \mathbf{b} $$
我们首先求解 $ L\mathbf{y} = \mathbf{b} $，然后求解 $ U\mathbf{x} = \mathbf{y} $。这两个三角系统的求解分别使用前向代入（Forward Substitution）和回代法（Back Substitution）即可。
前向代入算法如下：
algorithm FwdSubstitution(L, b):
for i from 0 to n-1:
sum = 0
for j from 0 to i-1:
sum += y[j] * L[i][j]
y[i] = (b[i] - sum) / L[i][i]
return y
3.1. LU 分解算法
设 $ A $ 是一个 $ n \times n $ 的非奇异方阵，我们可以将其写成块形式：
$$ \begin{bmatrix} a_{11} & \mathbf{a_{12}} \ \mathbf{a_{21}} & A_{22} \end{bmatrix} = \begin{bmatrix} 1 & \mathbf{0} \ \mathbf{l_{21}} & L_{22} \end{bmatrix} \begin{bmatrix} u_{11} & \mathbf{u_{12}} \ \mathbf{0} & U_{22} \end{bmatrix} $$
通过比较两边，可以得到：
$$ \begin{align*} u_{11} &= a_{11} \ \mathbf{u_{12}} &= \mathbf{a_{12}} \ \mathbf{l_{21}} &= \frac{1}{u_{11}} \mathbf{a_{21}} \ L_{22}U_{22} &= A_{22} - \mathbf{l_{21}} \cdot \mathbf{u_{12}} \end{align*} $$
我们可以递归地解决这个子问题，从而得到完整的 LU 分解。
伪代码如下：
algorithm LUdcmp(A):
// INPUT
// A = 非奇异方阵
// OUTPUT
// L = 单位下三角矩阵
// U = 上三角矩阵
N = A 的维度
if N == 1:
return L = [1], U = A
a11 = A[0][0]
a12 = A[0][1..N-1]
a21 = A[1..N-1][0]
A22 = A[1..N-1][1..N-1]
u11 = a11
u12 = a12
l21 = a21 / u11
S22 = A22 - l21 * u12
L22, U22 = LUdcmp(S22)
构造 L 和 U 矩阵并返回
4. QR 分解
除了 LU 分解，还有一种常用的矩阵分解方法是 QR 分解。它在某些场景下特别有用。
4.1. Gram-Schmidt 正交化过程
假设我们有一组基向量 $ \mathcal{B} = { \mathbf{w}_1, \ldots, \mathbf{w}_n } $，我们希望从中构造一个正交基 $ { \mathbf{v}_1, \ldots, \mathbf{v}_n } $。
Gram-Schmidt 过程如下：
$$ \mathbf{v}_k = \mathbf{w}k - \sum{j=1}^{k-1} \frac{\mathbf{w}_k \cdot \mathbf{v}_j}{|\mathbf{v}_j|^2} \mathbf{v}_j $$
通过这个过程，我们可以逐步构造出正交基。
4.2. Gram-Schmidt 的改进
我们可以通过正交化后归一化的方式直接构造标准正交基 $ { \mathbf{u}_1, \ldots, \mathbf{u}_n } $，并表示原基向量为：
$$ \mathbf{w}j = \sum{i=1}^j r_{ij} \mathbf{u}_i $$
其中：
$$ r_{ij} = \mathbf{w}j \cdot \mathbf{u}i, \quad i < j $$ $$ r{jj} = | \mathbf{w}j | - \sum{i=1}^{j-1} r{ij}^2 $$
4.3. 正交矩阵
如果一个矩阵的列向量构成了 $ \mathbf{R}^n $ 的标准正交基，则称该矩阵为正交矩阵。正交矩阵 $ Q $ 满足：
$$ Q^T Q = Q Q^T = I $$
因此，其逆矩阵为：
$$ Q^{-1} = Q^T $$
4.4. QR 分解算法
QR 分解可以将矩阵 $ A $ 分解为一个正交矩阵 $ Q $ 和一个上三角矩阵 $ R $：
$$ A = Q R $$
实现 QR 分解的伪代码如下：
algorithm QRFact(A):
// INPUT
// A = 非奇异方阵
// OUTPUT
// Q = 正交矩阵
// R = 上三角矩阵
N = A 的维度
u = 零矩阵
R = 零矩阵
for j from 0 to N-1:
w_j = A[:, j]
// 计算 w_j 的模平方
norm2 = sum(w_j[k]^2 for k in 0..N-1)
// 计算 R[i][j]
sum_r = 0
for i from 0 to j-1:
R[i][j] = dot(w_j, u[i])
sum_r += R[i][j]^2
R[j][j] = sqrt(norm2 - sum_r)
// 计算 u[j]
sum_vec = 0
for i from 0 to j-1:
sum_vec += R[i][j] * u[i]
u[j] = (w_j - sum_vec) / R[j][j]
Q = transpose(u)
return Q, R
4.5. 使用 QR 分解求解线性方程组
一旦我们有了 QR 分解，就可以高效地求解线性方程组：
$$ A\mathbf{x} = \mathbf{b} \Rightarrow QR\mathbf{x} = \mathbf{b} \Rightarrow R\mathbf{x} = Q^T\mathbf{b} $$
由于 $ R $ 是上三角矩阵，所以可以使用回代法快速求解。
5. 总结
本文回顾了高斯消去法这一中学阶段就接触的解线性方程组的方法，并介绍了 LU 和 QR 分解等更高效的矩阵分解技术。这些分解方法是更复杂算法的基础。
如今，像 MATLAB、Eigen、LAPACK、Linpack、NAG 等数值计算库都已广泛支持 $ A = LU $、$ A = QR $ 等分解方法，用于求解高达 $ 10^6 \times 10^6 $ 的大型线性系统。