统计计算一|非线性方程的求解

系列文章目录

常见优化器对比：梯度下降法、带动量的梯度下降法、Adagrad、RMSProp、Adam

---

一、最大似然估计，频率派和贝叶斯学派

（一）似然函数

给定输出$X$时，关于参数$\Theta$的似然函数$L(\Theta |X)$在数值上等于给定参数$\Theta$后变量$X$的概率。如果 X = { x 1 , . . . , x n } X=\{x_1,...,x_n\} X={x1​,...,xn​}是独立同分布（i.i.d）的，每个变量都服从密度函数$f(x|\Theta )$，该密度函数依赖于$p$个未知参数组成的向量 Θ = { θ 1 , . . . , θ p } \Theta=\{\theta_1,...,\theta_p\} Θ={θ1​,...,θp​}，那么联合似然函数就为：
$$L(\Theta |X)=\prod _{i=1}^{n}f(x_i|\Theta )$$
当数据不是独立同分布时，联合似然函数仍然表示关于$\Theta$的联合密度：$f(x_1,...,x_n|\Theta )$。

（二）频率学派

1、参数估计

对频率学派来说，参数是固定的，通常使用MLE的方法来估计参数，寻找以较高概率产生观察数据的参数值。似然函数在MLE $\hat{\Theta }$处达到最大值：
$$\hat{\Theta }=\underset{\Theta }{\mathrm{argmax}}L(\Theta |X)=\underset{\Theta }{\mathrm{argmax}}l(\Theta |X)$$
其中$l(\Theta |X)$是对数似然函数。在相当普遍的条件下，当$n\to \infty$时，$\hat{\Theta }$是渐进无偏的。

2、得分函数

求解MLE，即寻找$l^{\prime }(\theta |x)=0$的根。其中偏导数$l^{\prime }(\theta |x)$就是得分函数

• 得分函数关于$x$在$\theta$固定下的条件期望为0：
  $$E[\frac{\partial }{\partial \theta }logf(x;\theta )|\theta ]={\int }_R\frac{\frac{\partial }{\partial \theta }f(x;\theta )}{f(x;\theta )}f(x;\theta )\mathrm{d}x=\frac{\partial }{\partial \theta }{\int }_{R}f(x;\theta )\mathrm{d}x=0$$

其中，${\int }_{R}f(x;\theta )\mathrm{d}x=1$。

• 得分函数的方差等于Fisher信息：
  $$Var(U(\theta ))=I(\theta )=E[(\frac{\partial l(\theta )}{\partial \theta }{)}^2]$$

3、Fisher信息

Fisher 信息是统计学和信息理论中的一个重要概念，用于量化一个参数估计问题中数据对参数信息的敏感程度。具体来说，Fisher 信息度量了参数变化对似然函数变化的影响，被定义为得分函数关于 $x$ 的方差：

$$I(\theta )=E[(\frac{\partial }{\partial \theta }logf(x;\theta ){)}^2|\theta ]={\int }_R(\frac{\partial }{\partial \theta }logf(x;\theta ){)}^{2}f(x;\theta )\mathrm{d}x$$
如果$logf(x;\theta )$对于$\theta$具有两次可微性，并且在一定的正则条件下，Fisher信息也可以写成：
$$I(\theta )=-E[\frac{{\partial }^2}{\partial {\theta }^2}logf(x;\theta )|\theta ]$$

证明过程：

4、多元Fisher信息

多元情形下的 Fisher 信息度量多个参数对似然函数的敏感性，并且涉及 Fisher 信息矩阵：
I ( θ ) : = E { l ′ ( θ ) l ′ ( θ ) T } = − E { l ′ ′ ( θ ) } = − E { Δ 2 l ( θ ) } I(\theta):=E\{l'(\theta)l'(\theta)^T\}=-E\{l''(\theta)\}=-E\{\Delta^2l(\theta)\} I(θ):=E{l′(θ)l′(θ)T}=−E{l′′(θ)}=−E{Δ2l(θ)}
其中，$l^{\prime \prime }(\theta )$表示一个$p\times p$的矩阵，矩阵的第$(i,j)$个元素为${\partial }^{2}l(\theta )/(\partial {\theta }_i\partial {\theta }_j)$。

• 期望Fisher信息：$I(\theta )$
• 观测Fisher信息：$-l^{\prime \prime }(\theta )$

在正则条件下，MLE $\hat{\theta }$的渐进方差-协方差矩阵是$I({\theta }^{\ast }{)}^{-1}$，其中${\theta }^{\ast }$表示$\theta$的真实值。当$n\to \infty$时，$\hat{\theta }$的极限分布是：
$$N_p({\theta }^{\ast },I({\theta }^{\ast }{)}^{-1})$$

5、正态分布未知参数的MLE估计

有些似然函数太复杂，无法通过分析求解（如涉及多元正态分布的情况），寻找最大似然估计（MLE）→ 优化问题 → 解非线性方程。

（三）贝叶斯学派

贝叶斯统计是统计学的一个分支，强调通过结合数据和先验信息来更新对参数的认知。贝叶斯学派基于贝叶斯定理，该定理为在观察数据后更新信念提供了一种正式的方法。贝叶斯学派假设参数是随机的。

• 先验分布$\pi (\Theta )$：关于参数的一些先验信息（例如先前的知识、专家建议）
• 似然函数$L(\Theta |X)$：与频率统计中的定义相同
• 后验分布$f(\Theta |X)$：在观察数据后更新对参数的信息

$$f(\Theta |X)=\frac{f(\Theta ,X)}{f(x)}=\frac{\pi (\Theta )f(X|\Theta )}{\int f(\tau ,X)d\tau }\propto \pi (\Theta )f(X|\Theta )$$

在贝叶斯统计中，如果后验分布$f(\Theta |X)$ 与先验分布 $\pi (\Theta )$ 属于同类，则先验分布与后验分布被称为共轭分布，而先验分布被称为似然函数的共轭先验。共轭先验的使用简化了贝叶斯更新的计算过程，因为它确保了后验分布与先验分布具有相同的分布形式。

贝叶斯学派下求参数有多种方法，如通过 MCMC 采样出后验分布下的参数样本，另外也可以用最大化后验分布的形式，求得参数的点估计，即后验分布的 mode。

二、非线性方程的求解

（一）二分法（Bisection Method）

如果$g^{\prime }$在$[a_0,b_0]$上连续且满足$g^{\prime }(a_0)g^{\prime })b_0)\le 0$，那么中值定理表明在$[a_0,b_0]$上至少存在一个$x^{\ast }\in [a_0,b_0]$，使得$g^{\prime }(x)=0$，因此$x^{\ast }$是$g$的一个局部极值点。

二分法（Bisection Method）是一种简单而有效的数值方法，用于求解连续函数在给定区间内的根。该方法利用连续函数的性质，通过逐步缩小包含根的区间，逐步逼近根的值。令$x^{(0)}=(a_0+b_0)/2$为初始值，则更新公式为：
$$[a_{t+1},b_{t+1}]=\{\begin{array}{ll}[a_t,x^{(t)}],& \text{if g′(at)g′(x(t)≤0}\\ [x^{(t)},b_t],& \text{if g′(at)g′(x(t)>0}\end{array}$$
和
$$x^{(t+1)}=\frac{1}{2}(a_{t+1}+b_{t+1})$$

• 二分法是一种相对较慢的方法：它需要较多的迭代次数才能达到所需的精度。
• 如果 $g\prime$ 在 $[a_0,b_0]$ 上连续，无论 $g\prime \prime$ 是否存在、求导的难易性如何，都可以找到一个根。

（二）不动点迭代法（Fixed-Point Iteration）

1、不动点

一个函数 $f(x)$ 的不动点是指满足 $f(x)=x$ 的 $x$ 值。换句话说，如果将 $x$ 代入函数 $f(x)$ 后，结果等于 $x$，那么 $x$ 就是函数 $f(x)$ 的不动点。在数学中，不动点通常表示函数的稳定状态或者函数在某种意义下的平衡状态。以$f(x)=cos(x)$为例。$f(x)$ 的不动点就是满足 $cos(x)=x$ 的 $x$ 值。用图示来看，不动点就是函数 $f(x)$ 的图像与直线 $y=x$ 的交点。

2、不动点迭代法

不动点迭代法也称为函数迭代法（Functional Iteration），是一种用于求解非线性方程的数值方法。当函数形式较复杂时，求解器不容易直接求得。利用不动点的性质，可以将函数的求解问题转化为不动点方程，从而使用有限次迭代求得方程的解。

基本的不动点迭代算法步骤：

• 选择一个初始猜测值$x_0$，作为迭代的起点
• 迭代更新：使用一个迭代函数$g(x)$，通过不断迭代$x_{n+1}=g(x_n)$来逼近不动点
• 收敛判断：当迭代值$x_n$足够接近不动点$x^{\ast }$时停止迭代。

寻找根的不动点策略是确定一个函数$G$，满足当且仅当$g^{\prime }(x)=0$时$G(x)=x$。最简单方法是使用更新方程$x^{(t+1)}=G(x^{(t)})$。

$G$的选择不是固定的，但最明显的选择是$G(x)=g^{\prime }(x)+x$：
$$x^{(t+1)}=x^{(t)}+g^{\prime }(x^{(t)})$$

例如：求解$g^{\prime }(x)=cosx-x$的根，设置$G(x)=cosx,x_0=0$：

3、不动点迭代法的收敛性

该算法的收敛性取决于函数$G$在$[a,b]$上是否是收缩的：

• $\forall x\in [a,b]$，有$G(x)\in [a,b]$
• $\exists \lambda \in [0,1)$，使得对任意$x_1,x_2\in [a,b]$，有：（Lipschitz条件）
  $$|G(x_1)-G(x_2)|\le \lambda |x_1-x_2|$$

此时就有：

• $G(x)$有一个唯一的固定点$x_{\infty }\in [a,b]$
• $x_n=G(x_{n-1})\to x_{\infty },\forall x_0\in [a,b]$
• $|x_n-x_{\infty }|\le \frac{{\lambda }^n}{1-\lambda }|x_1-x_0|$

证明：（在一个柯西序列中，随着序列的项越来越靠后，序列的项彼此之间的距离会变得越来越小。）

如果对于$[a,b]$ 中的所有 $x$，都满足 $|G\prime (x)|\le \lambda <1$，则满足了Lipschitz 条件。

Lipschitz 连续函数限制了函数变化的速度。具体而言，满足 Lipschitz 条件的函数的斜率必须受到一个实数的限制，该实数被称为 Lipschitz 常数。

4、更新公式的调整

假设函数$G(x)=g^{\prime }(x)+x$在区间$[a,b]$上不满足收敛条件，即$|g^{\prime \prime }(x)+1|<1$不成立，可以选择$G(x)=\alpha g^{\prime }(x)+x$和$|\alpha g^{\prime \prime }(x)+1|<1$来重新缩放非收敛问题，通过引入缩放因子$\alpha$来调整更新公式，使其满足收敛条件。
$$x_{n+1}=x_n+\alpha g^{\prime }(x_n)$$
不动点迭代的效果高度依赖于所选择的$G$的形式。如找到$g^{\prime }(x)=logx+x$的根时，不同$G$的收敛速度为：
$$(x+e^{-x})/2>x^{-x}>-logx$$

（三）牛顿法（Newton’s method）

1、牛顿法算法

牛顿法，也称为牛顿-拉夫森法（Newton-Raphson method），是一种用于求解非线性方程 $f(x)=0$ 的迭代方法。它基于函数的泰勒级数展开，通过反复迭代逼近方程的根。

牛顿法算法步骤：

• 选择一个初始猜测值$x_0$，作为迭代的起点
• 迭代更新：对于$n=0,1,2,...$，利用以下迭代公式更新$x_n$：
  $$x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)}$$

2、牛顿法在寻找函数导数 $g\prime (x)$ 的根的迭代过程

假设我们要找到函数 $g\prime (x)$ 的根，即满足$g^{\prime }(x)=0$的$x$。假设$g^{\prime }$是连续可微的，并且$g^{\prime \prime }(x)\ne 0$。我们可以使用牛顿法来逼近这个根。

在第$t$次迭代中：

• 泰勒级数展开：在点$x^{(t)}$处对$g^{\prime }(x)$进行一次泰勒级数展开：
  $$g^{\prime }(x)\approx g^{\prime }(x^{(t)})+(x-x^{(t)})g^{\prime \prime }(x^{(t)})$$
• 线性近似：设$g^{\prime }(x)=0$求根：
  $$0\approx g^{\prime }(x^{(t)})+(x-x^{(t)})g^{\prime \prime }(x^{(t)})$$
• 解方程：解上式得到$x^{\ast }$的近似值：
  $$0=g^{\prime }(x^{(t)})+(x-x^{(t)})g^{\prime \prime }(x^{(t)})\Rightarrow x-x^{(t)}=-\frac{g^{\prime }(x^{(t)})}{g^{\prime \prime }(x^{(t)})}$$
• 进而得到更新公式：
  $$x^{\ast }=x^{(t)}-\frac{g^{\prime }(x^{(t)})}{g^{\prime \prime }(x^{(t)})}=x^{(t)}+h^{(t)}$$
  这个方程描述了一个依赖于当前猜测 $x^{(t)}$ 和修正项 $h^{(t)}$ 的 $x^{\ast }$ 的近似值。迭代这个策略给出了牛顿法的更新方程：
  $$x^{(t+1)}=x^{(t)}+h^{(t)}$$

3、牛顿法的收敛性

牛顿法的收敛性取决于函数 g 的形状和初始值的选择。

• 根附近的收敛性定理：如果$g^{\prime \prime \prime }$连续且$g^{\prime }(x^{\ast })=0,g^{\prime \prime }(x^{\ast })\ne 0$，那么存在一个$x^{\ast }$的邻域，在该邻域内从任意$x_0\in (x^{\ast }-\delta ,x^{\ast }+\delta )$开始的牛顿法迭代：
  $$x_n=G(x_{n-1})=x_{n-1}-\frac{g^{\prime }(x_{n-1})}{g^{\prime \prime }(x_{n-1})}$$
  会收敛到$x^{\ast }$

证明：

• 全局收敛性：当 $g\prime$ 具有二阶连续可导性、是凸函数并且根存在的时候，牛顿法可以从任意起始点收敛到该根。

凸函数的定义：
如果一个实值函数的图像上任意两点间的线段位于这两点之间的图像上方，则该函数被称为凸函数。对如果对任意的$x,y\in S$及$0\le \theta \le 1$，有：
$$f(\theta x+(1-\theta )y)\le \theta f(x)+(1-\theta )f(y)$$
则称$f$为$S$上的凸函数。
一个单变量的二次可微函数是凸函数，当且仅当它的二阶导数在整个定义域上非负。

• 区间上的收敛性：当从区间$[a,b]$ 中的某个点开始时，满足以下条件则牛顿法将从区间内的任意$x_0$收敛
  • $g^{\prime \prime }(x)\ne 0$在$[a,b]$上
  • $g^{\prime \prime \prime }(x)$在$[a,b]$不变符号
  • $g^{\prime }(a)g^{\prime }(b)<0$并且
  • $|g^{\prime }(a)/g^{\prime \prime }(a)|<b-a$并且$|g^{\prime }(b)/g^{\prime \prime }(b)|<b-a$

4、牛顿法的收敛速度

求根方法的速度通常通过其收敛阶数来衡量。如果一个方法具有收敛阶数$\beta$，则当$n$趋向无穷时，满足${\lim }_{n\to \infty }|x_n-x^{\ast }|=0$以及
$$\underset{n\to \infty }{\lim }\frac{|x_n-x^{\ast }|}{|x_{n-1}-x^{\ast }{|}^{\beta }}=c$$

• 二分法是线性收敛的，即$\beta =1$
• 牛顿法具有二次收敛，即$\beta =2$：
  

5、多维变量的情况

在变量多维时，对于牛顿法，再次使用一阶泰勒级数展开来近似$g^{\prime }(X^{\ast })$：
$$g^{\prime }(X^{\ast })\approx g^{\prime }(X^{(t)})+(X^{\ast }-X^{(t)}{)}^T{g}^{\prime \prime }(X^{(t)})$$
对应的更新公式为：
$$X^{(t+1)}=X^{(t)}+g^{\prime \prime }(X^{(t)}{)}^{-1}{g}^{\prime }(X^{(t)})$$
在最大似然估计中，可以用$I(\Theta (t))$，即在$\Theta (t)$处的期望Fisher信息，来替代$\Theta (t)$处的观测信息：
$${\Theta }^{(t+1)}={\Theta }^{(t)}+I({\theta }^{(t)}{)}^{-1}{l}^{\prime }({\Theta }^t)$$

（四）Fisher Scoring方法

Fisher Scoring 是一种用于优化模型参数的迭代方法，特别常用于统计学和机器学习中的最大似然估计。它与牛顿-拉夫森法类似，但用期望的（Fisher 信息矩阵）代替观测的海森矩阵。Fisher Scoring 主要应用于广义线性模型（GLM）和其他参数估计问题。

由于$I(\theta )$可以近似为$-l^{\prime \prime }(\theta )$，因此Fisher Scoring方法的更新公式即为：
$${\theta }^{(t+1)}={\theta }^{(t)}+l^{\prime }({\theta }^{(t)})I({\theta }^{(t)}{)}^{-1}$$

（五）割线法（Secant Method）

除牛顿法以外，单变量的求根问题也可以使用割线法。割线法是微积分中的一种数值逼近方法，用于求解函数的根（即方程的解）。它基于一个简单的思想：通过一条割线逼近函数的根。牛顿法基于切线逼近函数的根，割线法则是通过两个接近根的点的割线来逼近根。表现在更新方程上就是将$g^{\prime \prime }(x)$替换为$[g^{\prime }(x^{(t)})-g^{\prime }(x^{(t-1)})]/(x^{(t)}-x^{(t-1)})$：
$$x^{(t+1)}=x^{(t)}-g^{\prime }(x^{(t)})\frac{x^{(t)}-x^{(t-1)}}{g^{\prime }(x^{(t)})-g^{\prime }(x^{(t-1)})}$$

（六）基于牛顿法思想的数值求根方法

牛顿法并不能确保每次迭代目标函数都是上升的，而且计算成本大。合理的选择更新形式，比如最速上升算法，有可能会得到更好的表现，以及更小的计算量。一些非常有效的方法都依赖于以下形式的更新方程：
$$X^{(t+1)}=X^{(t)}-(M^{(t)}{)}^{-1}\nabla g(X^{(t)})$$
其中$M^{(t)}$是一个近似海森矩阵${\nabla }^{2}g(X^{(t)})$的$p\times p$矩阵。=

• Fisher Scoring方法:$M^{(t)}=-I(\Theta (t))$
• 最速上升算法：$M^{(t)}=-I$
  $$X^{(t+1)}=X^{(t)}+\nabla g(X^{(t)})$$
• 沿着最陡上升方向进行缩放步长：
  $$X^{(t+1)}=X^{(t)}+{\alpha }^{(t)}\nabla g(X^{(t)})$$

在做上述优化算法的时候，可迫使它的每一步，也就是每一次更新，目标函数都是上升的：

• 对于任何固定的 $x^{(t)}$和负定 $M(t)$，对于${\alpha }^{(t)}\to 0$ 有：
  $$g(x^{(t+1)})-g(x^{(t)})=g(x^{(t)}+h^{(t)})-g(x^{(t)})=-{\alpha }^{(t)}g\prime (x^{(t)}{)}^T(M^{(t)}{)}^{-1}g\prime (x^{(t)})+o({\alpha }^{(t)})$$
  其中 $h^{(t)}=-{\alpha }^{(t)}[M^{(t)}{]}^{-1}\nabla g(X^{(t)})$，第二个等式来自于 Taylor 展开$g(x^{(t)}+h^{(t)}))=g(x^{(t)})+g\prime (x^{(t)}{)}^T{h}^{(t)}+o({\alpha }^{(t)})$。

因此，如果 $-M(t)$ 是正定的，通过选择足够小的 $\alpha (t)$可以确保上升，从而得到 $g(x(t+1))-g(x(t))>0$。

（七）深度学习中的优化算法

考虑最小化平均损失函数：$g(x)=\frac{1}{n}{\sum }_{i=1}^n{l}_i(x)$

• GD (Gradient Descent): $x_t=x_{t-1}-\alpha \nabla g(x_{t-1})$

• SGD (Stochastic Gradient Descent)：$x_t=x_{t-1}-\alpha \nabla {\ell }_{t_i}(x_{t-1})$

• BGD (Batch Gradient Descent)： $x_t=x_{t-1}-\frac{\alpha }{|B|}{\sum }_{i\in B_t}\nabla {\ell }_i(x_{t-1})$

• SGDM（Stochastic Gradient Descent with Momentum）是对随机梯度下降（SGD）的扩展，它引入动量来加速收敛并减小震荡。
  $$g_t=\sum _{i\in B_t}\nabla {\ell }_i(x_{t-1})$$
  $$v_t=\beta v_{t-1}+g_t,x_t=x_{t-1}-\alpha v_t$$
  梯度的指数加权移动平均：
  $$v_t=g_t+\beta g_{t-1}+{\beta }^2{g}_{t-2}+{\beta }^3{g}_{t-3}+...$$
  $\beta$ 常用的值为： [0.5, 0.9, 0.95, 0.99].

• AdaGrad 是一种自适应学习率优化算法，它根据历史梯度调整每个参数的学习率。类似标准化。（其中 · 表示逐元素相乘。）
  $$s_t=s_{t-1}+g_t\cdot g_t,x_t=x_{t-1}-\frac{\alpha }{\sqrt{s_t}+ϵ}\cdot \nabla g_t$$

• RMSProp 是一种自适应学习率优化算法，通过引入衰减因子来解决AdaGrad 的一些局限性。
  $$s_t=\gamma s_{t-1}+(1-\gamma )g_t\cdot g_t,x_t=x_{t-1}-\frac{\alpha }{\sqrt{s_t}+ϵ}\cdot g_t$$

• Adam 是一种自适应学习率优化算法，结合了 Momentum 和 RMSProp的优点。
  $$v_t={\beta }_1{v}_{t-1}+(1-{\beta }_1)g_t,s_t={\beta }_2{s}_{t-1}+(1-{\beta }_2)g_t\cdot g_t$$
  常见的值是 ${\beta }_1=0.9，{\beta }_2=0.999$。由于当 $t$ 受限时，指数加权和不等于 1，因此进行偏差修正：
  $${\hat{v}}_t=\frac{v_t}{1-{\beta }_t^1},{\hat{s}}_t=\frac{s_t}{1-{\beta }_t^2}$$
  最终更新公式为：
  $$x_t=x_{t-1}-\alpha \frac{{\hat{v}}_t}{\sqrt{\widehat{s_t}+ϵ}}$$