[学习] 矩阵求逆常用算法（含实现代码）

矩阵求逆常用算法

---

引言

矩阵求逆是求解线性方程组 $Ax=b$ 的核心操作，在优化问题（如牛顿法）、统计建模（协方差矩阵求逆）及控制系统设计中广泛用。本文系统梳理六大类算法，从数学原理到代码实现，并给出场景选择指南。

---

1 矩阵求逆的基本概念

数学条件：矩阵 $A$ 可逆当且仅当：

• $\text{rank}(A)=n$（满秩）
• $\det (A)\ne 0$
• $A$ 的行/列向量线性无关

核心性质：

1. $(A^{-1}{)}^{-1}=A$
2. $(AB{)}^{-1}=B^{-1}{A}^{-1}$
3. $(A^T{)}^{-1}=(A^{-1}{)}^T$

实际意义：若 $A$ 可逆，则方程组 $Ax=b$ 有唯一解 $x=A^{-1}b$。

import numpy as np

# 基本测试案例
A = np.array([[2, 1], [1, 3]])
A_inv_gt = np.linalg.inv(A)  # 标准解用于验证
print("标准逆矩阵:\n", A_inv_gt)

---

2 高斯-约旦消元法（Gauss-Jordan Elimination）

原理：通过初等行变换将 $[A|I]$ 转化为 $[I|A^{-1}]$。

推导步骤：

1. 构造增广矩阵 $M=[A\mid I_n]$
2. 前向消元：对 $k=0$ 到 $n-1$：
   • 选主元 $m_{kk}$（避免除零，增强稳定性）
   • 对 $i=k+1$ 到 $n$：
     $${\text{row}}_i\leftarrow {\text{row}}_i-\frac{m_{ik}}{m_{kk}}\times {\text{row}}_k$$
3. 后向消元：对 $k=n-1$ 到 $0$：
   • ${\text{row}}_k\leftarrow \frac{{\text{row}}_k}{m_{kk}}$
   • 对 $i=0$ 到 $k-1$：
     $${\text{row}}_i\leftarrow {\text{row}}_i-m_{ik}\times {\text{row}}_k$$
4. 提取 $M$ 右半部分得 $A^{-1}$

时间复杂度：$O(n^3)$
适用场景：$n<1000$ 的稠密矩阵。

def gauss_jordan_inv(A):
    n = A.shape
    M = np.hstack((A, np.eye(n)))
    for k in range(n):
        # 列主元选择
        pivot = np.argmax(np.abs(M[k:, k])) + k
        M[[k, pivot]] = M[[pivot, k]]
        
        # 归一化主元行
        M[k] = M[k] / M[k, k]
        
        # 消元
        for i in range(n):
            if i != k:
                M[i] -= M[i, k] * M[k]
    return M[:, n:]

# 测试
A_inv_gj = gauss_jordan_inv(A)
print("高斯-约旦法结果:\n", A_inv_gj)
print("误差:", np.linalg.norm(A_inv_gj - A_inv_gt))

---

3 LU分解法

原理：分解 $A=LU$（$L$ 下三角，$U$ 上三角），则 $A^{-1}=U^{-1}{L}^{-1}$。

分解过程（Doolittle算法）：
$$u_{kj}=a_{kj}-\sum _{m=1}^{k-1}{l}_{km}{u}_{mj},l_{ik}=\frac{1}{u_{kk}}\left(a_{ik}-\sum _{m=1}^{k-1}{l}_{im}{u}_{mk}\right)$$

求逆步骤：

1. 解 $LY=I$ 得 $Y=L^{-1}$（前向替换）
2. 解 $UX=Y$ 得 $X=U^{-1}$（后向替换）
3. $A^{-1}=U^{-1}{L}^{-1}$

时间复杂度：分解 $O(n^3)$，求逆 $O(n^2)$（若需多次求逆，分解只需一次）
适用场景：需重复求解 $Ax=b_k$ 的问题。

def lu_inv(A):
    n = A.shape
    L = np.eye(n)
    U = np.zeros((n, n))
    # LU分解
    for k in range(n):
        U[k, k:] = A[k, k:] - L[k, :k] @ U[:k, k:]
        for i in range(k+1, n):
            L[i, k] = (A[i, k] - L[i, :k] @ U[:k, k]) / U[k, k]
    
    # 求L^{-1}和U^{-1}
    L_inv = np.linalg.inv(L)  # 实际用前向替换更高效
    U_inv = np.linalg.inv(U)  # 实际用后向替换
    return U_inv @ L_inv

# 测试
A_inv_lu = lu_inv(A)
print("LU分解法误差:", np.linalg.norm(A_inv_lu - A_inv_gt))

---

4 伴随矩阵法（Adjugate Method）

原理：$A^{-1}=\frac{1}{\det (A)}\text{adj}(A)$，其中 $\text{adj}(A)$ 是代数余子式矩阵的转置。

代数余子式：
$$\text{adj}(A{)}_{ij}=(-1{)}^{i+j}{M}_{ji}$$
$M_{ji}$ 是删除第 $j$ 行第 $i$ 列的子矩阵行列式。

时间复杂度：计算行列式需 $O(n!)$，仅适用于 $n\le 3$。

def adjugate_inv(A):
    n = A.shape
    if n == 1:
        return np.array([[1/A[0,0]]])
    det_A = np.linalg.det(A)
    adj = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            # 生成余子式矩阵
            minor = np.delete(np.delete(A, i, axis=0), j, axis=1)
            adj[i, j] = ((-1)**(i+j)) * np.linalg.det(minor)
    return adj.T / det_A

# 测试（仅用于小矩阵）
A_small = np.array([[1, 2], [3, 4]])
A_inv_adj = adjugate_inv(A_small)
print("伴随矩阵法结果:\n", A_inv_adj)

---

5 迭代法：牛顿迭代（Newton-Schulz）

原理：构造迭代序列 $X_{k+1}=2X_k-X_{k}A{X}_k$，收敛到 $A^{-1}$。收敛性证明：
定义误差矩阵 $E_k=I-A{X}_k$，则：
$$E_{k+1}=I-A(2X_k-X_{k}A{X}_k)=(I-A{X}_k{)}^2=E_k^2$$
若 $\Vert E_0\Vert <1$，则 $\Vert E_k\Vert \to 0$（二次收敛）。

适用场景：大规模稀疏矩阵（$n>10^4$），可并行化。

def newton_inv(A, max_iter=50, tol=1e-6):
    X = A.T / (np.linalg.norm(A, 1) * np.linalg.norm(A, np.inf))  # 初始估计
    I = np.eye(A.shape)
    for _ in range(max_iter):
        X_new = 2*X - X @ A @ X
        if np.linalg.norm(I - A @ X_new) < tol:
            break
        X = X_new
    return X_new

# 测试（对称正定矩阵收敛性好）
A_large = np.random.rand(100, 100) 
A_large = A_large @ A_large.T + 100*np.eye(100)  # 确保可逆
A_inv_newton = newton_inv(A_large)
print("牛顿迭代误差:", np.linalg.norm(A_large @ A_inv_newton - np.eye(100)))

---

6 分块矩阵求逆法

原理：对分块矩阵 $A=\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]$，其逆为：
$$A^{-1}=\left[\begin{array}{cc}{S}^{-1}& -S^{-1}{A}_{12}{A}_{22}^{-1}\\ -A_{22}^{-1}{A}_{21}{S}^{-1}& A_{22}^{-1}+A_{22}^{-1}{A}_{21}{S}^{-1}{A}_{12}{A}_{22}^{-1}\end{array}\right]$$
其中 $S=A_{11}-A_{12}{A}_{22}^{-1}{A}_{21}$（Schur补）。

适用场景：分块对角占优矩阵、GPU并行计算。

def block_inv(A, block_size):
    n = A.shape
    if n <= block_size:
        return np.linalg.inv(A)
    # 分块
    A11 = A[:block_size, :block_size]
    A12 = A[:block_size, block_size:]
    A21 = A[block_size:, :block_size]
    A22 = A[block_size:, block_size:]
    
    # 递归求逆
    A22_inv = block_inv(A22, block_size)
    S = A11 - A12 @ A22_inv @ A21
    S_inv = block_inv(S, block_size)
    
    # 组合结果
    B11 = S_inv
    B12 = -S_inv @ A12 @ A22_inv
    B21 = -A22_inv @ A21 @ S_inv
    B22 = A22_inv + A22_inv @ A21 @ S_inv @ A12 @ A22_inv
    return np.block([[B11, B12], [B21, B22]])

# 测试（分块大小需为2的幂）
A_block = np.random.rand(8, 8)
A_inv_block = block_inv(A_block, 4)
print("分块法误差:", np.linalg.norm(A_inv_block - np.linalg.inv(A_block)))

---

7 特殊矩阵的优化算法

矩阵类型	求逆公式	时间复杂度
对角矩阵	$(D^{-1}{)}_{ii}=1/d_{ii}$	$O(n)$
正交矩阵	$Q^{-1}=Q^T$	$O(1)$
三对角矩阵	追赶法（Thomas算法）	$O(n)$

# 三对角矩阵求逆示例（追赶法）
def tridiag_inv(diag, sub, super):
    n = len(diag)
    inv = np.zeros((n, n))
    # 此处实现追赶法求逆（略）
    return inv

---

8 数值稳定性与优化

关键问题：

• 病态矩阵：条件数 $\kappa (A)=\Vert A\Vert \Vert A^{-1}\Vert$ 过大时，微小扰动导致解失真
• 稳定性措施：
  1. 高斯法中使用列主元消去（Partial Pivoting）
  2. 添加正则化项：$ (A + \lambda I)^{-1} $（Tikhonov正则化）
  3. 用SVD分解替代：$A^{-1}=V{\Sigma }^{-1}{U}^T$（稳定但昂贵）

硬件加速：

• GPU并行：CuBLAS库的cublasDgetrfBatched加速LU分解
• 分布式计算：Spark MLlib的BlockMatrix分块求逆

---

9 总结与对比

算法	时间复杂度	空间复杂度	适用场景
高斯-约旦	$O(n^3)$	$O(n^2)$	中小稠密矩阵
LU分解	$O(n^3)$	$O(n^2)$	需多次求逆的方程组
伴随矩阵	$O(n!)$	$O(n^2)$	$n\le 3$ 的小矩阵
牛顿迭代	$O(k{n}^2)$	$O(n^2)$	大规模稀疏矩阵
分块法	$O(n^{2.8})$	$O(n^2)$	分块并行计算
特殊矩阵优化	$O(n)$	$O(n)$	对角/三对角/正交矩阵

选择建议：

• $n<1000$：优先选高斯-约旦法（实现简单）
• 多次求解：用LU分解（分解结果复用）
• 大规模稀疏：牛顿迭代法（内存友好）
• 病态矩阵：SVD分解（牺牲速度换稳定性）

未来方向：量子算法（HHL算法）可在 $O(\log n)$ 时间求逆，但需量子硬件支持。

---

研究学习不易，点赞易。
工作生活不易，收藏易，点收藏不迷茫 ：）

---