NumPy 的线性代数与随机数

周二下午两点，某上海私募的股票池经理把你叫到工位前：要 600519.SH 对沪深300 ETF（510300.SH）的市场 β，日简单收益（daily simple return），近252个交易日窗口，今晚9点前要见。教科书答案一行就能解决：beta = Cov(r_stock, r_mkt) / Var(r_mkt)。工程答案稍长：把 [1, r_mkt] 当设计矩阵，跑一遍普通最小二乘（ordinary least squares, OLS）回归。两条路径都要会，因为只有第二条能扩展到多因子（multi-factor）情形。再加一个 Monte-Carlo 自检：六行 NumPy 模拟几何布朗运动（geometric Brownian motion, GBM）的价格路径，构造一只「已知 β=0.85」的合成股票，让 OLS 把它反推回来。这一整套都跑在 NumPy 最常用的两个子包上：np.linalg 与 np.random。

矩阵乘法与形状规则

2026 年的 NumPy 里，​​矩阵乘法​​（matrix multiplication）的写法就是 PEP 465 引入的中缀算子 @：A @ B 等价于 np.matmul(A, B)。二维情形的形状规则记住一条：

import numpy as np
A = np.array([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]])    # shape (3, 2)
B = np.array([[1.0, 0.0, 1.0], [0.0, 1.0, 1.0]])      # shape (2, 3)
C = A @ B                                              # shape (3, 3)

(m, k) @ (k, n) -> (m, n)，中间维度 k 必须吻合并被吸收。把它和逐元素乘 A * B 分清楚：后者要求两个数组形状相同（或可广播），是位置对齐的乘法——这是新手最常把矩阵乘写成 * 的源头，也是 code review 里最高频的 NumPy bug 类型之一。np.dot 是更老的 API：二维下 np.dot(A, B) 与 A @ B 完全一致；张量秩高于 2 时两者的广播规则会发散，新代码统一写 @。底层一句话：@ 走 BLAS，单核 GEMM 已经接近物理上限，不需要你再去手工调优。

解线性方程组：solve 还是 inv？

要解 $A x = b$，你有两条路：先算​​矩阵逆​​（matrix inverse）再点乘 np.linalg.inv(A) @ b，或者直接 np.linalg.solve(A, b) 一步出结果。后者内部走 LU 分解（LU factorisation），只做一次分解、不显式构造逆矩阵——对病态（ill-conditioned）的 $A$ 浮点误差放大要小得多，工程上还大约快三倍。

A = np.array([[3.0, 1.0], [1.0, 2.0]])
b = np.array([9.0, 8.0])
x = np.linalg.solve(A, b)
assert np.allclose(A @ x, b)

经验法则：除非你真的需要 $A^{-1}$ 作为一个独立对象（量化里很少见），都用 solve。把残差 A @ x - b 用 np.allclose 验一下，是上线前最便宜的一道闸。

OLS 的三种解法

回到开头的 β 估计。设 r_mkt 与 r_stock 是两条长度 252 的 1-D 数组，分别为 510300.SH 与 600519.SH 的日简单收益。设计矩阵第一列为常数 1（截距 / intercept），第二列为市场收益；目标向量 y = r_stock。三种解法殊途同归：

# r_mkt: 510300.SH daily simple returns over 252 trading days
# r_stock: 600519.SH daily simple returns over the same window
# (Both are 1-D arrays of length 252; load from your data layer in production.)
X = np.column_stack([np.ones(252), r_mkt])
y = r_stock
coef_lstsq, *_ = np.linalg.lstsq(X, y, rcond=None)
coef_normal = np.linalg.solve(X.T @ X, X.T @ y)
L = np.linalg.cholesky(X.T @ X)
z = np.linalg.solve(L, X.T @ y)
coef_chol = np.linalg.solve(L.T, z)
assert np.allclose(coef_lstsq, coef_normal) and np.allclose(coef_lstsq, coef_chol)

coef[0] 是 alpha，coef[1] 是 beta；某次跑出来约 0.87 就是茅台对沪深300 的样本 β（仅作演示，不对应任何特定年度）。三个方法各有应用场景：lstsq 对任意 X 都安全，包括秩亏（rank-deficient）情形；​​正规方程​​（normal equations）X.T @ X 法在条件数良好时跑得最快，但条件数一上去就先吃亏；​​乔莱斯基分解​​（Cholesky decomposition）适合「同一个 Gram 矩阵反复求解」的场景，例如交叉验证时设计矩阵不变、目标 y 切换——X.T @ X 在 X 满列秩时是正定的，正好满足 Cholesky 的前提。rcond=None 是 lstsq 的现代写法，用来切换到新版奇异值阈值约定并抑制弃用警告。三种结果用 np.allclose 互相对账，是检验你的实现没写错的最小自检。

np.linalg 的其余工具与边界

np.linalg.norm(v) 默认给 L2 范数（L2 norm），传 ord=1 或 ord=np.inf 切换其他范数。np.linalg.det(A) 算行列式（determinant），实际工作中很少直接用，更多是教学价值。np.linalg.eig(A) 返回特征值（eigenvalue）与特征向量（eigenvector）；金融里它的经典用途是对​​协方差矩阵​​（covariance matrix）做主成分分析（principal component analysis, PCA），但在工程实践里更稳定的做法是对中心化后的数据矩阵做奇异值分解（singular value decomposition, SVD）np.linalg.svd(M)。面向生产级数值线性代数（例如对称矩阵的 eigh、可指定 LAPACK 驱动的 lstsq、广义特征值问题）你会转向 scipy.linalg，在 3.2.3 课讲授。

现代随机数生成器（Generator API）

np.random 在 NumPy 1.17（2019 年）引入了新的 Generator 接口。新代码只用 np.random.default_rng(seed)，不要碰旧的 np.random.seed / np.random.rand / np.random.randn 这一族裸函数——它们共享一个全局 RandomState，不是线程安全的（thread-safe），早已不推荐用于新开发。

rng = np.random.default_rng(seed=42)
z = rng.normal(loc=0.0, scale=1.0, size=10)
u = rng.uniform(low=0.0, high=1.0, size=10)
pick = rng.choice(np.arange(100), size=5, replace=False)
perm = rng.permutation(10)

种子契约（seed contract）：同一个 seed 在同一个 NumPy 版本下产出字节级一致（byte-identical）的样本，跨机器、跨进程都成立——这正是「可复现回测」（reproducible backtest）的物理基础。每个独立任务各建一个 rng 实例，避免共享可变全局状态。

几何布朗运动模拟

GBM 的离散更新写成：

$S_{t+1} = S_t \cdot \exp\left( \left(\mu - \tfrac{1}{2}\sigma^2\right)\Delta t + \sigma\sqrt{\Delta t}\, Z_t \right), \quad Z_t \sim \mathcal{N}(0, 1).$

其中 $\mu$ 为年化漂移率（annualised drift），$\sigma$ 为年化波动率（annualised volatility），$\Delta t$ 为一个时间步（time step）的长度，$Z_t$ 为独立标准正态。连续时间 SDE 的推导与显式解的证明留给 2.7.1 布朗运动课；这一课只用它的离散写法当作向量化练习。一次性把 (T, N) 张标准正态全部抽出，沿时间轴累加对数收益，再 vstack 上初值：

mu = 0.06
sigma = 0.20
dt = 1/252
T = 252
N = 5
S0 = 5.0
rng = np.random.default_rng(seed=42)
Z = rng.normal(size=(T, N))
increments = (mu - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * Z
log_returns = increments
S = S0 * np.exp(np.cumsum(log_returns, axis=0))
S = np.vstack([np.full((1, N), S0), S])

参数取沪深300 ETF 量级的年化漂移 mu=0.06、年化波动 sigma=0.20、初值 S0=5.0、T=252 天、N=5 条路径。整段模拟里没有一个 Python 层的循环，全部算力都落在底层的 ufunc 与 BLAS 调用上。回到开头那道「Monte-Carlo 自检」：用同一个 rng 抽一份 r_mkt，令 r_stock = 0.85 * r_mkt + 残差噪声，再把上一节的 OLS 代码原样跑一遍——coef[1] 应当收敛在 0.85 周围的 Monte-Carlo 误差带内。「先模拟、再让估计器把真值反推出来」（simulate-then-estimate）是任何统计代码上线前应当做的最小一道闸：估计器对自己生成的真值都恢复不了，对真实数据就更不必信。

练习

到这里你已经能在不写 Python 循环的前提下，用 np.linalg 跑出一条 β，用 np.random.default_rng 生成一段可复现的 GBM 价格路径，并能在脑中分清 solve 与 inv、@ 与 *、新版 Generator 与旧版 RandomState 的四道分水岭。下一课进入 Pandas（3.2.2）：当你的收益矩阵需要一个日期索引、一列股票代码，并要和一张 A 股财报日历做 asof 合并时，你就已经溢出了原始 NumPy 的能力边界——DataFrame 是为这些「带标签的二维表」准备的下一层抽象。