工程优化:梯度下降(SGD & BGD)、牛顿法的概念与实现

工程优化:梯度下降(SGD & BGD)、牛顿法的概念与实现

九月 20, 2019
当前设备屏幕尺寸过小,推荐使用PC模式浏览。

引用

Wikipedia:Gradient descent
Wikipedia:Newton’s method


Preliminaries

先要了解一下术语及前导知识。

目标/代价/损失/误差函数

  1. 我们称需要最大化或最小化的函数为目标函数(Objective Function)或准则(Criterion)。特别地,当进行最小化搜索时,也将其称为代价函数(Cost Function)、损失函数(Loss Function)、误差函数(Error Function)。

minxf(x)\min \limits_x f(x)

在工程优化中,常见称为目标函数地情况多些。在机器学习中,由于优化目标常常是目标真实值和预测值间的偏差,所以称呼为代价函数、损失函数的情况多一些。

  1. 在数域XX上搜索得到一个x^\hat x,使得目标函数目标函数取得最小/大值的过程称为参数搜索
  2. 我们将搜索得到的能使目标函数取得最小/大值的参数估计值xx称为最优参数估计值,记为xx^*
  3. 对于大部分函数及优化、搜索算法,并不能保证得到一个全局最优的参数估计值,此时称之为极大/极小点

局部最优、全局最优与驻点、鞍点

  1. 在搜索过程中获得一个最优值x0x_0时,使得当f(x0)f(x_0)取得全局最小值时,该解x0x_0全局最优解;反之,其为局部最优解

  2. 凸函数的局部最优值就是全局最优值,局部最优解属于全局最优解集。

  3. 对于一个nn维的目标函数f(x)f(x),对于x0Xnx_0 \in X^n当此处的偏导数f(x0)=0f'(x_0)=0,且此时f(x0)f(x_0)为一个局部极大/极小小值,称该点为一个驻点临界点;若此时f(x0)f(x_0)在某个维度上稳定(极大/极小),在其它至少一个维度上不稳定(仍有下降/上升空间),则其为鞍点

驻点与鞍点的存在,都意味着对目标函数的参数搜索陷入了某种算法困境;因为无论朝任意哪个法向看去,其导数方向都为0,不可能再通过梯度下降等方式寻得更优值。对于此类困境,退火算法等智能算法进行了一定上的解决。

导数与偏导数、梯度

  1. 对于一维函数f(x)f(x),其导数定义为:

f(x)=limΔx0f(x0+Δx)f(x0)Δxf'(x)=\lim \limits_{\Delta x \rightarrow 0} \frac{f(x_0+\small{\Delta} x)-f(x_0)}{\small{\Delta} x}

  1. 对于多维函数f(x1,...,xn)f(x_1,...,x_n),对xix_i求导数dfdxi\frac{df}{dx_i},将其记为偏导数fxI\frac{\partial f}{\partial x_I}。特别的,记录梯度f(x)\triangledown f(x)或简记为f\triangledown f为对xix_i求偏导后的列向量:

f(x)=(fx1,fx2,...,fx1)T\triangledown f(x)=(\frac{\partial f}{\partial x_1}, \frac{\partial f}{\partial x_2},..., \frac{\partial f}{\partial x_1})^T

常见函数的梯度与海森矩阵

给定向量aabbxx,矩阵AAHf(x)H_{f(x)}是函数f(x)f(x)的海森矩阵:

g(x)=aTx+bg(x)=aHg(x)=0g(x)=xTxg(x)=2xHg(x)=2Hf(x)g(x)=xTAxg(x)=(A+AT)xHg(x)=(A+AT)Hf(x)g(x)=xTAx+aTx+bg(x)=(A+AT)x+aHg(x)=(A+AT)Hf(x)\begin{aligned} g(x)=a^Tx+b \rightarrow \triangledown & g(x)=a\\ & H_{g(x)}=0 \\ g(x)=x^Tx \rightarrow \triangledown & g(x)=2x\\ & H_{g(x)}=2H_{f(x)} \\ g(x)=x^TAx \rightarrow \triangledown & g(x)=(A+A^T)x\\ & H_{g(x)}=(A+A^T)H_{f(x)} \\ g(x)=x^TAx+a^Tx+b \rightarrow \triangledown & g(x)=(A+A^T)x+a \\ & H_{g(x)}=(A+A^T)H_{f(x)} \end{aligned}

方向导数与下降方向

  • 设有 单位向量h=(h1,h2,...,hn)TRnh=(h_1,h_2,...,h_n)^T\in \R^n,表示为nn维空间中的一个方向,对于nn元函数f(x)f(x),若limΔx0f(x0+αh)f(x0)α\lim \limits_{\Delta x \rightarrow 0} \frac{f(x_0+\alpha h)-f(x_0)}{\alpha}存在,则称该极限值为f(x)f(x)x0x_0处沿方向hh方向导数,记作f(x0)h\frac{\partial f(x_0)}{\partial h}。如图所示:

1.jpg

根据导数的定义不难得到:

f(x0)h=limΔx0f(x0+αh)f(x0)α=limΔx0f(x0)T(αh)+o(αh)α=f(x0)Th=f(x0)cos(f(x0),h)\begin{aligned} \frac{\partial f(x_0)}{\partial h} &= \lim \limits_{\Delta x \rightarrow 0} \frac{f(x_0+\alpha h)-f(x_0)}{\alpha} \\ &= \lim \limits_{\Delta x \rightarrow 0} \frac{\triangledown f(x_0)^T (\alpha h)+o(\parallel \alpha h \parallel)}{\alpha} \\ &= \triangledown f(x_0)^T h \\ &= \parallel \triangledown f(x_0) \parallel \cos(\triangledown f(x_0), h) \end{aligned}

  • f(x0)h<0\frac{\partial f(x_0)}{\partial h}<0hhf(x)f(x)x0x_0处的一个下降方向;反之称为上升方向
  • 使得f(x0)h<0\frac{\partial f(x_0)}{\partial h}<0最大负值时,称其为最速下降方向
  • 显然,目标函数f(x)f(x)的负梯度方向f(x)\triangledown f(x)就是目标函数的最速下降方向。

证明:在x0x_0处即为使得f(x0)h\frac{\partial f(x_0)}{\partial h}取得最大负值f(x0)-\parallel \triangledown f(x_0) \parallel,此时 h=f(x0)f(x0)h=-\frac {\triangledown f(x_0)}{\parallel f(x_0) \parallel},即 cos(f(x0),h)=1\cos(\triangledown f(x_0), h)=-1

>返回目录

梯度之上的凸函数预测

海森矩阵(Hessian matrix)

  1. 若存在f:RnRf:\R^n \rightarrow \R,即一个多维输入xRx\in\R到一维输出的映射,若其在任意维度上都二阶可导,则定义其海森矩阵:

H=[f2x12f2x1x2f2x1xnf2x2x1f2x22f2x2xnf2xnx1f2xnx2f2xn2]H = \begin{bmatrix} \frac{\partial f^2}{\partial x_1^2} & \frac{\partial f^2}{\partial x_1\partial x_2} & \cdots &\frac{\partial f^2}{\partial x_1 \partial x_n} \\ \frac{\partial f^2}{\partial x_2 \partial x_1} & \frac{\partial f^2}{\partial x_2^2} & \cdots &\frac{\partial f^2}{\partial x_2 \partial x_n} \\ \vdots & \vdots & \ddots& \vdots \\ \frac{\partial f^2}{\partial x_n \partial x_1} & \frac{\partial f^2}{\partial x_n\partial x_2} & \cdots &\frac{\partial f^2}{ \partial x_n^2} \end{bmatrix}

显然 ,HT=HH^T = HHH 的尺寸为 n×nn\times n

雅可比矩阵(Jacobian matrix)

若存在f:RnRmf:\R^n \rightarrow \R^m,即一个多维输入xRnx\in\R^n到多维输出f(x)Rmf(x)\in\R^m的映射,则ff的雅可比矩阵:

J=[fx1fx2fxn]=[f1x1f1xnfmx1fmxn]\begin{aligned} J &= [ \frac{\partial f}{\partial x_1} \frac{\partial f}{\partial x_2} \cdots \frac{\partial f}{\partial x_n}] \\ &= \begin{bmatrix} \frac{\partial f_1}{\partial x_1} & \cdots & \frac{\partial f_1}{\partial x_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial f_m}{\partial x_1} & \cdots & \frac{\partial f_m}{\partial x_n} \\ \end{bmatrix} \end{aligned}

显然,雅可比矩阵的尺寸为n×mn \times mJij=fixjJ_{ij} = \frac{\partial f_i}{\partial x_j}

曲率预测

针对一个下降优化问题,当目标函数ff在某个维度djd_j上在x0x_0处的二阶偏导取得负值,说明其一阶导在此处呈下降趋势,即根据一阶导预测的函数下降趋势实际上要比原预测的要大;反之亦然。

由二阶导的此类性质可以更好地拟合目标函数在局部地变化趋势,引用工程优化中地理论,就是比“一阶导”多向前看了一步。优化中常用地牛顿法(由牛顿求根法得到的优化方法)就是建立在曲率预测地基础上的。

微分学三大中值定理

  • 三大中值定理的关系:罗尔中值定理 to\mathop \rightarrow \limits^{to} 拉格朗日中值定理 to\mathop \rightarrow \limits^{to} 柯西中值定理。

罗尔(中值)定理

Wikipedia:Rolle’s theorem

罗尔(中值)定理有多种表述形式,其常见标准形式为:若函数映射ff[a,b][a,b]上连续、在开区间(a,b)(a,b)上可微且f(a)=f(b)f(a)=f(b)则在区间(a,b)(a,b)中至少存在一个值cc使得f(x)=0f'(x)=0

2.jpg

罗尔(中值)定理是泰勒定理和拉格朗日中值定理的基础,是三大微分中值理论的基础。

拉格朗日中值定理

Wikipedia:Mean value theorem

标准形式:在罗尔(中值)定理的假设基础上,开区间(a,b)(a,b)中至少存在一点cc,使得:

f(c)=f(b)f(a)baf'(c)=\frac{f(b)-f(a)}{b-a}

3.jpg

若对标准形式,取c=x+θΔxc=x+\theta \Delta xΔy=f(b)f(a)\Delta y=f(b)-f(a)Δx=ba\Delta x= b-a,则可得到增量公式:

Δy=f(x+θΔx)Δx\Delta y= f'(x+\theta \Delta x)\Delta x

拉格朗日中值定理是罗尔(中值)定理的推广,是泰勒公式的一阶展开。

柯西中值定理

Wikipedia:Cauchy mean value theorem

标准形式:在罗尔(中值)定理的假设基础上,开区间(a,b)(a,b)中至少存在一点cc,使得:

f(c)g(c)=f(a)f(b)g(b)g(a)\frac{f'(c)}{g'(c)}=\frac{f(a)-f(b)}{g(b)-g(a)}

4.jpg

证明过程:
构建辅助函数:F(x)=[f(b)f(a)][g(x)g(a)][g(b)g(a)][f(x)f(a)]F(x)=[f(b)-f(a)][g(x)-g(a)]-[g(b)-g(a)][f(x)-f(a)]在区间上连续可导且F(a)=F(b)=0,有罗尔定理可得其cc点导数为00。求导即可证明。

  • 该方法同样适用于对拉格朗日中值定理的求证。

柯西中值定理是拉格朗日中值定理的推广。

泰勒定理

Wikipedia:Taylor’s theorem

泰勒定理的核心思想是通过若干线性构造函数θ(x)\theta(x)去拟合原函数(线性逼近),从而由一个点开始推出原函数的分布:

f(x)=f(a)+f(a)(xa)+θ(x)(xa),limxaθ(x)=0f(x)=f(a)+f'(a)(x-a)+\theta(x)(x-a),\lim \limits_{x\rightarrow a}\theta(x)=0

它的思想很类似牛顿迭代法(牛顿-拉弗森求根法,Newton’s method)。由拉格朗日中值定理,很容易得到它的正确性。

  1. 假设此处做一次一阶逼近函数P1(x)=f(a)+f(a)(xa)P_1(x)=f(a)+f'(a)(x-a),则P1(x)P_1(x)f(x)f(x)的逼近误差为:R1(x)=f(x)P1(x)=θ(x)(xa)R_1(x)=f(x)-P_1(x)=\theta(x)(x-a)。不难看出,P1(x)P_1(x)就是f(x)f(x)aa处的切线。
  2. 若取P2(x)=f(a)+f(a)(xa)+f(x)2(xa)2P_2(x)=f(a)+f'(a)(x-a)+\frac{f''(x)}{2}(x-a)^2,误差R2(x)=θ(x)(xa)2R_2(x)=\theta(x)(x-a)^2
  3. 以此类推,当取Pn(x)=f(a)+f(a)(xa)++f(n)(x)n!(xa)nP_n(x)=f(a)+f'(a)(x-a)+\cdots+\frac{f^{(n)}(x)}{n!}(x-a)^n时,误差Rn(x)=θ(x)(xa)nR_n(x)=\theta(x)(x-a)^n
  4. limxaθ(x)=0\lim \limits_{x\rightarrow a}\theta(x)=0,则RnR_n要比Rn1R_{n-1}更快地逼近00,即误差更小。

由此,给出泰勒公式:

f(x)=f(a)+f(a)(xa)++f(n)(x)n!(xa)n+Rnf(x) = f(a)+f'(a)(x-a)+\cdots+\frac{f^{(n)}(x)}{n!}(x-a)^n+R_n

RnR_n在此处被称为皮亚诺余项(Peano form of the remainder)。经由拉格朗日推导和柯西推导,还有拉格朗日余项和柯西余项等,但它们大部分都是等价的。

>返回目录


梯度下降(Gradient Descent)

梯度下降是一种一阶优化算法。

  • 优化算法的阶数将在下一节提及牛顿方向与牛顿法小节中提到。

下降更新

根据上文Preliminaries中关于梯度和最速下降方向的知识,梯度下降法建议的更新点为:

xt+1=xtλf(xt)x_{t+1}={x_t}-\lambda\triangledown f(x_t)

其中,λ\lambda下降步长xtx_ttt时刻的搜索结果。

在机器学习领域中, λ\lambda有时也被写作 ϵ\epsilonα\alpha,称为学习率(learning rate)。

  • 从直觉上来说,此时f(xt+1)<f(xt)f(x_{t+1})<f(x_t),然而事实并非如此。因此我们还要谈谈步长关于下降效果的影响。

下降步长对梯度下降的影响

一个常见的例子:
给定目标函数f(x)=x2f(x)=x^2,初始点x0=1x_0=1,其梯度(导数)为[2x][2x],给定步长λ=1\lambda=1,根据上文的更新公式计算x1=x0λ×2x0=1x_1=x_0-\lambda \times 2x_0=-1,发现f(x1)=f(x0)f(x_1)=f(x_0);继续进行下降更新,发现x2=1=x0x_2=1=x_0…;如此反复,搜寻结果会反复在[1,1][-1, 1]之间横跳,且目标函数的值并不下降。

可见,基于最速下降方向的梯度下降虽然可以找到一个合适的下降方向,但并不能确定”下降的量“,对于λ\lambda的搜索,我们有以下方法:

  1. 求得f(xt)\triangledown f(x_t),列出式子xt+1=xtλf(xt)x_{t+1}={x_t}-\lambda \triangledown f(x_t)
  2. 将该式子带入目标函数f(x)f(x),即t+1t+1时刻有ψ(λ)\psi(\lambda)
  3. ψ(λ)=0\psi'(\lambda)=0,解得最优下降步长λ\lambda^*(即,尝试求得沿该下降方向能得到的驻点)。
  • 理论上,该方法对一阶函数下降一步即可得到驻点(如果是凸函数则下降一步即可得最优解)。

扩展:对步长的搜索有很多中方法,例如AdaGradAdamMomentum等方法,详见:Wikipedia:SGD-Extensions and variants

这个方法看起来似乎非常棒,能够帮助我们很快地进行下降,但是在编程时,难免会遇到目标函数过于复杂,无法直接推得λ\lambda^*的表示式,或者沿某个方向函数持续下降(如f(x)=exf(x)=e^{-x}),导致步长的取值向上溢出。所以,更聪明、健壮的方法是进行一些更“暴力”的一维搜索,如接下来要说的“成功-失败法”。

基于"成功-失败法"的梯度下降步长搜索

Alg 基于成功-失败法的梯度下降:

  1. 给定初始点x0Rnx_0\in\R^n,搜索精度ϵ>0\epsilon>0,预设增量步长λ0>ϵ>0\lambda_0>\epsilon>0
  2. 对于tt时刻,计算梯度f(xt)\triangledown f(x_t),若f(xt)<ϵ\triangledown f(x_t)<\epsilon,算法结束转步骤7,否则转步骤3;
  3. (开始准备搜索最优步长)将搜索步长初始化为λ=0\lambda^*=0λ=λ0\triangledown \lambda=\lambda_0,>f(xt(λ))=f(xt)f(x_t^{(\lambda)})=f(x_t);
  4. 计算更新点xt+1(λ)=xt(λ+λ)f(xt)x_{t+1}^{(\lambda)}=x_t-(\lambda^*+\triangledown \lambda)\triangledown f(x_t),判断若>λ<ϵ|\triangledown >\lambda|<\epsilon,搜索结束,得到此时的最优步长搜索结果,即λ\lambda^*,转步骤6;
  5. 计算f(xt+1(λ))f(x_{t+1}^{(\lambda)}),若f(xt+1(λ))>f(xt(λ))f(x_{t+1}^{(\lambda)})>f(x_t^{(\lambda)})则搜索成功,更新>f(xt(λ))=f(xt+1(λ))f(x_t^{(\lambda)})=f(x_{t+1}^{(\lambda)})λ:=λ+λ\lambda^*:=\lambda^*+\triangledown \lambda、>λ:=2λ\triangledown\lambda:=2\triangledown\lambda;否则,认为搜索失败,需要缩小步长,即更新>λ:=0.25λ\triangledown\lambda:=-0.25\triangledown\lambda,转步骤4;
  6. xt+1=xtλf(xt)x_{t+1}=x_t-\lambda^*\triangledown f(x_t)xt:=xt+1x_t:=x_{t+1},跳回步骤2。
  7. 得到一个(局部)最优解x:=xtx^*:=x_t

该算法的python伪代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# 带有成功失败法的梯度下降
epsilon = 1e-1
lambda0 = 0.1
while True:
delta_x = grad(x) # 得到梯度
if abs(delta_x) < epsilon :
break
t0_cost = cost(x)
# 成功失败发搜索最优步长
delta_lambda, lambda_star = lambda0, 0
while abs(delta_lambda) >= epsilon:
t1_cost = cost(x-(lambda_star+delta_lambda))
if t1_cost > t0_cost: # 搜索失败
delta_lambda *= -0.25
else: # 搜索成功
t0_cost = t1_cost
lambda_star = delta_lambda+lambda_star
delta_lambda *= 2
x -= lambda_star * delta_x # 根据最优步长更新解集
print("x* = ", x)

作为对比,将应用了成功失败法和定小步长的梯度下降法做了对比,可以看到前者经过十个点得搜索就完成了梯度下降(>代码见此<):

5.jpg

点此下载代码

另外,基于成功搜索法的梯度下降可以给出一个非常大的步长的同时保证收敛,而定步长法若步长太大,则很可能不收敛而左右摇摆。

  • 在下文的牛顿法中,还将介绍一种基于二阶导数测试,即泰勒公式的最优步长确定方式。

重要概念区分及术语

随机/(小)批量梯度下降(SGD/BGD/MBGD)

Wikipedia:Stochastic gradient descent

在机器学习中,往往存在着大量样本,例如随本文附上的>代码<,就是对f:xR2yR=ax02+sin(bx1)f: x \in \R^2 \rightarrow y \in \R =ax_0^2+sin(bx_1)的代价函数 cost(y,x)=12Ni=0N(yf(x))2cost(y,x) = \frac{1}{2N}\sum_{i=0}^N(y-f(x))^2的优化。
NN取得很大的值时,即样本数量很大时,代价函数的数值容易溢出,因此往往只随机选择一个或若干个样本参与运算。该做法也有加速过程训练过程的作用。(SGD本质上是在做梯度的无偏估计)

我们称之为随机梯度下降(Stochastic gradient descent);而针对整个数据集的梯度下降方法,称为批量梯度下降(Batch Gradient Descent)。

“只随机抽取一个样本的SGD”和“BGD”相比,前者的优化过程是震荡的,因为对于单个样本来说,其有可能是一个噪音,并不能指向全局最优点。所以,作为折中方案,小批量梯度下降(Mini-Batch Gradient Descent, MBGD)被提了出来:样本被随机划分为nn个子集,每次取一个子集进行训练。本质上MBGD就是SGD。

显然,"只取一个样本的SGD"实际上是"MBGD"的一个特例。

基于套袋法的思想,我们给定epoch的概念:每一个epoch中都随机划分子集做MBGD算法,重复epochs次。

  • 台大李宏毅教授的课程24:00左右处提到,基于Batch的方法在运算量和单样本SGD的运算量是一致的,但基于GPU的并行计算技术,适当的Batch可以提升训练速度,因此我们通常所使用的"SGD"方法实际上是基于"GPU+MBGD+Epoch"的方法,在各类教程上也常常利用"SGD"来做简单指代,但在严格意义上来说,SGD是该方法的超集,BGD和SGD是同一级别的概念,MBGD是SGD的子概念。这是最严谨最正式的概念区分

  • 但是,随着GPU技术的发展,我们已经很少只用单样本来做梯度下降了。不同的文献对随机梯度下降有不同的指代方法,例如有时把此处提到的抽取mm个样本(m>1m>1)的方法称为小批量梯度下降(MBGD),而只随机抽取一个样本的梯度下降称为SGD。我个人偏向以“花书”《深度学习》为准,即提到SGD时直接默认指代MBGD,。

梯度下降(Gradient Descent)与"最速下降法"(method of steepest descent)的歧义性

一个很有趣的一点是,梯度下降(一阶迭代优化算法)与最速下降法(method of steepest descent)似乎并不是同一种东西。之所以用上了“似乎”二字,是因为我即在中英文文献中读到了认为二者是同一算法的观点,也读到了认为二者不同的观点。在Wikipedia中做出了很显眼的特殊注释:

Gradient descent is also known as steepest descent. However, gradient descent should not be confused with the method of steepest descent for approximating integrals.

大概确实是太多人弄混了(包括我),在中文文献中似乎区分不那么明显,但我在看某些英文文献时,看到文中提到了method of steepest descent下意识地认为这是Gradient Descent,然而继续看下去发现好像并不是,当时就觉得非常困惑,直到看到了wiki上的相应解释。

  • 但也有一些英文文献提到"最速下降法"(也许此处应该换一个中文译法)method of steepest descent时,指代的就是Gradient Descent
  • 在"花书"《深度学习》第55页明确指出method of steepest descent就是gradient descent

因此,为了严谨考虑,此处还是只把优化中的Gradient Descent称为梯度下降,把分析中的method of steepest descent最速下降法称为最速下降法。

梯度下降的缺点及改进

梯度下降的最大问题在于锯齿现象:当沿着某个方向下降时已经找到了某个维度的极小点,而下一次下降又会破坏该维度的极小点,从而使得搜索路径来回震荡,不断处于“找到某一维度极小点-破坏该维度极小点”的恶性病态循环过程中。最终使得共轭梯度法的收敛速度较慢。

  • 为此,提出了共轭方向法来解决这一问题。

参考:[工程优化]共轭方向法(Conjugate direction method)的共轭梯度法(Conjugate gradient method)实现【附python代码】

>返回目录


牛顿法(Newton’s method)

数值分析中的牛顿-拉弗森求根法

完整理解牛顿法,应当从牛顿求根法开始。
牛顿迭代求根法:

  1. 给定曲线f(x)f(x)和初始点x0x_0,求导数f(x0)f'(x_0),根据泰勒公式前二项,做一次逼近,即:x1=x0f(x0)f(x0)x_1=x_0-\frac{f(x_0)}{f'(x_0)}
  2. 依次重复这个以上过程,第nn次逼近之后,xn=xn1f(xn1)f(xn1)x_n=x_{n-1}-\frac{f(x_{n-1})}{f'(x_{n-1})},当f(xn)0<ϵ|f(x_n)-0|<\epsilon时算法结束,即得到此时方程的根。

事实上,xn=xn1f(xn1)f(xn1)x_n=x_{n-1}-\frac{f(x_{n-1})}{f'(x_{n-1})}是式子f(xn)=f(xn1)+f(xn1)(xnxn1)f(x_n)=f(x_{n-1})+f'(x_{n-1})(x_n-x_{n-1})关于xx轴的交点。而后者其实就是其在xn1x_{n-1}处的切线。

给出下面的python代码,对f(x)=x2+sin(πx)f(x)=x^2+sin(\pi x)求根:

1
2
3
4
5
6
7
8
import numpy as np
f = lambda x: np.power(x, 2) + np.sin(np.pi*x)
df = lambda x: 2*x + np.pi*np.cos(np.pi*x)
def newton(xt, m_lambda=1.0):
# 这里给定了一个默认步长1,与上文中的梯度下降中的步长意义相同
while abs(f(xt) - 0) > 1e-1:
xt = xt - m_lambda * f(xt) / df(xt)
return xt

由数学知识我们可以很容易得知方程在(1,1)(-1,1)上有唯一的两个根,其中有一个根显然为0:

6.jpg

运行上面的算法,分别从10和-10处逼近,我们希望分别得到这两个根:

1
2
3
4
5
6
7
8
9
10
11
12
import numpy as np
import matplotlib.pyplot as plt
xt1 = newton(10)
xt2 = newton(-10)
# print("xt1=", xt1)
# print("xt2=", xt2)
plt.plot(np.arange(-1, 1, 0.1), f(np.arange(-1, 1, 0.1)))
plt.plot(np.arange(-1, 1, 0.1), [0 for _ in np.arange(-1, 1, 0.1)], '--')
plt.ylim((-1, 1))
plt.scatter(xt1, f(xt1), s=20, marker='o', c='r') # 第一个根标记为红点
plt.scatter(xt2, f(xt2), s=20, marker='o', c='g') # 第二个根标记为绿点
plt.show()

函数在很长一段时间内不收敛,经过努力和一些调整画图结果如下:

7.jpg

我们发现两个根居然统一逼近到了0这个实根,并没有我们预期的双根搜索结果。
经过上文对于梯度下降的理解,我们尝试加入一个较小的步长:

1
2
xt1 = newton(10, 0.01)  # 也可尝试从-0.1逼近0根
xt2 = newton(-10, 0.01)

其求根的结果如图所示(根据计算精度的不同会略有不同):

8.jpg

可以知道,牛顿求根法实际上就是梯度下降的一种应用,当算法的终止条件变为 f(x)<ϵ|f'(x)|<\epsilon 时,就成了求曲线的驻点。

二阶导数测试(Second derivative test)

牛顿求根法不断用切线修正、拟合原函数。由泰勒公式我们可以知道,由高次项对函数进行拟合测试,能更快的从一点x0x_0出发拟合整个函数,所以此处展开二阶泰勒公式:

f(x)f(x0)+f(x0)(xx0)+12f(x0)(xx0)2f(x)\approx f(x_0)+f'(x_0)(x-x_0)+\frac{1}{2}f''(x_0)(x-x_0)^2

转化为nn阶函数,其向量形式为:

f(x)f(x0)+(xx0)Tg(x0)+12(xx0)TH(x0)(xx0)f(x)\approx f(x_0)+(x-x_0)^Tg_{(x_0)}+\frac{1}{2}(x-x_0)^TH_{(x_0)}(x-x_0)

其中,H(x0)H_{(x_0)}是函数在x0x_0处的海森矩阵(Hessian matrix),g(x0)g_{(x_0)}是函数在x0x_0处的梯度f(x0)\triangledown f(x_0)

>回忆梯度下降<的做法,我们给出更新点x1=x0λf(x0)x_1=x_0-\lambda \triangledown f(x_0),带入原始,并尝试使ψ(λ)=0\psi'(\lambda)=0求得λ\lambda^*

ψ(λ)=f(x1)f(x0)λg(x0)Tg(x0)+λ22g(x0)TH(x0)g(x0)ψ(λ)=g(x0)Tg(x0)+λg(x0)TH(x0)g(x0)=0λ=g(x0)Tg(x0)g(x0)TH(x0)g(x0)\begin{aligned} \psi (\lambda) &= f(x_1) \approx f(x_0)- \lambda g_{(x_0)}^Tg_{(x_0)} + \frac{\lambda^2}{2}g_{(x_0)}^TH_{(x_0)}g_{(x_0)} \\ \psi' (\lambda) &= -g_{(x_0)}^Tg_{(x_0)} + \lambda g_{(x_0)}^TH_{(x_0)}g_{(x_0)} = 0 \\ \lambda^* &=\frac{g_{(x_0)}^Tg_{(x_0)}}{g_{(x_0)}^TH_{(x_0)}g_{(x_0)}} \end{aligned}

根据二阶导数测试,我们还可以得知一个一阶导数为0的究竟是鞍点还是驻点:
f(x0)=0f'(x_0)=0处,若f(x0)>0f''(x_0)>0则意味着一阶导数的微分增量为正,即随着xx的增大f(x)f'(x)增大;若f(x)f'(x)增大,说明f(x)f(x)将会增大,即其为一个驻点或局部极小点(取决于函数的凹凸性)。

  • f(x0)=0f'(x_0)=0f(x0)>0f''(x_0)>0,其为一个驻点或局部极小点。
  • f(x0)=0f'(x_0)=0f(x0)<0f''(x_0)<0,其为一个鞍点或局部极大点。

利用上文中“梯度下降”小节中的>代码<,加入二次导数测试(SDT)求得的λ\lambda^*,可以得到SDT方法与”成功-失败法“的结果对比图:

9.jpg

我们注意到基于二次导数测试与利用成功失败法求得最优步长的方法具有近似的下降性能,同时由于免去了对步长的测试,它在事实上进行了更少的迭代。
但是,在代码中,我们输出了几种方法的用时,发现反而是“成功-失败法”用时较短。这是由于由于每一次进行二阶导数测试都要求取一次海森矩阵并进行相应矩阵运算,这样复杂的运算延长了整体的执行时间。针对这个问题,拟牛顿法进行了优化。

牛顿方向与牛顿法

由二阶导数测试,我们知道二阶拟合可以更快的对一个优化问题进行目标下降。基于这个思想,回忆牛顿求根法中迭代求根的步骤,我们继续来展开泰勒公式的二阶形式:

f(x)f(x0)+(xx0)Tg(x0)+12(xx0)TH(x0)(xx0)f(x)\approx f(x_0)+(x-x_0)^Tg_{(x_0)}+\frac{1}{2}(x-x_0)^TH_{(x_0)}(x-x_0)

对其求导数(>回忆Preliminaries中常见函数的二阶求导<):

f(x)g(x0)+12(H(x0)+H(x0)T)(xx0)g(x0)+H(x0)(xx0)\begin{aligned} f'(x) &\approx g_{(x_0)}+\frac{1}{2}(H_{(x_0)}+H_{(x_0)}^T)(x-x_0)\\ &\approx g_{(x_0)}+H_{(x_0)}(x-x_0) \end{aligned}

令其等于0即可得到下降根:

x1=x0H(x0)1g(x0)x_1 = x_0 - H_{(x_0)}^{-1}g_{(x_0)}

即牛顿下降公式:xt+11=xtH(xt)1g(xt)x_{t+11} = x_t - H_{(x_t)}^{-1}g_{(x_t)}

这里,H(xt)1g(xt)H_{(x_t)}^{-1}g_{(x_t)} 被称为牛顿方向。

同样,使用>代码<测试,让它与基于“成功-失败法”、“二阶导数测试”的梯度下降法进行比较:

10.jpg

我们惊喜地发现,它的性能直逼SDT,并且如果你运行了代码,还会发发现在时间开销上它是最小的,与SDT方法相比快了将近15倍:

11.jpg

优化算法的阶数

我们将牛顿法这种优化算法称为二阶优化算法(Second-order optimization algorithms),将梯度下降算法称为一阶优化算法(First-order optimization algorithms)。

尽管牛顿法看上去比一阶优化算法要完美许多,然而由于其每次下降都要求取一次海森矩阵,同时海森矩阵未必存在逆矩阵(HH正定),因此并不那么实用。涉及到f(x)=1xf(x)=\frac{1}{x}这类的非连续、存在大跳变的函数时,如果没有设定较好的步长,还会导致不收敛的病态条件出现。

因此,拟牛顿法就是为了解决这个问题而提出的。

参考:[工程优化]牛顿法的缺陷及拟牛顿法(Newton’s method):DFP\BFGS\L-BFGS【附python代码实现】