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

矩阵求逆常用算法

文章目录

矩阵求逆常用算法引言1 矩阵求逆的基本概念2 高斯-约旦消元法（Gauss-Jordan Elimination）3 LU分解法4 伴随矩阵法（Adjugate Method）5 迭代法：牛顿迭代（Newton-Schulz）6 分块矩阵求逆法7 特殊矩阵的优化算法8 数值稳定性与优化9 总结与对比

引言

矩阵求逆是求解线性方程组

A

x

=

b

Ax = b

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

1 矩阵求逆的基本概念

数学条件：矩阵

A

A

A 可逆当且仅当：

rank

(

A

)

=

n

\text{rank}(A) = n

rank(A)=n（满秩）

det

⁡

(

A

)

≠

0

\det(A) \neq 0

det(A)=0

A

A

A 的行/列向量线性无关

核心性质：

(

A

−

1

)

−

1

=

A

(A^{-1})^{-1} = A

(A−1)−1=A

(

A

B

)

−

1

=

B

−

1

A

−

1

(AB)^{-1} = B^{-1}A^{-1}

(AB)−1=B−1A−1

(

A

T

)

−

1

=

(

A

−

1

)

T

(A^T)^{-1} = (A^{-1})^T

(AT)−1=(A−1)T

实际意义：若

A

A

A 可逆，则方程组

A

x

=

b

Ax = b

Ax=b 有唯一解

x

=

A

−

1

b

x = A^{-1}b

x=A−1b。

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

]

[A | I]

[A∣I] 转化为

[

I

∣

A

−

1

]

[I | A^{-1}]

[I∣A−1]。

推导步骤：

构造增广矩阵

M

=

[

A

∣

I

n

]

M = [A \mid I_n]

M=[A∣In​]前向消元：对

k

=

0

k=0

k=0 到

n

−

1

n-1

n−1：

选主元

m

k

k

m_{kk}

mkk​（避免除零，增强稳定性）对

i

=

k

+

1

i = k+1

i=k+1 到

n

n

n：

row

i

←

row

i

−

m

i

k

m

k

k

×

row

k

\text{row}_i \leftarrow \text{row}_i - \frac{m_{ik}}{m_{kk}} \times \text{row}_k

rowi​←rowi​−mkk​mik​​×rowk​ 后向消元：对

k

=

n

−

1

k=n-1

k=n−1 到

0

0

0：

row

k

←

row

k

m

k

k

\text{row}_k \leftarrow \frac{\text{row}_k}{m_{kk}}

rowk​←mkk​rowk​​对

i

=

0

i = 0

i=0 到

k

−

1

k-1

k−1：

row

i

←

row

i

−

m

i

k

×

row

k

\text{row}_i \leftarrow \text{row}_i - m_{ik} \times \text{row}_k

rowi​←rowi​−mik​×rowk​ 提取

M

M

M 右半部分得

A

−

1

A^{-1}

A−1

时间复杂度：

O

(

n

3

)

O(n^3)

O(n3) 适用场景：

n

<

1000

n < 1000

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

=

L

U

A = LU

A=LU（

L

L

L 下三角，

U

U

U 上三角），则

A

−

1

=

U

−

1

L

−

1

A^{-1} = U^{-1}L^{-1}

A−1=U−1L−1。

分解过程（Doolittle算法）：

u

k

j

=

a

k

j

−

∑

m

=

1

k

−

1

l

k

m

u

m

j

,

l

i

k

=

1

u

k

k

(

a

i

k

−

∑

m

=

1

k

−

1

l

i

m

u

m

k

)

u_{kj} = a_{kj} - \sum_{m=1}^{k-1} l_{km}u_{mj}, \quad l_{ik} = \frac{1}{u_{kk}} \left( a_{ik} - \sum_{m=1}^{k-1} l_{im}u_{mk} \right)

ukj​=akj​−m=1∑k−1​lkm​umj​,lik​=ukk​1​(aik​−m=1∑k−1​lim​umk​)

求逆步骤：

解

L

Y

=

I

LY = I

LY=I 得

Y

=

L

−

1

Y = L^{-1}

Y=L−1（前向替换）解

U

X

=

Y

UX = Y

UX=Y 得

X

=

U

−

1

X = U^{-1}

X=U−1（后向替换）

A

−

1

=

U

−

1

L

−

1

A^{-1} = U^{-1}L^{-1}

A−1=U−1L−1

时间复杂度：分解

O

(

n

3

)

O(n^3)

O(n3)，求逆

O

(

n

2

)

O(n^2)

O(n2)（若需多次求逆，分解只需一次） 适用场景：需重复求解

A

x

=

b

k

Ax=b_k

Ax=bk​ 的问题。

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

=

1

det

⁡

(

A

)

adj

(

A

)

A^{-1} = \frac{1}{\det(A)} \text{adj}(A)

A−1=det(A)1​adj(A)，其中

adj

(

A

)

\text{adj}(A)

adj(A) 是代数余子式矩阵的转置。

代数余子式：

adj

(

A

)

i

j

=

(

−

1

)

i

+

j

M

j

i

\text{adj}(A)_{ij} = (-1)^{i+j} M_{ji}

adj(A)ij​=(−1)i+jMji​

M

j

i

M_{ji}

Mji​ 是删除第

j

j

j 行第

i

i

i 列的子矩阵行列式。

时间复杂度：计算行列式需

O

(

n

!

)

O(n!)

O(n!)，仅适用于

n

≤

3

n \leq 3

n≤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

=

2

X

k

−

X

k

A

X

k

X_{k+1} = 2X_k - X_k A X_k

Xk+1​=2Xk​−Xk​AXk​，收敛到

A

−

1

A^{-1}

A−1。收敛性证明： 定义误差矩阵

E

k

=

I

−

A

X

k

E_k = I - A X_k

Ek​=I−AXk​，则：

E

k

+

1

=

I

−

A

(

2

X

k

−

X

k

A

X

k

)

=

(

I

−

A

X

k

)

2

=

E

k

2

E_{k+1} = I - A(2X_k - X_k A X_k) = (I - A X_k)^2 = E_k^2

Ek+1​=I−A(2Xk​−Xk​AXk​)=(I−AXk​)2=Ek2​ 若

∥

E

0

∥

<

1

\|E_0\| < 1

∥E0​∥<1，则

∥

E

k

∥

→

0

\|E_k\| \to 0

∥Ek​∥→0（二次收敛）。

适用场景：大规模稀疏矩阵（

n

>

10

4

n > 10^4

n>104），可并行化。

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

=

[

A

11

A

12

A

21

A

22

]

A = \begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{bmatrix}

A=[A11​A21​​A12​A22​​]，其逆为：

A

−

1

=

[

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

]

A^{-1} = \begin{bmatrix} 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{bmatrix}

A−1=[S−1−A22−1​A21​S−1​−S−1A12​A22−1​A22−1​+A22−1​A21​S−1A12​A22−1​​] 其中

S

=

A

11

−

A

12

A

22

−

1

A

21

S = A_{11} - A_{12}A_{22}^{-1}A_{21}

S=A11​−A12​A22−1​A21​（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

)

i

i

=

1

/

d

i

i

(D^{-1})_{ii} = 1/d_{ii}

(D−1)ii​=1/dii​

O

(

n

)

O(n)

O(n)正交矩阵

Q

−

1

=

Q

T

Q^{-1} = Q^T

Q−1=QT

O

(

1

)

O(1)

O(1)三对角矩阵追赶法（Thomas算法）

O

(

n

)

O(n)

O(n)

# 三对角矩阵求逆示例（追赶法）

def tridiag_inv(diag, sub, super):

n = len(diag)

inv = np.zeros((n, n))

# 此处实现追赶法求逆（略）

return inv

8 数值稳定性与优化

关键问题：

病态矩阵：条件数

κ

(

A

)

=

∥

A

∥

∥

A

−

1

∥

\kappa(A) = \|A\| \|A^{-1}\|

κ(A)=∥A∥∥A−1∥ 过大时，微小扰动导致解失真稳定性措施：

高斯法中使用列主元消去（Partial Pivoting）添加正则化项：$ (A + \lambda I)^{-1} $（Tikhonov正则化）用SVD分解替代：

A

−

1

=

V

Σ

−

1

U

T

A^{-1} = V \Sigma^{-1} U^T

A−1=VΣ−1UT（稳定但昂贵）

硬件加速：

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

9 总结与对比

算法时间复杂度空间复杂度适用场景高斯-约旦

O

(

n

3

)

O(n^3)

O(n3)

O

(

n

2

)

O(n^2)

O(n2)中小稠密矩阵LU分解

O

(

n

3

)

O(n^3)

O(n3)

O

(

n

2

)

O(n^2)

O(n2)需多次求逆的方程组伴随矩阵

O

(

n

!

)

O(n!)

O(n!)

O

(

n

2

)

O(n^2)

O(n2)

n

≤

3

n \leq 3

n≤3 的小矩阵牛顿迭代

O

(

k

n

2

)

O(kn^2)

O(kn2)

O

(

n

2

)

O(n^2)

O(n2)大规模稀疏矩阵分块法

O

(

n

2.8

)

O(n^{2.8})

O(n2.8)

O

(

n

2

)

O(n^2)

O(n2)分块并行计算特殊矩阵优化

O

(

n

)

O(n)

O(n)

O

(

n

)

O(n)

O(n)对角/三对角/正交矩阵

选择建议：

n

<

1000

n < 1000

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

未来方向：量子算法（HHL算法）可在

O

(

log

⁡

n

)

O(\log n)

O(logn) 时间求逆，但需量子硬件支持。

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