1. 阶乘不只是"从1乘到n"：把它拆成三层才看得懂
很多人第一次在机器学习教材里看到阶乘，反应都是"这不是高中排列组合的东西吗"。我刚开始啃这块的时候也这么想，直到在推导泰勒展开的余项、算泊松回归的对数似然、甚至排查一个二项式系数溢出成inf的 bug 时，才发现阶乘和阶数这两个概念是被自己严重低估的。阶乘在微积分与机器学习里的角色，远不是一个"算排列数"的工具，它更像一把标尺：衡量增长速度的标尺、给泰勒展开定分母的标尺、给数值计算划安全边界的标尺。而"阶数"这个词在数学语境里至少有四层不同的含义，混用会直接导致对公式的理解偏差。
这一篇想干的事情很具体：把阶乘从"n! = 1·2·3···n"这个定义里拽出来，讲清楚它在实数域上怎么被接续、增长有多快、在机器学习里到底出现在哪几个真实位置；再把"阶数"在导数、无穷小、多项式近似、计算复杂度、收敛速度这五个场景下的不同含义逐个拆开。适合两类人看：一类是正在补微积分基础、看公式时对n!出现得莫名其妙感到困惑的同学；另一类是已经能跑模型，但在写损失函数、做数值计算时被溢出和精度问题咬过的从业者。
提示：阶乘和阶数是两个独立概念，"阶乘"是一个函数，"阶数"是一个度量标签。后面所有容易误解的地方，几乎都源于把这两个词混为一谈。
1.1 0! = 1 不是拍脑袋的规定，而是唯一自洽的选择
教科书通常直接写"规定 0! = 1"，然后翻页。这种讲法让人记不住，因为它看起来像人为约定。实际上只要从阶乘的递推关系出发，0! 的值是被逼出来的，没得选。
阶乘最本质的性质是递推：n! = n × (n-1)!。这个等式对所有正整数成立。现在把 n = 1 代进去，得到1! = 1 × 0!。既然 1! 必须等于 1（只有一个元素的排列只有一种），那么 0! 只能等于 1。如果强行让 0! = 0，整个递推链条在 n=1 处就断了，组合数公式C(n,k) = n!/(k!(n-k)!)在 k=0 时也会除零崩掉——而"从 n 个里选 0 个"显然只有一种选法。
这个推导本身没多难，但它的意义在于：给了你一个判断"某个数学约定是否合理"的模板。凡是遇到看起来像硬性规定的东西，先去找它必须满足的结构性约束，往往能反推出唯一解。后面讲 Gamma 函数时，用的完全是同一个思路。
我自己带新人的时候会让他们做一个小练习：不用任何资料，只用递推关系和"排列数的物理意义"，把 0! 推出来。推得出来的人，后面理解 Gamma 函数的延拓会顺畅很多，因为逻辑是同一套。
1.2 Gamma 函数：把阶乘从整数点"接"到实轴上
阶乘只在非负整数上有定义，但微积分要求函数在连续区间上可导。这里就出现了一个典型矛盾：机器学习里做连续优化，你没法对一个只定义在整数点上的函数求导。解决办法是把阶乘延拓成 Gamma 函数。
Gamma 函数的积分定义是Γ(z) = ∫₀^∞ t^(z-1) e^(-t) dt，对Re(z) > 0收敛。它满足和阶乘一模一样的递推关系Γ(z+1) = z·Γ(z)，并且Γ(1) = 1。于是对正整数 n 有Γ(n+1) = n!。注意这个偏移一位的关系，是初学者最容易写错的地方——代码里调math.gamma(n)想算n!，必须写成math.gamma(n+1)，写成math.gamma(n)得到的是(n-1)!。
Gamma 函数真正的价值不仅是"能算小数阶乘"，而是它揭示了一个事实：阶乘的增长本质是指数级乘以幂级的混合体。这个性质在Γ(1/2) = √π这个特殊值上体现得最明显——半整数点的阶乘会出现 π，说明它和正态分布、和积分∫e^(-x²)dx是同源的。后面讲泊松分布、讲高斯积分的时候，这个 √π 会反复冒出来。
顺便说一个实操中的坑：Gamma 函数在负整数点有极点，math.gamma(-1)会直接抛异常，scipy.special.gamma(-1)返回inf或nan。有些数值库在不同版本间行为不一致，写代码时不要依赖异常来判负，自己先做入参校验更稳。
1.3 增长到底多快：Stirling 近似给出的实用边界
n!的增长速度是个必须建立直觉的东西，因为它是判断"这段代码会不会炸"的第一道防线。Stirling 近似给出n! ≈ √(2πn) · (n/e)^n，相对误差在 n ≥ 10 时已经小到 1% 以内。把它拆开看：(n/e)^n是指数级的底数在增长，√(2πn)是多项式修正。所以阶乘不是单纯的指数增长，它的底数本身随 n 增加，这比任何固定底数的指数函数都猛。
这里可以做一个粗糙的手算，帮助建立量级感：取 n = 20，(20/e)^20 ≈ 7.36^20，而7.36^10 ≈ 4.6×10^8，平方一下约2.1×10^17，再乘√(2π·20) ≈ 11.2，得到约2.4×10^18。而 20! 的真实值是2432902008176640000，量级完全对得上。
这个数量级意味着什么？64 位整数能表示到约9.2×10^18，也就是 20! 勉强塞得下，21! 就溢出了。用双精度浮点数（最大约1.8×10^308），大概能撑到 170!。超过这个范围，任何直接计算阶乘的代码都会返回inf，然后你的损失函数变成nan，训练直接崩。
| 计算方式 | 能表示的阶乘上限 | 溢出后的表现 | 典型使用场景 |
|---|---|---|---|
| 32 位整数 | 12! | 静默回绕成错误值 | 几乎不该用 |
| 64 位整数 | 20! | 静默回绕或抛异常 | 小规模组合计数 |
| 双精度浮点 | 约 170! | 变成inf | 一般科学计算 |
| 对数域（lgamma） | 理论上无上限 | 精度随数值增大缓慢下降 | 概率模型、似然计算 |
| 大整数（任意精度） | 受内存限制 | 不会溢出，但极慢 | 精确组合计数 |
我踩过的第一个坑就在这里：早期写一个组合数相关的评估指标，直接用math.factorial(50) // (math.factorial(25) * math.factorial(25))，结果算出来的数值完全不对——Python 的 int 是任意精度的，30! 以上虽然不会溢出，但这个方案在 n=1000 时慢得离谱，因为中间过程要算巨大的整数。后来改成对数域计算才恢复正常。
2. "阶数"这个词在微积分里有五张面孔
阶数是个被严重滥用的词。同一次讨论里，"二阶"可能指的是二阶导数，也可能指的是二阶无穷小，还可能指的是二阶优化算法，听众如果不在同一个语境下，沟通会彻底跑偏。我经历过一次评审会，有人讲"这个模型的阶数太高了"，结果三个人分别理解成了多项式次数太高、求导阶数太高、模型复杂度太高，争执了二十分钟才发现大家说的根本不是一件事。所以这一节把常见的几种含义列清楚，是省时间的做法。
2.1 导数的阶：二阶导才带"曲率"这个信息
一阶导给的是方向，二阶导给的是形状。这个区别在优化里是决定性的：梯度下降只看一阶导，所以它必须靠一个足够小的步长来避免冲过谷底；牛顿法看二阶导，它知道当前点处损失的弯曲程度，可以自动调整步长。
用一个具体的函数感受一下。取f(x) = x²，一阶导是2x，二阶导是常数2。在x = 10处，梯度是 20，如果你用学习率 0.1，更新量是 2，一下子从 10 跳到 8，再跳到 6.4……收敛得挺稳。但如果你用f(x) = x⁴，在x = 10处一阶导是 4000，学习率 0.1 会让参数跳到-390，直接发散。而二阶导在x = 10处是12x² = 1200，牛顿法的步长是f'/f'' = 4000/1200 ≈ 3.33，一步只走到 6.67 附近，非常稳。
这就是"阶"的实际含义：它决定了你能从函数曲面上读到多少局部信息，进而决定了步长该取多大。经验法则是我自己在调参时常用的：学习率的上界大约等于2/L，其中 L 是损失函数关于参数的二阶导（Hessian）的最大特征值。这个2/L不是玄学，把更新公式θ ← θ - η·f'(θ)在极小值附近做一次泰勒展开就能推出来，要求|1 - η·L| < 1，解出η < 2/L。多分类任务里 L 往往很大，这就是为什么学习率得设得那么小。
2.2 无穷小的阶：等价替换什么时候会让你算错
微积分里判断两个无穷小谁"更小"，用的是阶的概念：如果lim f(x)/g(x) = 0，就说 f 是比 g 高阶的无穷小；如果极限是非零常数，就说二者同阶。等价无穷小替换（比如x→0时sin x ~ x）就是基于这个判断。
这个工具的威力在于把复杂极限化简，但它有个明确的失效条件：在加减法中被替换的部分不能相互抵消。举个容易翻车的例子，求lim (x - sin x)/x³。如果你直接把sin x换成x，分子变成x - x = 0，得到答案 0。但正确答案是1/6。原因就是sin x和x的一阶项刚好抵消了，真正起作用的是三阶项。换句话说，你保留的"阶"不够高，把决定结果的那一项给扔掉了。
| 替换操作 | 是否安全 | 原因 |
|---|---|---|
| 乘除法中的因子替换 | 安全 | 极限的四则运算性质保证 |
| 加减法中各项不抵消 | 一般安全 | 主导项未被消去 |
| 加减法中最高阶项抵消 | 危险 | 需保留更高阶项 |
| 替换后出现 0/0 或 ∞/∞ | 需重新判断 | 阶数估计不足 |
这条规则在机器学习里也有对应物：做损失函数的近似简化时（比如把小量ε的高阶项丢掉），必须确认你丢掉的项不会和保留项同阶抵消。我见过有人在小样本修正公式里丢掉了 O(1/n²) 的项，结果修正量本身也是 O(1/n²)，符号反了，偏差反而变大。
2.3 多项式与泰勒近似的阶：阶数越高不一定越好
多项式的阶数就是最高次项的次数，这个定义很清楚。但"用几阶多项式去近似一个函数"这个问题，答案并不是"越高越好"。
理论上，泰勒展开的阶数越高，在展开点附近误差越小。但实际计算时有两个反向力量。第一是收敛半径限制：1/(1-x)的展开只在|x| < 1收敛，你在x = 2处展开到 100 阶也没用，误差反而爆炸。第二是数值误差累积：高阶项系数往往是阶乘分之一这种极小的数，在浮点运算里会损失有效位。
一个具体的对比：用不同阶数的多项式近似e^x在x = 1处的值。0 阶（常数 1）误差约 1.718；1 阶（1+x）误差约 0.718；2 阶误差约 0.218；5 阶误差降到约 0.0016；10 阶误差已经到 1e-8 量级。看起来是单调变好，但如果你把x换成 10，10 阶泰勒多项式的值会严重偏离真实值e^10 ≈ 22026，因为 10 阶展开在x = 10处远未收敛。
我的实际经验是：在需要全局近似时，优先用切比雪夫多项式或者直接上分段近似；泰勒展开只适合局部、小范围的线性化分析。深度学习里的梯度推导、误差分析基本都属于后者，所以泰勒展开够用；但如果有人拿泰勒展开去做大范围的函数拟合，那基本是在浪费时间。
2.4 复杂度的阶：O(n!) 是一条不可跨越的红线
算法复杂度里的"阶"是另一套体系，讲的是输入规模趋向无穷时运算量的增长趋势。这一套和微积分里的阶在思想上完全一致——都是比较增长速度——但应用场景完全不同。
工程上值得记住的几条分界线：O(n log n)是可接受的上限（排序、大部分树结构操作）；O(n²)在 n 上千时开始难受（朴素矩阵乘法、两两相似度计算）；O(2^n)只能撑到 n 二三十（穷举子集）；O(n!)则连 n = 15 都撑不住（15! 约 1.3×10^12，按每秒十亿次运算也要跑二十多分钟）。
这个数字很有用：当你看到一个问题需要枚举所有排列时，先算一下 n! 的数量级，再决定要不要换思路。旅行商问题、任务排序、排班调度这类问题表面上看就是"找个最优排列"，但 n 稍大就必须转向启发式或者动态规划。我在做一个小规模的调度功能时，最初方案是枚举所有排列找最优，n = 9 的时候跑得挺欢（36 万种），上线后遇到 n = 12 的输入就直接超时（4.8 亿种）。
2.5 收敛阶：判断一个迭代算法好不好用的硬指标
收敛阶的定义是：如果迭代序列满足|x_{k+1} - x*| / |x_k - x*|^p → C（C 为非零常数），就说这个算法是 p 阶收敛的。p = 1 叫线性收敛，p = 2 叫二次收敛。
这个指标为什么重要？因为线性收敛意味着每步误差缩小固定倍数，要达到 1e-10 的精度可能需要几十步；而二次收敛意味着每步有效位数翻倍，五六步就能搞定。牛顿法在极值点附近是二次收敛的，这就是它在理论上远优于梯度下降的原因。
| 算法 | 收敛阶 | 每步代价 | 实际表现 |
|---|---|---|---|
| 梯度下降 | 线性 | 一阶梯度 | 稳定但慢，需要精细调学习率 |
| 动量法 | 线性（更快） | 一阶梯度 + 历史 | 在狭长谷底明显加速 |
| 牛顿法 | 二次 | 二阶 Hessian + 求逆 | 极快但单步代价高，Hessian 可能非正定 |
| 拟牛顿（BFGS 类） | 超线性 | 一阶梯度 + 近似 Hessian | 折中方案，中等规模问题常用 |
| 一阶矩自适应方法 | 线性 | 一阶梯度 + 两个指数滑动平均 | 鲁棒性好，超参不敏感 |
不过要注意：收敛阶是局部性质，它只告诉你"在离最优解足够近的地方"收敛速度如何，不保证从任意初值出发都能收敛。很多教程把牛顿法的二次收敛吹得很神，但实际用起来最头疼的是初值选择——离得远的时候牛顿法可能直接跳到更差的点。这也是为什么实际工程里大家都用一阶方法加自适应，稳字当头。
3. 泰勒展开里那个 n!，究竟是被谁"塞"进去的
几乎所有人在第一次看到e^x = 1 + x + x²/2! + x³/3! + ...的时候，都会下意识问一句：分母上为什么是阶乘？为什么不是 n，不是 n²，偏偏是 n!？教材往往一句"可以证明"就跳过，但这个"为什么"其实非常值得展开，因为搞懂它之后，你对"阶乘是一种标尺"这件事会有切身体会。
3.1 用待定系数法反推，n! 是被逼出来的唯一解
假设我们想用一个多项式P(x) = a₀ + a₁x + a₂x² + a₃x³ + ...在x = 0附近逼近f(x)，并且要求"在 0 点处，P 的 k 阶导数等于 f 的 k 阶导数"，对所有 k 都成立。
先看x = 0：P(0) = a₀，所以a₀ = f(0)。
求一阶导：P'(x) = a₁ + 2a₂x + 3a₃x² + ...，代入x = 0得P'(0) = a₁，所以a₁ = f'(0)。
求二阶导：P''(x) = 2a₂ + 6a₃x + 12a₄x² + ...，代入x = 0得P''(0) = 2a₂，所以a₂ = f''(0)/2。
求三阶导：P'''(0) = 6a₃，所以a₃ = f'''(0)/6。
看出来了吗？每求一次导，x^k这一项就会掉下来一个系数，连续掉 k 次，掉下来的正是k×(k-1)×...×1 = k!。所以a_k = f^(k)(0) / k!。
阶乘在里面扮演的角色是"导数掉下来的那个倍数"。这就是它的全部秘密。你希望多项式的 k 阶导在原点匹配目标函数，那x^k项就必须先除掉它自带的那 k! 倍系数，否则匹配不上。
这个推导的价值不只是"知道公式怎么来的"。它给了一个通用的方法论：凡是看到公式里冒出阶乘，先去找那个"连续下降 k 次"的过程。排列组合里的 n! 是"每次选一个、可选数递减 1"；泰勒展开里的 n! 是"求导 k 次、每次降一阶"；Gamma 函数里的 n! 是"递推 n 次、每次乘一个数"。本质上是同一件事。
3.2 阶乘决定了展开的"信用额度"：余项与收敛半径
泰勒展开在有限阶截断时，会留下一个余项。拉格朗日余项的形式是R_n(x) = f^(n+1)(ξ) · x^(n+1) / (n+1)!，其中 ξ 介于 0 和 x 之间。
这个式子里阶乘在分母上，这意味着随着 n 增大，分母爆炸式增长，余项被强力压制。这是泰勒展开好用的根本原因：阶乘的分母保证了高阶项的贡献迅速衰减。反过来说，如果你的函数的高阶导数增长得比阶乘还快，展开就不收敛。
具体算一下e^x在x = 1处的余项。取 n = 5，R_5 ≈ e^ξ · 1^6/6! ≤ e/720 ≈ 0.0038。实测 5 阶近似值是 2.71667，真值 2.71828，误差 0.0016，在估计范围内。取 n = 10，1/11! ≈ 2.5×10^-8，误差直接掉到 1e-8 量级。这就是"每加一阶，误差除以 n"的效果。
但同样的展开在x = 10处就废了。10^11/11! ≈ 10^11/4×10^7 ≈ 2500，余项估计都大于函数值本身了。所以这里有个非常实用的结论：泰勒展开的可用范围由 x 的大小和阶数共同决定，粗略判断标准是|x| < n/2左右，超过这个范围就该换方法。
3.3 机器学习里最常用的三个展开，以及它们各自在哪一步用到阶乘
第一个是e^x的展开。这是 logistic 回归和 softmax 推导的起点。sigmoid 函数σ(x) = 1/(1+e^(-x))在x = 0附近展开，一阶项给了1/2 + x/4，这是"为什么初始化时输出接近 0.5"的数学解释。softmax 在 logits 差异不大时近似均匀分布，也是同一个来源。
第二个是log(1+x)的展开：x - x²/2 + x³/3 - ...，注意这里分母是 n 而不是 n!。这个区别很重要——对数函数的高阶导数衰减得比指数函数慢，所以它的展开系数只有 1/n 而不是 1/n!。实际后果是log(1+x)在|x|接近 1 时收敛极慢，用它做近似要非常谨慎。交叉熵损失里经常出现log(1-p)这类项，p 接近 1 时数值非常不稳定，根源就在这里。
第三个是(1+x)^α的广义二项式展开，系数是C(α, k) = α(α-1)...(α-k+1)/k!，阶乘又出现在分母上。这个展开在推导 focal loss、解释样本重加权机制时会用到。
import math def taylor_exp(x, n): """n 阶泰勒多项式近似 e^x""" return sum(x**k / math.factorial(k) for k in range(n + 1)) def taylor_log1p(x, n): """n 阶泰勒多项式近似 log(1+x)，注意分母是 k 不是 k!""" return sum((-1)**(k+1) * x**k / k for k in range(1, n + 1)) # 实测对比 for n in [3, 5, 8, 12]: approx = taylor_exp(1.0, n) print(f"n={n:2d} exp近似={approx:.10f} 误差={abs(approx - math.e):.2e}") # log1p 在接近收敛边界时的表现 for x in [0.1, 0.5, 0.9, 0.99]: approx = taylor_log1p(x, 50) print(f"x={x:.2f} 50阶近似={approx:.6f} 真值={math.log(1+x):.6f}")跑这段代码会发现，exp的 8 阶展开误差已经在 1e-6 以下，而log1p在x = 0.99时即使展开 50 阶，误差仍然肉眼可见。这个对比很能说明问题：分母上有没有阶乘，决定了这个展开"好不好用"。
4. 阶乘在机器学习代码里真正出现的位置
讲完理论，回到工程。阶乘在机器学习代码里出现的频率其实不低，只是常常藏在你看不到的地方——要么被约掉了，要么被移到对数域里了。搞清楚它出现在哪、为什么常常"消失"，能省下大量排查nan的时间。
4.1 交叉熵里的阶乘：它确实存在，但对梯度没贡献
先看一个容易忽略的事实：softmax 加交叉熵，严格来说是从多项式分布的最大似然推出来的。多项式分布的概率质量函数是P(x₁,...,x_k) = n!/(x₁!···x_k!) · p₁^x₁···p_k^x_k。取对数之后，log(n!/(x₁!···x_k!))这一项里全是阶乘。
然后对参数求梯度的时候，这一项就消失了——因为它只和观测到的计数有关，跟模型参数p无关，导数为零。这就是为什么你在几乎所有教程里看到的交叉熵都是-Σ y log(ŷ)，阶乘项直接不见踪影。
这个"常数项可以丢"的操作本身没问题，但有两个实践上的注意点。第一，如果你在实现一个需要比较不同样本的似然值的功能（比如异常检测、样本加权），常数项不能随便丢，因为它和样本有关，会影响排序。第二，如果你在算对数似然的绝对值（比如做模型选择、算 BIC/AIC），必须把常数项补回去，否则数值完全对不上。
我做样本重加权时踩过一次：损失函数里按惯例省掉了组合数项，结果重加权后的样本间相对权重算错了，因为不同样本的log(n!/(x₁!...x_k!))不一样。改成用math.lgamma把这一项显式加回去之后才对上。
4.2 泊松回归与计数建模：阶乘直接坐在分母上
如果任务本身是预测"某段时间内发生了多少次事件"（点击次数、调用次数、故障次数），常见的建模选择是泊松回归。它的概率质量函数是P(k|λ) = λ^k · e^(-λ) / k!。
取对数：log P = k·log λ - λ - log(k!)。前两项对 λ 求导很好处理，第三项log(k!)是常数，对参数更新没影响。但它在计算单个样本的对数似然值时是必须的，而且 k 稍微大一点（比如 k = 200），k!就已经超出双精度浮点的表示范围。
正确写法是用lgamma(k+1)代替log(factorial(k))。这两个在数学上等价，但前者在 k 很大时依然稳定，后者在 k > 170 时直接返回inf，log(inf) = inf，整个似然变成负无穷。
import math def poisson_logpdf_bad(k, lam): # 不要这样写：k > 170 时 factorial 溢出 return k * math.log(lam) - lam - math.log(math.factorial(k)) def poisson_logpdf_good(k, lam): # 正确写法：全程在对数域 return k * math.log(lam) - lam - math.lgamma(k + 1) for k in [5, 50, 200, 1000]: try: bad = poisson_logpdf_bad(k, 10.0) except OverflowError: bad = "溢出" good = poisson_logpdf_good(k, 10.0) print(f"k={k:5d} bad={bad} good={good:.4f}")这段跑下来，k = 200 时上面的写法会直接抛OverflowError（因为math.factorial返回 Python 大整数，转 float 时才溢出，行为依实现而定），下面的写法老老实实输出一个有限的负数。这个对比值得每个人自己跑一遍，印象会非常深。
4.3 组合数爆炸：采样与搜索里最容易被低估的成本
排列数和组合数都和阶乘直接相关：P(n,k) = n!/(n-k)!，C(n,k) = n!/(k!(n-k)!)。这两个公式在采样、搜索、评估指标里大量出现。
一个具体场景：评估排序质量时，如果要穷举所有可能的排序来计算某个指标的期望值，n 个物品就有 n! 种排序。n = 10 时是 362 万，还能跑；n = 13 时是 62 亿，跑不动了。这时候要么换成采样式估计，要么用解析解。
另一个场景：负采样。如果要从 n 个候选中选 k 个，组合数是C(n,k)。当 n = 10000、k = 5 时，组合数约8.3×10^17，虽然算得出来，但如果你想遍历所有组合做精确评估，同样不现实。
| 场景 | 涉及的阶乘规模 | n 的可用上限 | 替代方案 |
|---|---|---|---|
| 全排列枚举 | n! | 约 11-12 | 随机采样、beam search |
| 全组合枚举 | C(n,k) | 视 k 而定，通常 n < 25 | 分层采样、重要性采样 |
| 精确排列数计算 | n! | 任意（用大整数或对数） | lgamma或任意精度整数 |
| 二项分布概率 | C(n,k)p^k(1-p)^(n-k) | 任意（对数域） | logaddexp组合 |
| 多项分布采样 | 多项式系数 | 任意 | 逐项采样，避免算系数 |
注意：真正需要"精确计算"组合数本身（例如输出给用户看）和"在概率计算中用到"组合数，是两件不同的事。前者要精确整数，后者必须在对数域做。混用会同时踩精度和性能两个坑。
5. 数值实现：三条让阶乘不炸的习惯
理论讲完了，这一节全是实操。我把自己在代码里反复用的几条规则整理出来，每一条都对应一个真实踩过的坑。
5.1 任何进入概率计算的阶乘，一律换成 lgamma
这条是硬性规矩。理由很简单：log(n!)几乎总是出现在概率的对数里，而n!本身从来不需要单独算。既然目标是取对数，那就一步到位用lgamma(n+1)，跳过中间那个巨大的数。
Python 里有两个入口：math.lgamma(x)返回log(Γ(x))，scipy.special.gammaln(x)行为一致但支持数组。做批量计算时用后者，向量化能快出几十倍。
需要提醒的是lgamma的入参是n+1而不是n。这个偏移我见过太多人写错，而且错了之后结果不会报错，只是数值完全不对，非常隐蔽。建议在项目里封装一个语义明确的函数：
from scipy.special import gammaln import numpy as np def log_factorial(n): """返回 log(n!)，接受标量或 numpy 数组""" return gammaln(np.asarray(n, dtype=np.float64) + 1.0) # 验证 for n in [0, 1, 5, 20, 100, 1000]: direct = np.log(float(np.math.factorial(int(n)))) if n <= 170 else None lg = log_factorial(n) print(f"n={n:5d} log(n!)={lg:.6f}" + (f" 直接法={direct:.6f}" if direct else " 直接法=溢出"))第三个坑来自负数和零的边界。gammaln(0)返回inf，gammaln(-1)返回inf，不会抛异常。所以在做累加的时候，一个inf会把整批结果污染掉。入参校验必须做在前面。
5.2 二项式系数：用对数差，不要用阶乘相除
C(n,k) = n!/(k!(n-k)!)，如果直接用三个阶乘算，n = 100 就炸了。改写成对数域：log C(n,k) = lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1)，然后按需取指数。
但取指数这一步也有讲究：如果 k 和 n 都很大，log C(n,k)可能达到几百，exp回来直接溢出。所以更好的做法是不要在中间取指数，让整个计算链都留在对数域，最后如果需要概率值再统一做归一化（用 log-sum-exp 技巧）。
from scipy.special import gammaln import numpy as np def log_binom(n, k): return gammaln(n + 1) - gammaln(k + 1) - gammaln(n - k + 1) def stable_binom_pmf(k, n, p): """二项分布 pmf 的稳定写法，全程对数域""" log_pmf = log_binom(n, k) + k * np.log(p) + (n - k) * np.log1p(-p) return np.exp(log_pmf) # 大 n 场景：直接算法会挂，对数域算法稳 n, k, p = 5000, 2500, 0.5 print(f"log C({n},{k}) = {log_binom(n, k):.4f}") print(f"pmf = {stable_binom_pmf(k, n, p):.6e}")这里用np.log1p(-p)而不是np.log(1-p)。当 p 很小时，1-p在浮点里会损失精度，log1p专门解决这个问题。这是数值计算里的标准操作，凡是log(1+x)形式都该用log1p。同理还有expm1处理exp(x)-1。
5.3 Stirling 近似什么时候可以用，什么时候不能用
Stirling 近似在算法分析、复杂度估计、大样本理论推导里用得很广，但它是渐近近似，在小 n 时误差巨大。n = 1 时 Stirling 给出√(2π)·(1/e) ≈ 0.922，而真值是 1，误差 8%；n = 2 时误差还有 4% 左右。
实际操作中我通常按这个标准判断：
| n 范围 | 是否可用 Stirling | 建议做法 |
|---|---|---|
| n < 5 | 不可用 | 直接算 |
| 5 ≤ n < 20 | 谨慎，误差 1%-3% | 用修正版√(2πn)(n/e)^n · (1 + 1/(12n)) |
| n ≥ 20 | 可用，误差 < 0.5% | 标准 Stirling 足够 |
| n ≥ 100 | 完全可用 | 对数形式n log n - n + 0.5 log(2πn) |
带一阶修正的版本精度提升很明显：n = 5 时误差从 1.6% 降到 0.03%。这个修正项1 + 1/(12n)的推导来自对 log 版的进一步展开，实际用的时候直接用结论就行。
真正需要用到 Stirling 的场合，往往是推导而不是计算。比如分析某个采样算法的复杂度、估计大 n 下二项分布的渐近行为、推导信息准则里的惩罚项，这些场景里你需要的是一个可解析的表达式，而不是一个精确数值。
6. 从"阶"看优化器：一阶和二阶到底差在哪
写到这里，"阶乘"这条线基本讲完了。最后想顺着"阶数"再聊一层，因为它直接关系到你选优化器时的判断依据，也是我在实际项目里被问得最多的问题之一。
6.1 梯度下降的"阶数盲区"
一阶方法只知道当前点往哪个方向下降最快，但不知道这个方向能走多远。想象你站在一个山谷里，闭着眼睛只能感受脚下的坡度——坡度陡的地方你不敢迈大步，坡度缓的地方你也不知道前面是不是悬崖。这就是一阶方法的处境。
数学上，一阶方法隐含的假设是"函数在局部近似线性"。如果函数曲率很大（二阶导大），这个线性近似只在极小的邻域内成立，步长必须很小。这就是为什么深度学习的训练需要那么多学习率调参技巧——本质上是在用各种手段弥补"看不到曲率"这个缺陷。
一个非常具体的数字感受：考虑f(x) = 100x²，在x = 1处梯度是 200。学习率 0.001 时，一步走到0.8，还行；学习率 0.01 时，一步走到-1，来回震荡；学习率 0.011 时，直接发散。而这里的临界值2/L = 2/200 = 0.01，和前面的理论完全对上。这个2/L的推导值得每个人自己动手算一遍，比背结论有用得多。
6.2 二阶信息的获取成本与替代方案
牛顿法直接用 Hessian 矩阵，参数量是 P 的时候 Hessian 有 P² 个元素，求逆是 O(P³)。P = 100 万时完全不可行。所以实际方案全是"用便宜的方式近似二阶信息"。
第一类是拟牛顿法，用历史梯度差去构造 Hessian 的近似，存储是 O(P²)，适合中等规模。第二类是共轭梯度，只存方向向量，但需要多次函数求值。第三类也是深度学习里最主流的，是用一阶矩和二阶矩的指数滑动平均来模拟二阶信息的效应。这一类方法对每个参数单独维护两个统计量，代价是 O(P)，非常便宜。
| 方法族 | 需要的二阶信息 | 单步内存 | 单步计算 | 适用规模 |
|---|---|---|---|---|
| 梯度下降 | 无 | O(P) | O(P) | 任意 |
| 动量法 | 无（历史一阶） | O(P) | O(P) | 任意 |
| 双矩自适应 | 二阶矩的滑动估计 | O(P) | O(P) | 任意 |
| 拟牛顿类 | Hessian 近似 | O(P²) | O(P²) | P < 10⁴ |
| 牛顿法 | 完整 Hessian | O(P²) | O(P³) | P < 10³ |
这张表里最值得琢磨的一列是"单步内存"。对于一个 10 亿参数的模型，O(P²) 意味着 10^18 个浮点数，物理上不可能。所以规模一大，二阶方法必然退场，这不是精度问题，是存储和算力的硬约束。
6.3 曲率感知在实践中的三个落地方式
既然不能显式计算二阶导，那就想办法"间接感知"曲率。我总结下来常用三种。
第一种是自适应学习率。参数更新幅度在梯度的历史尺度上做归一化，梯度长期很大的参数步子收小，梯度长期很小的参数步子放大。这本质上是在用梯度的历史方差去估计局部的 L 值。我在一个梯度尺度差异极大的多任务场景里对比过，用固定的全局学习率怎么调都有任务被牺牲，换到自适应之后才平衡。
第二种是梯度裁剪。这招看起来很粗暴——直接把梯度范数截到一个阈值——但它确实在抑制"曲率突然变大导致的爆炸"上非常有效。判定阈值的一个实用方法是：先跑一小段不加裁剪的训练，记录梯度范数的分位数，把阈值设在 95 分位到 99 分位之间。
第三种是对参数做小范围的二阶差分探测。这个做法比较小众但很有用：在某个参数上做一次小扰动ε，观察损失变化，用(f(x+ε) - 2f(x) + f(x-ε))/ε²估计该方向的二阶导。代价是三次前向，适合诊断阶段用，不适合每步都算。
import numpy as np def estimate_curvature(loss_fn, params, idx, eps=1e-4): """沿第 idx 个参数方向做二阶差分，估计局部曲率""" p_plus = params.copy(); p_plus[idx] += eps p_minus = params.copy(); p_minus[idx] -= eps f0 = loss_fn(params) f_plus = loss_fn(p_plus) f_minus = loss_fn(p_minus) return (f_plus - 2 * f0 + f_minus) / (eps ** 2) # 用一个简单二次函数验证：f(x) = 3x^2，理论上所有方向曲率都是 6 def quad_loss(x): return float(3 * np.sum(x ** 2)) x = np.array([1.0, 2.0, 3.0]) for i in range(3): c = estimate_curvature(quad_loss, x, i) print(f"参数 {i} 方向的曲率估计: {c:.4f}")这段代码跑出来三个方向都接近 6，验证了估计的准确性。eps的选取是这里的关键：太大时高阶项污染结果，太小时浮点误差淹没差分。经验值是eps ≈ 1e-4（对已经归一化到 O(1) 的参数），如果参数尺度差异大，应该按参数自身的尺度取相对值。
我在实际排查"某个层梯度一直爆炸"的问题时用过这个手段，发现是那一层的输入尺度比预期大了一个数量级，导致该方向的曲率特别大。后来在前面加了一层标准化就解决了。这个诊断过程比盲目调学习率高效得多，因为你能直接看到"哪个方向特别陡"。
6.4 一个容易被忽略的细节：阶与精度的关系
最后提一个我踩过的坑。做二阶差分估计曲率的时候，如果损失函数本身量级很大（比如几万），而eps太小，f(x+ε) - f(x)这个差值会落在浮点数的有效位之外，直接变成 0。这时候算出来的曲率是 0，你会误以为那个方向很平。
判断方法很简单：算完之后检查f_plus - f_minus是否明显大于该数量级的浮点精度（约1e-16 × |f|）。如果两者接近，说明eps需要放大。这个检查我建议做成一个断言，写在调试代码里，比事后追查便宜得多。
| 现象 | 可能原因 | 处理方式 |
|---|---|---|
| 曲率估计恒为 0 | eps 太小，差值被浮点误差吞掉 | 放大 eps 到差值与 1e-12· |
| 曲率估计剧烈跳变 | eps 太大，高阶项污染 | 缩小 eps，或改用中心差分 |
| 曲率估计为负 | 该方向非凸，或 eps 过大 | 先用小 eps 复测确认 |
| 估计值与理论差一个数量级 | 参数尺度未归一化 | 按参数尺度取相对 eps |
这些东西教科书上基本不写，但在真实调参现场天天遇到。我个人在实际操作中的体会是：数值计算的坑，八成不是算法错了，而是尺度没对齐、边界没处理、精度没检查。阶乘这个看起来最简单的数学对象，恰恰是这些坑最密集的地方——因为它增长得太快，快到任何一点疏漏都会被指数级放大。
如果这条线还想继续往下走，有两个自然的延伸方向。一是把阶乘和 Stirling 近似用到大样本统计推断里，看信息准则那些惩罚项是怎么从似然里长出来的；二是把组合爆炸的问题系统地过一遍，看看在什么规模下该从精确算法切换到近似算法，那个切换点在哪里。这两个方向都比单纯刷题有意思得多，也更接近真实工程里需要做的判断。
