1. 项目概述从理论到实践的期望上界计算在概率论和算法分析里我们常常会遇到一些“长相”复杂的随机变量。它们可能由一个递归方程定义比如S 1 max{U * S‘, (1-U) * S‘’}其中U是均匀分布S‘和S‘’是它的独立副本。这种结构在分析快速选择算法QuickSelect的比较次数、某些分支过程的灭绝时间或是极值统计模型中并不少见。面对这样的S一个很自然的问题就是它的平均值也就是数学期望E[S]到底有多大理论上我们知道它存在且有限但给出一个具体、紧致的数值上界往往非常困难因为其分布函数F(x) P(S ≤ x)通常没有漂亮的闭式解。直接暴力数值模拟可以给出一个估计但那不是证明。我们需要的是一种可验证的、严格的上界计算方法。这就引出了本次分享的核心一套结合了单调分布函数迭代与保守网格离散化的方案。它的核心思想很直观既然我们无法精确求出F(x)那就构造一个函数序列L_n(x)确保每一项都是F(x)的下界即L_n(x) ≤ F(x)对所有x成立并且这个序列是单调不减的。那么当我们用1 - L_n(x)去替代1 - F(x)来计算期望的积分公式时得到的就是E[S]的一个上界。随着迭代进行L_n(x)越来越接近F(x)我们得到的上界也就越来越紧。这套方法的精妙之处在于它将一个无限维的函数不动点问题F KF通过精心设计的网格离散化转化为了一个有限维的向量迭代计算。同时在积分区间的尾部x A我们利用一个解析推导出的指数型上界exp(-x log(x/2) x)进行控制从而将无穷积分截断为有限积分加上一个可控的余项。整个流程就像是在用一套严格保守的“脚手架”从下方一点点搭建并逼近真实分布函数的形状最终为我们关心的期望值圈定一个可靠的范围。2. 核心原理分布函数、不动点方程与上界框架要理解整个方案我们需要先打好三个理论基础随机变量期望的积分表示、分布函数满足的不动点方程以及如何利用下界函数来构造期望上界。2.1 期望的尾部积分表示与问题拆分对于一个支撑集在[2, ∞)的非负随机变量S其期望有一个非常实用的表达式E[S] 2 ∫_{2}^{∞} (1 - F(x)) dx其中F(x) P(S ≤ x)是其分布函数。这个公式来源于分部积分直观理解就是S超过2的部分的期望等于其尾部概率的积分。我们的目标就是控制这个积分。直接计算不可能因为F(x)未知。策略是分而治之选取一个足够大的截断点A 2将积分拆分为两部分E[S] 2 ∫_{2}^{A} (1 - F(x)) dx ∫_{A}^{∞} (1 - F(x)) dx对于第一部分∫_{2}^{A} (1 - F(x)) dx我们计划用一个函数L(x)来逼近F(x)。如果我们能找到一个函数L(x)使得在区间[2, A)上恒有L(x) ≤ F(x)那么就有1 - F(x) ≤ 1 - L(x)从而∫_{2}^{A} (1 - F(x)) dx ≤ ∫_{2}^{A} (1 - L(x)) dx这样我们就用L(x)给出了第一部分积分的一个上界。对于第二部分尾部积分∫_{A}^{∞} (1 - F(x)) dx我们需要一个关于P(S x)的解析上界。原文中通过切尔诺夫界Chernoff bound等技巧推导出了一个强有力且形式简洁的上界对于t 2有P(S t) ≤ exp(-t log(t/2) t)这个上界衰减得比指数函数还快因为t log t项确保了尾部积分是收敛的。利用这个上界我们可以得到∫_{A}^{∞} (1 - F(x)) dx ∫_{A}^{∞} P(S x) dx ≤ ∫_{A}^{∞} exp(-x log(x/2) x) dx而这个积分可以通过一个引理原文中的 Lemma A.3进一步被一个显式表达式控制∫_{A}^{∞} exp(-x log(x/2) x) dx ≤ exp(-A log(A/2) A) / log(A/2)至此我们得到了计算E[S]上界的一个完整框架E[S] ≤ 2 ∫_{2}^{A} (1 - L(x)) dx exp(-A log(A/2) A) / log(A/2)其中L(x)是F(x)在[2, A)上的一个下界。接下来的所有工作都围绕着如何系统地、可计算地构造出一个尽可能紧致的L(x)。2.2 分布函数的不动点方程与单调算子为什么我们可以构造出F(x)的下界这源于S的定义所导致的分布函数F(x)必须满足的一个积分方程。根据S的递归定义S 1 max{U * S‘, (1-U) * S‘’}经过条件概率和独立性的推导可以得到如下不动点方程 对于所有实数xF(x) ∫_{0}^{1} F((x-1)/u) * F((x-1)/(1-u)) du这个方程的意义在于它将F在x点的值与F在另外两个通常更大的点(x-1)/u和(x-1)/(1-u)的值联系了起来。由于u在 (0,1) 内变化这实际上是一个非线性泛函方程。为了处理这个方程我们定义一个算子K它对任意一个取值在[0,1]的可测函数G作用如下(KG)(x) ∫_{0}^{1} G((x-1)/u) * G((x-1)/(1-u)) du那么上述不动点方程就可以简洁地写成F KF。这个算子K有两个至关重要的性质是后续所有迭代方案的基石保序性如果两个函数G和H满足G(x) ≤ H(x)对所有x成立那么也有(KG)(x) ≤ (KH)(x)对所有x成立。这是因为被积函数是非负的且每一项都满足不等式。保单调性如果G是非减函数这是分布函数的天然性质那么KG也是非减函数。这两个性质意味着如果我们从一个非减的、且是F的下界函数L_0出发即L_0(x) ≤ F(x)且L_0非减那么通过反复应用算子K得到的序列L_{n1} K L_n将始终保持是非减的并且始终是F的下界。更妙的是如果我们的初始猜测L_0还满足L_0 ≤ K L_0 L_1那么这个序列还是单调递增的L_0 ≤ L_1 ≤ L_2 ≤ ... ≤ F。这就为我们提供了一个从下方单调逼近真实分布函数F的完美机制。2.3 保守网格离散化将无限维问题有限化理论很美好但K算子涉及对一个连续变量u在 (0,1) 区间上的积分以及函数G在任意实数点上的求值这在计算机上是无法直接实现的。我们必须对其进行离散化。这里的关键词是“保守”意思是我们的离散化必须保证离散后得到的近似算子K̃其作用结果仍然是原算子K结果的一个下界。只有这样由K̃迭代产生的函数序列才能保证仍然是F的下界从而确保最终计算出的期望值是一个可靠的上界。具体做法如下对自变量x离散化在区间[2, A)上我们引入一个网格点2 x_0 x_1 ... x_N A。对于任意x属于[x_k, x_{k1})我们用x_k点的函数值来代表整个小区间的函数值。由于我们构造的函数都是非减的这种阶梯状近似是合理的。对积分变量u离散化在积分区间[0,1]上引入另一个网格点0 u_0 u_1 ... u_M 1。构造保守的积分近似对于固定的x_k和u属于子区间[u_r, u_{r1}]由于被积函数L((x_k-1)/u) * L((x_k-1)/(1-u))中的L是非减函数而(x_k-1)/u关于u是递减的(x_k-1)/(1-u)关于u是递增的我们可以找到一个保守的下界估计∫_{u_r}^{u_{r1}} L((x_k-1)/u) * L((x_k-1)/(1-u)) du ≥ (u_{r1} - u_r) * L((x_k-1)/u_{r1}) * L((x_k-1)/(1-u_r))直观理解在小区间[u_r, u_{r1}]上为了让乘积最小我们让第一个因子在u最大即u_{r1}时取值因为分母u越大(x_k-1)/u越小而L非减所以函数值也越小让第二个因子在u最小即u_r时取值因为1-u越大(x_k-1)/(1-u)越小函数值也越小。这个下界是容易计算的因为它只涉及在网格点u_{r1}和u_r处的函数值。将所有这些小区间的保守下界加起来我们就得到了原积分(KL)(x_k)的一个可计算的保守下界(K̃L)(x_k)(K̃L)(x_k) Σ_{r0}^{M-1} (u_{r1} - u_r) * L((x_k-1)/u_{r1}) * L((x_k-1)/(1-u_r))对于x不在网格点上的情况我们沿用阶梯近似若x ∈ [x_k, x_{k1})则定义(K̃L)(x) (K̃L)(x_k)。对于x 2定义(K̃L)(x) 0对于x ≥ A则与初始函数L_0(x)取最大值以保证迭代的稳定性。这样定义的离散算子K̃具有一个关键性质对于任何非减函数L在[2, A)上恒有(K̃L)(x) ≤ (KL)(x)。因此如果我们从下界函数L_0出发定义迭代L_{m1} K̃ L_m那么由数学归纳法可以证明每一个L_m都是F的非减下界。这就将无限的函数迭代问题完美地转化为了一个在有限网格点{x_0, ..., x_{N-1}}上进行的向量迭代问题。3. 算法实现从数学公式到可执行代码理论框架搭建好后实现就变得清晰而直接。整个算法可以看作一个“函数迭代器”只不过这个函数由它在N个网格点上的值所完全表征。下面我们一步步拆解实现细节。3.1 数据结构与初始化算法的输入是几个关键参数和初始函数截断点A例如A10。它需要大于2其选择平衡了计算量和精度A越大需要计算的网格区间[2, A)就越长但尾部余项exp(-A log(A/2) A) / log(A/2)会越小。网格数量N和M分别对应x网格和u网格的划分数。通常取N M以简化。N越大离散化越精细结果越接近理论值但计算量和内存消耗以O(N^2)增长。初始向量v0一个长度为N的数组代表初始下界函数L_0在x网格点x_0, x_1, ..., x_{N-1}上的取值。L_0必须是一个F(x)的保守下界。一个简单可靠的选择是利用已知的尾部上界L_0(x) max(0, 1 - exp(-x log(x/2) x))对于x ≥ 2。 因为P(S x) ≤ exp(...)所以P(S ≤ x) 1 - P(S x) ≥ 1 - exp(...)。注意当x较小时1 - exp(...)可能为负所以要与0取最大值以保证非负性。在网格点x_k上计算这个值就得到了v0[k]。此外我们还需要预计算两个重要的映射数组以加速核心迭代索引映射表idx_u和idx_1mu对于每一对(k, r)我们需要计算L((x_k-1)/u_{r1})和L((x_k-1)/(1-u_r))。由于(x_k-1)/u_{r1}和(x_k-1)/(1-u_r)很可能不在x网格点上我们需要找到它们所落入的网格区间并用该区间左端点的函数值来近似因为L是阶梯函数且非减左端点值是最保守的下界估计。我们可以预先计算好这些参数对应的x网格索引。具体来说对于每个(k, r)计算target1 (x_k - 1) / u_{r1}target2 (x_k - 1) / (1 - u_r)然后在x网格中寻找最大的索引i使得x_i ≤ target1。这个i就是target1对应的索引。如果target1 x_0 2则索引记为-1对应函数值为0。target2同理。将这些索引存储在两个N x M的数组idx_u和idx_1mu中可以避免在每次迭代中进行耗时的二分查找。3.2 核心迭代步骤假设在第m轮迭代我们有一个向量v_old代表L_m在x网格点上的值。我们的目标是计算v_new代表L_{m1} K̃ L_m在网格点上的值。对于每一个网格点索引k(0 ≤ k N)我们需要计算v_new[k] Σ_{r0}^{M-1} (u_{r1} - u_r) * v_old[ idx_u[k, r] ] * v_old[ idx_1mu[k, r] ]这里idx_u[k, r]和idx_1mu[k, r]就是上一步预计算的索引。如果索引为-1则对应的v_old值取0。这个过程本质上是一个巨大的求和非常适合用向量化操作或并行计算来加速。在 Python 中利用 NumPy 可以高效地实现import numpy as np def mesh_iteration(v_old, du_grid, idx_u, idx_1mu): 执行一次保守网格迭代。 v_old: 当前迭代函数值向量形状 (N,) du_grid: u网格的宽度数组形状 (M,)du_grid[r] u_{r1} - u_r idx_u, idx_1mu: 预计算的索引数组形状 (N, M) # 将索引为-1的位置对应的值置为0 val_u np.where(idx_u 0, v_old[idx_u], 0.0) val_1mu np.where(idx_1mu 0, v_old[idx_1mu], 0.0) # 计算每个r对应的乘积并沿r方向axis1与du_grid加权求和 # 注意这里du_grid需要扩展维度以进行广播 v_new np.sum(du_grid[np.newaxis, :] * val_u * val_1mu, axis1) return v_new3.3 上界计算与迭代终止每次迭代得到新的向量v_new后我们可以立即计算当前对应的期望上界U_mU_m 2 Σ_{k0}^{N-1} (x_{k1} - x_k) * (1 - v_new[k]) exp(-A * log(A/2) A) / log(A/2)其中Σ_{k0}^{N-1} (x_{k1} - x_k) * (1 - v_new[k])是对积分∫_{2}^{A} (1 - L_{m1}(x)) dx的矩形法近似。由于v_new是L_{m1}在左端点的值且1 - L_{m1}(x)是非增的这个矩形和是积分的一个上界与我们整体的“上界”策略一致。迭代何时停止理论上序列{U_m}是单调递减的因为{L_m}单调递增并且有下界E[S]。在实践中我们可以设定一个最大迭代次数例如50次或者当相邻两次迭代的上界值U_m和U_{m-1}的相对变化小于某个容差例如1e-12时停止。由于离散化误差的存在迭代最终会收敛到离散问题的一个不动点这个不动点对应的上界U_*就是我们能得到的最佳估计。3.4 一个完整的实现示例以下是一个简化的、用于演示算法流程的 Python 代码框架。在实际严谨的应用中需要考虑高精度计算和误差累积问题。import numpy as np from math import log, exp def compute_expectation_upper_bound(A10.0, N200, M200, max_iter50, tol1e-12): 计算期望E[S]上界的主函数。 # 1. 创建网格 x_grid np.linspace(2.0, A, N1) # N1个点包括端点A u_grid np.linspace(0.0, 1.0, M1) # 网格宽度 dx x_grid[1] - x_grid[0] du u_grid[1] - u_grid[0] # x的采样点左端点共N个 x_samples x_grid[:-1] # u的采样点左端点和宽度数组 u_samples u_grid[:-1] du_array np.full(M, du) # 均匀网格所有宽度相同 # 2. 初始化 L0 def tail_bound(t): if t 2: return 1.0 # 实际上对于t2, P(St) 1这里初始下界取0更简单。 return exp(-t * log(t/2.0) t) # 初始下界函数 L0(x) max(0, 1 - tail_bound(x)) v_current np.maximum(0.0, 1.0 - np.array([tail_bound(x) for x in x_samples])) # 3. 预计算索引映射 (这里简化处理采用逐元素计算实际应优化) # 为清晰起见这里用循环表示逻辑。实际应用应使用向量化优化。 idx_u np.full((N, M), -1, dtypeint) idx_1mu np.full((N, M), -1, dtypeint) for k in range(N): xk x_samples[k] for r in range(M): u_r_plus_1 u_grid[r1] u_r u_grid[r] target1 (xk - 1.0) / u_r_plus_1 if u_r_plus_1 1e-12 else float(inf) target2 (xk - 1.0) / (1.0 - u_r) if (1.0 - u_r) 1e-12 else float(inf) # 在x_grid中寻找左端点索引 if target1 2.0: # 找到最大的i使得 x_grid[i] target1 i np.searchsorted(x_grid, target1, sideright) - 1 if i N: # 确保索引在v_current范围内 idx_u[k, r] i if target2 2.0: j np.searchsorted(x_grid, target2, sideright) - 1 if j N: idx_1mu[k, r] j # 4. 迭代计算 bounds [] # 记录每次迭代的上界 for it in range(max_iter): # 核心迭代 val_from_u np.where(idx_u 0, v_current[idx_u], 0.0) val_from_1mu np.where(idx_1mu 0, v_current[idx_1mu], 0.0) v_next np.sum(du_array[np.newaxis, :] * val_from_u * val_from_1mu, axis1) # 确保函数值在[0,1]区间内理论上应自动满足数值计算可加裁剪 v_next np.clip(v_next, 0.0, 1.0) # 计算当前上界 integral_approx np.sum(dx * (1.0 - v_next)) tail_term exp(-A * log(A/2.0) A) / log(A/2.0) current_bound 2.0 integral_approx tail_term bounds.append(current_bound) # 检查收敛 if it 0 and abs(bounds[-1] - bounds[-2]) / (abs(bounds[-2]) 1e-15) tol: print(f在 {it1} 次迭代后收敛。) break # 更新 v_current v_next else: print(f达到最大迭代次数 {max_iter}。) return bounds, v_current # 运行示例 bounds, final_L compute_expectation_upper_bound(A10.0, N100, M100, max_iter30) print(f最终上界估计: {bounds[-1]:.6f}) print(f上界序列 (前几次): {bounds[:5]})注意以上代码是概念演示未进行严格的数值误差控制。在真正的“计算机辅助证明”中需要使用区间算术Interval Arithmetic库来保证每一步计算产生的舍入误差都被严格包含在某个区间内从而使得最终结果U_m即使在考虑所有数值误差后仍然是E[S]的一个数学上严格成立的上界。Python 中的mpmath库或专门的pyinterval库可以用于此类目的。4. 数值实验解析与参数影响分析根据原文的浮动小数点实验我们能看到这套方法的实际表现。以A10, NM100为例迭代大约在10步以内就快速收敛得到的期望上界U_m稳定在4.34附近。当把网格加密到NM10000时上界进一步下降到4.09左右。这个数值为我们理解随机变量S的尺度提供了宝贵的定量信息。4.1 网格密度与计算精度的权衡网格数量N和M是影响结果精度和计算成本的核心参数。精度更密的网格更大的N, M能更精确地逼近原积分算子K从而得到更紧致更小的上界U_*。从100网格到10000网格上界从4.34降至4.09改善显著。成本计算成本主要来自核心迭代步骤中的双重循环求和。其时间复杂度为O(N * M)空间复杂度则至少需要存储idx_u和idx_1mu这两个N x M的索引数组即O(N * M)。当NM10000时仅索引数组就需要存储2 * 10^8个整数约800MB迭代计算量也极大。因此在实际应用中需要在可接受的计算资源内选择尽可能大的网格。实操建议可以从中等网格如NM500开始测试观察上界随迭代的变化趋势和收敛值。然后逐步增加网格密度直到连续两次加密网格得到的结果差异小于你的精度要求。也可以考虑非均匀网格在函数变化剧烈的区域如x靠近2或u靠近0、1的区域使用更密的划分。4.2 截断点A的选择策略截断点A的选择是一个权衡艺术A太小尾部余项exp(-A log(A/2) A) / log(A/2)会比较大从而直接抬高了上界即使[2, A)区间内的积分估计再精确总上界也可能不理想。A太大主要计算区间[2, A)变长为了保持相同的离散化精度需要按比例增加网格点数量N否则dx会变大导致积分近似误差增大。同时当x很大时(x-1)/u会变得极大可能超出我们预设的x网格范围我们只定义到A在算法中这部分会被截断处理例如用L(A)或L0(x)的值可能引入误差。一个实用的方法是先根据对S量级的粗略估计比如通过模拟选择一个适中的A例如使P(S A)非常小。运行算法后检查尾部余项在总上界U_m中的占比。如果占比过高比如超过1%则应增大A重新计算如果占比极低则可以尝试略微减小A以缩减计算区间。原文中选择A10是合理的因为对于这个具体的S其期望值大约在4左右A10已经覆盖了分布的主要部分。4.3 初始函数L_0的选取与迭代收敛初始函数L_0需要满足两个条件1) 是F的下界2) 是非减函数。使用尾部上界导出的L_0(x) max(0, 1 - exp(-x log(x/2) x))是一个安全且有效的选择。它虽然不是最优的下界但保证了迭代的起点是“保守”的。从图3(a)可以看到迭代序列{L_m}从初始的L_0一条从0开始快速上升的曲线开始随着迭代次数的增加曲线整体向上移动形状也逐渐趋于稳定。这直观地展示了算子K的平滑效应和迭代的收敛性。通常前几次迭代带来的改进最大之后逐渐趋于平稳。收敛性判断除了监视上界序列{U_m}的变化还可以观察向量v_m的差异例如计算连续两次迭代向量的无穷范数差||v_{m1} - v_m||_∞。当这个差小于预设容差时可以认为离散的不动点已经找到。4.4 从数值实验到严格证明需要极度警惕的是文中提到的浮动小数点计算如在 Google Colab 中用 NumPy 完成只是一个“试点计算”或“数值实验”。它为我们提供了上界可能在哪里的强烈证据但不是数学证明。原因在于浮点运算存在舍入误差这些误差在迭代过程中可能会被放大破坏“保守性”。我们无法保证浮点计算出的v_next严格满足v_next ≤ K(v_current)在离散意义上。要获得一个计算机辅助的严格证明必须使用区间算术。区间算术的基本思想是用一个区间[a, b]来表示一个实数这个区间 guaranteed 包含了该数的真实值。所有算术运算加、减、乘、除都定义在区间上并且结果区间 guaranteed 包含了所有可能真实结果。在本文的算法中所有网格点x_k,u_r用区间表示例如精确的有理数或高精度浮点区间。初始向量v0的每个分量都是一个保守的区间下界例如[0, some_upper_bound]但左端点需是真实下界。在计算v_new[k]的求和时每一项(u_{r1} - u_r) * v_old[i] * v_old[j]都使用区间乘法求和使用区间加法。最终计算出的v_new的每个分量都是一个区间并且可以数学证明真实的不动点函数值位于该区间内。同样用区间算术计算积分∫ (1 - L(x)) dx和尾部余项最终得到的U_m是一个区间并且可以证明真实的E[S]一定小于等于这个区间的右端点。只有这样我们才能说“E[S] ≤ 4.09”是一个被严格证明的定理而不仅仅是一个数值观察。这方面的工具包括MPFI(基于MPFR的区间算术库)、Julia语言的IntervalArithmetic.jl包等。5. 方案总结与扩展应用思考回顾整个方案其核心贡献在于提供了一条系统化的、可计算的路径来获取复杂随机变量期望的严格上界。它将概率分析推导尾部上界和不动点方程、数值分析保守网格离散化和计算机科学算法实现与可验证计算紧密结合。核心优势严格性整个框架建立在严格的不等式基础上。只要初始条件L_0 ≤ F成立且离散化是保守的那么最终得到的上界在数学上就是成立的。可自动化整个过程可以编码实现一旦模型即S的递归定义和尾部上界形式确定计算上界的过程几乎是机械化的。可优化性通过加密网格、优化初始函数、调整截断点A可以系统地改进降低上界。局限性与挑战计算复杂度O(N*M)的复杂度对于高精度要求是个挑战。可能需要借助更高效的数值积分方法如高斯积分或自适应网格来优化。维度灾难本方法针对一维分布函数。对于高维随机向量其联合分布函数的逼近会变得异常复杂。模型依赖性方法的成功实施强烈依赖于能够推导出类似F KF的不动点方程以及一个形式良好的尾部概率上界。对于更复杂的递归关系推导可能非常困难。扩展应用 这套单调分布函数方案的思想并不局限于原文中的特定随机变量S。它可以推广到任何满足类似分布递归等式的随机变量例如其他快速选择算法变体分析QuickSelect在不同 pivot 选择策略下的比较次数。分支随机游走某些分支过程的生存时间或最大位移。递归算法复杂度满足特定分治递归的算法运行时间。极值统计某些关联模型中最大值的分布。在实际操作中最关键的一步是为你所研究的随机变量Y建立起它的分布函数F_Y(x)所满足的方程通常通过全概率公式或条件期望并找到一个合适的保守积分算子K。一旦做到了这一点剩下的网格迭代和上界计算就是一套相对固定的流程了。最后分享一个我个人的调试心得在实现这类保守迭代算法时可视化是极其强大的工具。就像原文中的图2和图3绘制出每次迭代的上界值U_m和函数曲线L_m(x)能让你直观地看到收敛过程快速发现参数设置如A,N是否合理以及算法实现中可能存在的错误比如曲线出现非单调性这一定意味着代码有 bug。将数学的严谨与计算的直观相结合是解决这类问题的有效法门。
基于单调迭代与保守离散化的随机变量期望上界计算
发布时间：2026/6/2 4:56:22
1. 项目概述从理论到实践的期望上界计算在概率论和算法分析里我们常常会遇到一些“长相”复杂的随机变量。它们可能由一个递归方程定义比如S 1 max{U * S‘, (1-U) * S‘’}其中U是均匀分布S‘和S‘’是它的独立副本。这种结构在分析快速选择算法QuickSelect的比较次数、某些分支过程的灭绝时间或是极值统计模型中并不少见。面对这样的S一个很自然的问题就是它的平均值也就是数学期望E[S]到底有多大理论上我们知道它存在且有限但给出一个具体、紧致的数值上界往往非常困难因为其分布函数F(x) P(S ≤ x)通常没有漂亮的闭式解。直接暴力数值模拟可以给出一个估计但那不是证明。我们需要的是一种可验证的、严格的上界计算方法。这就引出了本次分享的核心一套结合了单调分布函数迭代与保守网格离散化的方案。它的核心思想很直观既然我们无法精确求出F(x)那就构造一个函数序列L_n(x)确保每一项都是F(x)的下界即L_n(x) ≤ F(x)对所有x成立并且这个序列是单调不减的。那么当我们用1 - L_n(x)去替代1 - F(x)来计算期望的积分公式时得到的就是E[S]的一个上界。随着迭代进行L_n(x)越来越接近F(x)我们得到的上界也就越来越紧。这套方法的精妙之处在于它将一个无限维的函数不动点问题F KF通过精心设计的网格离散化转化为了一个有限维的向量迭代计算。同时在积分区间的尾部x A我们利用一个解析推导出的指数型上界exp(-x log(x/2) x)进行控制从而将无穷积分截断为有限积分加上一个可控的余项。整个流程就像是在用一套严格保守的“脚手架”从下方一点点搭建并逼近真实分布函数的形状最终为我们关心的期望值圈定一个可靠的范围。2. 核心原理分布函数、不动点方程与上界框架要理解整个方案我们需要先打好三个理论基础随机变量期望的积分表示、分布函数满足的不动点方程以及如何利用下界函数来构造期望上界。2.1 期望的尾部积分表示与问题拆分对于一个支撑集在[2, ∞)的非负随机变量S其期望有一个非常实用的表达式E[S] 2 ∫_{2}^{∞} (1 - F(x)) dx其中F(x) P(S ≤ x)是其分布函数。这个公式来源于分部积分直观理解就是S超过2的部分的期望等于其尾部概率的积分。我们的目标就是控制这个积分。直接计算不可能因为F(x)未知。策略是分而治之选取一个足够大的截断点A 2将积分拆分为两部分E[S] 2 ∫_{2}^{A} (1 - F(x)) dx ∫_{A}^{∞} (1 - F(x)) dx对于第一部分∫_{2}^{A} (1 - F(x)) dx我们计划用一个函数L(x)来逼近F(x)。如果我们能找到一个函数L(x)使得在区间[2, A)上恒有L(x) ≤ F(x)那么就有1 - F(x) ≤ 1 - L(x)从而∫_{2}^{A} (1 - F(x)) dx ≤ ∫_{2}^{A} (1 - L(x)) dx这样我们就用L(x)给出了第一部分积分的一个上界。对于第二部分尾部积分∫_{A}^{∞} (1 - F(x)) dx我们需要一个关于P(S x)的解析上界。原文中通过切尔诺夫界Chernoff bound等技巧推导出了一个强有力且形式简洁的上界对于t 2有P(S t) ≤ exp(-t log(t/2) t)这个上界衰减得比指数函数还快因为t log t项确保了尾部积分是收敛的。利用这个上界我们可以得到∫_{A}^{∞} (1 - F(x)) dx ∫_{A}^{∞} P(S x) dx ≤ ∫_{A}^{∞} exp(-x log(x/2) x) dx而这个积分可以通过一个引理原文中的 Lemma A.3进一步被一个显式表达式控制∫_{A}^{∞} exp(-x log(x/2) x) dx ≤ exp(-A log(A/2) A) / log(A/2)至此我们得到了计算E[S]上界的一个完整框架E[S] ≤ 2 ∫_{2}^{A} (1 - L(x)) dx exp(-A log(A/2) A) / log(A/2)其中L(x)是F(x)在[2, A)上的一个下界。接下来的所有工作都围绕着如何系统地、可计算地构造出一个尽可能紧致的L(x)。2.2 分布函数的不动点方程与单调算子为什么我们可以构造出F(x)的下界这源于S的定义所导致的分布函数F(x)必须满足的一个积分方程。根据S的递归定义S 1 max{U * S‘, (1-U) * S‘’}经过条件概率和独立性的推导可以得到如下不动点方程 对于所有实数xF(x) ∫_{0}^{1} F((x-1)/u) * F((x-1)/(1-u)) du这个方程的意义在于它将F在x点的值与F在另外两个通常更大的点(x-1)/u和(x-1)/(1-u)的值联系了起来。由于u在 (0,1) 内变化这实际上是一个非线性泛函方程。为了处理这个方程我们定义一个算子K它对任意一个取值在[0,1]的可测函数G作用如下(KG)(x) ∫_{0}^{1} G((x-1)/u) * G((x-1)/(1-u)) du那么上述不动点方程就可以简洁地写成F KF。这个算子K有两个至关重要的性质是后续所有迭代方案的基石保序性如果两个函数G和H满足G(x) ≤ H(x)对所有x成立那么也有(KG)(x) ≤ (KH)(x)对所有x成立。这是因为被积函数是非负的且每一项都满足不等式。保单调性如果G是非减函数这是分布函数的天然性质那么KG也是非减函数。这两个性质意味着如果我们从一个非减的、且是F的下界函数L_0出发即L_0(x) ≤ F(x)且L_0非减那么通过反复应用算子K得到的序列L_{n1} K L_n将始终保持是非减的并且始终是F的下界。更妙的是如果我们的初始猜测L_0还满足L_0 ≤ K L_0 L_1那么这个序列还是单调递增的L_0 ≤ L_1 ≤ L_2 ≤ ... ≤ F。这就为我们提供了一个从下方单调逼近真实分布函数F的完美机制。2.3 保守网格离散化将无限维问题有限化理论很美好但K算子涉及对一个连续变量u在 (0,1) 区间上的积分以及函数G在任意实数点上的求值这在计算机上是无法直接实现的。我们必须对其进行离散化。这里的关键词是“保守”意思是我们的离散化必须保证离散后得到的近似算子K̃其作用结果仍然是原算子K结果的一个下界。只有这样由K̃迭代产生的函数序列才能保证仍然是F的下界从而确保最终计算出的期望值是一个可靠的上界。具体做法如下对自变量x离散化在区间[2, A)上我们引入一个网格点2 x_0 x_1 ... x_N A。对于任意x属于[x_k, x_{k1})我们用x_k点的函数值来代表整个小区间的函数值。由于我们构造的函数都是非减的这种阶梯状近似是合理的。对积分变量u离散化在积分区间[0,1]上引入另一个网格点0 u_0 u_1 ... u_M 1。构造保守的积分近似对于固定的x_k和u属于子区间[u_r, u_{r1}]由于被积函数L((x_k-1)/u) * L((x_k-1)/(1-u))中的L是非减函数而(x_k-1)/u关于u是递减的(x_k-1)/(1-u)关于u是递增的我们可以找到一个保守的下界估计∫_{u_r}^{u_{r1}} L((x_k-1)/u) * L((x_k-1)/(1-u)) du ≥ (u_{r1} - u_r) * L((x_k-1)/u_{r1}) * L((x_k-1)/(1-u_r))直观理解在小区间[u_r, u_{r1}]上为了让乘积最小我们让第一个因子在u最大即u_{r1}时取值因为分母u越大(x_k-1)/u越小而L非减所以函数值也越小让第二个因子在u最小即u_r时取值因为1-u越大(x_k-1)/(1-u)越小函数值也越小。这个下界是容易计算的因为它只涉及在网格点u_{r1}和u_r处的函数值。将所有这些小区间的保守下界加起来我们就得到了原积分(KL)(x_k)的一个可计算的保守下界(K̃L)(x_k)(K̃L)(x_k) Σ_{r0}^{M-1} (u_{r1} - u_r) * L((x_k-1)/u_{r1}) * L((x_k-1)/(1-u_r))对于x不在网格点上的情况我们沿用阶梯近似若x ∈ [x_k, x_{k1})则定义(K̃L)(x) (K̃L)(x_k)。对于x 2定义(K̃L)(x) 0对于x ≥ A则与初始函数L_0(x)取最大值以保证迭代的稳定性。这样定义的离散算子K̃具有一个关键性质对于任何非减函数L在[2, A)上恒有(K̃L)(x) ≤ (KL)(x)。因此如果我们从下界函数L_0出发定义迭代L_{m1} K̃ L_m那么由数学归纳法可以证明每一个L_m都是F的非减下界。这就将无限的函数迭代问题完美地转化为了一个在有限网格点{x_0, ..., x_{N-1}}上进行的向量迭代问题。3. 算法实现从数学公式到可执行代码理论框架搭建好后实现就变得清晰而直接。整个算法可以看作一个“函数迭代器”只不过这个函数由它在N个网格点上的值所完全表征。下面我们一步步拆解实现细节。3.1 数据结构与初始化算法的输入是几个关键参数和初始函数截断点A例如A10。它需要大于2其选择平衡了计算量和精度A越大需要计算的网格区间[2, A)就越长但尾部余项exp(-A log(A/2) A) / log(A/2)会越小。网格数量N和M分别对应x网格和u网格的划分数。通常取N M以简化。N越大离散化越精细结果越接近理论值但计算量和内存消耗以O(N^2)增长。初始向量v0一个长度为N的数组代表初始下界函数L_0在x网格点x_0, x_1, ..., x_{N-1}上的取值。L_0必须是一个F(x)的保守下界。一个简单可靠的选择是利用已知的尾部上界L_0(x) max(0, 1 - exp(-x log(x/2) x))对于x ≥ 2。 因为P(S x) ≤ exp(...)所以P(S ≤ x) 1 - P(S x) ≥ 1 - exp(...)。注意当x较小时1 - exp(...)可能为负所以要与0取最大值以保证非负性。在网格点x_k上计算这个值就得到了v0[k]。此外我们还需要预计算两个重要的映射数组以加速核心迭代索引映射表idx_u和idx_1mu对于每一对(k, r)我们需要计算L((x_k-1)/u_{r1})和L((x_k-1)/(1-u_r))。由于(x_k-1)/u_{r1}和(x_k-1)/(1-u_r)很可能不在x网格点上我们需要找到它们所落入的网格区间并用该区间左端点的函数值来近似因为L是阶梯函数且非减左端点值是最保守的下界估计。我们可以预先计算好这些参数对应的x网格索引。具体来说对于每个(k, r)计算target1 (x_k - 1) / u_{r1}target2 (x_k - 1) / (1 - u_r)然后在x网格中寻找最大的索引i使得x_i ≤ target1。这个i就是target1对应的索引。如果target1 x_0 2则索引记为-1对应函数值为0。target2同理。将这些索引存储在两个N x M的数组idx_u和idx_1mu中可以避免在每次迭代中进行耗时的二分查找。3.2 核心迭代步骤假设在第m轮迭代我们有一个向量v_old代表L_m在x网格点上的值。我们的目标是计算v_new代表L_{m1} K̃ L_m在网格点上的值。对于每一个网格点索引k(0 ≤ k N)我们需要计算v_new[k] Σ_{r0}^{M-1} (u_{r1} - u_r) * v_old[ idx_u[k, r] ] * v_old[ idx_1mu[k, r] ]这里idx_u[k, r]和idx_1mu[k, r]就是上一步预计算的索引。如果索引为-1则对应的v_old值取0。这个过程本质上是一个巨大的求和非常适合用向量化操作或并行计算来加速。在 Python 中利用 NumPy 可以高效地实现import numpy as np def mesh_iteration(v_old, du_grid, idx_u, idx_1mu): 执行一次保守网格迭代。 v_old: 当前迭代函数值向量形状 (N,) du_grid: u网格的宽度数组形状 (M,)du_grid[r] u_{r1} - u_r idx_u, idx_1mu: 预计算的索引数组形状 (N, M) # 将索引为-1的位置对应的值置为0 val_u np.where(idx_u 0, v_old[idx_u], 0.0) val_1mu np.where(idx_1mu 0, v_old[idx_1mu], 0.0) # 计算每个r对应的乘积并沿r方向axis1与du_grid加权求和 # 注意这里du_grid需要扩展维度以进行广播 v_new np.sum(du_grid[np.newaxis, :] * val_u * val_1mu, axis1) return v_new3.3 上界计算与迭代终止每次迭代得到新的向量v_new后我们可以立即计算当前对应的期望上界U_mU_m 2 Σ_{k0}^{N-1} (x_{k1} - x_k) * (1 - v_new[k]) exp(-A * log(A/2) A) / log(A/2)其中Σ_{k0}^{N-1} (x_{k1} - x_k) * (1 - v_new[k])是对积分∫_{2}^{A} (1 - L_{m1}(x)) dx的矩形法近似。由于v_new是L_{m1}在左端点的值且1 - L_{m1}(x)是非增的这个矩形和是积分的一个上界与我们整体的“上界”策略一致。迭代何时停止理论上序列{U_m}是单调递减的因为{L_m}单调递增并且有下界E[S]。在实践中我们可以设定一个最大迭代次数例如50次或者当相邻两次迭代的上界值U_m和U_{m-1}的相对变化小于某个容差例如1e-12时停止。由于离散化误差的存在迭代最终会收敛到离散问题的一个不动点这个不动点对应的上界U_*就是我们能得到的最佳估计。3.4 一个完整的实现示例以下是一个简化的、用于演示算法流程的 Python 代码框架。在实际严谨的应用中需要考虑高精度计算和误差累积问题。import numpy as np from math import log, exp def compute_expectation_upper_bound(A10.0, N200, M200, max_iter50, tol1e-12): 计算期望E[S]上界的主函数。 # 1. 创建网格 x_grid np.linspace(2.0, A, N1) # N1个点包括端点A u_grid np.linspace(0.0, 1.0, M1) # 网格宽度 dx x_grid[1] - x_grid[0] du u_grid[1] - u_grid[0] # x的采样点左端点共N个 x_samples x_grid[:-1] # u的采样点左端点和宽度数组 u_samples u_grid[:-1] du_array np.full(M, du) # 均匀网格所有宽度相同 # 2. 初始化 L0 def tail_bound(t): if t 2: return 1.0 # 实际上对于t2, P(St) 1这里初始下界取0更简单。 return exp(-t * log(t/2.0) t) # 初始下界函数 L0(x) max(0, 1 - tail_bound(x)) v_current np.maximum(0.0, 1.0 - np.array([tail_bound(x) for x in x_samples])) # 3. 预计算索引映射 (这里简化处理采用逐元素计算实际应优化) # 为清晰起见这里用循环表示逻辑。实际应用应使用向量化优化。 idx_u np.full((N, M), -1, dtypeint) idx_1mu np.full((N, M), -1, dtypeint) for k in range(N): xk x_samples[k] for r in range(M): u_r_plus_1 u_grid[r1] u_r u_grid[r] target1 (xk - 1.0) / u_r_plus_1 if u_r_plus_1 1e-12 else float(inf) target2 (xk - 1.0) / (1.0 - u_r) if (1.0 - u_r) 1e-12 else float(inf) # 在x_grid中寻找左端点索引 if target1 2.0: # 找到最大的i使得 x_grid[i] target1 i np.searchsorted(x_grid, target1, sideright) - 1 if i N: # 确保索引在v_current范围内 idx_u[k, r] i if target2 2.0: j np.searchsorted(x_grid, target2, sideright) - 1 if j N: idx_1mu[k, r] j # 4. 迭代计算 bounds [] # 记录每次迭代的上界 for it in range(max_iter): # 核心迭代 val_from_u np.where(idx_u 0, v_current[idx_u], 0.0) val_from_1mu np.where(idx_1mu 0, v_current[idx_1mu], 0.0) v_next np.sum(du_array[np.newaxis, :] * val_from_u * val_from_1mu, axis1) # 确保函数值在[0,1]区间内理论上应自动满足数值计算可加裁剪 v_next np.clip(v_next, 0.0, 1.0) # 计算当前上界 integral_approx np.sum(dx * (1.0 - v_next)) tail_term exp(-A * log(A/2.0) A) / log(A/2.0) current_bound 2.0 integral_approx tail_term bounds.append(current_bound) # 检查收敛 if it 0 and abs(bounds[-1] - bounds[-2]) / (abs(bounds[-2]) 1e-15) tol: print(f在 {it1} 次迭代后收敛。) break # 更新 v_current v_next else: print(f达到最大迭代次数 {max_iter}。) return bounds, v_current # 运行示例 bounds, final_L compute_expectation_upper_bound(A10.0, N100, M100, max_iter30) print(f最终上界估计: {bounds[-1]:.6f}) print(f上界序列 (前几次): {bounds[:5]})注意以上代码是概念演示未进行严格的数值误差控制。在真正的“计算机辅助证明”中需要使用区间算术Interval Arithmetic库来保证每一步计算产生的舍入误差都被严格包含在某个区间内从而使得最终结果U_m即使在考虑所有数值误差后仍然是E[S]的一个数学上严格成立的上界。Python 中的mpmath库或专门的pyinterval库可以用于此类目的。4. 数值实验解析与参数影响分析根据原文的浮动小数点实验我们能看到这套方法的实际表现。以A10, NM100为例迭代大约在10步以内就快速收敛得到的期望上界U_m稳定在4.34附近。当把网格加密到NM10000时上界进一步下降到4.09左右。这个数值为我们理解随机变量S的尺度提供了宝贵的定量信息。4.1 网格密度与计算精度的权衡网格数量N和M是影响结果精度和计算成本的核心参数。精度更密的网格更大的N, M能更精确地逼近原积分算子K从而得到更紧致更小的上界U_*。从100网格到10000网格上界从4.34降至4.09改善显著。成本计算成本主要来自核心迭代步骤中的双重循环求和。其时间复杂度为O(N * M)空间复杂度则至少需要存储idx_u和idx_1mu这两个N x M的索引数组即O(N * M)。当NM10000时仅索引数组就需要存储2 * 10^8个整数约800MB迭代计算量也极大。因此在实际应用中需要在可接受的计算资源内选择尽可能大的网格。实操建议可以从中等网格如NM500开始测试观察上界随迭代的变化趋势和收敛值。然后逐步增加网格密度直到连续两次加密网格得到的结果差异小于你的精度要求。也可以考虑非均匀网格在函数变化剧烈的区域如x靠近2或u靠近0、1的区域使用更密的划分。4.2 截断点A的选择策略截断点A的选择是一个权衡艺术A太小尾部余项exp(-A log(A/2) A) / log(A/2)会比较大从而直接抬高了上界即使[2, A)区间内的积分估计再精确总上界也可能不理想。A太大主要计算区间[2, A)变长为了保持相同的离散化精度需要按比例增加网格点数量N否则dx会变大导致积分近似误差增大。同时当x很大时(x-1)/u会变得极大可能超出我们预设的x网格范围我们只定义到A在算法中这部分会被截断处理例如用L(A)或L0(x)的值可能引入误差。一个实用的方法是先根据对S量级的粗略估计比如通过模拟选择一个适中的A例如使P(S A)非常小。运行算法后检查尾部余项在总上界U_m中的占比。如果占比过高比如超过1%则应增大A重新计算如果占比极低则可以尝试略微减小A以缩减计算区间。原文中选择A10是合理的因为对于这个具体的S其期望值大约在4左右A10已经覆盖了分布的主要部分。4.3 初始函数L_0的选取与迭代收敛初始函数L_0需要满足两个条件1) 是F的下界2) 是非减函数。使用尾部上界导出的L_0(x) max(0, 1 - exp(-x log(x/2) x))是一个安全且有效的选择。它虽然不是最优的下界但保证了迭代的起点是“保守”的。从图3(a)可以看到迭代序列{L_m}从初始的L_0一条从0开始快速上升的曲线开始随着迭代次数的增加曲线整体向上移动形状也逐渐趋于稳定。这直观地展示了算子K的平滑效应和迭代的收敛性。通常前几次迭代带来的改进最大之后逐渐趋于平稳。收敛性判断除了监视上界序列{U_m}的变化还可以观察向量v_m的差异例如计算连续两次迭代向量的无穷范数差||v_{m1} - v_m||_∞。当这个差小于预设容差时可以认为离散的不动点已经找到。4.4 从数值实验到严格证明需要极度警惕的是文中提到的浮动小数点计算如在 Google Colab 中用 NumPy 完成只是一个“试点计算”或“数值实验”。它为我们提供了上界可能在哪里的强烈证据但不是数学证明。原因在于浮点运算存在舍入误差这些误差在迭代过程中可能会被放大破坏“保守性”。我们无法保证浮点计算出的v_next严格满足v_next ≤ K(v_current)在离散意义上。要获得一个计算机辅助的严格证明必须使用区间算术。区间算术的基本思想是用一个区间[a, b]来表示一个实数这个区间 guaranteed 包含了该数的真实值。所有算术运算加、减、乘、除都定义在区间上并且结果区间 guaranteed 包含了所有可能真实结果。在本文的算法中所有网格点x_k,u_r用区间表示例如精确的有理数或高精度浮点区间。初始向量v0的每个分量都是一个保守的区间下界例如[0, some_upper_bound]但左端点需是真实下界。在计算v_new[k]的求和时每一项(u_{r1} - u_r) * v_old[i] * v_old[j]都使用区间乘法求和使用区间加法。最终计算出的v_new的每个分量都是一个区间并且可以数学证明真实的不动点函数值位于该区间内。同样用区间算术计算积分∫ (1 - L(x)) dx和尾部余项最终得到的U_m是一个区间并且可以证明真实的E[S]一定小于等于这个区间的右端点。只有这样我们才能说“E[S] ≤ 4.09”是一个被严格证明的定理而不仅仅是一个数值观察。这方面的工具包括MPFI(基于MPFR的区间算术库)、Julia语言的IntervalArithmetic.jl包等。5. 方案总结与扩展应用思考回顾整个方案其核心贡献在于提供了一条系统化的、可计算的路径来获取复杂随机变量期望的严格上界。它将概率分析推导尾部上界和不动点方程、数值分析保守网格离散化和计算机科学算法实现与可验证计算紧密结合。核心优势严格性整个框架建立在严格的不等式基础上。只要初始条件L_0 ≤ F成立且离散化是保守的那么最终得到的上界在数学上就是成立的。可自动化整个过程可以编码实现一旦模型即S的递归定义和尾部上界形式确定计算上界的过程几乎是机械化的。可优化性通过加密网格、优化初始函数、调整截断点A可以系统地改进降低上界。局限性与挑战计算复杂度O(N*M)的复杂度对于高精度要求是个挑战。可能需要借助更高效的数值积分方法如高斯积分或自适应网格来优化。维度灾难本方法针对一维分布函数。对于高维随机向量其联合分布函数的逼近会变得异常复杂。模型依赖性方法的成功实施强烈依赖于能够推导出类似F KF的不动点方程以及一个形式良好的尾部概率上界。对于更复杂的递归关系推导可能非常困难。扩展应用 这套单调分布函数方案的思想并不局限于原文中的特定随机变量S。它可以推广到任何满足类似分布递归等式的随机变量例如其他快速选择算法变体分析QuickSelect在不同 pivot 选择策略下的比较次数。分支随机游走某些分支过程的生存时间或最大位移。递归算法复杂度满足特定分治递归的算法运行时间。极值统计某些关联模型中最大值的分布。在实际操作中最关键的一步是为你所研究的随机变量Y建立起它的分布函数F_Y(x)所满足的方程通常通过全概率公式或条件期望并找到一个合适的保守积分算子K。一旦做到了这一点剩下的网格迭代和上界计算就是一套相对固定的流程了。最后分享一个我个人的调试心得在实现这类保守迭代算法时可视化是极其强大的工具。就像原文中的图2和图3绘制出每次迭代的上界值U_m和函数曲线L_m(x)能让你直观地看到收敛过程快速发现参数设置如A,N是否合理以及算法实现中可能存在的错误比如曲线出现非单调性这一定意味着代码有 bug。将数学的严谨与计算的直观相结合是解决这类问题的有效法门。
