面向非平滑随机复杂度与流形采样的精确Schur-Sylvester降维

# 非光滑随机复杂度和流形采样的精确 Schur-Sylvester 降维方法
来源: https://arxiv.org/html/2606.23867

###### 摘要

对于正则非光滑估计量（例如 Lasso），精确计算归一化最大似然（NML）编码长度历来受限于流形约束投影和体积积分中的立方缩放壁垒。在几何 Propose-and-Project Metropolis–Hastings (PPMH) 采样器的每一步中，评估投影算子需要求逆一个 \((N+k) \times (N+k)\) 的广义 KKT 矩阵，而计算体积因子需要求一个 \((N-k) \times (N-k)\) Gram 矩阵的行列式。本文提出了一种精确的、数学上等价的公式，通过利用块 Schur 补和 Sylvester 行列式恒等式规避了这两个瓶颈。我们证明，这两种操作的计算复杂度从 \(\mathcal{O}(N^3)\) 降低到每步 \(\mathcal{O}(k^3 + N^2 k)\)。我们将此降维方法推广到稀疏支持向量机 (SVM)、弹性网络和 Group Lasso。最后，我们提供了严格的数值稳定性分析，并使用每秒有效样本量 (ESS) 评估了采样器的效率。我们在高维数据集上的实证基准测试证实，在保持双精度数值等价性的前提下，持续加速超过 14,100 倍，使得大规模统计推断中的精确非光滑 NML 估计变得高度可行。

## 1 引言

最小描述长度 (MDL) 原则认为，最优统计模型是使模型和压缩数据的联合编码长度最小化的模型 [1, 2, 3]。理论最优值由归一化最大似然 (NML) 分布实现，其归一化因子被称为随机复杂度 [4, 5]:

\[ C(\mu_\lambda) = \int_{\mathcal{X}} p(\mathbf{x} \mid \hat{\boldsymbol{\theta}}_\lambda(\mathbf{x})) \, d\mathcal{L}^N(\mathbf{x}), \tag{1} \]

其中 \(\hat{\boldsymbol{\theta}}_\lambda(\mathbf{x})\) 是在正则化参数 \(\lambda\) 下的最大似然估计量 (MLE)，且 \(\mathcal{X} \subseteq \mathbb{R}^N\) 是连续数据空间。对于现代机器学习中普遍存在的非光滑估计量（如 Lasso、稀疏 SVM 和低秩奇异值约束），评估这个积分所需的经典光滑性假设不再成立。

在我们先前的工作 [6] 中，通过将余面积公式推广到正则路径可微 Lipschitz (PDL) 估计量，我们建立了非光滑 NML 的测度理论基础。这将对随机复杂度的公式描述为对估计量的非可微水平集（流形）的积分：

\[ C(\mu_\lambda) = \int_{\Theta} \left[ \int_{\hat{\boldsymbol{\theta}}_\lambda^{-1}(\boldsymbol{\theta}')} \frac{p(\mathbf{x} \mid \boldsymbol{\theta}')}{J_{\text{cons}}(\hat{\boldsymbol{\theta}}_\lambda(\mathbf{x}))} \, d\mathcal{H}^{N-k}(\mathbf{x}) \right] d\mathcal{L}^k(\boldsymbol{\theta}'), \tag{2} \]

其中 \(J_{\text{cons}}(\hat{\boldsymbol{\theta}}_\lambda(\mathbf{x}))\) 表示水平集上的约束雅可比矩阵。将余面积公式严格应用于样本分布理论在统计学中有着悠久历史，特别是在刻画随机变量的代数和 Lipschitz 变换的密度 [7] 以及定义矩阵流形上的非光滑概率分布 [8] 方面。

为了精确计算公式 (2) 中的水平集积分，我们引入了 Propose-and-Project Metropolis–Hastings (PDL-PPMH) 采样器，它直接遍历这些水平集。虽然可以使用运动学公式，通过将代数簇与随机线性空间相交来实现在代数流形上的全局独立同分布采样 [9]，但此类方法在结构上仅限于多项式方程，并且可能随着代数簇的次数增加而扩展性变差。相比之下，局部 MCMC 提议方案会沿切空间构造提议步，然后投影回目标流形 [10, 11, 12]，或者部署针对约束域设计的惩罚 Langevin 更新 [13]。然而，PPMH 采样器利用位置依赖的提议协方差来跟随流形的局部几何。如 Livingstone [14] 所分析的，位置依赖的提议协方差可能会改变 Metropolis–Hastings 链的各态历经性质，需要在尾部进行步长控制以防止接受率趋近于 1。在隐式神经流形上，通常使用迭代共轭梯度 (CG) 例程来处理这些投影步骤，以绕过显式雅可比矩阵的构建 [15]。然而，迭代近似会引入非平凡的误差传播。

此外，当在流形约束下进行采样时，从数学上讲，需要一个由曲率引起的确定性修正项，以确保提议严格保持在目标流形上，而不会引入系统性偏差或维度损失 [16, 12]。为了绕过对曲率信息或二阶导数的需求，可以采用特定的正交投影方法，这些方法中雅可比因子可以相互抵消，从而保持细致平衡 [12]。PDL-PPMH 算法在其投影和曲率修正循环中面临两个显著的计算瓶颈：

1. 1. **投影壁垒**：将提议点投影回水平集 \(\hat{\boldsymbol{\theta}}_\lambda^{-1}(\boldsymbol{\theta}')\) 需要求解一个约束 KKT 系统。对该映射求微分以计算 Radon-Nikodym 修正因子，需要直接求逆一个大小为 \((N+k) \times (N+k)\) 的广义 KKT 矩阵，每步开销为 \(\mathcal{O}((N+k)^3)\) 次运算。这个零空间投影问题类似于冗余机器人系统中遇到的问题，在该系统中，使用雅可比引导的零空间遍历来生成隐式流形表示 [17]。
2. 2. **体积因子壁垒**：计算水平集体积因子 \(J_{\text{cons}}(\hat{\boldsymbol{\theta}}_\lambda(\mathbf{x}))\) 需要构建大小为 \((N-k) \times (N-k)\) 的投影切基 Gram 矩阵并评估其行列式，每步开销为 \(\mathcal{O}((N-k)^3)\) 次运算。

此外，如果底层约束水平集分解为多个不连通的拓扑分量，则局部提议机制从根本上会被困住。解决这一拓扑障碍通常需要并行回火或副本交换框架来促进跨越不连通区域的转换 [18]。

为了克服上述计算瓶颈，本文提出了一种精确的、数学上等价的公式，将两个缩放壁垒都降低到每步 \(\mathcal{O}(k^3 + N^2 k)\) 次运算，从而将计算成本与环境数据维度完全解耦。在第 2 节中，我们首先介绍本工作的信息理论基础。然后，我们在第 3 节中描述我们提出的 Schur-Sylvester 降维公式，并展示如何将其用于 Lasso 及其他正则非光滑模型。在第 4 节中，我们进一步建立评估我们提出的 Schur-Sylvester 采样器误差传播和收敛分析的数学基础。在第 5 节中，我们给出数值实验结果以评估我们提出的方法。我们在第 6 节中总结工作并讨论未来可能的研究方向。

## 2 信息理论基础与非渐近遗憾

为了将本贡献置于信息论框架内，回想一下，NML 分布在最坏情况参数序列下达到了通用编码中的极小极大遗憾 [5]。对于数据序列 \(\mathbf{x}^N \in \mathcal{X}^N\)，相对于正则化最大似然模型的精确遗憾由下式给出：

\[ R(\mathbf{x}^N) = \ln p(\mathbf{x}^N \mid \hat{\boldsymbol{\theta}}_\lambda(\mathbf{x}^N)) - \ln C(\mu_\lambda). \tag{3} \]

经典的渐近 MDL 近似（Rissanen 公式）将随机复杂度简化为 [5]：

\[ \ln C(\mu_\lambda) \approx \frac{k}{2} \ln\left(\frac{N}{2\pi}\right) + \ln \int_{\Theta} \sqrt{\det\left(\mathbf{I}(\boldsymbol{\theta})\right)} \, d\boldsymbol{\theta} + o(1), \tag{4} \]

其中 \(\mathbf{I}(\boldsymbol{\theta})\) 是 \(k \times k\) 的 Fisher 信息矩阵。该展开依赖于严格的渐近假设：(i) 样本量 \(N \to \infty\) 而活跃模型维度 \(k\) 保持不变，(ii) MLE \(\hat{\boldsymbol{\theta}}\) 严格位于光滑参数流形的内部，(iii) Fisher 信息是正定且连续的。

对于有限样本或高维 (\(D \ge N\)) 情况下的正则非光滑模型（如 Lasso 或稀疏 SVM），这些假设不再成立。参数向量恰好位于非光滑边界上，这些边界代表 \(L_1\) 球的奇异角，Fisher 信息在这些角处变得退化或无法良定义。在这些配置中使用渐近展开会产生不准确的复杂度估计，导致过度惩罚和模型选择不稳定性的。相比之下，通过我们的余面积水平集框架（公式 (2)）计算精确随机复杂度，保证了有限样本下的数学正确性和通用编码最优性，绕过了边界奇异问题。

## 3 提出的 Schur-Sylvester 降维公式

在本节中，我们描述我们提出的 Schur-Sylvester 降维公式。为简化讨论，我们首先为 Lasso 估计量推导精确的维度约简，并展示高维投影和体积操作可以完全映射到低维活跃参数子空间上。然后我们将该方法推广到其他正则非光滑模型。

### 3.1 通过 Schur 补实现 KKT 求逆

为了具体阐述，我们首先使用经典的 Lasso 估计量作为主要基线，详细说明我们的 Schur 补降维方法。对 Lasso 的水平集约束求微分需要在每个投影步骤求解一个约束优化系统。相关的广义 KKT 矩阵具有如下块结构：

\[ \mathbf{V} = \begin{bmatrix} \mathbf{I}_N & \mathbf{G}^\top \\ \mathbf{G} & \mathbf{0} \end{bmatrix}, \tag{5} \]

其中 \(N\) 表示样本量，\(\mathbf{I}_N\) 表示 \(N \times N\) 的单位矩阵，\(\mathbf{G} = \mathbf{X}_A^\top \in \mathbb{R}^{k \times N}\) 是活跃特征矩阵 \(\mathbf{X}_A \in \mathbb{R}^{N \times k}\) 的转置，且 \(k \ll N\) 是活跃集大小。

根据块矩阵求逆定理，由于 \(\mathbf{I}_N\) 显然是可逆的，\(\mathbf{I}_N\) 的 Schur 补是 [19]：

\[ \mathbf{S} = \mathbf{0} - \mathbf{G} \mathbf{I}_N^{-1} \mathbf{G}^\top = -\mathbf{X}_A^\top \mathbf{X}_A \in \mathbb{R}^{k \times k}. \tag{6} \]

\(\mathbf{V}^{-1}\) 的左上块（定义我们的投影算子 \(\mathbf{D}_{\text{proj}}\)）可以精确地写成：

\[ \mathbf{D}_{\text{proj}} = \mathbf{I}_N - \mathbf{G}^\top (\mathbf{G} \mathbf{G}^\top)^{-1} \mathbf{G} = \mathbf{I}_N - \mathbf{X}_A (\mathbf{X}_A^\top \mathbf{X}_A)^{-1} \mathbf{X}_A^\top. \tag{7} \]

评估 \(\mathbf{X}_A^\top \mathbf{X}_A\) 需要 \(\mathcal{O}(N k^2)\) 次运算，其求逆需要 \(\mathcal{O}(k^3)\) 次运算。这完全避免了构建和求逆庞大的 \((N+k) \times (N+k)\) KKT 矩阵 \(\mathbf{V}\) 的需要。

### 3.2 通过 Sylvester 恒等式实现体积因子降维

水平集体积因子由投影切基的 Gram 矩阵 \(\mathbf{\Gamma}\) 计算得到：

\[ \mathbf{\Gamma} = \mathbf{B}^\top \mathbf{D}_{\text{proj}}^\top \mathbf{D}_{\text{proj}} \mathbf{B} \in \mathbb{R}^{(N-k) \times (N-k)}, \tag{8} \]

其中 \(\mathbf{B} \in \mathbb{R}^{N \times (N-k)}\) 是表示 Clarke 切锥的标准正交基矩阵。根据定义，\(\mathbf{B}\) 具有满足 \(\mathbf{B}^\top \mathbf{B} = \mathbf{I}_{N-k}\) 的标准正交列，其中 \(\mathbf{I}_{N-k}\) 是 \((N-k) \times (N-k)\) 的单位矩阵。由于投影算子 \(\mathbf{D}_{\text{proj}}\) 是对称且幂等的（\(\mathbf{D}_{\text{proj}}^\top \mathbf{D}_{\text{proj}} = \mathbf{D}_{\text{proj}}\)），我们将 Lasso 投影算子的 Schur 表示（公式 (7)）代入公式 (8) 得到：

\[ \mathbf{\Gamma} = \mathbf{B}^\top \left( \mathbf{I}_N - \mathbf{X}_A (\mathbf{X}_A^\top \mathbf{X}_A)^{-1} \mathbf{X}_A^\top \right) \mathbf{B} = \mathbf{I}_{N-k} - \mathbf{U}^\top (\mathbf{X}_A^\top \mathbf{X}_A)^{-1} \mathbf{U}, \tag{9} \]

其中 \(\mathbf{U} = \mathbf{X}_A^\top \mathbf{B} \in \mathbb{R}^{k \times (N-k)}\) 表示活跃列约束到切子空间的投影，且 \(\mathbf{U}^\top \in \mathbb{R}^{(N-k) \times k}\) 是其转置。关键的是，在 PPMH 框架中，投影算子 \(\mathbf{D}_{\text{proj}}\) 是在提议点 \(\mathbf{y}\) 处评估的（得到活跃集 \(\mathbf{X}_{A(\mathbf{y})}\)），而切基 \(\mathbf{B}\) 是在当前点 \(\mathbf{x}\) 处定义的。这种不匹配确保了 \(\mathbf{U} = \mathbf{X}_{A(\mathbf{y})}^\top \mathbf{B}(\mathbf{x}) \neq \mathbf{0}\)，使得体积修正 n