问题
minimize f(x) 其中, 函数 f:Rn→R 是一个实值函数。 Ω=Rn
subject to x∈Ω
函数梯度迭代示意图
梯度方法
- 思想:令x0作为初始搜索点,沿着梯度负方向构造一个新点x0−α▽f(x0),如果▽f(x0)=0,那么α>0足够小时,有f(x0−α▽f(x0))<f(x0).
- 迭代公式:xk+1=xk−αk▽f(xk).
- 最速下降法:选择合适的步长αk,使得目标函数能够得到最大程度的减小。
αk=arg mina≥0 f(xk−α▽f(xk))
最速下降法
- 定理: 最速下降法搜索函数 f:R2→R 的极小点, 迭代过程产生的序列为 {x(k)}k=0∞, 那么, x(k+1)−x(k) 与 x(k+2)−x(k+1) 正交对于所有的 k⩾0 都成立。
- 下降性质:利用梯度的本质,只要▽f(xk)=0,就能保证每次迭代得到的新数值比原数值小。
- 实际计算中的停止搜索准则:定义最小变化量。
一阶必要条件: ∇f(x(k))<ε
最速下降法例题
最速下降法求 f(x1,x2,x3)=(x1−4)4+(x2−3)2+4(x3+5)4 的极小点。初始搜索点为 x(0)=[4,2−1]⊤, 开展 3 次迭代。
α0x(1)=α⩾0argminf(x(0)−α∇f(x(0)))=α⩾0argmin(0+(2+2α−3)2+4(−1−1024α+5)4)=3.967×10−3=x(0)−α0∇f(x(0))=[4.000,2.008,−5.062]⊤
Step2: α1=α⩾0argmin(0+(2.008+1.984α−3)2+4(−5.062+0.003875α+5)4)=0.5000 x(2)=x(1)−α1∇f(x(1))=[4.000,3.000,−5.060]⊤
Step3: α2=α⩾0argmin(0.000+0.000+4(−5.060+0.003525α+5)4)=16.29 x(3)=[4.000,3.000,−5.002]⊤
理论上:函数 f 的极小点就是 [4,3,−5]⊤
梯度方法特性以及收敛性分析
- 特性:
单调下降算法;全局收敛
- 收敛性分析:
将目标函数设定为二次函数:
V(x)=f(x)+21x∗⊤Qx∗▽f(x)=21(x−x∗)⊤Q(x−x∗)=Qx−b
其中, Q=Q⊤>0 x∗=Q−1b
λmin(Q)∥x∥2⩽x⊤Qx⩽λmax(Q)∥x∥2λmin(Q−1)∥x∥2⩽x⊤Q−1x⩽λmax(Q−1)∥x∥2
- 收敛性:
对于最速下降法,对于任意的初始点 x(0) ,都有 x(k)→x∗ 。
固定步长梯度法
对于所有步长 αk=α∈R,迭代公式为:x(k+1)=x(k)−αg(k)
0<α<λmax(Q)2
时, x(k)→x∗ 。
收敛率
- 定义8.1(收敛阶):
存在一个序列 {x(k)}, 能够收敛到 x∗, 即 limk→∞x(k)−x∗=0 。 如果
0<k→∞lim∥x(k)−x∗∥p∥x(k+1)−x∗∥<∞
则序列 {x(k)} 的收敛阶数为 p, 其中, p∈R。
如果对任意的 p>0 ,上式极限都为 0 ,那么称收敛阶数为 ∞ 。
- 定理8.6:
最速下降法在求解目标函数 f 的极小点时,产生一个收敛的迭代点 {x(k)}, 该序列在最坏情况下的收敛阶数为 1 。也就是说, 存在一个目标函数 f 和某 始点 x(0), 能够使得 {x(k)} 的收敛阶数为 1 。
牛顿法
-
思路:最速下降法只用到了目标函数的一阶导数,如果能够在迭代方法中引入高阶导数,其效率可能优于最速下降法。
- 首先构造一个二次型函数,其与目标函数在该点的一阶和二阶导数相等,以此可以作为目标函数的近似表达式。
- 求改二次型函数的极小点,以此作为下一次迭代的起始点。
- 重复上述过程,直到满足迭代条件后退出。
-
迭代操作
目标函数 f:Rn→R 二阶连续可微, 将函数 f 在点 x(k) 处进行泰勒展开, 可得到二次型近似函数:
f(x)≈f(x(k))+(x−x(k))⊤g(k)+21(x−x(k))⊤ˉF(x(k))(x−x(k))≜q(x)
其中 g(k)=∇f(x(k)),局部极小点的一阶必要条件
0=∇q(x)=g(k)+F(x(k))(x−x(k))
如果 F(x(k))>0, 函数 q 的极小点为x(k+1)=x(k)−F(x(k))−1g(k)
牛顿法收敛性分析
- 如果初始点靠近极大(小)点,那么牛顿法将具有非常好的收敛性,如果初始点离极大(小)点较远,牛顿法并不一定收敛。
修正牛顿法
- 带步长牛顿法:x(k+1)=x(k)−αkF(x(k))−1g(k) , αk=argminf(x(k)−αF(x(k))−1g(k))
牛顿法求解非线性最小二乘
- 非线性最小二乘问题 minimize∑i=1m(ri(x))2 其中, ri:Rn→R,i=1,⋯,m 为给定的函数。
- 分析:令 r=[r1,⋯,rm]⊤, 可将目标函数写为 f(x)=r(x)⊤r(x) 。 计算函数 f 的梯度和黑塞矩阵
(∇f(x))j=∂xj∂f(x)=2i=1∑mri(x)∂xj∂ri(x)
r 的雅可比矩阵为:
J(x)=∂x1∂r1(x)⋮∂x1∂rm(x)⋯⋯∂xn∂r1(x)⋮∂xn∂rm(x)
函数 f 的梯度可以表示为 ∇f(x)=2J(x)⊤r(x)
函数 f 的黑塞矩阵F(x)=2(J(x)⊤J(x)+S(x))
S(x)k,j=∑i=1mri(x)∂xk∂xj∂2ri(x)
牛顿法求解非线性最小二乘迭代公式:
x(k+1)=x(k)−(J(x)⊤J(x)+S(x))−1J(x)⊤r(x)
矩阵 S(x) 包含函数 r 的二阶导数, 其中的元素都很小
高斯一牛顿法求解非线性最小二乘迭代公式:即略掉 S(x) 项
x(k+1)=x(k)−(J(x)⊤J(x))−1J(x)⊤r(x)
共轭类方法
方法特性:
- 对于n维二次型问题,能够在n步之内得到结果。
- 共轭梯度法不需要计算黑塞矩阵
- 不需要储存n·n矩阵,也不需要对矩阵进行求逆
- 计算速度介于最速下降法和牛顿法之间。
定义10.1(共轭):
Q 为 n×n 的对称实矩阵, 对于方向 d(0),d(1),d(2),⋯,d(m), 如果对 于所有的 i=j, 有 d(i)Qd(j)=0, 则称它们是关于 Q 共轭的
引理10.1:
Q 为 n×n 的对称正定矩阵, 如果方向 d(0),d(1),⋯,d(k)∈Rn,k⩽n− 1 非零, 且是关于 Q 共轭的, 那么它们是线性无关的。
共轭方向法
针对 n 维二次型函数的最小化问题:
f(x)=21x⊤Qx−x⊤b
其中, Q=Q⊤>0,x∈Rn 。注意, 由于 Q>0, 因此函数 f 有一个全局极小点, 可通过求解 Qx=b 得到。
基本的共轭方向算法
给定初始点 x(0) 和一组关于 Q 共轭的方向 d(0),d(1),⋯, d(n−1)
g(k)αkx(k+1)=∇f(x(k))=Qx(k)−b=−d(k)⊤Qd(k)g(k)⊤d(k)=x(k)+αkd(k)
定理10.1:
对于任意的初始点 x(0),基本的共轭方向算法都能在 n 次迭代之内收敛到唯一的全局极小点 x∗ ,即 x(n)=x∗ 。
引理10.2(精确步长):
在共轭方向算法中, 对于所有 k,0⩽k⩽n−1,0⩽i⩽k, 都有
g(k+1)⊤d(i)=0
扩展子空间定理
共轭方向法满足 f(x(k+1))=minαf(x(k)+αd(k)),
而且还能满足 f(x(k+1))=mina0,⋯,akf(x(0)+∑i=0kaid(i))
记 Vk=x(0)+span[d(0),d(1),⋯,d(k)] 则有 f(x(k+1))=minx∈νkf(x)
共轭梯度法
特性:
共轭方向法虽然计算效率高,但是需要提供一组 Q 共轭方向,共轭梯度法不需要预先给定 Q 共轭方向,而是随着迭代的进行不断产生 Q共轭方向。利用上一个搜索方向和目标函数在当前已经产生的搜索方向组成 Q 共轭方向。
计算步骤:
- 令 k=0; 选择初始值 x(0) 。
- 计算 g(0)=∇f(x(0)), 如果 g(0)=0, 停止迭代; 否则, 令 d(0)=−g(0) 。
- 计算 αk=−d(k)⊤Qd(k)g(k)⊤d(k) 。
- 计算 x(k+1)=x(k)+αkd(k) 。
- 计算 g(k+1)=∇f(x(k+1)), 如果 g(k+1)=0, 停止迭代。
- 计算 βk=d(k)⊤Qd(k)g(k+1)⊤Qd(k) 。
- 计算 d(k+1)=−g(k+1)+βkd(k) 。
- 令 k:=k+1, 回到第 3 步。
共轭梯度法实例
例: f(x1,x2,x3)=23x12+2x22+23x32+x1x3+2x2x3−3x1−x3 用共轭梯度法求其极小点, 初始点为 x(0)=[0,0,0]⊤ 。
Step 1: Step 2: Step 3: g(0)=[−3,0,−1]⊤d(0)=−g(0)α0=−d(0)⊤Qd(0)g(0)⊤d(0)=3610x(1)=x(0)+α0d(0)=[0.8333,0,0.2778]⊤g(1)=∇f(x(1))=[−0.2222,0.5556,0.6667]⊤β0=d(0)⊤Qd(0)g(1)⊤Qd(0)=0.08025d(1)=−g(1)+β0d(0)=[0.4630,−0.5556,−0.5864]⊤α1=−d(1)⊤Qd(1)g(1)⊤d(1)=0.2187x(2)=x(1)+α1d(1)=[0.9346,−0.1215,0.1495]⊤g(2)=∇f(x(2))=[−0.04673,−0.1869,0.1402]⊤β1=d(1)⊤Qd(1)g(2)⊤Qd(1)=0.07075d(2)=−g(2)+β1d(1)=[0.07948,0.1476,−0.1817]⊤α2=−d(2)⊤Qd(2)g(2)⊤d(2)=0.8231x(3)=x(2)+α2d(2)=[1.000,0.000,0.000]⊤g(3)=∇f(x(3))=0 即 x∗=x(3)
非线性共轭梯度法
理论上把线性共轭梯度法中的矩阵 Q 换成黑 塞矩阵。但是,对一般的非线性函数,每次迭代都必须重新计算黑 塞矩阵,这需要非常大的运算量
所以需要对其进行修改,消除每次迭代中进行求黑塞矩阵的环节。
βk=d(k)⊤[g(k+1)−g(k)]g(k+1)⊤[g(k+1)−g(k)]
- Polak-Ribiere公式:
将 Hestenes−Stiefel 公式的分母部分展开。
βk=g(k)⊤g(k)g(k+1)⊤[g(k+1)−g(k)]
- Fletcher_reeves公式:
将 Polak−Ribiere 公式的分子部分展开。
βk=g(k)⊤g(k)g(k+1)⊤g(k+1)
对于二次型问题,这三个公式是等价的;但是。当目标函数为一般的非线性函数时,这三个公式并不一致。
拟牛顿法
- 当目标函数为一般性的非线性函数时,牛顿法不能保证能够从任意起始点 x(0) 收敛到函数的极小点。
- 必须计算黑塞矩阵 F(x(k)) 和求解方程 F(x(k))d(k)=−g(k) 上进行一维搜索。
- 黑塞矩阵是非正定或奇异的。
- 拟牛顿法的思路:
设计黑塞矩阵的近似矩阵来代替黑塞矩阵。此时需要用到目标函数值和梯度。令 H0,H1,H2,… 表示对应黑塞矩阵的逆的一系列近似矩阵。
则 Hk 满足:
- 对称正定。(保证函数值迭代后下降)
- 从Hk−1 到 Hk 计算量较小。
- Hk 与对应的黑塞矩阵近似。
分析: 黑塞矩阵性质 (即二阶导):
1维:近似计算 f′′(x(k)):x(k)−x(k−1)f′(x(k))−f′(x(k−1)) f′(x(k−1))≈f′(x(k))+f′′(x(k))(x(k−1)−x(k))
n维: ∇f(x(k−1))≈∇f(x(k))+∇2f(x(k))(x(k−1)−x(k)) 记 Δg(k−1)=∇f(x(k))−∇f(x(k−1)),Δx(k−1)=x(k)−x(k−1)
简化为: Δg(k−1)≈∇2f(x(k))Δx(k−1) 或 [∇2f(x(k))]−1Δg(k−1)≈Δx(k−1)
因此:HkΔg(k−1)=Δx(k−1)
d(k)αkx(k+1)=−Hkg(k)=α⩾0argminf(x(k)+αd(k))=x(k)+αkd(k)
其中, 矩阵 H0,H1,H2,⋯ 是n×n 对称实矩阵,是对黑塞矩阵的逆的近似。目标函数为二次型函数时, 它们必须满足
Hk+1Δg(i)=Δx(i),0⩽i⩽k
其中, Δx(i)=x(i+1)−x(i)=αid(i),Δg(i)=g(i+1)−g(i)=QΔx(i) 。实际上, 拟牛顿法也是 一种共轭方向法, 接下来给出证明。
- 定理 11. 1:
将拟牛顿法应用到二次型问题中, 黑塞矩阵为 Q=Q⊤, 对于 0⩽k<n−1, 有
Hk+1Δg(i)=Δx(i),0⩽i⩽k
其中, Hk+1=Hk+1⊤ 。 如果 αi=0,0⩽i⩽k, 那么 d(0),d(1),⋯,d(k+1) 是 Q 共轭的。
从矩阵 Hk 必须满足的方程来看,Hk 并不能唯一确定,常见的矩阵Hk+1 是通过在矩阵 Hk 上增加一个修正项来得到。(秩1修正、DFP,BFGS)
秩1修正公式
Hk+1=Hk+akz(k)z(k)⊤
Hk+1Δg(k)=Δx(k)
在给定 Hk,Δg(k),Δx(k) 后, 确定 ak 和 z(k)
Hk+1=Hk+Δg(k)⊤(Δx(k)−HkΔg(k))(Δx(k)−HkΔg(k))(Δx(k)−HkΔg(k))⊤
- 令 k=0; 选择初始点 x(0), 任选一个对称正定实矩阵 H0 。
- 如果 g(k)=0, 停止迭代; 否则, 令 d(k)=−Hkg(k) 。
- 计算 αk=α⩾0argminf(x(k)+αd(k))x(k+1)=x(k)+αkd(k)
- 计算 Δx(k)=αkd(k)
Δg(k)Hk+1=g(k+1)−g(k)=Hk+Δg(k)⊤(Δx(k)−HkΔg(k))(Δx(k)−HkΔg(k))(Δx(k)−HkΔg(k))⊤
- 令 k:=k+1, 回到第 2 步。
秩1校正公式实例
例:用秩 1 算法求函数 f 极小点。 f(x1,x2)=x12+21x22+3 初始值为 x(0)=[1,2]⊤,H0=I2
解:
g(k)=[2001]x(k)
Step 1:
d(0)=−g(0)=[−2,−2]⊤α0=α⩾0argminf(x(0)+αd(0))=−d(0)⊤Qd(0)g(0)⊤d(0)=32x(1)=x(0)+α0d(0)=[−31,32]⊤
Step 2:
Δx(0)=α0d(0)=[−34,−34]⊤g(1)=Qx(1)=[−32,32]⊤Δg(0)=g(1)−g(0)=[−38,−34]⊤H1=H0+Δg(0)⊤(Δx(0)−H0Δg(0))(Δx(0)−H0Δg(0))(Δx(0)−H0Δg(0))⊤=[21001]d(1)=−H1g(1)=[31,−32]⊤α1=−d(1)⊤Qd(1)g(1)⊤d(1)=1x(2)=x(1)+α1d(1)=[0,0]⊤g(2)=0, 也就是说 x(2)=x∗
DFP算法
- 秩1算法产生的矩阵 Hk+1 可能非正定, 将导致 d(k+1) 可能不是下降方向(即使是二次型问题)。
- 秩 1 公式的分母如果接近 0 ,会出现计算困难。
寻找 F(x(k))−1 的近似
具体公式:
Δx(k)=αkd(k)Δg(k)=g(k+1)−g(k)Hk+1=Hk+Δx(k)⊤Δg(k)Δx(k)Δx(k)⊤−Δg(k)⊤HkΔg(k)[HkΔg(k)][HkΔg(k)]⊤
使用秩一算法或者 DFP 算法求解二次型问题时,黑塞矩阵为 Q=Q⊤ 有 Hk+1Δg(i=Δx(i) ,0≤i≤k
- DFP算法特性:
- 在 DFP 算法中,只要矩阵 Hk 正定,Hk+1 就一定是正定的。
- 当矩阵 Hk 接近为奇异矩阵时,迭代有时无法展开。
DFP公式实例
例:用 DFP 算法求函数 f(x)=21x⊤[4222]x−x⊤[−11],x∈R2 的极小点。初始点为 x(0)=[0,0]⊤,H0=I2 。
step 1:
g(0)=[1,−1]⊤d(0)=−H0g(0)=[−11]α0=α⩾0argminf(x(0)+αd(0))=1x(1)=x(0)+α0d(0)=[−1,1]⊤g(k)=[4222]x(k)−[−11]
step 2:
g(1)=[4222][−11]−[−11]=[−1−1]Δg(0)=g(1)−g(0)=[−2,0]⊤Δx(0)=x(1)−x(0)=[−1,1]⊤H1=H0+Δx(0)⊤Δg(0)Δx(0)Δx(0)⊤−Δg(0)⊤H0Δg(0)(H0Δg(0))(H0Δg(0))⊤=[21−21−2123]d(1)=−H1g(1)=[0,1]⊤α1=α⩾0argminf(x(1)+αd(1))=21x(2)=x(1)+α1d(1)=[−1,3/2]⊤
函数 f 为两变量二次型函数, x(2) 就是极小点 x∗。可以验证 d(0) 和 d(1) 是 Q 共轭方向。
BFSGS公式
- 旧思路:根据 Δg(i)=QΔx(i),0⩽i⩽k
黑塞矩阵逆矩阵的近似矩阵需要满足以下条件:
Hk+1Δg(i)=Δx(i),0⩽i⩽k
基于上述等式, 可以构造黑塞矩阵逆矩阵 Q−1 的近似矩阵的更新公式。
如秩1公式和DFP公式都是据此而来的。
- 新思路: 除了构造 Q−1 的近似矩阵, 还可以构造矩阵 Q 的近似矩阵
令矩阵 Bk 表示在第 k 次迭代中关于矩阵 Q 的估计
则 Bk+1 应该满足 Δg(i)=Bk+1Δx(i),0⩽i⩽k
由DFP公式: Hk+1DFP=Hk+Δx(k)⊤Δg(k)Δx(k)Δx(k)⊤−Δg(k)⊤HkΔg(k)HkΔg(k)Δg(k)⊤Hk
用对称性,得BFGS公式
Bk+1=Bk+Δg(k)⊤Δx(k)Δg(k)Δg(k)⊤−Δx(k)⊤BkΔx(k)BkΔx(k)Δx(k)⊤Bk
注: DFP公式与BFGS公式称为互补的(或对偶的)
(Bk+1)−1=(Bk+Δg(k)⊤Δx(k)Δg(k)Δg(k)T−Δx(k)⊤BkΔx(k)BkΔx(k)Δx(k)⊤Bk)−1
原理:应用两次 Sherman-Morison公式 即可