Python科学计算中的线性代数实战：NumPy与SciPy矩阵分解、特征值与稀疏矩阵完全指南

线性代数是科学计算的基石，无论是机器学习中的参数求解、信号处理中的频域变换，还是物理仿真中的偏微分方程离散化，都离不开矩阵运算。Python 生态中，NumPy 和 SciPy 提供了从基础矩阵操作到高级分解算法的完整工具链。本文将系统讲解如何使用这两个库完成实际工作中的线性代数任务，重点覆盖矩阵分解、特征值问题、奇异值分解（SVD）以及稀疏矩阵的高效处理。

一、NumPy 矩阵基础与高级操作

NumPy 的

ndarray

是一切数值计算的起点。虽然 NumPy 曾提供

np.matrix

类（已在 1.24 中弃用），但现代写法统一使用二维数组配合

运算符进行矩阵乘法。

1.1 矩阵创建与基本运算


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
import numpy as np



# 创建常见矩阵

A = np.array([[1, 2, 3],

              [4, 5, 6],

              [7, 8, 9]], dtype=np.float64)



# 单位矩阵、零矩阵、对角矩阵

I = np.eye(3)

Z = np.zeros((3, 3))

D = np.diag([1, 2, 3])



# 随机矩阵（可复现）

rng = np.random.default_rng(42)

R = rng.standard_normal((3, 3))



# 矩阵乘法：使用 @ 运算符

C = A @ D  # 等价于 np.matmul(A, D)



# 逐元素乘法

E = A * D  # 注意区分 @ 和 *



# 矩阵转置与共轭转置

print(A.T)           # 转置

print(A.conj().T)    # 共轭转置（对复矩阵有意义）

一个常见的陷阱是混淆

（矩阵乘法）和

<em>

（逐元素乘法）。在科学计算中，绝大多数场景需要的是

。如果你从 MATLAB 转过来，务必注意

</em>

在 NumPy 中的语义完全不同。

1.2 矩阵求逆与行列式


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# 矩阵求逆

A_inv = np.linalg.inv(A)



# 验证：A @ A_inv ≈ I

print(np.allclose(A @ A_inv, np.eye(3)))



# 行列式

det_A = np.linalg.det(A)

print(f"行列式: {det_A:.4f}")



# 迹（对角线元素之和）

trace_A = np.trace(A)

print(f"迹: {trace_A:.4f}")



# 秩

rank_A = np.linalg.matrix_rank(A)

print(f"秩: {rank_A}")

对于病态矩阵（条件数极大），直接求逆在数值上是不稳定的。实际应用中，应优先使用分解方法（如 LU 分解）来求解线性方程组，而非显式计算逆矩阵。NumPy 的

1	np.linalg.solve()

内部就采用了这种策略。

1.3 线性方程组求解


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# 求解 Ax = b

A = np.array([[3, 1], [1, 2]], dtype=np.float64)

b = np.array([9, 8], dtype=np.float64)



x = np.linalg.solve(A, b)

print(f"解: {x}")  # [2. 3.]



# 验证

print(np.allclose(A @ x, b))  # True



# 最小二乘解（超定系统）

A_over = np.array([[1, 1], [1, 2], [1, 3]], dtype=np.float64)

b_over = np.array([1, 2, 2], dtype=np.float64)

x_lstsq, residuals, rank, sv = np.linalg.lstsq(A_over, b_over, rcond=None)

print(f"最小二乘解: {x_lstsq}")

1	np.linalg.solve

要求矩阵是方阵且非奇异，而

1	np.linalg.lstsq

可以处理任意形状的矩阵，在数据拟合和回归分析中极为常用。

二、矩阵分解：LU、QR 与 Cholesky

矩阵分解是数值线性代数的核心工具。通过将矩阵拆解为具有特殊性质的矩阵乘积，不仅能加速线性方程组求解，还能揭示矩阵的结构信息。SciPy 的

1	scipy.linalg

模块提供了比 NumPy 更丰富的分解算法。

2.1 LU 分解

LU 分解将矩阵

A = PLU

，其中 P 是置换矩阵，L 是下三角矩阵，U 是上三角矩阵。它是求解线性方程组和高斯消元法的标准化实现。


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
from scipy.linalg import lu, lu_factor, lu_solve



A = np.array([[2, 1, 1],

              [4, 3, 3],

              [8, 7, 9]], dtype=np.float64)



# 完整LU分解

P, L, U = lu(A)

print("P:

", P)

print("L:

", L)

print("U:

", U)

print("验证 PA=LU:", np.allclose(P @ A, L @ U))



# 实际求解：lu_factor + lu_solve 更高效

# 适用于同一 A、多个 b 的场景

b1 = np.array([1, 2, 3], dtype=np.float64)

b2 = np.array([4, 5, 6], dtype=np.float64)



lu_piv = lu_factor(A)          # 只分解一次

x1 = lu_solve(lu_piv, b1)     # 求解第一个

x2 = lu_solve(lu_piv, b2)     # 求解第二个

print(f"x1: {x1}, x2: {x2}")

当需要用同一个矩阵 A 求解多个不同的 b 时，

lu_factor

lu_solve

的两步法比每次调用

1	np.linalg.solve

高效得多，因为消元过程只执行一次。

2.2 QR 分解

QR 分解将

A = QR

，其中 Q 是正交矩阵，R 是上三角矩阵。它是求解最小二乘问题、计算特征值（QR 算法）的基础。


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
from scipy.linalg import qr



A = np.array([[1, 2, 3],

              [4, 5, 6],

              [7, 8, 10],

              [1, 1, 1]], dtype=np.float64)



# 完整QR分解

Q, R = qr(A, mode='full')

print(f"Q shape: {Q.shape}")  # (4, 4)

print(f"R shape: {R.shape}")  # (4, 3)

print("Q正交性:", np.allclose(Q.T @ Q, np.eye(4)))



# 紧凑QR分解（经济型）

Q_eco, R_eco = qr(A, mode='economic')

print(f"Q_eco shape: {Q_eco.shape}")  # (4, 3)



# 用QR分解求解最小二乘问题

b = np.array([1, 2, 3, 4], dtype=np.float64)

# Rx = Q^T b，其中 x 只取前3个方程

x_qr = np.linalg.solve(R_eco, Q_eco.T @ b)

print(f"QR最小二乘解: {x_qr}")

2.3 Cholesky 分解

对于对称正定矩阵，Cholesky 分解

1	A = L * L^T

比 LU 分解快约两倍，且数值稳定性更好。在贝叶斯统计、协方差矩阵处理和优化算法中极为常用。


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
from scipy.linalg import cholesky, cho_solve



# 构造对称正定矩阵

A = np.array([[4, 2, 1],

              [2, 5, 3],

              [1, 3, 6]], dtype=np.float64)



# 验证正定性

eigvals = np.linalg.eigvalsh(A)

print("特征值:", eigvals)  # 全部为正



# Cholesky 分解（默认返回上三角）

U = cholesky(A)

print("U:

", U)

print("验证 U^T U = A:", np.allclose(U.T @ U, A))



# 下三角形式

L = cholesky(A, lower=True)

print("验证 L L^T = A:", np.allclose(L @ L.T, A))



# 高效求解 Ax = b

b = np.array([7, 10, 10], dtype=np.float64)

c_and_lower = (L, True)  # 传入L和lower标志

x = cho_solve(c_and_lower, b)

print(f"解: {x}")

实际项目中，当你知道矩阵是对称正定的（如协方差矩阵、Gram 矩阵），务必使用 Cholesky 分解而非通用的 LU 分解——不仅速度快，精度也更高。

三、特征值与奇异值分解

特征值分解和奇异值分解（SVD）是理解矩阵结构的两把钥匙。前者只适用于方阵，后者适用于任意矩阵。在主成分分析（PCA）、推荐系统、图像压缩等场景中，SVD 尤其不可或缺。

3.1 特征值分解


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# 对称矩阵的特征值分解（保证实特征值和正交特征向量）

A_sym = np.array([[2, -1, 0],

                  [-1, 2, -1],

                  [0, -1, 2]], dtype=np.float64)



eigenvalues, eigenvectors = np.linalg.eigh(A_sym)

print(f"特征值: {eigenvalues}")

print(f"特征向量矩阵正交性: {np.allclose(eigenvectors.T @ eigenvectors, np.eye(3))}")



# 一般方阵的特征值分解

A_gen = np.array([[1, 2, 3],

                  [0, 4, 5],

                  [0, 0, 6]], dtype=np.float64)



vals, vecs = np.linalg.eig(A_gen)

print(f"特征值: {vals}")  # 1, 4, 6



# 验证 Av = λv

for i in range(len(vals)):

    v = vecs[:, i]

    print(f"λ={vals[i]:.2f}: Av ≈ λv? {np.allclose(A_gen @ v, vals[i] * v)}")

1	np.linalg.eigh

专门用于对称/埃尔米特矩阵，结果保证是实数且特征向量正交。

1	np.linalg.eig

用于一般方阵，可能返回复数特征值。如果矩阵已知对称，务必使用

eigh

，它更快且更稳定。

3.2 奇异值分解（SVD）


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# SVD: A = U Σ V^T

A = np.array([[1, 2, 3],

              [4, 5, 6],

              [7, 8, 9],

              [10, 11, 12]], dtype=np.float64)



U, s, Vt = np.linalg.svd(A, full_matrices=False)

print(f"U shape: {U.shape}")    # (4, 3)

print(f"s shape: {s.shape}")    # (3,)

print(f"Vt shape: {Vt.shape}")  # (3, 3)



# 验证

Sigma = np.diag(s)

print("重构:", np.allclose(U @ Sigma @ Vt, A))



# 低秩近似：用前k个奇异值重构

k = 2

A_approx = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]

print(f"原矩阵F范数: {np.linalg.norm(A):.4f}")

print(f"近似误差: {np.linalg.norm(A - A_approx):.4f}")

print(f"信息保留率: {1 - np.linalg.norm(A - A_approx)/np.linalg.norm(A):.2%}")

3.3 SVD 实战：图像压缩

低秩近似最直观的应用就是图像压缩。保留前 k 个奇异值，可以在损失一定精度的前提下大幅减少存储量。


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
from PIL import Image

import numpy as np



# 读取灰度图像

img = Image.open('photo.jpg').convert('L')

A = np.array(img, dtype=np.float64)



U, s, Vt = np.linalg.svd(A, full_matrices=False)



# 不同压缩率的对比

for k in [10, 50, 100, 200]:

    A_k = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]

    A_k = np.clip(A_k, 0, 255).astype(np.uint8)

    ratio = (k * (A.shape[0] + A.shape[1] + 1)) / (A.shape[0] * A.shape[1])

    print(f"k={k}: 压缩比={ratio:.2%}, 存储量={ratio:.2%}")

    # Image.fromarray(A_k).save(f'compressed_k{k}.jpg')

对于一张 1000×1000 的灰度图，k=50 时存储量仅为原图的 10.05%，但视觉上仍然保留了主要轮廓和结构。这就是 SVD 在数据压缩和降维中的威力。

四、稀疏矩阵：大规模问题的内存与计算优化

在有限元分析、社交网络图计算、自然语言处理的词-文档矩阵等场景中，矩阵通常99%以上都是零元素。如果用密集矩阵存储，内存和计算时间都会严重浪费。SciPy 的

1	scipy.sparse

模块提供了多种稀疏存储格式和对应的运算。

4.1 稀疏矩阵存储格式选择

格式	全称	适用场景	修改效率	运算效率
COO	Coordinate	矩阵构建、FEM组装	高	低
CSR	Compressed Sparse Row	行切片、矩阵乘法Ax	低	高
CSC	Compressed Sparse Column	列切片、矩阵乘法x^TA	低	高
DIA	Diagonal	对角线结构（如三对角矩阵）	低	极高
LIL	List of Lists	逐元素修改、构建稀疏模式	极高	低


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
from scipy import sparse



# 从坐标创建 COO 矩阵

row = np.array([0, 1, 2, 3, 0])

col = np.array([0, 1, 2, 3, 4])

data = np.array([1, 2, 3, 4, 5])

A_coo = sparse.coo_matrix((data, (row, col)), shape=(5, 5))

print(A_coo)

print(f"非零元素数: {A_coo.nnz}")

print(f"密集存储: {A_coo.toarray()}")



# 转换为 CSR 用于运算

A_csr = A_coo.tocsr()



# 对角线矩阵（极其高效）

diag_data = [1, 2, 3, 4, 5]

A_dia = sparse.diags(diag_data, offsets=0, shape=(5, 5), format='csr')



# 三对角矩阵

A_tridiag = sparse.diags(

    [np.ones(4), 2*np.ones(5), np.ones(4)],

    offsets=[-1, 0, 1],

    shape=(5, 5),

    format='csr'

)

print("三对角矩阵:

", A_tridiag.toarray())

4.2 稀疏矩阵运算


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
from scipy.sparse import csr_matrix, random as sparse_random



# 生成大型稀疏矩阵

n = 10000

A = sparse_random(n, n, density=0.001, format='csr', dtype=np.float64)

b = rng.standard_normal(n)



# 稀疏矩阵-向量乘法

y = A @ b  # 自动识别稀疏格式



# 稀疏矩阵乘法

C = A @ A.T  # 结果仍是稀疏矩阵



# 内存对比

dense_size = n * n * 8 / (1024**3)  # GB

sparse_size = A.data.nbytes / (1024**2)  # MB

print(f"密集存储需要: {dense_size:.2f} GB")

print(f"稀疏存储仅需: {sparse_size:.2f} MB")

对于 10000×10000、密度 0.1% 的矩阵，密集存储需要约 745 GB，而稀疏存储只需不到 1 MB。这种差距使得原本无法放入内存的问题变得可解。

4.3 稀疏线性方程组求解

SciPy 提供了专门的稀疏求解器，分为直接法和迭代法两类：


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
from scipy.sparse.linalg import spsolve, cg, gmres, LinearOperator



# 构造稀疏正定系统

n = 5000

diag_main = 4 * np.ones(n)

diag_off = -1 * np.ones(n - 1)

A_sp = sparse.diags([diag_off, diag_main, diag_off], [-1, 0, 1], format='csr')

b = np.ones(n)



# 直接法（适合中小规模）

x_direct = spsolve(A_sp, b)

print(f"直接法残差: {np.linalg.norm(A_sp @ x_direct - b):.2e}")



# 共轭梯度法（适合大规模对称正定系统）

x_cg, info = cg(A_sp, b)

print(f"CG收敛: {info == 0}, 残差: {np.linalg.norm(A_sp @ x_cg - b):.2e}")



# GMRES（适合一般非对称系统）

x_gmres, info = gmres(A_sp, b)

print(f"GMRES收敛: {info == 0}, 残差: {np.linalg.norm(A_sp @ x_gmres - b):.2e}")

直接法（spsolve）：基于 SuperLU，适合维度在几万以内、需要高精度的场景
CG（共轭梯度）：只适用于对称正定矩阵，收敛快，内存占用低
GMRES：通用迭代求解器，适合一般非对称系统，可配合预条件子加速收敛

五、特征值问题的进阶应用

5.1 大规模稀疏特征值问题

当矩阵维度达到十万甚至百万量级时，

1	np.linalg.eigh

无法使用（需要密集存储）。SciPy 提供了 ARPACK 接口，可以高效地只计算前 k 个最大或最小的特征值。


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
from scipy.sparse.linalg import eigsh, eigs



# 大型稀疏对称矩阵

n = 10000

A_sp = sparse.diags(

    [-np.ones(n-1), 2*np.ones(n), -np.ones(n-1)],

    [-1, 0, 1], format='csr'

)



# 只求前5个最小特征值（用于振动分析等）

eigvals, eigvecs = eigsh(A_sp, k=5, which='SM')

print(f"5个最小特征值: {eigvals}")



# 只求前5个最大特征值（用于PCA等）

eigvals_max, eigvecs_max = eigsh(A_sp, k=5, which='LM')

print(f"5个最大特征值: {eigvals_max}")



# 非对称矩阵用 eigs

eigvals_gen, _ = eigs(A_sp.astype(complex), k=3)

5.2 广义特征值问题

在结构动力学中，广义特征值问题

Ax = λBx

非常常见，其中 A 是刚度矩阵，B 是质量矩阵。


1
2
3
4
5
6
7
8
# 广义特征值问题：Kx = λMx

n = 100

K = sparse.diags([-np.ones(n-1), 2*np.ones(n), -np.ones(n-1)], [-1, 0, 1], format='csr')

M = sparse.eye(n, format='csr')



eigvals, eigvecs = eigsh(K, k=5, M=M, which='SM')

print(f"固有频率(平方): {eigvals}")

print(f"固有频率: {np.sqrt(np.abs(eigvals))}")

六、性能优化与最佳实践

6.1 条件数与数值稳定性


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# 条件数衡量矩阵对扰动的敏感度

A_well = np.array([[1, 0], [0, 1]], dtype=np.float64)

A_ill = np.array([[1, 1], [1, 1.0001]], dtype=np.float64)



print(f"良态矩阵条件数: {np.linalg.cond(A_well):.2f}")

print(f"病态矩阵条件数: {np.linalg.cond(A_ill):.2f}")



# 病态系统：微小扰动导致解的大幅变化

b = np.array([1, 1], dtype=np.float64)

x1 = np.linalg.solve(A_ill, b)

b_perturbed = np.array([1, 1.0001], dtype=np.float64)

x2 = np.linalg.solve(A_ill, b_perturbed)

print(f"原解: {x1}")

print(f"扰动解: {x2}")

print(f"相对误差: {np.linalg.norm(x2 - x1) / np.linalg.norm(x1):.2%}")

条件数超过 10^15 时，双精度浮点数的解已经不可信。此时应考虑正则化方法（如 Tikhonov 正则化）或使用 SVD 截断。

6.2 BLAS/LAPACK 后端选择

NumPy 和 SciPy 的底层线性代数运算依赖 BLAS/LAPACK 实现。不同后端的性能差异可能达到 10 倍以上：


1
2
3
4
5
6
# 查看当前使用的BLAS后端

np.show_config()



# OpenBLAS（conda默认）: 性能良好，通用

# MKL（pip安装的numpy-mkl）: Intel CPU上最优

# BLIS（新兴选择）: 兼容性好，性能接近MKL

建议在服务器部署时通过

1	conda install numpy scipy mkl

安装 MKL 后端，在个人开发环境中 OpenBLAS 即可满足需求。

6.3 批量运算避免循环


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# 不推荐：循环求解多个系统

n_systems = 1000

A = np.eye(3)

B = rng.standard_normal((n_systems, 3))



# 方法1: 循环

results_loop = [np.linalg.solve(A, b) for b in B]



# 方法2: 利用广播（A固定时更高效）

# 对于对角矩阵或特殊结构矩阵，手动展开

A_diag = np.diag([2.0, 3.0, 4.0])

x_batch = B / np.diag(A_diag)  # 对角矩阵的解直接除以对角元素

print(f"批量解 shape: {x_batch.shape}")



# 方法3: 多右端项求解

B_multi = rng.standard_normal((3, n_systems))

X_multi = np.linalg.solve(A, B_multi)  # 一次调用解所有

七、实战案例：使用 SVD 进行 PCA 降维

主成分分析（PCA）是 SVD 最经典的应用之一。下面用一个完整示例展示从数据标准化到降维可视化的全流程：


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
from sklearn.datasets import load_iris

from sklearn.preprocessing import StandardScaler



# 加载数据

data = load_iris()

X = data.data  # (150, 4)

y = data.target



# 标准化

scaler = StandardScaler()

X_scaled = scaler.fit_transform(X)



# 方法1: 用SVD直接做PCA

U, s, Vt = np.linalg.svd(X_scaled, full_matrices=False)



# 主成分 = Vt 的行

components = Vt  # (4, 4)



# 投影到前2个主成分

X_pca = U[:, :2] @ np.diag(s[:2])  # 或 X_scaled @ Vt[:2].T



# 方法2: 用特征值分解（等价）

cov_matrix = X_scaled.T @ X_scaled / (X_scaled.shape[0] - 1)

eigvals, eigvecs = np.linalg.eigh(cov_matrix)

# 特征值按升序，取最大的两个

idx = np.argsort(eigvals)[::-1]

X_pca2 = X_scaled @ eigvecs[:, idx[:2]]



print(f"方差解释比: {s**2 / np.sum(s**2)}")

print(f"前2个主成分解释: {np.sum(s[:2]**2) / np.sum(s**2):.2%}")

两种方法在数学上等价，但 SVD 方法更稳定——它不需要显式计算协方差矩阵，避免了矩阵乘法引入的额外数值误差。这也是 scikit-learn 内部 PCA 实现采用 SVD 的原因。

总结

Python 线性代数工具链的核心要点可以归纳为以下几点：

矩阵运算优先用
1
@

，不要用

1
*

（那是逐元素乘法）
求解线性系统用
1
solve

，不要显式求逆——数值更稳定、速度更快
对称正定矩阵用 Cholesky 分解，比 LU 分解快一倍且更稳定
对称矩阵的特征值用
1
eigh

，保证实特征值和正交特征向量
稀疏矩阵选对格式：构建用 COO/LIL，运算用 CSR/CSC
大规模稀疏特征值问题用 ARPACK（
1
eigsh

/

1
eigs

），只算需要的特征值
PCA 用 SVD，不计算协方差矩阵，数值更稳定
关注条件数，条件数过大时考虑正则化

掌握这些工具和原则，足以应对绝大多数科学计算和数据分析中的线性代数需求。NumPy 处理中小规模的密集矩阵，SciPy 接管大规模稀疏问题和高级分解——两者配合，构成了 Python 科学计算最坚实的底座。