Archimedean Copula Inference via Taylor-Mode AD¶

会议: NeurIPS 2026
arXiv: 2605.23134
OpenReview: https://openreview.net/forum?id=JHPRlItZ0c
代码: https://github.com/thisiscam/acopula
领域: 优化/理论
关键词: 阿基米德Copula, 嵌套Copula, Taylor模式自动微分, 逐变量删失, 贝尔多项式

一句话总结¶

提出了首个原生基于 JAX 的嵌套阿基米德 Copula 推理框架 acopula，通过 Taylor 模式自动微分的高阶泰勒系数多项式乘方卷积，消除了对各族解析贝尔多项式表的依赖，实现了在任意逐变量右删失模式与嵌套树拓扑下的多项式时间精确似然计算及自动参数梯度。

研究背景与动机¶

在现代生存分析、生物医学多组学和量化金融风险建模中，带有层级嵌套结构和不完全观测的多变量数据极为普遍。例如在重症监护室（ICU）场景下，患者各项生理生化实验室指标自然按临床谱系（如全血细胞、电解质、肝功能、凝血等）形成分组依赖，且因出院或转科导致每位患者每项生化指标呈现异构的逐变量独立右删失；在金融市场中，股票资产按行业板块聚类，在极端市场崩塌时表现出明显的尾部非对称联动。嵌套阿基米德 Copula（Nested Archimedean Copula, NAC）通过将生成元函数组织为层级树结构，能用高度可解释的参数精准刻画组内与组间的非对称相关性及尾部依赖，是处理此类层级相依数据的核心工具。

然而，现有的嵌套阿基米德 Copula 拟合工具在实际落地中面临难以逾越的三重推理壁垒。首先，高维逐变量删失难以统一表征：在包含数十个变量的实际场景中，每个样本的观测掩码各不相同，对应 Copula 累积分布函数的不同混合偏导数，传统基于符号微分或全排列角点容斥展开的方法随维度呈指数爆炸，导致前序工具大多局限在二元、低维（$d \le 10$）或单事件共享时间删失；其次，生成元族与嵌套拓扑严重受限，主流工具（如 R 语言的 copula 包）仅预先硬编码了 5 类经典生成元的手推斯特林数表，且仅支持两层嵌套，根本无法推广到新型神经生成元（如 Deep Archimedean Copula）或更深层次的任意树形；最后，缺乏精确参数梯度，高维拟合只能依赖无导数优化、有限差分或蒙特卡洛随机似然逼近（SMLE），其估计方差和有限样本偏差随维度剧增。

本文的核心 idea 是利用 Taylor 模式自动微分（Taylor-Mode AD）提取边生成元复合函数的截断泰勒级数，借助多项式幂卷积直接恢复 Faà di Bruno 高阶链式法则所需的偏贝尔多项式，将嵌套似然计算转化为统一的单阶段可微多项式代数运算，从而在多项式时间内精确求解任意删失掩码下的混合偏导数与参数梯度。

方法详解¶

整体框架¶

acopula 的核心机制是将嵌套 Copula 的概率密度与混合偏导数计算映射为树形自底向上的系数向量传播。对于包含 $d$ 个叶节点的嵌套树，模型首先自底向上计算各内部节点的逆生成元累加自变量；随后从叶节点向上递归：每个叶节点根据其删失状态初始化贡献多项式，各内部节点对其子树传入的多项式应用由生成元复合函数高阶导数构造的“贝尔变换”（Bell Transform），并通过多项式乘积（柯西乘积）聚合兄弟子节点；最后在根节点处结合根生成元的高阶泰勒系数向量完成整树似然装配。在整个前向传播流中，所有多项式系数运算与中间张量均在对数空间与有符号对数和指数（Signed Log-Sum-Exp）表征下执行，保证了极端高阶导数下的数值稳定性。

%%{init: {'flowchart': {'rankSpacing': 24, 'nodeSpacing': 28, 'padding': 6, 'wrappingWidth': 400}}}%%
flowchart TD
    A["观测样本 u 与逐变量删失掩码 δ"] --> B["自底向上参数求值<br/>计算各节点聚合自变量 t_v"]
    B --> C["叶节点状态初始化<br/>未删失 [0,1] vs 删失单位元 [1]"]
    C --> D["隐式泰勒方程求解<br/>获取边复合函数 h_vc 泰勒系数"]
    D --> E["多项式幂卷积与贝尔变换<br/>迭代卷积 p^k 提取偏贝尔多项式"]
    E --> F["兄弟节点柯西乘积<br/>多项式截断卷积聚合子树信息"]
    F --> G["根节点装配与对数似然输出<br/>根生成元 Jet 结合与端到端梯度"]

关键设计¶

1. 基于泰勒模式自动微分的多项式幂卷积贝尔变换

嵌套阿基米德 Copula 密度的核心数学瓶颈在于求解多层复合函数的高阶链式法则（Faà di Bruno 公式），该公式本质由偏贝尔多项式 $B_{n,k}$ 决定。传统方法要么穷举集合划分，要么针对特定族手推解析表。acopula 发掘了偏贝尔多项式与指数生成函数之间的基本恒等式：若令 $p_{vc}(\varepsilon) = \sum_{m \ge 1} \frac{h_{vc}^{(m)}(t_c)}{m!} \varepsilon^m$ 为边复合函数 $h_{vc} = \psi_v^{-1} \circ \psi_c$ 在节点自变量 $t_c$ 处的截断泰勒多项式，则有 $$B_{n,k}\bigl(h'_{vc}(t_c), h''_{vc}(t_c), \dots\bigr) = \frac{n!}{k!} \cdot \operatorname{Coeff}\bigl(p_{vc}(\varepsilon)^k;\, \varepsilon^n\bigr)$$ 利用 JAX 的 Taylor 模式自动微分（Jet 原语），算法可在单次前向传播中以 $\mathcal{O}(d_c^2)$ 复杂度直接获得 $p_{vc}$ 的全部泰勒系数向量。随后，通过对截断多项式进行 $k$ 次迭代卷积 $q \leftarrow q * p_{vc}$ 累积多项式幂 $p_{vc}^k$，各阶偏贝尔多项式与子树系数 $\beta^{(c)}$ 的加权和直接简化为一个内积： $$\alpha_k^{(v,c)} = \frac{1}{k!} \sum_{j=k}^{|\delta_c|} (j! \beta_j^{(c)}) \cdot (q_k)_j$$ 该机制使得任意用户自定义的解析生成元或神经网络隐式生成元均可即插即用，彻底消除了每族手工推导的硬编码限制。

2. 统一逐变量删失掩码的自底向上多项式传播

在生存分析中，观测样本包含任意指示向量 $\delta \in \{0, 1\}^d$，其中 $\delta_j = 1$ 表示第 $j$ 个变量被精确观测，$\delta_j = 0$ 表示该变量在观测时刻发生右删失。删失机制直接决定了需要对 Copula 联合分布函数求导的子集与阶数。acopula 将删失操作内化为卷积代数的多项式单位元：对于未删失叶节点 $\ell$，初始化贡献向量为 $\alpha^{(\ell)} = [0, 1]$（代表对该叶变量求 1 阶偏导，对应 $x^1$ 且系数为 1）；而对于删失叶节点，直接初始化为 $\alpha^{(\ell)} = [1]$（对应常数项多项式 1）。由于 $[1]$ 是多项式乘法（柯西乘积）的单位元，删失叶节点在与其兄弟节点进行多项式卷积时不会引入任何额外的求导阶数，同时通过静态形状填充和条件分支 lax.switch，确保所有样本的不同删失掩码共享单一 JIT 编译计算图，无需针对每种 $2^d$ 组合重新编译。此外，删失发生时有效求导阶数 $|\delta_v|$ 自底向上缩减，不仅不会引入额外开销，反而随删失率增加呈 $\mathcal{O}(|\delta_v|^3)$ 显著加速。

3. 避免下溢崩溃的隐式泰勒方程求解

在复合函数 $h_{vc}(t) = \psi_v^{-1}(\psi_c(t))$ 的求导过程中，朴素的 Taylor 模式自动微分直接对显式链式调用求导，但在高维或高相关性时，内层生成元 $\psi_c(t)$ 会在 float64 下迅速下溢为 0，导致外层逆生成元 $\psi_v^{-1}(0)$ 出现数值爆炸或返回 NaN。acopula 提出了隐式泰勒反解算法：利用恒等式 $\psi_v(h_{vc}(t)) = \psi_c(t)$，将 $h_{vc}(t_c + \varepsilon) = h_0 + Q(\varepsilon)$ 代入，利用两者的展开系数建立关于 $Q(\varepsilon)$ 未知系数 $q_k$ 的下三角非线性代数方程组： $$a_1 q_k + \sum_{m \ge 2} a_m [\varepsilon^k] Q(\varepsilon)^m = b_k$$ 其中 $a_m = \psi_v^{(m)}(h_0)/m!$，$b_m = \psi_c^{(m)}(t_c)/m!$ 分别仅从单变量生成元 $\psi_v$ 和 $\psi_c$ 的局部优良定义域中独立求出。算法按阶数 $k$ 递推解出 $q_k$，单次求解复杂度为严格的 $\mathcal{O}(d_c^3)$，完全消除了链式复合导致的数值截断与下溢。同时，对同族嵌套（如 Clayton/Clayton、Gumbel/Gumbel 等），框架支持注册代数简化闭式，将求导复杂度进一步降至 $\mathcal{O}(d_c^2)$。

损失函数 / 训练策略¶

模型拟合遵循两阶段半参数最大似然估计（Two-Stage Semi-Parametric MLE）：第一阶段通过非参数 Kaplan-Meier 估计量（或边缘经验分布函数）拟合边缘生存函数，将原始生存时间转换为伪 Copula 观测变量 $u_j^{(i)} = \hat{S}_j(t_j^{(i)})$；第二阶段通过最小化负对数联合似然函数拟合 Copula 生成元参数 $\boldsymbol{\theta}$： $$\mathcal{L}(\boldsymbol{\theta}) = - \sum_{i=1}^n \left( \log c_{\delta^{(i)}}(\mathbf{u}^{(i)}; \boldsymbol{\theta}) + \sum_{j:\,\delta_j^{(i)}=1} \log \hat{f}_j(t_j^{(i)}) \right)$$ 借助 JAX 原生自动微分 jax.grad 提取精确参数梯度，使用 L-BFGS-B 优化器进行端到端优化。为确保嵌套合法性（即各层复合函数满足完全单调导数，即 Bernstein 函数条件），同族嵌套采用 McNeil 增量再参数化约束 $\theta_{\text{inner}} \ge \theta_{\text{outer}}$；跨族或神经网络生成元则采用闭式约束或在边复合上利用泰勒系数符号交替条件进行无损合法性监控。

实验关键数据¶

主实验¶

论文在三个前序工具无法处理的真实高维场景展开了主实验：包含 85,229 例患者、53 项生化指标、存在高度异构逐变量删失的 MIMIC-IV ICU 队列；涵盖 98 只成分股、11 个行业板块的标普 500 日收益率层级相依结构；以及包含 10 个阿基米德生成元族（其中 5 个此前从未有删失实现）的糖尿病视网膜病变队列。

在 MIMIC-IV（$d=53$）数据集上，acopula 成功拟合了包含 7 个临床板块的两层嵌套经典 Copula 与嵌套神经 Copula（包含以 Softmax 凸组合指数级联参数化的 ACNet 和基于 MLP 潜变量拉普拉斯变换的 Gen-AC）。下表汇报了测试集上的负对数似然对比：

| 模型族 | 结构/参数量 $|\theta|$ | 训练集 $-\bar{\ell}_{\text{train}}$ | 测试集 $-\bar{\ell}_{\text{test}}$ | $\text{AIC}/n$ | |---|---|---|---|---| | Gumbel / Gumbel | 8 | 5.224 | 5.242 | 10.448 | | AMH / Clayton | 8 | 5.161 | 5.158 | 10.321 | | Frank / Frank | 8 | 5.066 | 5.080 | 10.132 | | Clayton / Clayton | 8 | 4.951 | 4.949 | 9.902 | | Gen-AC (HACSurv 配置) | 1,142 | 4.912 | 4.914 | 9.857 | | Gen-AC-Fourier | 260 | 4.828 | 4.828 | 9.663 | | Gen-AC-Fourier-$\int_{260}$ | 260 | 4.707 | 4.712 | 9.430 | | ACNet(4,4)/I4 (本文推荐) | 70 | 4.308 | 4.310 | 8.619 | | ACNet(8,8)/I4 | 130 | 4.306 | 4.307 | 8.616 | | ACNet(4,4,4)/I4 | 90 | 4.308 | 4.310 | 8.619 |

在标普 500（$d=98$）高维嵌套金融建模中，Nested Gumbel 取得了最佳经典 AIC（$-\bar{\ell} = -21.2$，$\text{AIC}/n = -42.5$），优于 Clayton（$-20.3$）与 Frank（$-20.6$）；而高参数量的 R-vine 虽然在无删失下 NLL 达 $-46.6$，但参数量高达 3,706 且在数学机制上无法扩展至高维逐变量删失。

消融实验与运行效率分析¶

论文针对计算开销、删失加速效应以及数值精度进行了深入的消融验证：

删失率对计算耗时的影响（Nested Clayton, $d=50, n=100$ 单核 CPU）：

删失比例	有效未删失维度 $d_{\text{eff}}$	运行时间 (ms)	加速比
0% (全观测)	50	32.0	$1.0\times$ (基准)
25%	38	11.8	$2.7\times$
50%	25	6.4	$5.0\times$
75%	13	2.4	$13.6\times$

多项式乘法方案消融（截断多项式幂扫描 $p^1, \dots, p^K$ 单线程 CPU）：在真实 Copula 常见的板块规模 $K \le 192$ 范围内，直接卷积（Direct Convolution）在运行速度和分量数值精度上全面碾压基于快速傅里叶变换的乘法（FFT）。在 $K=16$ 时，直接卷积耗时 0.020 ms，FFT 耗时 0.067 ms（直接法领先 $3.29\times$）；在 $K=64$ 时直接法领先 $5.47\times$（0.060 ms vs 0.331 ms）。更关键的是，FFT 仅满足范数意义的前向误差界，导致阶乘级动态范围下极小的低阶系数因浮点截断误差全部失真，在 $K \ge 64$ 时相对误差发散；而直接卷积完美保留了分量级精度。
与 R 语言 nacLL 的基准对比：在 $d=35$ 且每个内部节点包含 5 个叶节点的两层对称树上，acopula 的单次密度计算速度比 R 语言 copula::nacLL 快约 $650\times$；当维度达到 $d > 40$ 时，R 语言实现因组合数溢出和硬编码限制直接中止退出，而 acopula 可稳定扩展至 $d=8,000$ 且耗时完美贴合 $\mathcal{O}(d^2)$ 理论渐近线。

关键发现¶

删失反向赋能计算效率：在传统的容斥积分框架中，删失维度越多计算越慢；而 acopula 的代数设计将删失变量转化为单位元 $[1]$，直接裁剪了内部节点的求导多项式次数，使得 75% 删失率下的似然计算速度提升了 $13.6\times$。
神经生成元的表达力差异：在 MIMIC-IV 建模中，仅含 70 个参数的 ACNet 凭借凸组合的指数级原子混合，测试似然（4.310）显著优于参数量高达 1,142 的 Gen-AC（4.914）。这是因为标准小型 MLP 对潜变量的平滑先验限制了其捕捉高曲率极值相依的能力，引入傅里叶特征并采用高精度数值积分消除蒙特卡洛方差后，Gen-AC 性能方有明显回升。
医学机制的一致还原：在 MIMIC-IV 上拟合的 Kendall's $\tau$ 准确复现了医学领域的生理学耦合，如白细胞分类计数（$\tau = 0.72$）、凝血酶原时间与 INR（$\tau = 0.61$）以及血气分析平衡（$\tau = 0.60$），而板块间全局相关性仅为 $\tau_{\text{outer}} = 0.15$。

亮点与洞察¶

将复杂的组合学微分转化为纯粹的多项式代数：通过泰勒模式 AD 与偏贝尔多项式生成函数的对应关系，把困扰 Copula 领域数十年的高阶导数求解问题降维为多项式幂卷积，消除了针对每个生成元族手动编写解析导数或查表的历史包袱。
极致的数值稳定性设计：在底层的 JAX Jet 规则中深度结合了有符号对数域表征（LogSeries）与隐式方程反解技术，使得算法即使面对 Frank、Joe 等指数衰减型生成元在超高维度（$d=98$）下的极度欠拟合与浮点下溢，依然能够稳定输出精确似然和非 NaN 梯度。
兼顾通用性与特化加速的灵活架构：对外提供统一易用的 Python 装饰器（仅需定义单变量生成元即自动推导全局），对内为已知解析闭式的同族嵌套预留了自动分流直通路径，兼具学术研究的灵活性与工业落地的高性能。

局限与展望¶

区间删失的组合爆炸：当前框架原生支持任意逐变量右删失，但尚未解决区间删失（Interval Censoring）问题；处理区间删失需要对超矩形角点进行容斥展开，涉及 $2^m$ 次密度或 CDF 计算，在缺乏次指数级精确算法前难以直接扩展至高维。
树拓扑结构的先验依赖：嵌套树的层级结构（各变量从属于哪个板块、树的深度等）目前仍依赖领域专家知识先验或启发式聚类预先指定，无法像藤 Copula（Vine Copula）那样在拟合过程中端到端自适应学习最优树拓扑。
编译时常与显存考量：在首次针对极端复杂的高维树结构进行 JAX 静态编译时，冷启动耗时可达 30–60 秒，尽管有磁盘缓存机制，但在动态变结构建模时仍存在一定开销。

评分¶

新颖性: ⭐⭐⭐⭐⭐ [巧妙将 Taylor-Mode AD 与贝尔多项式生成函数结合，攻克了嵌套 Copula 领域长期以来的高阶导数推导壁垒。]
实验充分度: ⭐⭐⭐⭐⭐ [涵盖数学符号严格对照、模拟参数恢复、真实医学高维删失（MIMIC-IV）、金融高维分析（S&P 500）与 10 大生成元横向对比。]
写作质量: ⭐⭐⭐⭐⭐ [数学推导严密流畅，算法细节与工程数值稳定性机制阐述详尽，主文与附录结构紧凑。]
价值: ⭐⭐⭐⭐⭐ [开源了工业级 JAX 库 acopula，彻底解决了长期阻碍高维嵌套 Copula 在生存分析、可靠性工程与金融风险落地的计算瓶颈。]

删失比例	有效未删失维度 \(d_{\text{eff}}\)	运行时间 (ms)	加速比
0% (全观测)	50	32.0	\(1.0\times\) (基准)
25%	38	11.8	\(2.7\times\)
50%	25	6.4	\(5.0\times\)
75%	13	2.4	\(13.6\times\)