NumPy进阶：广播、线性代数、随机数与性能
NumPy进阶：广播、线性代数、随机数与性能
在掌握NumPy基础之后，本文将带你深入了解广播机制、线性代数运算、随机数生成以及性能优化等进阶主题，帮助你编写更高效、更优雅的数值计算代码。
广播机制
广播（Broadcasting）是NumPy中强大的特性，允许不同形状的数组进行算术运算，无需显式复制数据。规则如下：
- 两个数组的维度从右向左逐一比较
- 维度大小相等或其中一个为1时可以兼容
- 缺失的维度视为1
import numpy as np
# 标量与数组运算（自动广播）
arr = np.array([[1, 2, 3], [4, 5, 6]])
result = arr + 10
print(result)
# [[11 12 13]
# [14 15 16]]
# 行向量与矩阵
row = np.array([10, 20, 30])
print(arr + row)
# [[11 22 33]
# [14 25 36]]
# 列向量与矩阵
col = np.array([[100], [200]])
print(arr + col)
# [[101 102 103]
# [204 205 206]]
# 计算每个样本与均值的差值
data = np.array([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
mean = data.mean(axis=0) # 按列求均值 [4. 5. 6.]
centered = data - mean # 广播减去均值
print(centered)
# [[-3. -3. -3.]
# [ 0. 0. 0.]
# [ 3. 3. 3.]]
线性代数运算
NumPy的linalg模块提供了丰富的线性代数函数：
import numpy as np
# 矩阵乘法
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
# 三种等价方式
C = np.dot(A, B)
C = A @ B
C = np.matmul(A, B)
print(C)
# [[19 22]
# [43 50]]
# 转置
print(A.T)
# [[1 3]
# [2 4]]
# 行列式
det = np.linalg.det(A)
print(f"行列式: {det}") # -2.0
# 逆矩阵
inv = np.linalg.inv(A)
print(f"逆矩阵:\n{inv}")
# 特征值和特征向量
eigenvalues, eigenvectors = np.linalg.eig(A)
print(f"特征值: {eigenvalues}")
print(f"特征向量:\n{eigenvectors}")
# 求解线性方程组 Ax = b
b = np.array([5, 11])
x = np.linalg.solve(A, b)
print(f"解: {x}")
# 验证
print(f"验证: {A @ x}") # [5. 11.]
# SVD奇异值分解
U, S, Vt = np.linalg.svd(A)
print(f"奇异值: {S}")
# 矩阵范数
norm = np.linalg.norm(A)
print(f"Frobenius范数: {norm:.4f}")
随机数生成
NumPy的随机数模块支持多种概率分布，是模拟和统计分析的基础：
import numpy as np
# 基础随机数
rng = np.random.default_rng(seed=42)
# 均匀分布 [0, 1)
uniform = rng.random((3, 3))
print("均匀分布:\n", uniform)
# 正态分布
normal = rng.normal(loc=0, scale=1, size=(3, 3))
print("正态分布:\n", normal)
# 整数随机
integers = rng.integers(0, 100, size=(3, 4))
print("随机整数:\n", integers)
# 随机选择
choices = rng.choice([10, 20, 30, 40, 50], size=3, replace=False)
print("不重复随机选择:", choices)
# 打乱数组
arr = np.arange(10)
rng.shuffle(arr)
print("打乱后:", arr)
# 常用分布
beta = rng.beta(2, 5, size=100) # Beta分布
chi2 = rng.chisquare(3, size=100) # 卡方分布
poisson = rng.poisson(lam=5, size=100) # 泊松分布
# 设置种子保证可重复性
rng1 = np.random.default_rng(seed=123)
rng2 = np.random.default_rng(seed=123)
print(rng1.random(3) == rng2.random(3)) # True
性能优化技巧
NumPy的性能优势需要正确使用才能充分发挥：
import numpy as np
import time
# 避免Python循环，使用向量化
def slow_sum(n):
result = 0
for i in range(n):
result += i
return result
def fast_sum(n):
return np.sum(np.arange(n))
# 性能对比
n = 10_000_000
start = time.time()
slow_sum(n)
print(f"循环方式: {time.time()-start:.4f}秒")
start = time.time()
fast_sum(n)
print(f"向量化: {time.time()-start:.4f}秒")
# 预分配内存而非动态扩展
def bad_append(n):
arr = np.array([])
for i in range(n):
arr = np.append(arr, i)
return arr
def good_prealloc(n):
arr = np.empty(n)
for i in range(n):
arr[i] = i
return arr
# 使用视图而非复制
arr = np.arange(100)
view = arr[::2] # 视图，不复制数据
copy = arr[::2].copy() # 显式复制
# 就地操作减少内存分配
a = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
np.add(a, 1.0, out=a) # 就地加法
结构化数组
结构化数组允许创建类似数据库表的数据结构：
import numpy as np
# 定义结构化数据类型
dt = np.dtype([('name', 'U10'), ('age', 'i4'), ('salary', 'f8')])
# 创建结构化数组
employees = np.array([
('Alice', 30, 75000.0),
('Bob', 25, 65000.0),
('Charlie', 35, 85000.0),
], dtype=dt)
print(employees['name']) # ['Alice' 'Bob' 'Charlie']
print(employees['salary']) # [75000. 65000. 85000.]
# 按条件筛选
high_salary = employees[employees['salary'] > 70000]
print(f"高薪员工: {high_salary['name']}")
实战：蒙特卡洛模拟估算π
import numpy as np
def estimate_pi(n_samples):
rng = np.random.default_rng(seed=42)
x = rng.random(n_samples)
y = rng.random(n_samples)
# 判断点是否在单位圆内
inside = x**2 + y**2 <= 1
# π ≈ 4 * (圆内点数 / 总点数)
pi_estimate = 4 * np.sum(inside) / n_samples
return pi_estimate
# 不同样本量的估算
for n in [1000, 10000, 100000, 1000000]:
pi_est = estimate_pi(n)
error = abs(pi_est - np.pi)
print(f"样本数: {n:>10}, 估算π: {pi_est:.6f}, 误差: {error:.6f}")
总结
掌握NumPy的进阶特性——广播、线性代数、随机数和性能优化——能让你在数据科学工作中事半功倍。记住：优先使用向量化操作，避免不必要的循环和数据复制。