Empirical Bayes linear regression in high dimensions: Method of moments and sub-linear sample complexity

作者: Zhou Fan, Yandi Shen, Haoyu Wang, Yihong Wu
主题: 高维统计 / 随机矩阵
相关性: 7/10
链接: https://arxiv.org/abs/2608.16771

---

一、领域脉络与小综述

• 这个方向是什么：本子方向研究高维线性回归模型中的经验贝叶斯（Empirical Bayes, EB）先验估计。具体而言，在模型 y = Xβ + ε 中，回归系数 β 的 p 个分量被假设为独立同分布（i.i.d.）自一个未知的次高斯先验分布 g，目标是仅从观测数据 (X, y) 中估计出这个先验 g。这与经典的“序列模型”（y_i = θ_i + ε_i）形成对比：在序列模型中，观测值 y_i 是 i.i.d. 的，先验估计退化为一个标准的反卷积问题；而在线性模型中，设计矩阵 X 将 β 的分量耦合起来，使得观测值的边际分布不再是 g 与噪声的简单卷积，而是 g 与自身的卷积（自反卷积）。这使得问题在统计和计算上都更具挑战性。当前，该方向正从需要线性样本量 n = Ω(p) 的似然方法，向理论上最优的次线性样本复杂度 n = p^{1-o(1)} 迈进。

• 发展脉络（history）：

  1. 奠基工作：序列模型中的经验贝叶斯。Robbins [Rob51, Rob56] 开创了经验贝叶斯范式，其核心思想是从数据中估计先验，而非事先指定。在序列模型中，由于观测值 i.i.d.，先验估计等价于经典的反卷积问题，理论和计算都已成熟。Casella [Cas85]、Zhang [Zha03]、Efron [Efr12] 等对其进行了广泛研究。
  2. 线性模型中的早期尝试与似然方法。将 EB 推广到线性模型面临“自反卷积”的困难。早期工作 [NS86, GF00, YL05] 要么将先验限制在刚性参数族，要么计算成本高昂。近年来，大多数方法基于似然，其共同起点是边际似然的 Gibbs 变分表示，将最大似然转化为一个关于先验和后验的双变量优化问题。代表性工作包括：
     • EM 及其变体：起源于 [NH98]，通过交替更新先验和计算后验（通常用 MCMC 近似），在统计遗传学中广泛应用 [ZCS13, LJZS+19]。
     • 梯度流：Fan et al. [FGSW23, FKL+25a, FKL+25b] 分析了理想化的非参数最大似然估计（NPMLE），并提出一种双变量梯度流算法，同时更新先验并通过 Langevin 扩散从后验采样。作者指出，该工作证明了在随机次高斯设计下，当 n ≥ cp 时，理想化的 NPMLE 是一致的，但缺乏可证明的计算算法。
     • 变分推断：Mukherjee et al. [MSS23] 和 Lee & Deb [LD26] 研究了平均场变分近似下的 EB 估计，后者给出了 O(1/√p) 收敛率的充要条件。
  3. 矩方法（MoM）的复兴。矩方法历史悠久 [Pea94]，但在线性模型中的应用多限于参数化情形。例如，Henderson [Hen53] 和 Rao [Rao71a, Rao72] 的 MIVQUE 理论用于估计高斯先验的方差。在遗传学中，LD-score 回归 [BSLF+15] 和基于四阶矩的统计量 [OSH+19] 也被用于估计遗传力等参数。Wu & Yang [WY20a, WY20b] 的“去噪矩方法”（DMM）为从矩估计恢复分布提供了高效工具。
  4. 本文的位置：本文提出经验贝叶斯矩方法（EBMoM），首次在非参数先验设定下，为一般设计（包括一大类相关随机设计）提供了一个计算可行且统计最优（达到次线性样本复杂度 n = p^{1-o(1)}）的 EB 先验估计器。这填补了似然方法（需要 n = Ω(p)）与信息论下界之间的空白。

• 子线索聚类：

  • 似然方法：包括 EM 及其变体、梯度流、变分推断。共同点是计算后验或近似后验，理论分析困难，目前仅在线性样本量 n = Ω(p) 下有一致性保证。
  • 矩方法：包括经典的 MIVQUE、LD-score 回归、DMM 以及本文的 EBMoM。通常计算更高效，但传统上多用于参数模型或特定矩的估计。本文将其推广到非参数先验的完整估计。
  • 计算-统计权衡：虽然本文主要关注统计效率，但其 O(np²) 的计算复杂度（主要来自 Gram 矩阵）在 n ≪ p 时是可行的。这与似然方法（通常需要 MCMC 或复杂优化）形成对比。

• 这个方向在追问的核心问题（2-4 个）：

  1. 最优样本复杂度：对于非参数先验估计，需要多少样本（n 相对于 p）才能一致地估计出先验 g？本文回答了这个问题：n = p^{1-o(1)}。
  2. 计算可行性：能否设计一个计算高效的算法来达到这个最优样本复杂度？本文的 EBMoM 给出了肯定答案。
  3. 设计矩阵的普适性：这些结果对什么样的设计矩阵 X 成立？本文的条件（Assumption 2）相当温和，涵盖了高度相关的设计。
  4. 下游任务的影响：先验估计的误差如何影响后续的贝叶斯推断（如后验均值估计）？本文通过 AMP 实验初步探索了这一点，但理论上的“EB 遗憾”（regret）仍是开放问题。

• ⚠️ 作者的 framing（必须明确标注成"这是作者的说法"）：

  • 作者把缺口 frame 成什么：作者将缺口 frame 为“在一般设计下，是否存在一个计算高效的算法，能以最优的次线性样本复杂度一致地估计非参数先验？”他们声称，现有似然方法（如 [FGSW23]）虽然理论上暗示了次线性可能性，但缺乏可证明的算法，而 EBMoM 填补了这一空白。
  • 哪些竞争路线被他淡化或回避了：作者淡化了似然方法在实践中的广泛应用（如统计遗传学中的各种贝叶斯回归方法），并强调这些方法缺乏严格的次线性样本复杂度理论保证。他们回避了与这些方法在具体应用场景（如 GWAS）中的全面性能比较，仅与 EBflow [FGSW23] 进行了有限对比。
  • 什么明显该被引 / 该存在、却没出现在 intro 里？：这是一个值得研究者去查的问题。例如，是否存在其他基于矩的、但采用不同策略（如利用随机矩阵理论的谱方法）来估计先验的工作？或者，在“统计-计算权衡”领域，是否有工作从计算复杂度假说（如低度多项式障碍）的角度，为线性模型中的 EB 问题提供了下界？这些可能未被提及。

• 张力：未见明显对立引用。所有被引工作基本都承认线性模型中的 EB 问题比序列模型困难得多，且现有理论结果有限。本文与似然方法的主要“张力”在于：似然方法（如 [FGSW23]）的理论结果暗示了次线性样本复杂度的可能性，但本文声称是第一个提供可证明算法和匹配下界的。这更像是一种“填补空白”而非“推翻结论”。

二、最核心、最简单的例子 / 数学问题

第一步：把符号、模型、可观测数据交代清楚

• 符号：

  • y ∈ ℝⁿ：可观测的响应向量。
  • X ∈ ℝⁿˣᵖ：设计矩阵，可以是固定的或随机的。
  • β ∈ ℝᵖ：未知的回归系数向量。
  • ε ∈ ℝⁿ：噪声向量，ε ~ N(0, σ²Iₙ)，σ² 已知。
  • g：未知的先验分布，β₁, ..., βₚ ~ i.i.d. g。这是要估计的目标。
  • m_k(g) = E[β₁ᵏ]：先验 g 的 k 阶（原点）矩。
  • µ_k(g) = E[(β₁ - m₁)ᵏ]：先验 g 的 k 阶中心矩。
  • H = XᵀX ∈ ℝᵖˣᵖ：Gram 矩阵。
  • A_k = Tr(T⁽ᵏ⁾) = Σᵢⱼ Hᵏᵢⱼ：一个关键的标量，用于衡量估计 k 阶矩的信号强度。
  • T⁽ᵏ⁾ ∈ (ℝᵖ)⊗ᵏ：一个 k 阶张量，其元素为 T⁽ᵏ⁾_{s₁,...,sₖ} = Σⱼ H_{j,s₁} ... H_{j,sₖ}。
  • F_k(y)：一个基于 Hermite 多项式的统计量，其期望与 β 的 k 阶矩张量相关。
  • W₁(g, ĝ)：一维 Wasserstein-1 距离，用于衡量估计分布 ĝ 与真实分布 g 的差异。

• 模型：

  • 数据生成机制：y = Xβ + ε，其中 β ~ g⊗ᵖ（即 β 的分量独立同分布自 g），ε ~ N(0, σ²Iₙ)，且 β 与 ε 独立。
  • 统计模型：这是一个贝叶斯线性回归模型，但先验 g 是未知的、非参数的（仅假设为次高斯分布）。目标是从 (X, y) 中推断 g。
  • 已知量：设计矩阵 X，噪声方差 σ²。
  • 待估对象：先验分布 g（或其矩 m_k）。

• 可观测数据：

  • 可观测：设计矩阵 X 和响应向量 y。研究者拥有 (X, y) 这对数据。
  • 不可观测 / 潜在：回归系数 β 和噪声 ε 都是不可观测的潜在变量。β 的分布 g 是想要但观测不到的，只能通过 (X, y) 的联合分布来识别。

第二步：讲最小内核

本文的核心思路可以用一个最简特例来理解：假设先验 g 是零均值的（m₁=0），且我们只关心估计其方差 µ₂。

在这个特例下，模型为 y = Xβ + ε，β 分量 i.i.d. 自一个未知的零均值分布 g，方差为 µ₂。我们想估计 µ₂。

• 核心困难：在序列模型 (X=I) 中，yᵢ = βᵢ + εᵢ，yᵢ 是 i.i.d. 的，其方差为 µ₂ + σ²，因此 µ₂ 很容易估计。但在线性模型中，y 的协方差矩阵为 Var(y) = µ₂ XXᵀ + σ² Iₙ，y 的分量不再独立，其方差结构由 X 和 µ₂ 共同决定。

• 关键想法：利用一个二次型 yᵀ Q y 作为估计量。对于任意对称矩阵 Q，有： E[yᵀ Q y] = Tr(Q Var(y)) = µ₂ Tr(Q XXᵀ) + σ² Tr(Q)。 如果我们选择一个 Q 使得 Tr(Q XXᵀ) ≠ 0，就可以构造一个无偏估计量： µ̂₂ = (yᵀ Q y - σ² Tr(Q)) / Tr(Q XXᵀ)。

• 最小内核的体现：本文的方差估计量（公式 5）正是选择了 Q = XXᵀ（实际上 ∥Xᵀỹ∥² = ỹᵀ XXᵀ ỹ）。这个选择是自然的，因为它最大化信号（Tr(Q XXᵀ) = ∥XᵀX∥²_F）与噪声的某种比率。这个特例揭示了整篇论文的核心机制：通过精心构造的关于 y 的多项式统计量（这里是二次型），其期望可以写成关于先验矩的线性方程，从而可以“解出”该矩。对于高阶矩，这个线性方程变成了一个下三角方程组：k 阶矩的统计量的期望等于 A_k * µ_k 加上一个仅依赖于更低阶矩（µ₂, ..., µ_{k-1}）的项。因此，可以通过递归的方式，先用低阶矩的估计值“扣除”这些非线性项，然后解出 µ_k。这就是 EBMoM 的递归定义（公式 9）的数学本质。

三、这篇论文做了什么

• 三句话：

  1. 研究了什么问题：在高维线性回归 y = Xβ + ε 中，当回归系数 β 的分量独立同分布自一个未知的次高斯先验 g 时，如何从数据 (X, y) 中计算高效地估计出 g。
  2. 核心工具 / 方法：提出了经验贝叶斯矩方法（EBMoM），通过一个下三角的估计方程组递归地估计先验的中心矩，并利用去噪矩方法（DMM）或高斯混合模型（GMM）将矩估计转化为分布估计。
  3. 主要结论：在温和的设计条件下（包括一大类相关随机设计），EBMoM 在样本复杂度 n ≥ p^{1-o(1)} 时能一致地估计出先验 g（在 Wasserstein-1 距离下）。一个匹配的信息论下界表明，这个次线性样本复杂度对于非参数先验估计是最优的，改进了现有似然方法需要 n = Ω(p) 的结果。

• 关键设定与假设：

  • Assumption 1 (Prior)：β 的分量来自一个次高斯分布 g。这是一个标准的轻尾假设，保证了矩的存在和增长可控。
  • Assumption 2 (Design)：存在常数 c₀ > 0 和整数 k₀，使得：
    1. ∥X∥_op / σ ≥ c₀：保证模型有非退化的信噪比。
    2. A₁ = ∥X1∥² ≥ c₀ (n∧p) ∥X∥²_op 和 A₂ = ∥XᵀX∥²_F ≥ c₀ (n∧p) ∥X∥⁴_op：保证均值和方差的估计是可行的。
    3. 对于 3 ≤ k ≤ k₀，A_k = Σᵢⱼ Hᵏᵢⱼ ≥ c₀ᵏ ∥X∥^{2k}_op · p · (n/p ∧ 1)ᵏ：这是最关键的条件，它保证了估计 k 阶矩的信号强度 A_k 足够大。作者指出，这个条件非常温和，甚至允许像“复制正交设计”这样的极端相关设计。相比已有文献，这个条件比似然方法通常需要的 i.i.d. 随机设计假设要宽松得多。
  • Lemma 2：证明了对于一大类满足 Poincaré 不等式的随机设计（如高斯设计），Assumption 2 以高概率成立。这为理论提供了具体的设计实例。

• 主要结果：

  • Proposition 1 (均值和方差)：在 Assumption 2 的部分条件下，均值和方差的估计误差为 O(1/(n∧p))。这意味着只要 n 和 p 都趋于无穷，即使 n ≪ p，均值和方差也能被一致估计。
  • Theorem 3 (高阶矩的 MSE 界)：这是核心定理。它给出了 EBMoM 对 k 阶中心矩 µ_k 的截断均方误差（MSE）界： E[(µ̂_k - µ_k)² 1_{C_k}] ≤ Δ_k，其中 Δ_k = (Ck)^{5k²} * (1/(n∧p)) * (p/n ∨ 1)^{2k²}。
    • 直觉：这个界表明，当 n 相对于 p 很小时（n ≪ p），误差主要由 (p/n)^{2k²} 项主导，随着 k 增大而迅速恶化。为了估计 k 阶矩，需要 n ≥ p^{1 - 1/(1+2k²)}。因此，要估计的矩的阶数越高，所需的样本量就越接近线性。
    • 必要条件：Assumptions 1 和 2 必须成立。
    • 解决的技术难点：证明的核心在于处理高阶矩估计中，由设计矩阵耦合带来的复杂依赖关系。作者通过将估计量分解为“中心化模型下的 oracle 估计”和“由均值估计误差引起的扰动”两部分（Section 6.2），并分别用归纳法和精细的扰动分析来控制它们的误差。
  • Corollary 4 (先验估计的一致性)：将矩估计的误差转化为 Wasserstein-1 距离下的先验估计误差。结论是，如果 n ≥ p^{1-ε} 且 ε → 0，那么通过选择足够多的矩（k → ∞），可以使得 W₁(g, ĝ) → 0 以概率趋于 1。这正式确立了 n = p^{1-o(1)} 的充分性。
  • Theorem 6 (信息论下界)：对于一大类设计矩阵（满足条件 (30)），如果 n ≤ p^{1-δ}，则存在两个不同的 1-次高斯先验 g 和 g'，使得它们的观测分布 P_g(y) 和 P_{g'}(y) 的总变差距离小于 0.1。这意味着任何估计器都无法一致地区分它们，从而证明了 n = p^{1-o(1)} 的必要性。证明路线：利用一个巧妙的“数据加工不等式”（data processing inequality），将问题转化为对一维正态混合的区分，然后通过构造匹配前 L 个矩的先验对，并利用 Hermite 展开和 Berry-Esseen 类型的论证来界定额外的矩差异。

• 证明路线与技术技巧（理论型）：

  • 整体路线（以上界证明为例）：
    1. 中心化：先假设均值 m₁ 已知，将模型中心化，得到 oracle 估计量 µ̄_k。这一步简化了问题，因为中心化后 β 的均值为 0，矩张量的表达式更简洁。
    2. 归纳法证明 oracle 估计的 MSE：对 k 进行归纳。基础步骤 (k=2) 通过直接计算方差完成。归纳步骤假设对 2,...,k-1 阶矩的估计误差有界，然后利用 µ̄_k 的定义（公式 33）和 F_k(y) 的期望展开（Lemma 9），将 µ̄_k - µ_k 分解为：
       • 一个由 F_k(y) 的随机波动引起的方差项。
       • 一个由低阶矩估计误差 µ̄_d - µ_d 引起的偏差项。 方差项通过 Lemma 11 (Var(F_k) 的界) 和 Assumption 2 控制。偏差项通过归纳假设和 Lemma 10 (γ 系数的界) 控制。最终得到 E[(µ̄_k - µ_k)² 1_{E_{k-1}}] ≤ Δ_k。
    3. 扰动分析：将实际估计量 µ̂_k 与 oracle 估计量 µ̄_k 的差异归因于均值估计误差 m̂₁ - m₁。通过精细的代数展开和矩估计，证明 E[(µ̂_k - µ̄_k)² 1_{F_{k-1}}] ≤ Δ_k（Proposition 18）。
    4. 合并：通过三角不等式和事件包含关系，得到 µ̂_k 的最终 MSE 界。
  • 关键跳跃点：最吃功夫的部分是归纳步骤中偏差项的控制。µ̄_k 的定义中包含了低阶矩估计量的乘积 µ̄_{d₁} ... µ̄_{d_t}。当用归纳假设的 MSE 界来控制 E[(µ̄_{d₁} ... µ̄_{d_t} - µ_{d₁} ... µ_{d_t}) 1_{E_{k-1}}] 时，需要将乘积的差分解为一系列项的线性组合，每个项都包含一个 (µ̄_d - µ_d) 因子。然后利用 Cauchy-Schwarz 不等式和归纳假设，将误差累积起来。这个过程的系数控制（如公式 58）需要非常小心，最终导致了 MSE 界中 (Ck)^{5k²} 这样的因子。
  • 技术技巧点名：
    • Hermite 多项式展开：利用 Hermite 多项式 H_k 的性质（E[H_k(Z)] = µ^k for Z ~ N(µ,1)），将关于 y 的非线性统计量 F_k(y) 的期望与 β 的矩张量联系起来（公式 16）。
    • 张量代数与对角-非对角分解：将 k 阶矩张量 M⁽ᵏ⁾ 与 T⁽ᵏ⁾ 的内积，按照 β 分量索引的相等模式（即“对角”和“非对角”部分）展开，从而将 E[F_k(y)] 表示为 A_k µ_k 加上低阶矩的项（公式 17）。这是整个递归方法的基础。
    • 动态规划（Algorithm 1）：用于高效计算 γ 系数，将计算复杂度从 O(p^{k/2}) 降低到 O_k(p²)。这利用了组合恒等式和多项式展开的系数提取。
    • Poincaré 不等式：在 Lemma 2 中，用于证明随机设计下 A_k 的下界。Poincaré 不等式提供了函数方差的控制，从而可以证明 A_k 围绕其均值集中。
    • Le Cam 两点法：在 Theorem 6 的下界证明中，通过构造两个难以区分的先验，来证明样本复杂度的必要性。
    • 数据加工不等式：在下界证明中，将高维问题约化到一维，简化了分析。

• 真实例子与应用：

  • 数据 / 场景：使用模拟数据，包括独立高斯设计（X_{ij} ~ N(0, 1/p)）和一种高度相关的块状设计（X_{block, ρ=0.9}）。考虑了四种先验：Rademacher、稀疏、高斯、双峰。
  • 方法应用：
    1. 矩估计：首先用 EBMoM 估计先验的中心矩（对 Rademacher 和稀疏先验估计前 3 阶，对高斯和双峰先验估计前 4 阶）。
    2. 分布恢复：对离散先验（Rademacher, 稀疏），用去噪矩方法（DMM）将矩估计转化为离散分布。对连续先验（高斯, 双峰），用 Lindsay 算法拟合 2-成分同方差高斯混合模型，或用非参数样条拟合。
  • 结果：
    • Wasserstein-1 误差：在 n/p ∈ {0.5, 1, 2} 下，随着 p 增大，所有先验的估计误差都呈下降趋势，验证了理论。相关设计下的误差略高于独立设计。
    • 与 EBflow 对比：EBMoM 在多数设定下与 EBflow（一种似然方法）性能相当，甚至更好（尤其对稀疏和双峰先验，EBflow 可能陷入局部最优）。用 EBMoM 初始化 EBflow 能显著提升后者的性能。
    • 下游任务：后验均值估计：将 EBMoM 估计的先验用于 AMP 算法进行系数估计。结果显示，EBMoM 驱动的 AMP 在 MSE 上接近使用真实先验的 AMP，并优于 Ridge 和 Lasso 等经典方法。
  • 这个例子想说明什么：实验旨在验证 EBMoM 的有限样本性能，证明其理论优势在实践中是可实现的，并且其估计的先验足以支持高质量的下游推断（如系数估计）。与 EBflow 的对比则突出了 EBMoM 作为计算高效且性能稳健的替代方案的价值。

• 🔎 结论是否比证明窄：

  • Theorem 3 的 MSE 界：定理的结论是在一个高概率事件 C_k 上的截断 MSE。虽然这个事件的概率很高（1 - Δ_k），但严格来说，结论不是无条件的。作者通过 Markov 不等式将其转化为高概率误差界，这是理论分析中的标准做法。
  • Corollary 4 的 o(1) 因子：结论说 n ≥ p^{1-o(1)} 是充分的，其中 o(1) 可以任意慢地趋于 0。这对应于需要估计的矩的阶数 k 趋于无穷。证明中 k 的选取依赖于 ε（n = p^{1-ε}），且 k 的增长速度受限于 1/√ε。因此，对于任意固定的 ε > 0，只需要有限个矩就能达到一致性，但 ε 越小，需要的矩越多。这个“任意慢”的 o(1) 是理论上的最优，但实践中可能意味着需要非常多的矩。
  • 对相关随机设计的假设：Lemma 2 和 Lemma 7 证明了 Assumption 2 和条件 (30) 对一大类随机设计成立。但论文并未穷尽所有可能的设计。作者在 Section 2.2 中用一个“复制正交设计”的例子说明 Assumption 2 非常温和，但并未证明它对所有“合理”的设计都成立。因此，结论的普适性依赖于这些假设的广泛性，而这一点被作者通过例子和引理有力地论证了，但并未完全证明。

四、开放问题（点到为止，扎根具体语句）

1. 未知噪声方差 σ²：论文假设 σ² 已知。作者在 Section 5 中讨论，σ² 可以在更低的样本复杂度（n ≫ √p）下被一致估计，但并未给出与 EBMoM 联合使用的完整理论。扎根于：Section 5 "Unknown noise variance" 段落。一个开放问题是：当 σ² 未知时，EBMoM 的样本复杂度是否会退化？能否设计一个联合估计程序，同时达到最优的样本复杂度？

2. 参数先验的样本复杂度：论文关注非参数先验。对于由有限个矩（如 k 个）识别的参数族，作者推测最优样本复杂度可能是 n ≫ p^{1-1/k}，并指出 EBMoM 的当前分析（n ≫ p^{1-1/(1+2k²)}）和基于 CLT 的下界（n ≫ p^{1-2/k}）都不紧。扎根于：Section 5 "Sample complexity of estimating parametric priors" 段落。这是一个明确的理论缺口，且与研究者对 minimax 界的兴趣高度相关。

3. EB 遗憾（Regret）：论文区分了“先验估计”和“下游系数估计”。对于后者，一个核心概念是 EB 遗憾，即数据驱动估计器与最优贝叶斯估计器之间的风险差。作者指出，在序列模型中，遗憾可以是 poly-log(p) 量级，但在线性模型中，即使 NPMLE 理论上能达到 o(p) 的遗憾，其计算实现仍是挑战。扎根于：Section 5 "Regret in EB regression" 段落。开放问题是：能否设计一个计算高效的算法（如基于 EBMoM 的），在 n = Θ(p) 时达到次线性甚至多对数级别的 EB 遗憾？

4. EBMoM 的方差最优性：作者在 Section 2.1 和 Appendix A 中讨论了存在方差更小的无偏矩估计量（如 MIVQUE 的推广），但为了计算简便选择了当前形式。扎根于：Section 2.1 "Unbiased polynomial estimators" 段落和 Appendix A。一个开放问题是：对于高阶矩，能否设计出计算上可行且方差达到最优（或接近最优）的估计量？这可能需要解决一个高维张量优化问题，与研究者对统计-计算权衡的兴趣可能有关。

---

Maintained by 陈星宇 · Homepage · Source on GitHub

评论