玻尔兹曼分布详细推导、softmax 及 Energy Based Model

因为想申请 优快云 博客认证需要一定的粉丝量，而我写了五年博客才 700 多粉丝，本文开启关注才可阅读全文，很抱歉影响您的阅读体验

• 大多数人应该都是在利用神经网络解分类问题时第一次接触 softmax，对它的印象也只停留在把实数域的输出转换为合法的概率分布这里，本文介绍 softmax 的本源 —— 玻尔兹曼分布
• 参考
  1. 玻尔兹曼分布与Softmax的关联
  2. 近代物理导论-4.3麦克斯韦-玻尔兹曼分布

---

1. 理解玻尔兹曼分布

• 波尔兹曼分布是热力学与统计物理中的一个重要分布，它建立起 “宏观温度” 和 “微观粒子能量” 间的关系。具体而言，考虑一组理想状态下的粒子，要求

  1. 一个理想气体系统包含了巨大数量的以不同速度随机在各个方向上运动的全同气体分子
  2. 与分子间距相比，分子直径大小可以忽略不计
  3. 分子间除了相互间的完全弹性碰撞外别无其他相互作用
  4. 单个分子的运动遵循经典力学定律
  5. 在没有外势的情况下，当理想气体处在 平衡态 时气体分子均匀分布在系统中
  6. 在 平衡态，分子的速度在各方向上均匀分布，即理想气体的分子速度是随机的，并不存在一个特别偏好的方向

  将满足上述状态的理想气体（或理想状态的一组粒子）加热到一定温度（向系统内输入能量），各个微观粒子会以一定的速度进行随机运动和碰撞，在禁止一切能量交换和损失的情况下，各个粒子最终会稳定在一系列从低到高的不同能级上，即达到 平衡态，此时系统的熵最大

  除了混乱程度外，熵也可以表示有多少种到达此状态的方法。系统总能量是一定的，一些粒子占有的能量高，其他粒子占有的就少

  1. 只有少量的方法能将能量只分配给少量粒子（而且由于碰撞作用这些能量也会很快被传递给其他粒子）
  2. 只有少量方法能让所有粒子能量一致（同样很快会由于碰撞导致能量转移，从而离开该状态）

  因此最高熵时会呈现中等能级粒子多，高/低能级粒子少的状态，这就是 平衡态

  将各个能级的粒子数除以系统种粒子总数就得到一个分布形式，平衡态时对应的就是 玻尔兹曼分布/吉布斯分布（热力学中称最概然分布），这个分布是会随着系统温度变化的，形如下图所示
  
  注意几个点

  1. 理想气体几乎没有相互作用，那么它的能量（准确来说是哈密顿量）只有动能，因此可以用粒子速度表示，这里横轴是速度，其实是粒子能量
  2. 曲线越红，系统温度越高，粒子能达到的越高能量越高，玻尔兹曼分布曲线拉得越平，越接近均匀分布，系统熵越大

2. 推导玻尔兹曼分布

2.1 问题形式化

• 首先考虑粒子系统的表示形式，我们假设所有粒子只做平动不做转动，这样可以用xyz三个轴的位置 $r_x,r_y,r_z$ 和动量 $p_x,p_y,p_z$ 来刻画一个粒子，这六个维度 $(PP,rr)=(p_x,p_y,p_z;r_x,r_y,r_z)$ 构成了 μ相空间，每个粒子都可以表示为其中的一个点，这样一个 $N$ 粒子系统在某个时刻的微观状态就可以表示为 ( { P i } , { r i } ) (\{\pmb{P}_i\},\{\pmb{r}_i\}) ({PPi​},{rri​})

• 由于粒子太小太多，动量又是连续的物理量，我们不能分辨所有可能的状态，所以我们将相空间进行离散划分，只统计一个能级（一定能量区间）和一定空间区间内粒子的数目，同一个区间内的粒子我们认为是同质化的
  

  这个跟GDP的统计一样，我们不考虑张三年薪多少、李四年薪多少，只计数年薪20万上下的人数有多少、年薪30万上下的人数有多少。这里也可以看出，年薪是连续变化的，我们统计的时候要切分区间，这就是粗粒化的思想

  1. 宏观来讲，此相空间体积元足够的小以至于可以在操作中认为是 “无穷小”
  2. 微观来看，此相空间体积元又必须大到能足够装下数量巨大的粒子
  3. 在第 $i$ 个格子中的粒子具有共同的位置 ${rr}_i$、动量 ${PP}_i$，及能量 ${ϵ}_i={PP}_i^2/2m$（只考虑动能，包含三个维度，$m$ 是粒子质量）

• 这时我们要求的目标玻尔兹曼分布（本来是连续的概率密度函数，离散化后称分布函数）$f({PP}_i,{rr}_i)$ 满足
  $$f({PP}_i,{rr}_i)△\tau =n_i$$ 其中 $△\tau =△p_x△p_y△p_z△r_x△r_y△r_z$ 是体积元大小，$n_i$ 是此格子包含的粒子数，$f({PP}_i,{rr}_i)$ 就表示粒子落入此格子的概率

• 此外还有两个约束：对于孤立系统，粒子数守恒且能量守恒，即
  $$\sum _i{n}_i=N\sum _i{ϵ}_i=E$$ 对于没有外势的理想气体，在热力学极限下（$N\to \infty$），需要满足（六重积分）
  $$\int f({PP}_i,{rr}_i)d^{3}p{d}^{3}r=N\int f({PP}_i,{rr}_i)\frac{{PP}^2}{2m}{d}^{3}p{d}^{3}r=E$$

2.2 宏观状态与微观状态

• 在六维的 μ相空间 种，我们认为同一个体积元内的所有粒子是同质化的，在排列组合的语境下，我们认为这是一个无顺序的组合问题，我们不关心一个体积元内部粒子的微观顺序（微观状态），只关心宏观上这个体积元内的粒子数（宏观状态）
  
  从这个角度看，平衡态 = 系统最大熵状态 = 到达方法最多的宏观状态 = 对应微观状态最多的宏观状态 = 最概然分布
• 等几率假设：为了能通过除以总粒子数把一个宏观状态变形为概率分布，要求在一个能量恒定的孤立理想气体系统中，所有的微观状态出现的几率是相等的，这就是等几率假设
• 最后考察宏观状态对应的微观状态数，假设 $N,K$ 分别表示系统中的总粒子数和相空间中的体积元总数，宏观状态（分布）
  ( n 1 , n 2 , . . . , n K ) ≡ { n i } s . t . ∑ i K n i = N (n_1,n_2,...,n_K) \equiv \{n_i\} \quad s.t. \quad \sum_i^K n_i = N (n1​,n2​,...,nK​)≡{ni​}s.t.i∑K​ni​=N 对应的微观状态数是一个多重集的全排列问题
  Ω ( { n i } ) = C N n 1 C N − n 1 n 2 . . . C N − n 1 − n 2 − . . . − n K − 1 n K = N ! n 1 ! n 2 ! . . . n K ! \Omega(\{n_i\}) =C_{N}^{n_1}C_{N-n_1}^{n_2}...C_{N-n_1-n_2-...-n_{K-1}}^{n_K} = \frac{N!}{n_1!n_2!...n_K!} Ω({ni​})=CNn1​​CN−n1​n2​​...CN−n1​−n2​−...−nK−1​nK​​=n1​!n2​!...nK​!N!​

2.3 推导玻尔兹曼分布

• 根据上述分析，玻尔兹曼分布可以建模为以下约束优化问题的解：给定系统总能量 $E$ 和总粒子数 $N$，相空间划分体积元总数为 $K$
  max ⁡ { n i } ln ⁡ Ω ( { n i } ) s.t. ∑ i K n i = N ∑ i K ϵ i n i = E \begin{aligned} &\max_{\{n_i\}} &&\ln\Omega(\{n_i\}) \\ & \text{s.t.} &&\sum_i^K n_i = N \\ &&&\sum_i^K \epsilon_i n_i= E \end{aligned} ​{ni​}max​s.t.​​lnΩ({ni​})i∑K​ni​=Ni∑K​ϵi​ni​=E​ 其中 max ⁡ { n i } ln ⁡ Ω ( { n i } ) \max_{\{n_i\}} \ln\Omega(\{n_i\}) max{ni​}​lnΩ({ni​}) 等价于 max ⁡ { n i } Ω ( { n i } ) \max_{\{n_i\}}\Omega(\{n_i\}) max{ni​}​Ω({ni​})，通过取对数把连乘变成连加，从而简化计算
• 这是一个等式约束优化问题，用拉格朗日乘子法求解即可

2.3.1 斯特林（Stirling）公式

• 进行求解前先补充一个用于化简的数学工具，斯特林公式
  $$\ln n!\approx n\ln n-n$$ 这个公式的简单推导过程可以用化累加为积分来理解
  $$\begin{array}{rl}\ln n!& =\sum _{x=1}^n\ln x\\ & \approx {\int }_1^n\ln xdx\\ & =n\ln n-n+1\\ & \approx n\ln n-n\end{array}$$

2.3.2 用拉格朗日乘子法求解

• 引入拉格朗日乘子 $\alpha ,\beta$，构造拉格朗日函数为
  L = ln ⁡ Ω ( { n i } ) − α ∑ i K n i − β ∑ i K ϵ i n i = ln ⁡ N ! − ∑ i K ln ⁡ n i ! − α ∑ i K n i − β ∑ i K ϵ i n i ≈ N ln ⁡ N − N − ∑ i K n i ln ⁡ n i + ∑ i K n i − α ∑ i K n i − β ∑ i K ϵ i n i = N ln ⁡ N − N − ∑ i K [ n i ln ⁡ n i + ( α − 1 ) n i + β ϵ i n i ] \begin{aligned} \mathcal{L} &= \ln \Omega(\{n_i\}) - \alpha\sum_i^K n_i - \beta \sum_i^K \epsilon_i n_i \\ &= \ln N! - \sum_i^K \ln n_i ! -\alpha\sum_i^K n_i - \beta \sum_i^K \epsilon_i n_i \\ &\approx N\ln N - N - \sum_i^K n_i\ln n_i +\sum_i^K n_i -\alpha\sum_i^K n_i - \beta \sum_i^K \epsilon_i n_i \\ & = N\ln N - N - \sum_i^K\big[n_i\ln n_i+(\alpha-1)n_i +\beta \epsilon_in_i\big] \end{aligned} L​=lnΩ({ni​})−αi∑K​ni​−βi∑K​ϵi​ni​=lnN!−i∑K​lnni​!−αi∑K​ni​−βi∑K​ϵi​ni​≈NlnN−N−i∑K​ni​lnni​+i∑K​ni​−αi∑K​ni​−βi∑K​ϵi​ni​=NlnN−N−i∑K​[ni​lnni​+(α−1)ni​+βϵi​ni​]​

• 然后令偏导为零 $\frac{\partial \mathcal{L}}{\partial n_i}=\ln n_i+\alpha +\beta {ϵ}_i=0$ 解出
  $$n_i=e^{-\alpha }{e}^{-\beta {ϵ}_i}\text{\hspace{1em}(1)}$$ 这时可以直接和两个约束联立解出 $\alpha ,\beta$，但是那样比较复杂。我们可以如下计算
  $$N=\sum _i^K{n}_i=e^{-\alpha }\sum _i^K{e}^{-\beta {ϵ}_i}\Rightarrow e^{-\alpha }=\frac{N}{{\sum }_i{e}^{-\beta {ϵ}_i}}$$ 带入 (1) 式得到
  $$n_i=e^{-\alpha }{e}^{-\beta {ϵ}_i}=N\frac{e^{-\beta {ϵ}_i}}{{\sum }_j{e}^{-\beta {ϵ}_j}}\text{\hspace{1em}(2)}$$

• 现在考察一个平动自由度上的能量（以 $x$ 轴为例，注意这是相空间六个维度中的一维）
  $${\overline{ϵ}}_{i_x}=\frac{{\sum }_i{n}_i{ϵ}_{i_x}}{{\sum }_j{n}_j}=\frac{{\sum }_i{n}_i{ϵ}_{i_x}}{N}=\frac{{\sum }_i{ϵ}_{i_x}{e}^{-\beta {ϵ}_{i_x}}}{{\sum }_j{e}^{-\beta {ϵ}_{j_x}}}\text{\hspace{1em}(3)}$$ 把能量展开为动能，再用动量表示
  $${ϵ}_{i_x}=\frac{1}{2}m{v}_{i_x}^2=\frac{(m{v}_{i_x}{)}^2}{2m}=\frac{p_{i_x}^2}{2m}$$ 回代到 (3) 式，得到
  $${\overline{ϵ}}_{i_x}=\frac{{\sum }_i\frac{p_{i_x}^2}{2m}\exp (-\beta \frac{p_{i_x}^2}{2m})}{{\sum }_j\exp (-\beta \frac{p_{j_x}^2}{2m})}$$ 动量是连续变化的，改成积分形式
  $${\overline{ϵ}}_{i_x}=\frac{\frac{1}{2m}{\int }_{-\infty }^{\infty }{p}_{i_x}^2\exp (-\beta \frac{p_{i_x}^2}{2m})d{p}_{i_x}}{{\int }_{-\infty }^{\infty }\exp (-\beta \frac{p_{j_x}^2}{2m})d{p}_{j_x}}\text{\hspace{1em}(4)}$$

  这里要用两个积分公式（可以先自己乘自己，然后转换为极坐标求解）
  $$\begin{array}{rl}& {\int }_{-\infty }^{+\infty }{e}^{-a{x}^2}dx={\left(\frac{\pi }{a}\right)}^{1/2}\\ & {\int }_{-\infty }^{+\infty }{x}^2\cdot e^{-a{x}^2}dx=\frac{1}{2}{\left(\frac{\pi }{a^3}\right)}^{1/2}\end{array}$$

  分别对 (4) 中分母和分子应用上面的两个公式，解得
  $${\overline{ϵ}}_{i_x}=\frac{1}{2\beta }$$ 根据气体分子动理论，一个粒子在一个（平动）自由度上的所分配的平均能量为 $\frac{1}{2}{K}_{B}T$（其中 $K_B$ 是 玻尔兹曼常数），于是有
  $${\overline{ϵ}}_{i_x}=\frac{1}{2\beta }=\frac{1}{2}{K}_{B}T\Rightarrow \beta =\frac{1}{K_{B}T}$$

• 把 $\beta$ 回代到 (2) 式，就解出了给定温度 $T$ 下，各个相空间内粒子数为
  $$n_i=N\frac{e^{-{ϵ}_i/K_{B}T}}{{\sum }_j{e}^{-{ϵ}_j/K_{B}T}}$$ 上式除以粒子总数 $N$ 转换为概率，即得到 玻尔兹曼分布
  $$P_i=\frac{e^{-{ϵ}_i/K_{B}T}}{{\sum }_j{e}^{-{ϵ}_j/K_{B}T}}\text{\hspace{1em}(5)}$$ 其中 $e^{-{ϵ}_i/K_{B}T}$ 称为 玻尔兹曼因子；分母部分 ${\sum }_j{e}^{-{ϵ}_j/K_{B}T}$ 称为 配分函数，它是个归一化参数

2.4 得到 softmax

• 考察上面得到玻尔兹曼分布式，如果定义 $f(i)=-\frac{{ϵ}_i}{K_B}$，那么玻尔兹曼分布就变为 softmax 函数
  $$P_i=\frac{\exp (\frac{f(i)}{T})}{{\sum }_j\exp (\frac{f(j)}{T})}$$ 这里保留了温度系数 $T$ 做进一步讨论
  1. $T\to \infty$ 时，类似于随机选择状态，即每个状态的概率相等，这就是所谓的最大熵状态，也就是最混乱的状态；
  2. $T\to 0$ 时，趋向于选择最好的状态，其概率是1，其它状态概率为0（$T$ 在分母上，因此不可能为0，物理学中绝对零度是不可能达到的。统计力学中温度反映所有粒子的总的动能，绝对零度意味着所有粒子全部处于静止状态，这是不可能的）

3. 基于能量的模型（Energy Based Model, EBM）

• 玻尔兹曼分布建立了从能量到概率的映射关系，而机器学习模型本质也是学一个映射，所以可以利用玻尔兹曼分布建立一种基于能量的机器学习模型。具体来说，对于监督学习中的一个训练样本 $(x,y)$，我们可以把样本 $x$ 和标记 $y$ 的匹配度表示为能量 $E(x,y)$，能量越小匹配度越高，这样就把我们所关心变量的各种组合和一个标量能量联系在一起，训练模型的过程就是不断改变标量能量的过程
• EBM 为许多概率和非概率的学习方法提供了一个统一的框架，特别是图模型和其他结构化模型的非概率训练，监督学习、生成对抗、自监督学习等等范式都能涵盖进去。基于能量的学习可以被看作是预测、分类或决策任务的概率估计的替代方法。由于不需要适当的归一化，基于能量的方法避免了概率模型中与估计归一化常数相关的问题。此外，由于没有标准化条件，在学习机器的设计中允许了更多的灵活性。大多数概率模型都可以看作是特殊类型的基于能量的模型，其中能量函数满足一定的归一化条件，损失函数通过学习优化，具有特定的形式
• 关于 Energy Based Model，可以进一步参考 综述：基于能量的模型