1 下降算法中的搜索方向
1.1 下降方向的判定
根据泰勒展开f ( x k + α k d k ) = f ( x k ) + α k g k T d k + o ( ∣ ∣ α k d k ∣ ∣ 2 ) f(x_k+\alpha_kd_k)=f(x_k)+\alpha_kg^T_kd_k+o(||\alpha_kd_k||^2)f(xk+αkdk)=f(xk)+αkgkTdk+o(∣∣αkdk∣∣2)，忽略极小项后，我们可以在x k x_kxk点处找到f ( x ) f(x)f(x)的一条切线s ( α ) = f ( x k ) + g k T d k α s(\alpha)=f(x_k)+g_k^Td_k \alphas(α)=f(xk)+gkTdkα，这条切线的斜率是g k T d k g_k^Td_kgkTdk。我们不难得出结论，如果g k T d k < 0 g_k^Td_k<0gkTdk<0，则该方向为下降方向。
1.2 下降算法的收敛性
前面我们给出过下降算法的收敛性的定义：对于迭代序列x ( k ) x^{(k)}x(k)，k趋近于无穷时一阶偏导向量的范数为0。之后我们在介绍精确搜索与非精确搜索时均强调了搜索方向必须是下降方向。事实上，如果步长由精确搜索或者Wolfe-Powell搜索产生，而每一步的搜索方向都是下降方向，则必定满足下降算法的收敛性。（证明见附录1）
因此接下来我们在研究计算搜索方向的算法时，最重要的前提就是计算出的搜索方向应当是下降方向。而步长计算则使用精确搜索、Wolfe-Powell搜索或强条件的Wolfe-Powell搜索。
2 最速梯度下降法
2.1 最速梯度下降法步骤
既然我们要寻找下降方向，我们首先想到的就是梯度的反方向。函数沿着梯度方向数值上升最快，那么沿着梯度的反方向数值下降也就最快。所以，对于每一步，将负梯度方向作为下降方向的方法，又叫最速梯度下降法。
步骤：
- 给定初始点x 0 ∈ R n x_0\in R^nx0∈Rn，精度ϵ > 0 \epsilon>0ϵ>0，令k = 0 k=0k=0
- 若∣ ∣ ∇ f ( x k ) ∣ ∣ ≤ ϵ ||\nabla f(x_k)||\leq\epsilon∣∣∇f(xk)∣∣≤ϵ，则得解x k x_kxk，算法终止
- 计算d k = − ∇ f ( x k ) d_k=-\nabla f(x_k)dk=−∇f(xk)
- 计算步长α k \alpha_kαk
- 令x k + 1 = x k + α k d k , k = k + 1 x_{k+1}=x_k+\alpha_kd_k,k=k+1xk+1=xk+αkdk,k=k+1，转步2
2.2 性能评估
- 收敛性：因为每一步均产生下降方向，所以必定收敛。
- 收敛速度：用正定二次函数逼近点附近的函数。若该正定二次函数黑森矩阵的所有特征值相等，则超线性收敛；其余时候线性收敛。（证明见附录2）
- 二次终止性：显然满足
- 计算量：小，只需要计算梯度
- 储存空间：小，只需要储存梯度
2.3 实战测试
对于本文集的第一篇文章深入浅出最优化(1) 最优化问题概念与基本知识中提出的最小二乘问题，x 1 , x 2 , x 3 , x 4 x_1,x_2,x_3,x_4x1,x2,x3,x4的初值均在[ − 2 , 2 ] [-2,2][−2,2]的范围内随机生成，总共生成100组起点。统计迭代成功（在1000步内得到最优解且单次步长搜索迭代次数不超过1000次）的样本的平均迭代步数、平均迭代时间和得到的最优解及残差平方和最小值。
| 平均迭代步数 | 平均迭代时间 | 最优解 | 残差平方和最小值 |
|---|---|---|---|
| 234.68 | 3.31s | x 1 = 0.1755 x 2 = 0.3717 x 3 = 0.0439 x 4 = 0.2290 x_1=0.1755~x_2=0.3717~x_3=0.0439~ x_4=0.2290x1=0.1755 x2=0.3717 x3=0.0439 x4=0.2290 | 2.3065 × 1 0 − 4 2.3065\times10^{-4}2.3065×10−4 |
3 阻尼牛顿法
3.1 阻尼牛顿法步骤
古典牛顿法的思想是用近似二次函数的极小点作为原问题的新的近似解。f ( x ) f(x)f(x)在x k x_kxk处二阶泰勒展开式为f ( x ) = f ( x k ) + ∇ f ( x k ) T ( x − x k ) + 1 2 ( x − x k ) T ∇ 2 f ( x k ) ( x − x k ) + o ( ∣ ∣ x − x k ∣ ∣ 2 ) f(x)=f(x_k)+\nabla f(x_k)^T(x-x_k)+\frac{1}{2}(x-x_k)^T\nabla^2f(x_k)(x-x_k)+o(||x-x_k||^2)f(x)=f(xk)+∇f(xk)T(x−xk)+21(x−xk)T∇2f(xk)(x−xk)+o(∣∣x−xk∣∣2)，二次近似函数Q ( x ) = f ( x k ) + ∇ f ( x k ) T ( x − x k ) + 1 2 ( x − x k ) T ∇ 2 f ( x k ) ( x − x k ) Q(x)=f(x_k)+\nabla f(x_k)^T(x-x_k)+\frac{1}{2}(x-x_k)^T\nabla^2f(x_k)(x-x_k)Q(x)=f(xk)+∇f(xk)T(x−xk)+21(x−xk)T∇2f(xk)(x−xk)，且∇ Q ( x ) = ∇ f ( x k ) + ∇ 2 f ( x k ) ( x − x k ) \nabla Q(x)=\nabla f(x_k)+\nabla^2f(x_k)(x-x_k)∇Q(x)=∇f(xk)+∇2f(xk)(x−xk)，若∇ 2 f ( x k ) \nabla^2f(x_k)∇2f(xk)正定，则Q ( x ) Q(x)Q(x)的极小点为∇ f ( x k ) + ∇ 2 f ( x k ) ( x − x k ) \nabla f(x_k)+\nabla^2f(x_k)(x-x_k)∇f(xk)+∇2f(xk)(x−xk)的解，把二次函数的极小点作为x k + 1 x_{k+1}xk+1，则x k + 1 = x k − ∇ 2 f ( x k ) − 1 ∇ f ( x k ) x_{k+1}=x_k-\nabla^2f(x_k)^{-1}\nabla f(x_k)xk+1=xk−∇2f(xk)−1∇f(xk)，称该迭代公式为古典牛顿法的迭代公式，其中d k = − ∇ 2 f ( x k ) − 1 ∇ ( x k ) d_k=-\nabla^2f(x_k)^{-1}\nabla(x_k)dk=−∇2f(xk)−1∇(xk)为x k x_kxk处的牛顿方向。
如果目标函数就是二次函数，则向着牛顿方向步长为1的搜索可以直接搜索到局部最优解。但如果二次函数仅仅是目标函数的近似，则步长需要使用Wolfe-Powell搜索来求取，这时候算法被称为阻尼牛顿法。
步骤：
- 定初始点x 0 ∈ R n x_0\in R^nx0∈Rn，精度ϵ > 0 \epsilon>0ϵ>0，令k = 0 k=0k=0
- 若∣ ∣ ∇ f ( x k ) ∣ ∣ ≤ ϵ ||\nabla f(x_k)||\leq\epsilon∣∣∇f(xk)∣∣≤ϵ，则得解x k x_kxk，算法终止
- d k = − ∇ 2 f ( x k ) − 1 ∇ ( x k ) d_k=-\nabla^2f(x_k)^{-1}\nabla(x_k)dk=−∇2f(xk)−1∇(xk)
- 计算步长α k \alpha_kαk
- 令x k + 1 = x k + α k d k , k = k + 1 x_{k+1}=x_k+\alpha_kd_k,k=k+1xk+1=xk+αkdk,k=k+1，转步2
3.2 性能评估
- 收敛性：当且仅当每一步都有目标函数的黑森矩阵，也就是近似二次函数的黑森矩阵正定时，牛顿方向才是下降方向，所以收敛性不能保证。
- 收敛速度：若在最优解附近二阶连续可微且最优解处梯度为0、黑森矩阵正定，则算法超线性收敛。特别地，若目标函数在最优解处二阶李普希兹连续，则算法二阶收敛。（证明见附录3）
- 二次终止性：显然满足
- 计算量：大，需要计算黑森矩阵，在变量数目多时计算量大
- 储存空间：大，需要储存黑森矩阵，在变量数目多时需要储存空间大
4 牛顿-梯度下降混合法
4.1 牛顿-梯度下降混合法步骤
由于我们不能保证牛顿法产生下降方向，但又希望能够借助牛顿法的二阶收敛性提高最速梯度下降法的收敛速度，我们可以将两者进行融合。对于牛顿法无法产生下降方向的时刻，使用最速梯度下降法来产生下降方向。
步骤：
- 定初始点x 0 ∈ R n x_0\in R^nx0∈Rn，精度ϵ > 0 \epsilon>0ϵ>0，令k = 0 k=0k=0
- 若∣ ∣ ∇ f ( x k ) ∣ ∣ ≤ ϵ ||\nabla f(x_k)||\leq\epsilon∣∣∇f(xk)∣∣≤ϵ，则得解x k x_kxk，算法终止
- d k = − ∇ 2 f ( x k ) − 1 ∇ ( x k ) d_k=-\nabla^2f(x_k)^{-1}\nabla(x_k)dk=−∇2f(xk)−1∇(xk)
- 若∇ f ( x k ) T d k ≥ 0 \nabla f(x_k)^Td_k\geq0∇f(xk)Tdk≥0，取d k = − ∇ f ( x k ) d_k=-\nabla f(x_k)dk=−∇f(xk)
- 计算步长α k \alpha_kαk
- 令x k + 1 = x k + α k d k , k = k + 1 x_{k+1}=x_k+\alpha_kd_k,k=k+1xk+1=xk+αkdk,k=k+1，转步2
4.2 实战测试
对于本文集的第一篇文章深入浅出最优化(1) 最优化问题概念与基本知识中提出的最小二乘问题，x 1 , x 2 , x 3 , x 4 x_1,x_2,x_3,x_4x1,x2,x3,x4的初值均在[ − 2 , 2 ] [-2,2][−2,2]的范围内随机生成，总共生成100组起点。统计迭代成功率（在1000步内得到最优解且单次步长搜索迭代次数不超过1000次）、平均迭代步数、平均迭代时间和得到的最优解及残差平方和最小值。
| 平均迭代步数 | 平均迭代时间 | 最优解 | 残差平方和最小值 |
|---|---|---|---|
| 56.0 | 1.02s | x 1 = 0.1926 x 2 = 0.1816 x 3 = 0.1158 x 4 = 0.1321 x_1=0.1926~x_2=0.1816~x_3=0.1158~ x_4=0.1321x1=0.1926 x2=0.1816 x3=0.1158 x4=0.1321 | 1.5397 × 1 0 − 4 1.5397\times10^{-4}1.5397×10−4 |
代码实现
本博客所有代码在https://github.com/HarmoniaLeo/optimization-in-a-nutshell开源，如果帮助到你，请点个star，谢谢这对我真的很重要！
你可以在上面的GitHub链接或本文集的第一篇文章深入浅出最优化(1) 最优化问题概念与基本知识中找到Function.py和lagb.py
最速下降法：
import numpy as np
from Function import Function #定义法求导工具
from lagb import * #线性代数工具库
n=2 #x的长度
def myFunc(x): #x是一个包含所有参数的列表
return x[0]**2 + 2*x[1]**2 + 2*x[0] - 6*x[1] +1 #目标方程
x=np.zeros(n) #初值点
rho=0.6
beta=1
sigma=0.4
e=0.001
k=0
tar=Function(myFunc)
while tar.norm(x)>e:
d=-tar.grad(x)
a=1
if not (tar.value(x+a*d)<=tar.value(x)+rho*a*dot(turn(tar.grad(x)),d) and dot(turn(tar.grad(x+a*d)),d)>=sigma*dot(turn(tar.grad(x)),d)):
a=beta
while tar.value(x+a*d)>tar.value(x)+rho*a*dot(turn(tar.grad(x)),d):
a*=rho
while dot(turn(tar.grad(x+a*d)),d)<sigma*dot(turn(tar.grad(x)),d):
a1=a/rho
da=a1-a
while tar.value(x+(a+da)*d)>tar.value(x)+rho*(a+da)*dot(turn(tar.grad(x)),d):
da*=rho
a+=da
x+=a*d
k+=1
print(k)
print(x)
牛顿-梯度下降混合法：
import numpy as np
from Function import Function #定义法求导工具
from lagb import * #线性代数工具库
from scipy import linalg
n=2 #x的长度
def myFunc(x): #x是一个包含所有参数的列表
return x[0]**2 + 2*x[1]**2 + 2*x[0] - 6*x[1] +1 #目标方程
x=np.zeros(n) #初值点
rho=0.6
beta=1
e=0.001
sigma=0.4
k=0
tar=Function(myFunc)
while tar.norm(x)>e:
try:
d=linalg.solve(tar.hesse(x),-tar.grad(x))
if tar.value(x)-tar.value(x+d)<0:
d=-tar.grad(x)
except Exception:
d=-tar.grad(x)
a=1
if not (tar.value(x+a*d)<=tar.value(x)+rho*a*dot(turn(tar.grad(x)),d) and dot(turn(tar.grad(x+a*d)),d)>=sigma*dot(turn(tar.grad(x)),d)):
a=beta
while tar.value(x+a*d)>tar.value(x)+rho*a*dot(turn(tar.grad(x)),d):
a*=rho
while dot(turn(tar.grad(x+a*d)),d)<sigma*dot(turn(tar.grad(x)),d):
a1=a/rho
da=a1-a
while tar.value(x+(a+da)*d)>tar.value(x)+rho*(a+da)*dot(turn(tar.grad(x)),d):
da*=rho
a+=da
x+=a*d
k+=1
print(k)
print(x)
附录
- 记向量d k d_kdk和− ∇ f ( x k ) -\nabla f(x_k)−∇f(xk)的夹角为θ k \theta_kθk，则有c o s = − ∇ f ( x k ) T d k ∣ ∣ ∇ f ( x k ) ∣ ∣ ∣ ∣ d k ∣ ∣ cos=\frac{-\nabla f(x_k)^Td_k}{||\nabla f(x_k)||~||d_k||}cos=∣∣∇f(xk)∣∣ ∣∣dk∣∣−∇f(xk)Tdk
给出下面的基本假设：f ( x ) f(x)f(x)连续可微且有下界，且∇ f ( x ) \nabla f(x)∇f(x)李普希兹连续，即存在常数L > 0 L>0L>0，使得∣ ∣ ∇ f ( x ) − ∇ f ( y ) ∣ ∣ ≤ L ∣ ∣ x − y ∣ ∣ ||\nabla f(x)-\nabla f(y)||\leq L||x-y||∣∣∇f(x)−∇f(y)∣∣≤L∣∣x−y∣∣
则有定理，若序列{ x k } \{x_k\}{xk}由下降算法产生，其中步长α k \alpha_kαk由精确搜索或Wolfe-Powell搜索产生，则∑ k = 0 ∞ ∣ ∣ ∇ f ( x k ) ∣ ∣ 2 c o s 2 θ k < + ∞ \displaystyle\sum_{k=0}^∞||\nabla f(x_k)||^2cos^2\theta_k<+∞k=0∑∞∣∣∇f(xk)∣∣2cos2θk<+∞。特别地，若存在常数δ > 0 \delta>0δ>0使得c o s θ k ≥ δ cos\theta_k\geq\deltacosθk≥δ，则lim k → ∞ ∣ ∣ ∇ f ( x k ) ∣ ∣ = 0 \displaystyle\lim_{k→∞}||\nabla f(x_k)||=0k→∞lim∣∣∇f(xk)∣∣=0。这个定理说明了产生的方向与负梯度方向夹角小于π 2 \frac{\pi}{2}2π时，可以保证算法收敛性。
下面根据假设证明该定理：由Wolfe-Powell搜索条件及假设，我们有− ( 1 − σ ) ∇ f ( x k ) T d k ≤ [ ∇ f ( x k + α k d k ) − ∇ f ( x k ) ] T d k ≤ α k L ∣ ∣ d k ∣ ∣ 2 -(1-\sigma)\nabla f(x_k)^Td_k\leq[\nabla f(x_k+\alpha_kd_k)-\nabla f(x_k)]^Td_k\leq\alpha_kL||d_k||^2−(1−σ)∇f(xk)Tdk≤[∇f(xk+αkdk)−∇f(xk)]Tdk≤αkL∣∣dk∣∣2，即可得α k ≥ − 1 − σ L ∇ f ( x k ) T d k ∣ ∣ d k ∣ ∣ 2 = − c 1 ∇ f ( x k ) T d k ∣ ∣ d k ∣ ∣ 2 \alpha_k\geq-\frac{1-\sigma}{L}\frac{\nabla f(x_k)^Td_k}{||d_k||^2}=-c_1\frac{\nabla f(x_k)^Td_k}{||d_k||^2}αk≥−L1−σ∣∣dk∣∣2∇f(xk)Tdk=−c1∣∣dk∣∣2∇f(xk)Tdk。这里c 1 = 1 − σ L c_1=\frac{1-\sigma}{L}c1=L1−σ，进一步由Wolfe-Powell搜索第一个条件，得f ( x k + α k d k ) − f ( x k ) ≤ − ρ c 1 [ ∇ f ( x k ) T d k ] 2 ∣ ∣ d k ∣ ∣ 2 = − ρ c 1 ∣ ∣ ∇ f ( x k ) ∣ ∣ 2 c o s 2 θ k f(x_k+\alpha_kd_k)-f(x_k)\leq-\rho c_1\frac{[\nabla f(x_k)^Td_k]^2}{||d_k||^2}=-\rho c_1||\nabla f(x_k)||^2cos^2\theta_kf(xk+αkdk)−f(xk)≤−ρc1∣∣dk∣∣2[∇f(xk)Tdk]2=−ρc1∣∣∇f(xk)∣∣2cos2θk
将上面的不等式左右两边从k = 0 k=0k=0到∞ ∞∞相加，并注意到f ( x ) f(x)f(x)有下界，可以得到定理第一条成立，由无穷级数收敛的必要条件可以得到第二条定理成立。
最速梯度下降法收敛速度证明：https://blog.csdn.net/weixin_43010548/article/details/97619095
牛顿法收敛速度证明：
提出定理：设f ff在x ∗ ∈ R n x^*\in R^nx∗∈Rn的某个邻域内二次连续可微且x ∗ x^*x∗满足∇ f ( x ∗ ) = 0 \nabla f(x^*)=0∇f(x∗)=0，∇ 2 f ( x ∗ ) \nabla^2 f(x^*)∇2f(x∗)正定，则存在常数δ > 0 \delta>0δ>0，使得当x 0 ∈ U δ ( x ∗ ) = { x ∣ ∣ ∣ x − x ∗ ∣ ∣ < δ } x_0\in U_\delta(x^*)=\{x|||x-x^*||<\delta\}x0∈Uδ(x∗)={x∣∣∣x−x∗∣∣<δ}时，由单位步长牛顿法x k + 1 = x k − ∇ 2 f ( x k ) − 1 ∇ f ( x k ) , k = 0 , 1 , 2 , . . . x_{k+1}=x_k-\nabla^2 f(x_k)^{-1}\nabla f(x_k),k=0,1,2,...xk+1=xk−∇2f(xk)−1∇f(xk),k=0,1,2,...产生的序列超线性收敛于x ∗ x^*x∗，此外，若∇ 2 f \nabla^2f∇2f在x ∗ x^*x∗李普希兹连续，则序列$ {x_k}二 次 收 敛 于 二次收敛于二次收敛于x^*$。
