阅读与运行方式

依据 Wang、Gao、Harvey、Liu、Tao 的 Machine Learning Meets Markowitz，NBER Working Paper 34861，2026年2月用户提供版本。页码均指论文印刷页码。代码是独立教学实现，不是作者源代码；不声称复现论文收益。

本文件将所有核心 Python 定义嵌入 R Markdown，各代码块按顺序共享同一 Python 会话，无须另行寻找 models.py。每节先列公式，再说明对应的代码对象和变量。Python 保留 PyTorch、CVXPYLayers 求解栈；R Markdown 通过 reticulate 执行 Python，并不是将数值实现重写为 R。

默认 run_python: false，可直接生成包含全部公式与代码的阅读版。执行方法：

install.packages(c("rmarkdown", "reticulate"))
# 在全新的 R 会话中，指定装好依赖的 Python；必须在初始化 Python 前设置。
Sys.setenv(RETICULATE_PYTHON = "/absolute/path/to/venv/bin/python")
rmarkdown::render("E2E_Markowitz_formula_code.Rmd",
                  params = list(run_python = TRUE))

Python 依赖：numpy scipy torch cvxpy cvxpylayers clarabel diffcp。默认关闭执行使没有 Python 依赖的 R 环境也可阅读。实际执行还需要 R、rmarkdown、Pandoc 和 reticulate。本次已核查 Python 代码与数值测试，但当前环境没有 Rscript，未实际 Knit HTML。

统一符号与计算对象

所有收益和费率使用十进制单位。共同可行域为

\[ \mathcal C(c)=\{w\in\mathbb R^N: \mathbf1^\top w=1,\ 0\le w_i\le c\},\qquad Nc\ge1.\tag{D1} \]

重要：式(15)的正则是未平方的 L2 范数；式(21)是标准差而非方差。 开启 zeta>0 用于 MV、TC、SD 时属于本报告的数值扩展，不是这些原式本身。

一般 E2E 与反向传播

预测层：Linear 与 NN1–NN5

\[ \mu_{i,t}=g_\theta(x_{i,t}),\qquad h^{(j)}=\operatorname{Dropout}_{.05}\!\left[ \operatorname{BN}(\operatorname{ReLU}(W_jh^{(j-1)}+b_j))\right].\tag{D6} \]

输入首先使用 LayerNorm；隐藏层宽度依次为32、16、8、4、2。depth=0使用线性模型，depth=1..5分别为NN1–NN5。这是同一优化框架的不同预测架构。

下面是后续各代码块共享的导入和参数定义。

from dataclasses import dataclass
import copy
import numpy as np
import cvxpy as cp
import torch
from torch import nn
from cvxpylayers.torch import CvxpyLayer
torch.set_default_dtype(torch.float64)

@dataclass
class Config:
    kind: str = 'sd'  # mv, ew, vw, tc, sd
    eta: float = 1.0
    gamma: float = 0.002
    zeta: float = 1e-4  # lower-level smoothing strength
    regularizer: str = "l2"  # paper Eq15: norm, NOT squared norm; alternative: squared
    cap: float = 0.1

class Predictor(nn.Module):
    def __init__(self, features, depth=3):
        super().__init__()
        if depth not in range(6):
            raise ValueError('depth must be 0 (linear) or 1..5')
        if depth == 0:
            self.net = nn.Linear(features, 1)
        else:
            layers = [nn.LayerNorm(features)]
            last = features
            for width in [32,16,8,4,2][:depth]:
                layers += [nn.Linear(last,width), nn.ReLU(), nn.BatchNorm1d(width), nn.Dropout(.05)]
                last = width
            self.net = nn.Sequential(*layers,nn.Linear(last,1))
    def forward(self,x):
        return self.net(x).squeeze(-1)

均值—方差模型：解析解到 QP

def mv_equality_solution(mu, Sigma, eta, zeta=0.0):
    n = len(mu)
    H = 2 * eta * Sigma + 2 * zeta * torch.eye(n, dtype=mu.dtype)
    e = torch.ones_like(mu)
    a = torch.linalg.solve(H, mu)
    b = torch.linalg.solve(H, e)
    nu = (e @ a - 1) / (e @ b)
    return a - nu * b  # 式(D7)

等权模型：LP、原文范数与平方范数

市值加权模型：两阶段导数

第一阶段得到选股权重 \(a\)，第二阶段用市值 \(v\) 转换：

\[ s=v^\top a>0,\qquad q_i=\frac{v_i a_i}{s}.\tag{D18} \]

\[ \frac{\partial q}{\partial a}=\frac{\operatorname{diag}(v)}s -\frac{(v\odot a)v^\top}{s^2},\qquad \frac{\partial(-y^\top q)}{\partial a_j} =-\frac{v_j}{s}(y_j-q^\top y).\tag{D19} \]

代码对应： w = w * market_cap; w = w / w.sum()。PyTorch 自动将式(D19)与优化层导数相乘。式(16)字面使用硬LP，因此仅增加归一化不能恢复其梯度；训练时需显式启用选股平滑。

最终 \(q\) 可能超过第一阶段的上限 \(c\)，不能把 \(a_i\le c\) 解释为 \(q_i\le c\)。若约束最终权重，需重新设计可行域。

交易成本模型：LP 与时序梯度

令 \(\kappa=\gamma/2\)。论文式(20)的下层及上层分别为

\[ w_t^*=\arg\min_{w\in\mathcal C(c)} -\mu_t^\top w+\kappa\|w-p_t\|_1,\quad p_t=w_{t-1}^*,\tag{D20} \]

\[ L_{TC}=\frac1T\sum_t\left[-y_t^\top w_t+\kappa\|w_t-p_t\|_1\right].\tag{D21} \]

引入 \(u\ge0\)，成本项可写为 \(\kappa\mathbf1^\top u\)，约束 \(w-p\le u\)、\(p-w\le u\)。原问题为LP；加平方正则得到QP；加未平方L2正则得到SOCP。这些正则是式(20)的可选扩展。

代码对应： 下层 cfg.gamma/2 * cp.norm1(w-prev)；上层 cfg.gamma/2 * (w-prev).abs().sum()。费率口径遵循论文的 \(\gamma/2\)，实际单边收费定义若不同，需要先统一。

令 \(w_t=h_t(\mu_t,w_{t-1})\)，\(A_t=\partial h_t/\partial\mu_t\)，\(B_t=\partial h_t/\partial w_{t-1}\)，则

\[ D_t=\frac{\partial w_t}{\partial\theta} =A_t\frac{\partial\mu_t}{\partial\theta}+B_tD_{t-1}.\tag{D22} \]

对单期损失 \(\ell_t(w_t,w_{t-1})\)，总反向伴随满足

\[ a_t=\partial_1\ell_t+\partial_2\ell_{t+1}+B_{t+1}^\top a_{t+1},\quad \nabla_\theta L=\frac1T\sum_t \left(\frac{\partial\mu_t}{\partial\theta}\right)^\top A_t^\top a_t.\tag{D23} \]

最后一期没有未来项。式(D23)对未平均的损失和定义伴随，最终统一除以 \(T\)。代码用 torch.stack(losses).mean() 自动完成尺度处理。

代码对应： rollout内部每期 prev = w，不调用 detach()，以保留式(D22)的跨期梯度。按日打乱或逐日detach会破坏这条链。连续8期分块并detach是截断反传近似，不是完整梯度。

若考虑价格漂移，可改为

\[ p_t=\frac{w_{t-1}\odot(1+r_t)}{1+w_{t-1}^\top r_t}, \]

同时保留这个变换的梯度；\(r_t\)必须在决策时已经实现。当前代码遵循论文简化，无资产池变化与价格漂移处理。

风险偏好模型：Mean–SD 与 SOCP

论文式(21)为

\[ w_t^*=\arg\min_{w\in\mathcal C(.1)} -\mu_t^\top w+\kappa\|w-p_t\|_1+\eta\sqrt{w^\top\Sigma_tw}.\tag{D24} \]

使用 \(R_t^\top R_t=\Sigma_t\)，风险项为 \(\eta\|R_tw\|_2\)，引入辅助变量 \(s\ge\|R_tw\|_2\) 即可表示为SOCP。上层用真实收益替换预测收益，其余经济惩罚保持相同。

对 \(\sigma=\sqrt{w^\top\Sigma w}>0\)，

\[ \nabla_w\sigma=\frac{\Sigma w}{\sigma},\qquad \nabla_w^2\sigma=\frac\Sigma\sigma- \frac{(\Sigma w)(\Sigma w)^\top}{\sigma^3}.\tag{D25} \]

因而 \(\nabla_w\ell=-y+\kappa\operatorname{sign}(w-p)+\eta\Sigma w/\sigma\)，在L1折点使用次梯度或锥求导。

代码对应： 下层 cp.norm(R @ w, 2)，上层 torch.linalg.vector_norm(R @ w)。不能写成方差。Mean–SD的Hessian在径向方向具有零曲率，不能据此承诺强凸或普遍更稳定。相同 \(\eta\) 也不意味着Mean–SD与Mean–Var持仓相同；在同一候选最优点匹配一阶条件时有 \(\eta_{SD}=2\eta_{MV}\sigma(w)\)。

通用优化层：把上述公式集中为可执行代码

以下完整定义覆盖式(D10)–(D24)的下层问题：

用 \(R\) 而非直接把可变协方差放入 quad_form，使参数化问题满足DPP要求。求解器使用Clarabel。

class Portfolio(nn.Module):
    def __init__(self,n,cfg):
        super().__init__()
        if cfg.kind not in {'mv','ew','vw','tc','sd'}:
            raise ValueError('unknown model')
        if cfg.cap*n < 1-1e-12 or not 0 < cfg.cap <= 1:
            raise ValueError('infeasible concentration cap')
        if min(cfg.eta,cfg.gamma,cfg.zeta)<0:
            raise ValueError('penalties must be nonnegative')
        self.cfg = cfg
        w=cp.Variable(n)
        mu=cp.Parameter(n)
        prev=cp.Parameter(n)
        # R.T @ R = Sigma. Passing the factor makes the problem DPP.
        R=cp.Parameter((n,n))
        if cfg.regularizer not in {"l2", "squared"}:
            raise ValueError("regularizer must be l2 or squared")
        penalty = cp.norm(w,2) if cfg.regularizer=="l2" else cp.sum_squares(w)
        obj=-mu@w+cfg.zeta*penalty
        params=[mu]
        if cfg.kind in {'tc','sd'}:
            obj += cfg.gamma/2*cp.norm1(w-prev)
            params.append(prev)
        if cfg.kind in {'mv','sd'}:
            obj += cfg.eta*(cp.sum_squares(R@w) if cfg.kind=='mv' else cp.norm(R@w,2))
            params.append(R)
        problem=cp.Problem(cp.Minimize(obj),[w>=0,w<=cfg.cap,cp.sum(w)==1])
        assert problem.is_dpp()
        self.layer=CvxpyLayer(problem,parameters=params,variables=[w])
    def forward(self,mu,prev,R,market_cap=None):
        args=[mu]
        if self.cfg.kind in {'tc','sd'}: args.append(prev)
        if self.cfg.kind in {'mv','sd'}: args.append(R)
        (w,)=self.layer(*args,solver_args={'solve_method':'Clarabel', 'tol_gap_abs':1e-10, 'tol_gap_rel':1e-10, 'tol_feas':1e-10})
        if self.cfg.kind=='vw':
            if market_cap is None or bool((market_cap<=0).any()):
                raise ValueError('strictly positive market caps required')
            w=w*market_cap
            w=w/w.sum()
        return w

def covariance_factor(history, ridge=1e-8, diagonal=True):
    """Use only returns observable at decision time. Output R.T R = Sigma."""
    cov=np.atleast_2d(np.cov(np.asarray(history),rowvar=False,ddof=1))
    if diagonal: cov=np.diag(np.diag(cov))
    cov=(cov+cov.T)/2+ridge*np.eye(cov.shape[0])
    return torch.tensor(np.linalg.cholesky(cov).T)

上层效用与时序展开

以下两段分别对应上层经济损失和式(D22)的顺序展开。realized_loss不包含下层平滑正则。VW传入的w已经是最终市值权重。

def realized_loss(w,y,prev,R,cfg):
    loss=-(w*y).sum()
    if cfg.kind in {'tc','sd'}: loss=loss+cfg.gamma/2*(w-prev).abs().sum()
    if cfg.kind=='mv': loss=loss+cfg.eta*(R@w).square().sum()
    if cfg.kind=='sd': loss=loss+cfg.eta*torch.linalg.vector_norm(R@w)
    return loss

def rollout(model,layer,X,Y,R,w0,market_caps=None):
    """No detach inside a sequence: full BPTT through previous holdings."""
    prev=w0
    losses=[]
    weights=[]
    for t in range(len(X)):
        w=layer(model(X[t]),prev,R[t],None if market_caps is None else market_caps[t])
        losses.append(realized_loss(w,Y[t],prev,R[t],layer.cfg))
        weights.append(w)
        prev=w
    return torch.stack(losses).mean(),torch.stack(weights)

预训练、Adam 与验证早停

预训练目标为

\[ L_{MSE}(\theta)=\frac1T\sum_t\frac1N\|g_\theta(X_t)-y_t\|_2^2.\tag{D26} \]

然后以 \(L_{E2E}\) 替换训练目标，保留网络参数作为warm start。MSE只用于初始化，不在本实现的E2E损失中额外混入。

fit每轮训练都从指定初始持仓重新推进；更新参数后，用当前参数重新计算训练段末持仓，再推进验证段。验证标签只用于损失，不进入下层决策。记录最佳验证参数并恢复。

本实现全序列反传，未逐字实现论文8个横截面的梯度累积安排；大规模应用需选择完整反传或明确标注的截断近似。学习率采用论文的一般设置，另加入梯度裁剪用于稳定数值。原文验证集为9个月且按季度滚动重训，本演示仅做合成时间切分。

def pretrain(model,X,Y,epochs=10,lr=1e-3):
    model.train()
    opt=torch.optim.Adam(model.parameters(),lr=lr)
    for _ in range(epochs):
        opt.zero_grad()
        loss=torch.stack([(model(x)-y).square().mean() for x,y in zip(X,Y)]).mean()
        loss.backward(); opt.step()

def fit(model,layer,train,valid,w0,epochs=50,patience=10):
    """train/valid=(X,Y,R,caps). Full-sequence reference training, not 8-day TBPTT.
    Validation holdings are recomputed from training under current fixed parameters.
    """
    import copy
    opt=torch.optim.Adam(model.parameters(),lr=1e-4 if layer.cfg.eta<1 else 1e-3)
    best=float('inf'); state=None; stale=0; log=[]
    for epoch in range(epochs):
        model.train(); opt.zero_grad()
        loss,_=rollout(model,layer,*train[:3],w0,train[3])
        loss.backward(); torch.nn.utils.clip_grad_norm_(model.parameters(),10); opt.step()
        model.eval()
        with torch.no_grad():
            _,train_w=rollout(model,layer,*train[:3],w0,train[3])
            val,_=rollout(model,layer,*valid[:3],train_w[-1],valid[3])
        score=float(val); log.append([epoch,float(loss.detach()),score])
        if score<best:
            best=score; state=copy.deepcopy(model.state_dict()); stale=0
        else: stale+=1
        if stale>=patience: break
    model.load_state_dict(state)
    return log

可执行合成数据示例

下列示例按顺序运行：构造数据、预训练、E2E训练、验证。为了观察小型横截面的稳定自由集，示例使用较强的平方正则替代模型。它不是论文默认超参数。\(N=20,c=.1\)最多单只10%，并不表示前10%选股；若要前10%，应设 \(c=10/N=.5\)。

真实数据应提供 X:[T,N,P]、Y:[T,N]、R:[T,N,N]、market_caps:[T,N]和初始持仓w0:[N]。初始均匀持仓意味着未计入首次由现金建仓的费用。

torch.manual_seed(17)
T, N, P = 12, 20, 5
X = torch.randn(T, N, P)
Y = .004 * X[:, :, 0] + .01 * torch.randn(T, N)
R = torch.diag(torch.linspace(.015, .03, N)).repeat(T, 1, 1)
market_caps = torch.exp(torch.randn(T, N))
w0 = torch.ones(N) / N

training_summary = {}
for kind in ["mv", "ew", "vw", "tc", "sd"]:
    cfg = Config(kind=kind, eta=1, gamma=.002, zeta=.02,
                 regularizer="squared", cap=.1)
    model = Predictor(P, depth=0)
    layer = Portfolio(N, cfg)
    pretrain(model, X[:8], Y[:8], epochs=2)
    train = (X[:8], Y[:8], R[:8], market_caps[:8])
    valid = (X[8:], Y[8:], R[8:], market_caps[8:])
    log = fit(model, layer, train, valid, w0, epochs=3)
    model.eval()
    with torch.no_grad():
        loss, W = rollout(model, layer, X, Y, R, w0, market_caps)
    assert torch.isfinite(W).all()
    assert torch.max(abs(W.sum(1) - 1)) < 1e-5
    assert W.min() > -1e-5
    if kind != "vw":
        assert W.max() < cfg.cap + 1e-5
    training_summary[kind] = {"loss": float(loss), "epochs": len(log)}
print(training_summary)

解析公式与自动微分的数值核查

norm_solution_errors = []
for d in [-.15, -.07, 0, .07, .15]:
    w_cp = cp.Variable(2)
    prob = cp.Problem(cp.Minimize(-np.array([d, 0]) @ w_cp + .1 * cp.norm(w_cp, 2)),
                      [w_cp >= 0, cp.sum(w_cp) == 1])
    prob.solve(solver="CLARABEL", tol_gap_abs=1e-10, tol_feas=1e-10)
    exact = (0 if d <= -.1 else 1 if d >= .1
             else .5 + d / (2 * np.sqrt(.02 - d*d)))
    norm_solution_errors.append(abs(float(w_cp.value[0]) - exact))

norm_layer = Portfolio(2, Config(kind="ew", cap=1, zeta=.1, regularizer="l2"))
mu = torch.tensor([.03, 0.], requires_grad=True)
w = norm_layer(mu, torch.ones(2)/2, torch.eye(2))
grad = torch.autograd.grad(w[0], mu)[0][0]
exact_grad = .01 / (.02 - .03**2)**1.5
assert max(norm_solution_errors) < 1e-4
assert abs(float(grad)-exact_grad) < 1e-3
print({"max_weight_error": max(norm_solution_errors),
       "derivative_error": abs(float(grad)-exact_grad)})

gradient_checks = {}
for kind in ["mv", "ew", "vw", "tc", "sd"]:
    layer_check = Portfolio(N, Config(kind=kind, eta=1, gamma=.002,
                            zeta=.02, cap=.1, regularizer="squared"))
    mu0 = torch.zeros(N, requires_grad=True)
    anchor = torch.zeros(N) if kind in {"tc", "sd"} else w0
    direction = torch.linspace(-1, 1, N)
    direction = direction / direction.norm()
    reward = torch.cos(torch.arange(N, dtype=torch.float64))
    def scalar(m):
        return (layer_check(m, anchor, R[0], market_caps[0]) * reward).sum()
    automatic = torch.autograd.grad(scalar(mu0), mu0)[0] @ direction
    h = 1e-3
    finite_diff = (scalar(mu0.detach()+h*direction)
                   - scalar(mu0.detach()-h*direction)) / (2*h)
    error = float(abs(automatic-finite_diff)/(1+abs(finite_diff)))
    assert error < 2e-3, (kind, error)
    gradient_checks[kind] = error
print(gradient_checks)

公式—代码定位总表

复现边界与扩展

1. MSE-Opt：先独立训练MSE预测器，冻结参数，再调用同一个优化层；不要反向传播经济损失。MSE-Top：冻结后用硬top-k；硬排序仅供评价，不能当成可训练梯度层。PPP直接输出权重，和本文E2E并非同一结构。
2. 原文EW/VW/TC上限为 \(10/N\)，Mean–SD上限为 \(.1\)，MV示例仅非负与预算。代码参数必须按所分析的模型选择。
3. 平滑训练后若切换硬top-k评价，须额外报告这一切换的效果。VW归一化后不承诺保留初始上限。
4. 真实研究需隔离训练、验证、测试，避免未来协方差与标签泄漏，并处理证券退出、停牌、可成交性、特征缺失和持有期重叠。这里固定资产池，不额外承诺可实盘成交。
5. 后续若采用连续时间块截断反传，应明确标为近似，并验证截断长度的敏感性。完整时序计算图内不能随意detach。

参考定位

• Wang et al. (2026), Machine Learning Meets Markowitz, NBER Working Paper 34861。
• §3.2，pp.19–22：一般E2E、均值—方差示例及链式法则。
• §3.3，pp.23–25：KKT与可微凸优化。
• §3.4–3.6，pp.26–29：等权、市值加权、跨期依赖。
• §4.3，式(20)，p.36：交易成本E2E。
• §4.4，式(21)，pp.40–42：均值—标准差、风险偏好。
• 附录A.2，pp.70–74：网络结构、训练程序与求解环境；A.3，pp.74–75：LP梯度退化与L2范数示例。