03-线性代数与优化

线性代数与优化

预计学习时间：2-3 小时

难度：⭐⭐⭐

核心问题：如何在多维空间中高效地表示、变换和优化投资组合？

---

概述

线性代数是量化投资的”语法”，优化理论则是”引擎”。本节涵盖：

• 矩阵运算：NumPy 高效矩阵操作
• 特征分解与 SVD：矩阵的深层结构
• PCA：因子降维与正交化
• 协方差矩阵：估计、收缩与指数加权
• 凸优化与 KKT：理解优化的理论基础
• 二次规划：均值方差组合优化
• 拉格朗日乘子法：约束优化的通用框架
• L1/L2 正则化：从优化视角理解过拟合
• 正交投影与最小二乘：回归的几何本质

---

一、矩阵运算基础

1.1 量化中的矩阵形式

概念	维度	量化含义
收益率矩阵	$(T\times n)$	T 日、n 只股票
因子暴露矩阵	$(n\times p)$	n 只股票、p 个因子
协方差矩阵	$(n\times n)$	股票间风险结构
权重向量	$(n\times 1)$	投资组合配置

1.2 NumPy 核心操作

import numpy as np
 
np.random.seed(42)
n_stocks, n_days, n_factors = 10, 252, 5
 
# 收益率矩阵 (T, n)
returns = np.random.normal(0.0005, 0.02, (n_days, n_stocks))
# 因子暴露 (n, p)
factor_exp = np.random.normal(0, 1, (n_stocks, n_factors))
# 权重向量 (n,)
weights = np.array([0.12, 0.08, 0.15, 0.10, 0.05, 0.10, 0.08, 0.12, 0.10, 0.10])
 
# 1. 组合日收益 = 收益 @ 权重
port_ret = returns @ weights
print(f"年化收益: {port_ret.mean()*252:.4f}")
print(f"年化波动: {port_ret.std()*np.sqrt(252):.4f}")
 
# 2. 协方差矩阵
cov = np.cov(returns, rowvar=False)
 
# 3. 组合方差 = w^T Σ w
port_var = weights @ cov @ weights
print(f"组合波动率: {np.sqrt(port_var):.4f}")
 
# 4. 因子收益 = 暴露^T @ 收益
fret = factor_exp @ returns.T
print(f"因子收益维度: {fret.shape}")
 
# 5. 组合因子暴露 = w^T F
port_fexp = weights @ factor_exp
print(f"组合因子暴露: {np.round(port_fexp, 3)}")
 
# 6. 边际风险贡献 MRC = (Σw) / σ_p
mrc = (cov @ weights) / np.sqrt(port_var)
print(f"MRC 之和（应=波动率）: {mrc.sum():.4f}")
 
# 7. 广播：标准化因子矩阵
fmean = factor_exp.mean(axis=0)
fstd = factor_exp.std(axis=0)
factor_z = (factor_exp - fmean) / fstd
print(f"标准化后均值: {np.round(factor_z.mean(axis=0), 6)}")

---

二、特征分解与 SVD

2.1 特征分解

对称矩阵（如协方差矩阵）的特征分解：$\Sigma =Q\Lambda Q^T$

• $Q$：正交矩阵（特征向量），$Q^{T}Q=I$
• $\Lambda$：对角矩阵（特征值 = 各方向上的方差）

2.2 SVD（适用于任意矩阵）

$X=U\Sigma V^T$

特性	特征分解	SVD
适用矩阵	对称矩阵	任意矩阵
金融对象	协方差矩阵	因子矩阵、收益矩阵
输出	特征值（可负）	奇异值（非负）

import numpy as np
 
np.random.seed(42)
 
# 特征分解：协方差矩阵
n = 20
returns = np.random.normal(0, 0.02, (200, n))
cov = np.cov(returns, rowvar=False)
 
eigenvalues, eigenvectors = np.linalg.eigh(cov)
sorted_idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[sorted_idx]
eigenvectors = eigenvectors[:, sorted_idx]
 
print("特征分解（协方差矩阵）:")
print(f"{'PC':<6} {'特征值':<10} {'方差比':<10} {'累计'}")
print("-" * 38)
cum = 0
for i, ev in enumerate(eigenvalues[:10]):
    r = ev / eigenvalues.sum()
    cum += r
    print(f"PC{i+1:<3} {ev:<10.6f} {r:<10.2%} {cum:.2%}")
 
# SVD：因子矩阵
n_stocks, n_factors = 100, 20
base = np.random.normal(0, 1, (n_stocks, 3))
noise = np.random.normal(0, 0.1, (n_stocks, n_factors - 3))
X = base @ np.random.normal(0, 1, (3, n_factors - 3)) + noise
 
U, S, Vt = np.linalg.svd(X, full_matrices=False)
explained = np.cumsum(S**2 / np.sum(S**2))
print(f"\nSVD: 前3个奇异值解释 {explained[2]:.1%} 方差")
 
# 低秩近似
k = 3
X_approx = U[:, :k] @ np.diag(S[:k]) @ Vt[:k, :]
print(f"前 {k} 个奇异值近似误差: {np.linalg.norm(X - X_approx)/np.linalg.norm(X):.4f}")

---

三、PCA（主成分分析）

3.1 原理

找投影方向使投影后方差最大：${\max }_{\mathbf{w}}{\mathbf{w}}^T\Sigma \mathbf{w}$, s.t. $\Vert \mathbf{w}\Vert =1$。解即协方差矩阵特征向量。

3.2 NumPy 手写 PCA + sklearn 对比

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
 
np.random.seed(42)
n_stocks, n_days = 50, 500
 
# 3 个隐性因子驱动
latent = np.random.normal(0, 1, (n_days, 3))
loadings = np.random.normal(0, 1, (3, n_stocks))
specific = np.random.normal(0, 0.5, (n_days, n_stocks))
returns = latent @ loadings + specific
 
# ---- 手写 PCA ----
def manual_pca(X, n_components=10):
    """手写 PCA（基于协方差矩阵特征分解）"""
    mean = X.mean(axis=0)
    X_c = X - mean
    cov = np.cov(X_c, rowvar=False)
    eigenvalues, eigenvectors = np.linalg.eigh(cov)
    idx = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]
    components = eigenvectors[:, :n_components].T
    Z = X_c @ components.T
    ratio = eigenvalues[:n_components] / eigenvalues.sum()
    return Z, ratio, components, mean
 
Z_manual, ratio_manual, _, _ = manual_pca(returns, 10)
 
# ---- sklearn PCA ----
pca_sk = PCA(n_components=10)
Z_sk = pca_sk.fit_transform(returns)
 
# 对比
print("手写 vs sklearn 前 5 个方差解释比:")
print(f"  手写:   {ratio_manual[:5].round(4)}")
print(f"  sklearn:{pca_sk.explained_ratio_[:5].round(4)}")
 
# 绘图
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
axes[0].bar(range(1, 11), ratio_manual, color='steelblue', alpha=0.7)
axes[0].set_title('PCA 碎石图'); axes[0].set_xlabel('主成分'); axes[0].set_ylabel('特征值')
 
cum = np.cumsum(ratio_manual)
axes[1].plot(range(1, 11), cum, 'o-', color='darkorange', lw=2)
axes[1].axhline(0.95, color='red', ls='--', label='95%')
axes[1].axhline(0.90, color='green', ls='--', label='90%')
axes[1].set_title('累计方差解释比'); axes[1].legend()
plt.tight_layout()
plt.show()
 
n_95 = np.argmax(cum >= 0.95) + 1
print(f"\n解释 95% 方差需 {n_95} 个主成分（原始 {n_stocks} 维）")

---

四、协方差矩阵估计

4.1 问题：大 n 小 T 时样本协方差病态

4.2 三种方法对比

import numpy as np
from sklearn.covariance import LedoitWolf
 
np.random.seed(42)
n_stocks, n_days = 50, 120
latent = np.random.normal(0, 1, (n_days, 3))
returns = latent @ np.random.normal(0, 0.5, (3, n_stocks)) + np.random.normal(0, 0.01, (n_days, n_stocks))
 
# 1. 样本协方差
sample_cov = np.cov(returns, rowvar=False)
 
# 2. Ledoit-Wolf 收缩
lw = LedoitWolf()
lw_cov = lw.fit(returns).covariance_
print(f"Ledoit-Wolf 收缩强度: {lw.shrinkage_:.4f}")
 
# 3. EWMA
def ewma_cov(ret, lam=0.94):
    """指数加权协方差（RiskMetrics）"""
    T, n = ret.shape
    w = np.array([(1-lam)*lam**(T-1-t) for t in range(T)])
    w /= w.sum()
    centered = ret - ret.mean(axis=0)
    return sum(w[t] * np.outer(centered[t], centered[t]) for t in range(T))
 
ewma = ewma_cov(returns, 0.94)
 
# 对比条件数（越小越好）
print(f"\n{'方法':<18} {'条件数':<12} {'最大特征值':<12} {'前5占总方差'}")
print("-" * 55)
for name, cov_mat in [("样本", sample_cov), ("LW", lw_cov), ("EWMA", ewma)]:
    ev = np.sort(np.linalg.eigvalsh(cov_mat))[::-1]
    cond = ev[-1] / ev[0]
    top5 = ev[:5].sum() / ev.sum()
    print(f"{name:<18} {cond:<12.2f} {ev[0]:<12.6f} {top5:.2%}")

衰减因子 $\lambda$	等效窗口	用途
0.94	~19 天	RiskMetrics 默认
0.97	~39 天	中期平滑
0.99	~100 天	长期稳定

---

五、凸优化与 KKT 条件

5.1 凸优化

凸函数的海森矩阵半正定，全局最优解唯一且高效可求。

5.2 KKT 条件

一般约束优化 $\min f(x)$, s.t. $g(x)\le 0,h(x)=0$ 的拉格朗日函数：

$\mathcal{L}=f(x)+\sum {\lambda }_i{g}_i(x)+\sum {\nu }_j{h}_j(x)$

条件	公式
平稳性	$\nabla f+\sum {\lambda }_i\nabla g_i+\sum {\nu }_j\nabla h_j=0$
原问题可行	$g_i\le 0,h_j=0$
对偶可行	${\lambda }_i\ge 0$
互补松弛	${\lambda }_i{g}_i=0$（约束松弛则乘子为零）

from scipy.optimize import minimize
 
np.random.seed(42)
n = 5
mu = np.array([0.08, 0.10, 0.06, 0.12, 0.07])
sigma = np.array([
    [0.040, 0.006, 0.002, 0.008, 0.003],
    [0.006, 0.090, 0.005, 0.010, 0.004],
    [0.002, 0.005, 0.030, 0.003, 0.002],
    [0.008, 0.010, 0.003, 0.100, 0.006],
    [0.003, 0.004, 0.002, 0.006, 0.050]
])
 
# 最小方差: min w^T Σ w, s.t. Σw=1, w>=0
w0 = np.ones(n) / n
result = minimize(
    fun=lambda w: w @ sigma @ w,
    x0=w0, method='SLSQP',
    bounds=[(0,1)]*n,
    constraints={'type': 'eq', 'fun': lambda w: w.sum() - 1}
)
print(f"最小方差组合波动率: {np.sqrt(result.fun):.4f}")
print(f"收益: {result.x @ mu:.4f}, 权重: {np.round(result.x, 3)}")

---

六、二次规划：均值方差优化

6.1 马科维茨模型

${\min }_{\mathbf{w}}\frac{1}{2}{\mathbf{w}}^T\Sigma \mathbf{w},\text{s.t.}{\mathbf{w}}^T\mathbf{\mu }\ge r_{\text{target}},1^T\mathbf{w}=1,\mathbf{w}\ge 0$

6.2 cvxpy 实现

import numpy as np
import cvxpy as cp
 
np.random.seed(42)
n = 10
mu = np.random.uniform(0.04, 0.15, n)
base = np.random.normal(0, 1, (n, n))
sigma = 0.01 * np.eye(n) + 0.001 * (base @ base.T)
 
# 1. 最小方差
w = cp.Variable(n)
prob = cp.Problem(cp.Minimize(cp.quad_form(w, sigma)),
                  [cp.sum(w) == 1, w >= 0])
prob.solve()
print(f"最小方差波动率: {np.sqrt(prob.value):.4f}")
 
# 2. 有效前沿
print(f"\n{'目标收益':<12} {'波动率':<10} {'夏普'}")
print("-" * 35)
targets = np.linspace(mu.min(), mu.max(), 20)
fvol, fret = [], []
for r in targets:
    w = cp.Variable(n)
    prob = cp.Problem(cp.Minimize(cp.quad_form(w, sigma)),
                      [cp.sum(w) == 1, w >= 0, mu @ w >= r])
    prob.solve(solver=cp.OSQP, warm_start=True)
    if prob.status == 'optimal':
        vol = np.sqrt(prob.value)
        fvol.append(vol); fret.append(w.value @ mu)
 
for i in range(0, len(fret), 4):
    sharpe = (fret[i] - 0.02) / fvol[i]
    print(f"{fret[i]:<12.4f} {fvol[i]:<10.4f} {sharpe:.4f}")
 
# 3. 解析解（允许做空时）
rf = 0.02
sigma_inv = np.linalg.inv(sigma)
w_anal = sigma_inv @ (mu - rf) / (np.ones(n) @ sigma_inv @ (mu - rf))
print(f"\n解析解（无约束）夏普: {((w_anal @ mu - rf) / np.sqrt(w_anal @ sigma @ w_anal)):.4f}")

---

七、拉格朗日乘子法

7.1 等式约束

$\min f(x)$, s.t. $h(x)=0\Rightarrow \mathcal{L}(x,\nu )=f(x)+\nu h(x)$，一阶条件 ${\nabla }_x\mathcal{L}=0$。

7.2 最小方差组合解析解

$\min \frac{1}{2}{\mathbf{w}}^T\Sigma \mathbf{w},\text{s.t.}{1}^T\mathbf{w}=1$

解：${\mathbf{w}}^{\ast }=\frac{{\Sigma }^{-1}1}{1^T{\Sigma }^{-1}1}$

import numpy as np
 
np.random.seed(42)
sigma_s = np.array([[0.040, 0.006, 0.002],
                     [0.006, 0.090, 0.005],
                     [0.002, 0.005, 0.030]])
ones = np.ones(3)
 
# 拉格朗日解析解
sigma_inv = np.linalg.inv(sigma_s)
w_lag = sigma_inv @ ones / (ones @ sigma_inv @ ones)
print("拉格朗日解: 最小方差组合")
print(f"  权重: {np.round(w_lag, 4)}, 之和: {w_lag.sum():.6f}")
print(f"  波动率: {np.sqrt(w_lag @ sigma_s @ w_lag):.4f}")

---

八、L1/L2 正则化

8.1 优化视角

正则化	惩罚项	特点
L2 (Ridge)	$\lambda \Vert \mathbf{w}{\Vert }_2^2$	权重缩小但不为零
L1 (Lasso)	$\lambda \Vert \mathbf{w}{\Vert }_1$	产生稀疏解
Elastic Net	混合	平衡稀疏和稳定

L2 约束（圆形）→ 最优解不易在坐标轴上 → 权重缩小
L1 约束（菱形）→ 最优解更可能在坐标轴上 → 产生稀疏解

8.2 正则化路径

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import Ridge, Lasso
from sklearn.preprocessing import StandardScaler
 
np.random.seed(42)
n_samples, n_features = 500, 20
true_w = np.zeros(n_features)
true_w[:3] = [0.5, 0.3, -0.4]
 
X = np.random.normal(0, 1, (n_samples, n_features))
y = X @ true_w + np.random.normal(0, 0.5, n_samples)
X_s = StandardScaler().fit_transform(X)
 
alphas = np.logspace(-4, 2, 50)
ridge_coefs = np.array([Ridge(alpha=a).fit(X_s, y).coef_ for a in alphas])
lasso_coefs = np.array([Lasso(alpha=a, max_iter=10000).fit(X_s, y).coef_ for a in alphas])
 
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
axes[0].plot(np.log10(alphas), ridge_coefs)
axes[0].axhline(0, color='black', lw=0.5)
axes[0].set_title('Ridge (L2) 路径')
 
axes[1].plot(np.log10(alphas), lasso_coefs)
axes[1].axhline(0, color='black', lw=0.5)
axes[1].set_title('Lasso (L1) 路径')
plt.tight_layout()
plt.show()
 
# 稀疏性对比
n_zero_l = np.sum(np.abs(lasso_coefs[15]) < 1e-3)
n_zero_r = np.sum(np.abs(ridge_coefs[15]) < 1e-3)
print(f"alpha 较大时零系数: Lasso={n_zero_l}/20, Ridge={n_zero_r}/20")

场景	推荐	理由
共线性强	Ridge	缩小共线因子权重
有效因子少	Lasso	自动选择
两者兼有	Elastic Net	平衡

---

九、正交投影与最小二乘

9.1 几何意义

线性回归 = 将 $\mathbf{y}$ 正交投影到 $\mathbf{X}$ 张成的空间：$\hat{\mathbf{y}}=\mathbf{Py}$，其中投影矩阵 $\mathbf{P}=\mathbf{X}({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T$。

残差 $\mathbf{e}=(\mathbf{I}-\mathbf{P})\mathbf{y}$ 与 $\hat{\mathbf{y}}$ 正交。

9.2 手写 OLS（解析、SVD、梯度下降三种方式）

import numpy as np
 
np.random.seed(42)
n_days, n_factors = 500, 5
X = np.random.normal(0, 1, (n_days, n_factors))
true_betas = np.array([0.5, 0.3, -0.2, 0.1, 0.4])
y = X @ true_betas + np.random.normal(0, 0.3, n_days)
 
# 方式 1: 解析解 (法方程)
X_aug = np.column_stack([np.ones(n_days), X])
w_analytical = np.linalg.solve(X_aug.T @ X_aug, X_aug.T @ y)
 
# 方式 2: SVD（数值更稳定）
U, S, Vt = np.linalg.svd(X_aug, full_matrices=False)
w_svd = Vt.T @ np.diag(1.0/S) @ U.T @ y
 
# 方式 3: 梯度下降
w_gd = np.zeros(n_factors + 1)
for _ in range(5000):
    grad = (2/n_days) * X_aug.T @ (X_aug @ w_gd - y)
    w_gd -= 0.001 * grad
 
print(f"{'参数':<10} {'解析':<8} {'SVD':<8} {'GD':<8} {'真实'}")
print("-" * 48)
labels = ['截距'] + [f'b{i}' for i in range(n_factors)]
true_vals = [0.0] + list(true_betas)
for i, lab in enumerate(labels):
    print(f"{lab:<10} {w_analytical[i]:<8.4f} {w_svd[i]:<8.4f} "
          f"{w_gd[i]:<8.4f} {true_vals[i]:.4f}")
 
print(f"\n解析 vs SVD 差异: {np.max(np.abs(w_analytical-w_svd)):.2e}")
print(f"解析 vs GD  差异: {np.max(np.abs(w_analytical-w_gd)):.2e}")

9.3 Gram-Schmidt 因子正交化

import numpy as np
 
def gram_schmidt(X):
    """Gram-Schmidt 正交化消除因子共线性"""
    n, p = X.shape
    Q = np.zeros_like(X)
    for j in range(p):
        v = X[:, j].copy()
        for i in range(j):
            proj = np.dot(Q[:, i], v) / np.dot(Q[:, i], Q[:, i])
            v -= proj * Q[:, i]
        Q[:, j] = v / np.linalg.norm(v)
    return Q
 
np.random.seed(42)
X = np.random.normal(0, 1, (200, 5))
X[:, 1] += 0.7*X[:, 0]; X[:, 2] += 0.5*X[:, 0] + 0.3*X[:, 1]
 
X_orth = gram_schmidt(X)
corr_before = np.corrcoef(X, rowvar=False)
corr_after = np.corrcoef(X_orth, rowvar=False)
 
print(f"正交化前最大非对角线: {np.max(np.abs(corr_before - np.eye(5))):.3f}")
print(f"正交化后最大非对角线: {np.max(np.abs(corr_after - np.eye(5))):.3f}")

---

总结

核心工具速查

工具	NumPy 函数	量化应用
矩阵乘法	@, np.dot	组合收益、因子暴露
矩阵求逆	np.linalg.inv, solve	最小二乘、马科维茨
特征分解	np.linalg.eigh	PCA、风险分解
SVD	np.linalg.svd	降维、去噪
二次规划	cvxpy.quad_form	均值方差优化
最小二乘	np.linalg.lstsq	因子回归

知识结构

线性代数                优化
├── 矩阵运算 → 组合计算    ├── 凸优化/KKT → 约束理解
├── 特征/SVD → PCA/去噪    ├── 二次规划 → 均值方差
├── PCA → 降维/正交化      ├── 拉格朗日 → 解析解
├── 协方差 → 风险建模      ├── 正则化 → 防过拟合
└── 投影 → 回归本质        └── 最小二乘 → 因子回归

常见陷阱

1. 协方差病态：大 n 小 T 时必须收缩或因子模型
2. 忽略约束：权重非负、行业约束影响很大
3. 非对称矩阵用 eigh：应使用 SVD
4. 正则化强度拍脑袋：需要交叉验证
5. 解析解不适用不等式约束：必须用数值优化器