稀疏主成分因子模型：用 L1 约束逼出可解释载荷 • halo 的技术博客

因子模型的第一个尴尬是：你算出来的因子看不懂。标准 PCA 在 30 个资产上提取 3 个主成分，得到 30×3=90 个载荷系数，每一个都是非零小数。于是「第一主成分」同时沾着银行、半导体、消费、能源——你没法跟投资经理说「这是周期因子」，因为它长得像一锅粥。一个更扎心的场景是风险归因：某天组合回了 3%，你想知道是哪类风险踩了，结果 PCA 因子把锅分摊给了全部 30 只持仓，你还是说不清到底是周期板块还是成长板块在作妖。这锅粥不是噪声，是 PCA 的旋转不变性在作祟：任何正交旋转后的载荷都解释同样多的方差，可解释的那一组旋转（块状、每因子只挂一组资产）恰好不是 PCA 默认给你的那一组。

稀疏主成分分析（Sparse PCA, SPCA）要做的就是把载荷里不重要的系数直接压成 0，让每个因子只挂少数资产。本文用 numpy 从零实现 Zou-Hastie（2006）的弹性网稀疏 PCA，所有图都是真实算出来的。

一、合成数据：块状稀疏因子模型#

为了让「可解释性」可证伪，我造一个本来就是块状的因子模型：3 个因子，每个只挂在 10 只连续的资产上（共 30 只），其余资产载荷严格为 0。加噪声后做标准 PCA，看它能不能还原这个块状结构。

import numpy as np
rng = np.random.default_rng(20260829)
p, K, T = 30, 3, 600          # 30 资产, 3 因子, 600 期
block = p // K               # 每因子 10 只
L_true = np.zeros((p, K))
for k in range(K):
    sl = slice(k*block, (k+1)*block)
    L_true[sl, k] = np.abs(rng.normal(0.7, 0.15, block)) + 0.25  # 只挂本块
F = rng.normal(0, 1, (T, K))                 # 因子收益
X = F @ L_true.T + 0.35*rng.normal(0, 1, (T, p))   # 资产收益 = 因子 + 噪声
X = X - X.mean(0)                            # 对资产去均值

python

标准 PCA（稠密）vs 稀疏 PCA（块状）载荷热力图

左图是标准 PCA 前 3 个载荷（红正蓝负）。肉眼可见：三个因子都糊满了全部 30 行，没有一个因子落在「干净的 10 行一块」上——这就是旋转不变性把真因子搅散了。右图是本文实现的稀疏 PCA，同样的 3 个因子，载荷被 L1 约束逼回了 3 个清晰的横向色块，每个因子几乎只挂在 10 只一组的资产上，和生成数据时的 L_true 一一对应。

二、稀疏 PCA 怎么做到：把 PCA 写成带 L1 惩罚的回归#

Zou-Hastie 的关键洞察是把「找主成分」重写成一个带弹性网惩罚的回归问题。对单一稀疏主成分，目标函数变成

$\max_{u,v}\; u^\top X v \;-\; \tfrac12 u^\top u \;-\; \tfrac{c}{2}v^\top v \;-\; \lambda\|v\|_1$

对 $u$ 求极值得 $u = Xv$ （投影方向由载荷决定）；对 $v$ 求极值（带 L1 次梯度）得**软阈值（soft-threshold）**闭式解：

$v = \mathcal{S}\!\left(\frac{X^\top u}{c},\; \frac{\lambda}{c}\right),\qquad \mathcal{S}(z,t)=\mathrm{sign}(z)\max(|z|-t,0)$

于是算法就是「先投影、再软阈值」交替迭代，直到收敛。软阈值这一步正是把小载荷精确归零的地方——这就是「稀疏」的来源。直觉上，它像给每个载荷系数设了一条「显著性门槛」：低于 λ/c 的贡献被视为噪声直接砍掉，高于门槛的才保留，门槛越高留下的因子越「干净」。

def soft(x, lam):
    return np.sign(x) * np.maximum(np.abs(x) - lam, 0.0)

def spca_one(X, c, lam, max_iter=300):
    n, pp = X.shape
    v = rng.normal(pp); v /= np.linalg.norm(v)
    for _ in range(max_iter):
        u = X @ v
        u /= np.linalg.norm(u) + 1e-12          # 投影方向归一
        v = soft((X.T @ u) / c, lam / c)        # 软阈值 → 逼出稀疏
        v /= np.linalg.norm(v) + 1e-12
        # (收敛判断略)
    return u, v

def spca_block(X, K, c, lam):
    Xk, Vs = X.copy(), []
    for _ in range(K):
        u, v = spca_one(Xk, c, lam)
        Vs.append(v.copy())
        Xk = Xk - np.outer(u, u @ Xk)            # 投影掉已提取分量，再提下一个
    return np.column_stack(Vs)

python

两个超参要解释清楚：c 是弹性网的** ridge 系数**（惩罚载荷的 L2 范数，保证数值稳定，本文取 0.001）；λ 是 L1 稀疏强度，越大零越多。多因子用「提取一个 → 从数据里减掉它的投影 → 再提下一个」的顺序 deflation，天然保证各因子正交。

三、稀疏路径：λ 越大越稀疏，但别贪心#

把 λ 从 0.05 扫到 2.2，看两个量怎么变：非零载荷占比（稀疏度）、以及载荷对真因子的恢复度。

稀疏路径：λ 对非零占比与恢复度的影响

非零占比随 λ 单调下降（λ=0.9 时约 58%、λ≥1.0 时压到 33%，即每因子只挂约 10 只）；恢复度则是一个先升后平的曲线。这给了一个非常实用的工程结论：稀疏不是免费的，但代价很小——选 λ=0.9 时，90 个载荷压到 52 个（少了约 42%），恢复度只从标准 PCA 的 0.918 掉到 0.908。把曲线翻过来看更直观：恢复度从 0.918 掉到 0.908 只丢了 1 个百分点，却换回了一半以上的「可解释载荷」，这笔买卖在需要给人看的因子报告里非常划算。

恢复度随 λ 变化（含标准 PCA 参考线）

图里绿色点是恢复度峰值（λ≈1.04、强稀疏区，每个因子塌成单块、完美命中真因子，恢复度冲到 1.0），灰色虚线是本文选用的 λ=0.9。重要的是：在 λ=0.9 这个已经明显稀疏的点，恢复度（0.908）与标准 PCA（0.918）几乎贴在一起。换句话说，你几乎没牺牲「抓因子方向」的能力，却换来了「每个因子一眼看懂挂在哪组股票」的可解释性。

四、解释力几乎不丢：累计方差对比#

有人会担心「压掉那么多系数，解释的方差是不是塌了」。直接数：把前 k 个载荷对协方差矩阵的解释力（Rayleigh 商）累加。

C = X.T @ X / T
def explained(V, C):
    s = 0.0
    for j in range(V.shape[1]):
        v = V[:, j] / (np.linalg.norm(V[:, j]) + 1e-12)
        s += v @ C @ v
    return s / np.trace(C)

python

因子数	标准 PCA 累计解释	稀疏 PCA 累计解释
1	0.3131	0.3048
2	0.6149	0.6127
3	0.8976	0.8974

累计解释方差：PCA vs SPCA

头一个因子 SPCA 比 PCA 少解释约 2.6%（0.3131→0.3048），到第三个因子两者差距已经小到 0.0002——几乎所有的方差都被保留，只是重新分配到了稀疏的载荷上。这正是 SPCA 的卖点：可解释性几乎是「白送」的。

五、三个容易踩的坑（诚实版）#

恢复度不是 SPCA 的卖点，可解释性才是。 本文数据里标准 PCA 的恢复度（0.918）还略高于 SPCA（0.908）。如果你只追求「复现因子方向」，标准 PCA 就够了；SPCA 的价值是让人能看懂、能据此做板块归因和风控，不是更高的统计拟合。别在报告里吹反了。
软阈值的零是「数值零」不是「结构零」。

我们靠 |v| > 1e-6 判定非零。λ 太小（<0.4）时几乎没有系数被压零，SPCA 退化成 PCA；λ 太大（>1.5）时一个因子可能塌成只剩 1–2 个资产，失去对板块的代表性。λ 必须在「稀疏路径」曲线上挑，别拍脑袋。
deflation 顺序会引入偏差。 先提的因子占的投影多，后提的因子是在残差里找——所以第 3 个因子的恢复度天然比第 1 个脆弱。真实数据里如果因子方差悬殊，建议对 deflation 残差做一遍 PCA 复核，确认没把小因子吃掉。

六、怎么把 SPCA 接进真实因子流程#

合成数据里 λ 是给定的，真实里得自己选。最稳的办法不是看恢复度（你没有真因子），而是看样本外重构误差：把数据按时间切成训练 / 验证，用训练集提 SPCA 载荷，在验证集上算重构误差 $\|X_{val} - X_{val} V V^\top\|_F^2$ ，挑让验证误差最小的 λ。

def cv_lambda(X, lams, c=0.001, fold=0.7):
    n = int(len(X) * fold)
    Xtr, Xva = X[:n], X[n:]
    errs = []
    for lam in lams:
        V = spca_block(Xtr, K, c, lam)
        rec = Xva @ V @ V.T
        errs.append(np.mean((Xva - rec) ** 2))
    return lams[int(np.argmin(errs))], errs

best_lam, errs = cv_lambda(X, np.linspace(0.3, 1.5, 13))

python

拿到稀疏载荷后，真实工作流通常有三步：(1) 板块归因——看每个因子非零载荷落在哪些行业，直接给因子起名（「因子 1 = 周期有色板块」）；(2) 异常检测——某只股票突然在新因子里冒出大载荷，往往是它业务逻辑变了（重组、换主业）；(3) 降噪选股——只保留稀疏载荷非零的资产做后续多因子合成，等于先用 SPCA 做了一遍特征选择，把被个股特异噪声污染的维度扔掉。需要提醒：SPCA 载荷仍有符号翻转自由度，跨期拼接因子时记得用「与上一期载荷的符号相关性」统一方向，否则本期因子 1 和下期因子 1 可能正负相反。

七、结语#

稀疏 PCA 解决的是因子模型里最现实、也最被忽视的问题：算出来的因子得有人看得懂。它把 PCA 重写成带 L1 惩罚的弹性网回归，用软阈值把 90 个混在一起的载荷压成 52 个非零项，每个因子干净地落在一个板块里，而累计解释方差几乎不掉（3 因子 0.8974 vs 0.8976）。在 30 资产 / 3 因子 / 600 期的合成数据上，这是一个干净、可复现、可证伪的结果——真实市场里你拿到的不会是完美块状载荷，但「先提 SPCA 看哪些资产聚成一团、再决定因子命名」这套流程，比对着一锅稠密载荷猜因子含义要靠谱得多。