设
n
n
n元实值函数
f
(
x
)
f(\boldsymbol{x})
f(x)连续可微,且在区域
Ω
\Omega
Ω内有唯一局部最小值点
x
0
\boldsymbol{x}_0
x0,根据极值点的必要条件知
f
′
(
x
0
)
=
0
f'(\boldsymbol{x}_0)=0
f′(x0)=0。给定初始点
x
1
\boldsymbol{x}_1
x1及容错误差
ε
\varepsilon
ε,利用迭代式
x
k
+
1
=
x
k
+
α
k
d
k
,
k
=
1
,
2
,
⋯
.
\boldsymbol{x}_{k+1}=\boldsymbol{x}_k+\alpha_k\boldsymbol{d}_k,k=1,2,\cdots.
xk+1=xk+αkdk,k=1,2,⋯.
其中,
d
k
∈
R
n
\boldsymbol{d}_k\in\text{ℝ}^n
dk∈Rn为方向向量,正数
α
k
>
0
\alpha_k>0
αk>0为搜索步长。每次迭代寻求合适的方向
d
k
\boldsymbol{d}_k
dk和步长
α
k
\alpha_k
αk,使得
f
(
x
k
+
1
)
=
f
(
x
k
+
α
k
d
k
)
<
f
(
x
k
)
,
f(\boldsymbol{x}_{k+1})=f(\boldsymbol{x}_k+\alpha_k\boldsymbol{d}_k)<f(\boldsymbol{x}_k),
f(xk+1)=f(xk+αkdk)<f(xk),
且
{
x
k
}
\{\boldsymbol{x}_k\}
{xk}收敛于
f
(
x
)
f(\boldsymbol{x})
f(x)的局部最优点
x
0
\boldsymbol{x}_0
x0。
可以证明,
d
k
=
−
∇
f
(
x
k
)
\boldsymbol{d}_k=-\nabla f(\boldsymbol{x_k})
dk=−∇f(xk)是使得
f
(
x
)
f(\boldsymbol{x})
f(x)在
x
k
\boldsymbol{x}_k
xk处下降最快的方向向量。定义一元函数
ϕ
k
(
α
)
=
f
(
x
k
+
α
d
k
)
,
α
∈
R
+
,
\phi_k(\alpha)=f(\boldsymbol{x}_k+\alpha\boldsymbol{d}_k),\alpha\in\text{ℝ}^+,
ϕk(α)=f(xk+αdk),α∈R+,
为使
f
(
x
k
+
α
k
d
k
)
=
ϕ
(
α
k
)
f(\boldsymbol{x}_k+\alpha_k\boldsymbol{d}_k)=\phi(\alpha_k)
f(xk+αkdk)=ϕ(αk)在
x
k
\boldsymbol{x}_k
xk处沿方向
d
k
\boldsymbol{d}_k
dk下降值最大,理想的是计算
α
k
=
arg
min
α
>
0
ϕ
k
(
α
)
.
\alpha_k=\arg\min_{\alpha>0}\phi_k(\alpha).
αk=argα>0minϕk(α).
用如此选取的方向向量
d
k
\boldsymbol{d}_k
dk和步长
α
k
\alpha_k
αk,从
k
=
1
k=1
k=1开始,每次迭代使由迭代式
x
k
+
1
=
x
k
+
α
k
d
k
\boldsymbol{x}_{k+1}=\boldsymbol{x}_k+\alpha_k\boldsymbol{d}_k
xk+1=xk+αkdk算得序列
{
x
k
}
\{\boldsymbol{x}_k\}
{xk}的函数值满足
f
(
x
k
+
1
)
<
f
(
x
k
)
f(\boldsymbol{x}_{k+1})<f(\boldsymbol{x}_k)
f(xk+1)<f(xk)。一旦得到
∣
f
′
(
x
k
)
∣
<
ε
|f'(\boldsymbol{x}_k)|<\varepsilon
∣f′(xk)∣<ε,则认为
x
k
\boldsymbol{x}_k
xk十分接近
f
(
x
)
f(\boldsymbol{x})
f(x)的驻点。由于可微函数
f
(
x
)
f(\boldsymbol{x})
f(x)有局部唯一最小值点,因此可以作为最优解
x
\boldsymbol{x}
x满足容错误差的近似值。这一算法因所选下降方向为
d
k
=
−
∇
f
(
x
k
)
\boldsymbol{d}_k=-\nabla f(\boldsymbol{x_k})
dk=−∇f(xk),故称为梯度下降法。下列代码实现梯度下降搜索算法。
from scipy.optimize import minimize_scalar,OptimizeResult
def gradientDescent(fun,x1,gtol,**options):
xk=x1
dk=-fprime1(fun,xk) #计算搜索方向
phi=lambda a:fun(xk+a*dk)
k=1
while np.linalg.norm(dk)>=gtol:
ak=(minimize_scalar(phi)).x #计算搜索步长
xk+=ak*dk #生成新迭代点
dk=-fprime1(fun,xk) #新的搜索方向
k=k+1
bestx=xk
besty=fun(bestx)
return OptimizeResult(fun=besty, x=bestx, nit=k)
程序的第2~14行定义函数gradientDescent,实现梯度下降算法。参数fun表示目标函数
f
(
x
)
f(\boldsymbol{x})
f(x),x1表示初始迭代点
x
1
\boldsymbol{x}_1
x1,gtol(自定义自由参数)表示容错误差
ε
\varepsilon
ε。
第3~6行进行初始化操作:第3行将表示迭代点
x
k
\boldsymbol{x}_k
xk的变量xk初始化为x1。第4行计算搜索方向
d
k
=
−
∇
f
(
x
k
)
\boldsymbol{d}_k=-\nabla f(\boldsymbol{x}_k)
dk=−∇f(xk)赋予dk。其中,调用计算多元函数梯度的函数grad(详见博文《最优化方法Python计算:n元函数梯度与Hesse阵的数值计算》)计算
f
(
x
)
f(\boldsymbol{x})
f(x)在当前迭代点
x
k
\boldsymbol{x}_k
xk处的梯度
∇
f
(
x
k
)
\nabla f(\boldsymbol{x}_k)
∇f(xk)。第5行定义函数
ϕ
(
α
)
=
f
(
x
k
+
α
d
k
)
\phi(\alpha)=f(\boldsymbol{x}_k+\alpha\boldsymbol{d}_k)
ϕ(α)=f(xk+αdk)
赋予phi。第6行将迭代次数k初始化为1。
第7~11行的while循环执行迭代操作:第8行调用minimize_scalar函数(第1行导入),传递phi,用缺省的brent方法计算搜索步长
α
k
=
arg
min
α
>
0
ϕ
(
α
)
\alpha_k=\arg\min\limits_{\alpha>0}\phi(\alpha)
αk=argα>0minϕ(α)赋予ak。第9行计算迭代
x
k
+
1
=
x
k
+
α
k
d
k
\boldsymbol{x}_{k+1}=\boldsymbol{x}_k+\alpha_k\boldsymbol{d}_k
xk+1=xk+αkdk更新xk。第10行调用grad计算
d
k
+
1
=
−
∇
f
(
x
k
+
1
)
\boldsymbol{d}_{k+1}=-\nabla f(\boldsymbol{x}_{k+1})
dk+1=−∇f(xk+1)更新dk。第11行将迭代次数k自增1。循环往复,直至条件
∥
∇
f
(
x
k
)
∥
<
ε
\lVert\nabla f(\boldsymbol{x}_k)\rVert<\varepsilon
∥∇f(xk)∥<ε满足为止。第12~14行用f(xk),xk和k构造OptimizeResult对象(第1行导入)返回。
例1:考虑Rosenbrock函数
f
(
x
)
=
100
(
x
2
−
x
1
2
)
2
+
(
1
−
x
1
)
2
f(\boldsymbol{x})=100(x_2-x_1^2)^2+(1-x_1)^2
f(x)=100(x2−x12)2+(1−x1)2,
x
=
(
x
1
x
2
)
∈
R
2
\boldsymbol{x}=\begin{pmatrix}x_1\\x_2\end{pmatrix}\in\text{ℝ}^2
x=(x1x2)∈R2。给定初始点
x
1
=
(
2
1
)
\boldsymbol{x}_1=\begin{pmatrix}2\\1\end{pmatrix}
x1=(21)。下列代码用gradientDescent分别以容错误差为
ε
=
1
0
−
3
\varepsilon=10^{-3}
ε=10−3和
ε
=
1
0
−
5
\varepsilon=10^{-5}
ε=10−5,计算
x
0
\boldsymbol{x}_0
x0的近似值。
import numpy as np #导入numpy
from scipy.optimize import minimize, rosen #导入minimize, rosen
x1=np.array([2,1]) #设置初始点
res=minimize(rosen,x1,method=gradientDescent,options={'gtol':1e-3}) #计算最优解
print(res)
res=minimize(rosen,x1,method=gradientDescent,options={'gtol':1e-5}) #计算最优解
print(res)
程序的第3行设置初始点 x 1 = ( 2 1 ) \boldsymbol{x}_1=\begin{pmatrix}2\\1\end{pmatrix} x1=(21)为x1。第4、6行调用minimize函数传递gradientDescent给参数method,并将表示容错误差的gtol参数分别设置为 1 0 − 3 10^{-3} 10−3和 1 0 − 5 10^{-5} 10−5,计算表示Rosenbrock函数的rosen(第2行导入)最优解。运行程序,输出
fun: 1.179532940214763e-06
nit: 1024
x: array([1.00108499, 1.00217598])
fun: 1.1791444533318637e-10
nit: 2068
x: array([1.00001085, 1.00002174])
即gradientDescent函数从
x
1
=
(
2
1
)
\boldsymbol{x}_1=\begin{pmatrix}2\\1\end{pmatrix}
x1=(21)开始,以
ε
=
1
0
−
3
\varepsilon=10^{-3}
ε=10−3为容错误差,迭代1024次,算得Rosenbrock函数的最优解
x
0
\boldsymbol{x}_0
x0的近似值
(
1.00108499
1.00217598
)
\begin{pmatrix}1.00108499\\1.00217598\end{pmatrix}
(1.001084991.00217598)。而以
ε
=
1
0
−
5
\varepsilon=10^{-5}
ε=10−5为容错误差,需迭代2068次,计算结果为
(
1.00001085
1.00002174
)
\begin{pmatrix}1.00001085\\1.00002174\end{pmatrix}
(1.000010851.00002174)。由此可见,不同的容错误差,对算法迭代次数的影响是比较大的,且算法的运行效率不是很高。
写博不易,敬请支持:
如果阅读本文于您有所获,敬请点赞、评论、收藏,谢谢大家的支持!