这一段时间学习优化理论中的一些经典方法,写下一些要点和体会,帮助自己和有需要的人理解最优化方法。本文介绍了如下优化方法的原理和计算步骤:
- 最速下降法。
- 牛顿法。
- 共轭梯度法。
- 变尺度法(变度量法)。
- 线性最小二乘法。
- 非线性最小二乘问题的高斯牛顿法和Levenberg-Marquardt法。
- 坐标轮换法。
1.基础知识
首先来看无约束最优化问题:
min
f
(
x
)
(1)
\min f(x) \tag{1}
minf(x)(1)
其中函数
f
:
R
n
→
R
f:R^n\rightarrow R
f:Rn→R.求解此问题的方法方法分为两大类:最优条件法和迭代法。
所谓的最优条件法,是指当函数存在解析形式,能够通过最优性条件求解出显式最优解。对于无约束最优化问题,如果
f
(
x
)
f(x)
f(x)在最优点
x
ˉ
\bar{x}
xˉ附近可微,那么
x
ˉ
\bar{x}
xˉ是局部极小点的必要条件为:
∇
f
(
x
ˉ
)
=
0
(2)
\nabla f(\bar{x})=0 \tag{2}
∇f(xˉ)=0(2)
我们常常就是通过这个必要条件去求取可能的极小值点,再验证这些点是否真的是极小值点。当上式方程可以求解的时候,无约束最优化问题基本就解决了。实际中,这个方程往往难以求解。这就引出了第二大类方法:迭代法。
迭代法,也称为“搜索”法,主要思想是通过简单的运算构造点列,逐渐逼近问题的最优解。这里说的“点”是多维空间中的点,也称为“向量”。还有少部分算法通过构造点的集合来逼近问题的最优解,如单纯形调优法。
用于求解无约束最优化问题的方法可以分为解析法和直接法两大类。解析法在构造迭代公式的过程中往往使用了泰勒展开来作近似或者推导,因此迭代步骤中含有梯度 ∇ f ( x ) \nabla f(x) ∇f(x)或黑塞(Hessian)矩阵 ∇ 2 f ( x ) \nabla^2f(x) ∇2f(x),在问题的解析形态较好的情况下使用往往能获得比较快的收敛速度。而直接法则从物理角度思考如何递推,不会用到梯度或者黑塞矩阵,它对问题的解析形态几乎没有要求,只要能计算出函数值即可。当然,它的收敛速度往往难于保证。
2.解析法
求解无约束最优化问题的解析法主要有:最速下降法、牛顿法、共轭梯度法(DFP法)和变尺度法(变度量法)。对于特殊的最小二乘问题,有最小二乘法。这些方法各有千秋,除了最小二乘法,后面的方法都针对前面方法的某个问题做了改进。这些方法的核心就是研究如何确定每一步迭代的方向和步长。
2.1最速下降法
最速下降法(Gradient Descent Method)的基本思想是:因为连续函数沿着负梯度方向的下降速度是最快的(这一结论由梯度和方向导数的定义可以推出),所以每次迭代我们都从当前点出发,沿着负梯度方向前进一个最优步长,可以期望能较快逼近函数的极值。
最速下降法仅有三个步骤:
- 设置初始值。设置迭代起点 x 0 ∈ R n x_0\in R^n x0∈Rn,允许误差 ϵ > 0 \epsilon>0 ϵ>0和迭代变量初值 k ← 0 k\leftarrow0 k←0。
- 检查终止条件。如果 ∣ ∣ ∇ f ( x k ) ∣ ∣ < ϵ ||\nabla f(x_k)||<\epsilon ∣∣∇f(xk)∣∣<ϵ,停止迭代输出 x k x_k xk作为近似最优解;否则转步骤3。
- 迭代,也就是通过一维搜索求下一个迭代点。取搜索方向为负梯度方向 d k = − ∇ f ( x ) d_k=-\nabla f(x) dk=−∇f(x),求 λ k \lambda_k λk使得 f ( x k + λ k d k ) = min λ ≥ 0 f ( x k + λ d k ) (3) f(x_k+\lambda_k d_k)=\min_{\lambda\geq0} f(x_k+\lambda d_k) \tag{3} f(xk+λkdk)=λ≥0minf(xk+λdk)(3)再令 x k + 1 = x k + λ d k (4) x_{k+1}=x_k+\lambda d_k \tag{4} xk+1=xk+λdk(4)转步骤2。
步骤中隐含了条件:函数
f
(
x
)
f(x)
f(x)必须可微,也就是说函数
f
(
x
)
f(x)
f(x)的梯度必须存在。这个步骤描述是最速下降发最抽象的形式。其中最关键的步骤是求解问题
min
λ
≥
0
f
(
x
k
+
λ
d
k
)
(5)
\min_{\lambda\geq0}f(x_k+\lambda d_k) \tag{5}
λ≥0minf(xk+λdk)(5)
我们写出它的最优性必要条件
d
f
(
x
k
+
λ
d
k
)
d
λ
=
∇
f
(
x
k
+
λ
d
k
)
T
d
k
=
0
(6)
\frac{d f(x_k+\lambda d_k)}{d\lambda}=\nabla f(x_k+\lambda d_k)^T d_k=0 \tag{6}
dλdf(xk+λdk)=∇f(xk+λdk)Tdk=0(6)
当
f
(
x
)
f(x)
f(x)的形式确定,我们可以通过求解这个一元方程来获得迭代步长
λ
\lambda
λ。当此方程形式复杂,解析解不存在,我们就需要使用“一维搜索”来求解
λ
\lambda
λ了。一维搜索是一些数值方法,有0.618法、Fibonacci法、抛物线法等等,这里不详细解释了。
在实际使用中,为了简便,也可以使用一个预定义的常数而不用一维搜索来确定步长 λ \lambda λ。这时步长的选择往往根据经验或者通过试算来确定。步长过小则收敛慢,步长过大可能震荡而不收敛。
最速下降法是最基本的迭代优化方法。它最简单,最基础,通常是收敛的。可以证明,最速下降法是一阶收敛的,往往需要多次迭代才能接近问题最优解。这是它的不足。
2.2牛顿法
牛顿法是从函数的二阶泰勒展开式推导而来,其思想就是利用目标函数二阶近似的解去逼近目标函数的解。将函数
f
f
f在迭代点
x
k
x_k
xk附近作二阶泰勒展开,有
f
(
x
)
≈
ϕ
(
x
)
=
f
(
x
k
)
+
∇
f
(
x
k
)
T
(
x
−
x
k
)
+
1
2
(
x
−
x
k
)
T
∇
2
f
(
x
k
)
(
x
−
x
k
)
(7)
f(x)\approx \phi(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) \tag{7}
f(x)≈ϕ(x)=f(xk)+∇f(xk)T(x−xk)+21(x−xk)T∇2f(xk)(x−xk)(7)
为了计算
ϕ
(
x
)
\phi(x)
ϕ(x)的极小值,令它的梯度为零,即
∇
ϕ
(
x
)
=
0
(8)
\nabla\phi(x)=0 \tag{8}
∇ϕ(x)=0(8)
有
∇
f
(
x
k
)
+
∇
2
f
(
x
k
)
(
x
−
x
k
)
=
0
(9)
\nabla f(x_k)+\nabla^2f(x_k)(x-x_k)=0 \tag{9}
∇f(xk)+∇2f(xk)(x−xk)=0(9)
从而推出
x
=
x
k
−
[
∇
2
f
(
x
k
)
]
−
1
∇
f
(
x
k
)
(10)
x = x_k - [\nabla^2f(x_k)]^{-1}\nabla f(x_k) \tag{10}
x=xk−[∇2f(xk)]−1∇f(xk)(10)
我们将这里的
x
x
x作为第
k
+
1
k+1
k+1次迭代的估值,就有了迭代公式
x
k
+
1
=
x
k
−
[
∇
2
(
f
x
k
)
]
−
1
∇
f
(
x
k
)
(11)
x_{k+1}=x_k-[\nabla^2(fx_k)]^{-1}\nabla f(x_k) \tag{11}
xk+1=xk−[∇2(fxk)]−1∇f(xk)(11)
牛顿法(Newton Method)利用了函数的二阶信息,即黑塞矩阵,来加速迭代收敛。具有二阶收敛速度是它的显著优势。牛顿法的步骤为:
- 设置初始值。给定迭代初值 x 0 ∈ R n x_0\in R^n x0∈Rn, ϵ > 0 \epsilon>0 ϵ>0,令 k ← 0 k \leftarrow 0 k←0。
- 检查终止条件。如果 ∣ ∣ ∇ f ( x ) ∣ ∣ < ϵ ||\nabla f(x)||<\epsilon ∣∣∇f(x)∣∣<ϵ,迭代终止, x k x_k xk为近似最优解;否则,转步骤3。
- 迭代计算。取迭代方向 d k = − [ ∇ 2 f ( x k ) ] − 1 ∇ f ( x k ) (12) d_k=-[\nabla^2f(x_k)]^{-1}\nabla f(x_k) \tag{12} dk=−[∇2f(xk)]−1∇f(xk)(12)令 x k + 1 = x k + d k k ← k + 1 (13) x_{k+1}=x_k+d_k \\ k \leftarrow k+1 \tag{13} xk+1=xk+dkk←k+1(13)转步骤2。
牛顿法要求初始点在最优点附近(泰勒展开的前提就是在邻域内),否则可能不收敛。为了使得牛顿法能够全局收敛,提出了阻尼牛顿法(Damped Newton Method)。阻尼牛顿法的改进在于每次的搜索步长不固定为1,而是通过一维搜索来确定步长。其步骤如下:
- 设置初始值。给定迭代初值 x 0 ∈ R n x_0\in R^n x0∈Rn, ϵ > 0 \epsilon>0 ϵ>0,令 k ← 0 k \leftarrow 0 k←0。
- 检查终止条件。如果 ∣ ∣ ∇ f ( x ) ∣ ∣ < ϵ ||\nabla f(x)||<\epsilon ∣∣∇f(x)∣∣<ϵ,迭代终止, x k x_k xk为近似最优解;否则,转步骤3。
- 取迭代方向 d k = − [ ∇ 2 f ( x k ) ] − 1 ∇ f ( x k ) (14) d_k=-[\nabla^2f(x_k)]^{-1}\nabla f(x_k) \tag{14} dk=−[∇2f(xk)]−1∇f(xk)(14)
- 进行一维搜索确定步长 λ k \lambda_k λk使得 f ( x k + λ k d k ) = min λ ≥ 0 f ( x k + λ d k ) (15) f(x_k+\lambda_k d_k)=\min_{\lambda\geq0}f(x_k+\lambda d_k) \tag{15} f(xk+λkdk)=λ≥0minf(xk+λdk)(15)令 x k + 1 = x k + λ d k k ← k + 1 (16) x_{k+1}=x_k+\lambda d_k\\k\leftarrow k+1 \tag{16} xk+1=xk+λdkk←k+1(16)转步骤2。
可以证明,当 f ( x ) f(x) f(x)具有二阶连续偏导数且 ∇ 2 f ( x ) \nabla^2f(x) ∇2f(x)正定,阻尼牛顿法是全局收敛的。
2.3共轭梯度法
最速下降法收敛速度慢,牛顿法虽然收敛速度快但是需要计算黑塞矩阵的逆矩阵,计算复杂度比较高。有的时候我们希望有一种收敛快且不用计算二阶导数的方法,共轭梯度法(Conjugate Gradient Method)应运而生了。当然,计算量减小的代价是理论推导更加繁琐。
要了解共轭梯度法,首先要了解共轭的概念。设有两个
n
n
n维向量
d
1
d_1
d1、
d
2
d_2
d2和一个
n
n
n阶正定矩阵
Q
Q
Q,如果它们满足
d
1
T
Q
d
2
=
0
(17)
d_1^T Q d_2=0 \tag{17}
d1TQd2=0(17)
则称向量
d
1
d_1
d1、
d
2
d_2
d2关于
Q
Q
Q共轭。当
Q
Q
Q为单位阵时,共轭就成了正交。可将共轭看成是正交概念的推广。共轭的概念还可以推广到含有m个向量的向量组,即
d
i
T
Q
d
j
=
0
,
∀
i
,
j
=
1
,
.
.
.
m
,
i
≠
j
(18)
d_i^TQd_j=0,\forall i,j=1,...m,i\neq j \tag{18}
diTQdj=0,∀i,j=1,...m,i=j(18)
我们用正定二次函数 f ( x ) = 1 2 x T Q x + b T x + c (19) f(x)=\frac{1}{2}x^TQx+b^Tx+c \tag{19} f(x)=21xTQx+bTx+c(19)来推导,最后推广到一般问题。注意这里的 Q ∈ R n × n Q\in R^{n\times n} Q∈Rn×n为正定矩阵, b ∈ R n b\in R^n b∈Rn, c ∈ R c\in R c∈R。
我们从
x
0
x_0
x0出发,沿着某个下降方向
d
0
d_0
d0作一维搜索得到下一个迭代点
x
1
x_1
x1,那么有
f
(
x
1
)
=
f
(
x
0
+
λ
0
d
0
)
=
min
λ
≥
0
f
(
x
0
+
λ
d
0
)
(20)
f(x_1)=f(x_0+\lambda_0 d_0)=\min_{\lambda\geq0}f(x_0+\lambda d_0) \tag{20}
f(x1)=f(x0+λ0d0)=λ≥0minf(x0+λd0)(20)
它的最优性条件是对
λ
\lambda
λ求导,于是有
∇
f
(
x
0
+
λ
0
d
0
)
T
d
0
=
0
(21)
\nabla f(x_0+\lambda_0d_0)^Td_0=0 \tag{21}
∇f(x0+λ0d0)Td0=0(21)
即
∇
f
(
x
1
)
T
d
0
=
0
(22)
\nabla f(x_1)^Td_0=0 \tag{22}
∇f(x1)Td0=0(22)
新的迭代方向
d
1
d_1
d1如果取负梯度方向,就成了最速下降法。负梯度方向虽然是当前函数值下降最快的方向,却未必是直指函数最优点的方向。设函数最优点为
x
ˉ
\bar{x}
xˉ,我们希望搜索方向
d
1
d_1
d1能直接指向最优值,即
x
ˉ
=
x
1
+
λ
1
d
1
(23)
\bar{x}=x_1+\lambda_1 d_1 \tag{23}
xˉ=x1+λ1d1(23)
而作为最优点,
x
ˉ
\bar{x}
xˉ应满足
0
=
∇
f
(
x
ˉ
)
=
Q
x
ˉ
+
b
=
Q
(
x
1
+
λ
1
d
1
)
+
b
=
(
Q
x
1
+
b
)
+
λ
1
Q
d
1
=
∇
f
(
x
1
)
+
λ
1
Q
d
1
(24)
\begin{align*} 0=\nabla f(\bar{x})&=Q\bar{x}+b \\ &=Q(x_1+\lambda_1 d_1)+b \\ &=(Qx_1+b)+\lambda_1 Q d_1 \\ &=\nabla f(x_1)+\lambda_1 Q d_1 \end{align*} \tag{24}
0=∇f(xˉ)=Qxˉ+b=Q(x1+λ1d1)+b=(Qx1+b)+λ1Qd1=∇f(x1)+λ1Qd1(24)
两边左乘
d
0
T
d_0^T
d0T并注意到
d
0
T
∇
f
(
x
1
)
=
0
d_0^T \nabla f(x_1)=0
d0T∇f(x1)=0就有
d
0
T
Q
d
1
=
0
(25)
d_0^TQd_1=0 \tag{25}
d0TQd1=0(25)
上式表明,搜索方向
d
1
d_1
d1与
d
0
d_0
d0是关于
Q
Q
Q共轭的。
可以证明,如果能构造出与 Q Q Q共轭的向量组,至多通过 n n n次迭代就可以求得正定二次函数最优化问题的最优解。由此产生一类称为共轭方向法的方法,当我们利用梯度信息来构造与 Q Q Q共轭的向量组,就产生了共轭梯度法。
下面推导如何利用梯度构造与
Q
Q
Q共轭的向量组。首先初始搜索方向取为负梯度方向,即
d
0
=
−
∇
f
(
x
0
)
(26)
d_0=-\nabla f(x_0) \tag{26}
d0=−∇f(x0)(26)
新的搜索方向由负梯度方向和上一次搜索方向的线型组合产生,即
d
k
+
1
=
−
∇
f
(
x
k
+
1
)
+
α
k
d
k
,
k
=
0
,
1
,
.
.
.
,
n
−
2
(27)
d_{k+1}=-\nabla f(x_{k+1})+\alpha_k d_k, k=0,1,...,n-2 \tag{27}
dk+1=−∇f(xk+1)+αkdk,k=0,1,...,n−2(27)
其中
α
k
\alpha_k
αk待定。我们利用
d
k
+
1
d_{k+1}
dk+1与
d
k
d_k
dk关于
Q
Q
Q共轭作为约束条件,有
0
=
d
k
+
1
T
Q
d
k
=
(
−
∇
f
(
x
k
+
1
)
T
+
α
k
d
k
T
)
Q
d
k
=
−
∇
f
(
x
k
+
1
)
T
Q
d
k
+
α
k
d
k
T
Q
d
k
(28)
\begin{align} 0&=d_{k+1}^TQd_k\\ &=(-\nabla f(x_{k+1})^T+\alpha_k d_k^T)Qd_k\\ &=-\nabla f(x_{k+1})^TQd_k+\alpha_k d_k^TQd_k \end{align} \tag{28}
0=dk+1TQdk=(−∇f(xk+1)T+αkdkT)Qdk=−∇f(xk+1)TQdk+αkdkTQdk(28)
从而解得
α
k
=
∇
f
(
x
k
+
1
)
T
Q
d
k
d
k
T
Q
d
k
(29)
\alpha_k=\frac{\nabla f(x_{k+1})^TQd_k}{d_k^TQd_k} \tag{29}
αk=dkTQdk∇f(xk+1)TQdk(29)
总结一下,我们得到了n个搜索方向的递推公式如下:
{
d
0
=
−
∇
f
(
x
0
)
d
k
+
1
=
−
∇
f
(
x
k
+
1
)
+
α
k
d
k
α
k
=
−
∇
f
(
x
k
+
1
)
T
Q
d
k
d
k
T
Q
d
k
(30)
\left\{ \begin{aligned} d_0 & =-\nabla f(x_0) \\ d_{k+1}&=-\nabla f(x_{k+1})+\alpha_k d_k \\ \alpha_k &=-\frac{\nabla f(x_{k+1})^TQd_k}{d_k^TQd_k} \end{aligned} \right. \tag{30}
⎩
⎨
⎧d0dk+1αk=−∇f(x0)=−∇f(xk+1)+αkdk=−dkTQdk∇f(xk+1)TQdk(30)
可以证明,上式得到的
d
k
{d_k}
dk是一组关于
Q
Q
Q共轭的方向向量。
将此公式向一般无约束最优化问题推广,必须消去
Q
Q
Q。这里直接给出三个消去
Q
Q
Q后
α
k
\alpha_k
αk的表达式:
(1)FR公式(Fletcher-Reeves)
α
k
=
∣
∣
∇
f
(
x
k
+
1
)
∣
∣
2
∣
∣
∇
f
(
x
k
)
∣
∣
2
(31)
\alpha_k=\frac{||\nabla f(x_{k+1})||^2}{||\nabla f(x_k)||^2} \tag{31}
αk=∣∣∇f(xk)∣∣2∣∣∇f(xk+1)∣∣2(31)
(2)DM公式(Dixon-Myers)
α
k
=
−
∣
∣
∇
f
(
x
k
+
1
)
∣
∣
2
d
k
T
∇
f
(
x
k
)
(32)
\alpha_k=-\frac{||\nabla f(x_{k+1})||^2}{d_k^T\nabla f(x_k)} \tag{32}
αk=−dkT∇f(xk)∣∣∇f(xk+1)∣∣2(32)
(3)PRP公式(Polak-Ribiere-Polyak)
α
k
=
∇
f
(
x
k
+
1
)
T
[
∇
f
(
x
k
+
1
)
−
∇
f
(
x
k
)
]
∣
∣
∇
f
(
x
k
)
∣
∣
2
=
∣
∣
∇
f
(
x
k
+
1
)
∣
∣
2
−
∇
f
(
x
k
+
1
)
T
∇
f
(
x
k
)
∣
∣
∇
f
(
x
k
)
∣
∣
2
(33)
\begin{align} \alpha_k&=\frac{\nabla f(x_{k+1})^T[\nabla f(x_{k+1})-\nabla f(x_k)]}{||\nabla f(x_k)||^2} \\ &=\frac{||\nabla f(x_{k+1})||^2 - \nabla f(x_{k+1})^T \nabla f(x_k)}{||\nabla f(x_k)||^2} \end{align} \tag{33}
αk=∣∣∇f(xk)∣∣2∇f(xk+1)T[∇f(xk+1)−∇f(xk)]=∣∣∇f(xk)∣∣2∣∣∇f(xk+1)∣∣2−∇f(xk+1)T∇f(xk)(33)
当问题为正定二次函数优化问题时,这三个公式是等价的。对非二次函数最优化问题,它们将产生不同的搜索方向。
共轭梯度法的步骤总结如下:
- 设置初始值。迭代初始点 x 0 x_0 x0和允许误差 ϵ \epsilon ϵ。
- 检查终止条件。如果 ∣ ∣ ∇ f ( x 0 ) ∣ ∣ < ϵ ||\nabla f(x_0)||<\epsilon ∣∣∇f(x0)∣∣<ϵ,迭代终止并输出 x 0 x_0 x0;否则转步骤3。
- 构造初始搜索方向。 d 0 = − ∇ f ( x 0 ) (34) d_0=-\nabla f(x_0) \tag{34} d0=−∇f(x0)(34) 并令 k = 0 k=0 k=0
- 进行一维搜索。求 λ k \lambda_k λk使得 f ( x k + λ k d k ) = min λ ≥ 0 f ( x k + λ d k ) (35) f(x_k+\lambda_k d_k)=\min_{\lambda\geq0}f(x_k+\lambda d_k) \tag{35} f(xk+λkdk)=λ≥0minf(xk+λdk)(35)并令 x k + 1 = x k + λ k d k (36) x_{k+1}=x_k+\lambda_k d_k \tag{36} xk+1=xk+λkdk(36)
- 检查终止条件。如果 ∣ ∣ ∇ f ( x k + 1 ) ∣ ∣ < ϵ ||\nabla f(x_{k+1})||<\epsilon ∣∣∇f(xk+1)∣∣<ϵ,迭代终止并输出 x k + 1 x_{k+1} xk+1;否则转步骤6。
- 检查迭代次数。如果 k + 1 = n k+1=n k+1=n,令 x 0 = x n (37) x_0=x_n \tag{37} x0=xn(37)转步骤3;否则转步骤7。
- 构造共轭方向。令
d
k
+
1
=
−
∇
f
(
x
k
+
1
)
+
α
k
d
k
α
k
=
∣
∣
∇
f
(
x
k
+
1
)
∣
∣
2
∣
∣
∇
f
(
x
k
)
∣
∣
2
(38)
\begin{align} d_{k+1}&=-\nabla f(x_{k+1})+\alpha_k d_k \\ \alpha_k &=\frac{||\nabla f(x_{k+1})||^2}{||\nabla f(x_k)||^2} \end{align} \tag{38}
dk+1αk=−∇f(xk+1)+αkdk=∣∣∇f(xk)∣∣2∣∣∇f(xk+1)∣∣2(38)
这里 α k \alpha_k αk的计算取了FR公式,也可以选取另外两个公式。再令 k = k + 1 k=k+1 k=k+1,转步骤4。
共轭梯度法在 n n n次迭代后将无法产生新的共轭方向(因为 R n R^n Rn中的共轭方向至多只能有 n n n个),故将第 n n n次迭代点 x n x_n xn作为新的起点,重新进行一轮共轭梯度迭代。共轭梯度法比最速下降法的收敛条件更弱。当 f f f具有一阶连续偏导数且有极值点时,共轭梯度法就是收敛的。对于二次函数最优化,理论上共轭梯度法 n n n次迭代即可收敛。在一定条件下,共轭梯度法具有二次收敛速度。
共轭梯度法的Matlab实现请参考我的另一篇文章:Matlab实现FR共轭梯度法。
2.4变尺度法
最速下降法和阻尼牛顿法的迭代公式可以写成
{
x
k
+
1
=
x
k
+
λ
k
d
k
d
k
=
−
G
k
∇
f
(
x
k
)
(39)
\left\{ \begin{aligned} x_{k+1}=x_k+\lambda_k d_k \\ d_k=-G_k \nabla f(x_k) \end{aligned} \tag{39} \right.
{xk+1=xk+λkdkdk=−Gk∇f(xk)(39)
当
G
k
G_k
Gk取单位阵时是最速下降法,当
G
k
=
[
∇
2
f
(
x
k
)
]
−
1
G_k=[\nabla^2 f(x_k)]^{-1}
Gk=[∇2f(xk)]−1时是阻尼牛顿法。事实上,我们把上式确定的迭代法统称为拟牛顿法(Quasi-Newton Method)。因此改进的思路可以让
G
k
G_k
Gk近似
[
∇
2
f
(
x
k
)
]
−
1
[\nabla^2 f(x_k)]^{-1}
[∇2f(xk)]−1,从而避免计算二阶导数和矩阵求逆,既简化计算又保持快速收敛性。
我们同样利用二阶泰勒展开来考虑如何近似。将
f
f
f在点
x
k
+
1
x_{k+1}
xk+1处作二阶泰勒展开,有
f
(
x
)
≈
f
(
x
k
+
1
)
+
∇
f
(
x
k
+
1
)
T
(
x
−
x
k
+
1
)
+
1
2
(
x
−
x
k
+
1
)
T
∇
2
f
(
x
k
+
1
)
(
x
−
x
k
+
1
)
(40)
\begin{aligned} f(x) \approx & f(x_{k+1})+\nabla f(x_{k+1})^T(x-x_{k+1})+ \\ & \frac{1}{2}(x-x_{k+1})^T\nabla^2 f(x_{k+1})(x-x_{k+1}) \end{aligned} \tag{40}
f(x)≈f(xk+1)+∇f(xk+1)T(x−xk+1)+21(x−xk+1)T∇2f(xk+1)(x−xk+1)(40)
对上式两边求梯度,有
∇
f
(
x
)
≈
∇
f
(
x
k
+
1
)
+
∇
2
f
(
x
k
+
1
)
(
x
−
x
k
+
1
)
(41)
\nabla f(x) \approx \nabla f(x_{k+1}) + \nabla^2f(x_{k+1})(x-x_{k+1}) \tag{41}
∇f(x)≈∇f(xk+1)+∇2f(xk+1)(x−xk+1)(41)
令
x
=
x
k
x=x_k
x=xk,有
∇
f
(
x
k
)
≈
∇
f
(
x
k
+
1
)
+
∇
2
f
(
x
k
+
1
)
(
x
k
−
x
k
+
1
)
(42)
\nabla f(x_k) \approx \nabla f(x_{k+1}) + \nabla^2f(x_{k+1})(x_k-x_{k+1}) \tag{42}
∇f(xk)≈∇f(xk+1)+∇2f(xk+1)(xk−xk+1)(42)
可得如下近似关系
[
∇
2
f
(
x
k
+
1
)
]
−
1
[
∇
f
(
x
k
+
1
)
−
∇
f
(
x
k
)
]
=
x
k
+
1
−
x
k
(43)
[\nabla^2 f(x_{k+1})]^{-1}[\nabla f(x_{k+1}) - \nabla f(x_k)]=x_{k+1}-x_k \tag{43}
[∇2f(xk+1)]−1[∇f(xk+1)−∇f(xk)]=xk+1−xk(43)
我们可以令
G
k
+
1
G_{k+1}
Gk+1满足此关系式,即
G
k
+
1
[
∇
f
(
x
k
+
1
)
−
∇
f
(
x
k
)
]
=
x
k
+
1
−
x
k
(44)
G_{k+1}[\nabla f(x_{k+1}) - \nabla f(x_k)]=x_{k+1}-x_k \tag{44}
Gk+1[∇f(xk+1)−∇f(xk)]=xk+1−xk(44)
从而近似阻尼牛顿法。为了简化后续推导,记
Δ
x
k
=
x
k
+
1
−
x
k
Δ
g
k
=
∇
f
(
x
k
+
1
)
−
∇
f
(
x
k
)
Δ
G
k
=
G
k
+
1
−
G
k
(45)
\begin{aligned} \Delta x_k &=x_{k+1}-x_k \\ \Delta g_k &= \nabla f(x_{k+1}) - \nabla f(x_k) \\ \Delta G_k &= G_{k+1}-G_k \end{aligned} \tag{45}
ΔxkΔgkΔGk=xk+1−xk=∇f(xk+1)−∇f(xk)=Gk+1−Gk(45)
则上式可以改写为
Δ
G
k
Δ
g
k
=
Δ
x
k
−
G
k
Δ
g
k
(46)
\Delta G_k \Delta g_k = \Delta x_k - G_k \Delta g_k \tag{46}
ΔGkΔgk=Δxk−GkΔgk(46)
为了求出
Δ
G
k
\Delta G_k
ΔGk的表达式,我们采用待定系数法,设
Δ
G
k
=
Δ
x
k
p
k
T
−
G
k
Δ
g
k
q
k
T
(47)
\Delta G_k = \Delta x_k p_k^T - G_k \Delta g_k q_k^T \tag{47}
ΔGk=ΔxkpkT−GkΔgkqkT(47)
其中
p
k
,
q
k
∈
R
n
p_k, q_k \in R^n
pk,qk∈Rn为待定向量。两边右乘
Δ
g
k
\Delta g_k
Δgk有
Δ
G
k
Δ
g
k
=
(
p
k
T
Δ
g
k
)
Δ
x
k
+
(
q
k
T
Δ
g
k
)
G
k
Δ
g
k
(48)
\Delta G_k \Delta g_k = (p_k^T \Delta g_k) \Delta x_k + (q_k^T \Delta g_k) G_k \Delta g_k \tag{48}
ΔGkΔgk=(pkTΔgk)Δxk+(qkTΔgk)GkΔgk(48)
与前面推导的近似式子对比,则
p
k
T
Δ
g
k
=
1
q
k
T
Δ
g
k
=
1
(49)
p_k^T \Delta g_k=1 \\ q_k^T \Delta g_k = 1 \tag{49}
pkTΔgk=1qkTΔgk=1(49)
考虑到
G
k
G_k
Gk近似于黑塞矩阵的逆(为对称矩阵),也应该为对阵矩阵,可构造性地设
p
k
=
α
k
Δ
x
k
q
k
=
β
k
G
k
Δ
g
k
(50)
\begin{aligned} p_k&=\alpha_k \Delta x_k \\ q_k&=\beta_kG_k \Delta g_k \end{aligned} \tag{50}
pkqk=αkΔxk=βkGkΔgk(50)
解得
α
k
=
1
Δ
x
k
T
Δ
g
k
β
k
=
1
Δ
g
k
G
k
T
Δ
g
k
(51)
\begin{aligned} \alpha_k&=\frac{1}{\Delta x_k^T\Delta g_k} \\ \beta_k&=\frac{1}{\Delta g_k G_k^T \Delta g_k} \end{aligned} \tag{51}
αkβk=ΔxkTΔgk1=ΔgkGkTΔgk1(51)
回代可得
Δ
G
k
=
Δ
x
k
Δ
x
k
T
Δ
x
k
T
Δ
g
k
−
G
k
Δ
g
k
Δ
g
k
T
G
k
Δ
g
k
T
G
k
Δ
g
k
(52)
\Delta G_k=\frac{\Delta x_k \Delta x_k^T}{\Delta x_k^T \Delta g_k}-\frac{G_k\Delta g_k \Delta g_k^T G_k}{\Delta g_k^T G_k \Delta g_k} \tag{52}
ΔGk=ΔxkTΔgkΔxkΔxkT−ΔgkTGkΔgkGkΔgkΔgkTGk(52)
上述公式由Davidon提出,被Fletcher和Powell改进得到,称为DFP公式。应用DFP公式的拟牛顿法称为DFP法或者变尺度法(变度量法)。
最后总结一下变尺度法的步骤:
- 选取初始值。初始点 x 0 ∈ R n x_0 \in R^n x0∈Rn,允许误差 ϵ > 0 \epsilon>0 ϵ>0和 G 0 = I n G_0=I_n G0=In为单位阵。
- 检查终止条件。如果 ∣ ∣ ∇ f ( x 0 ) ∣ ∣ < ϵ ||\nabla f(x_0)||<\epsilon ∣∣∇f(x0)∣∣<ϵ,终止迭代并输出 x 0 x_0 x0;否则转步骤6。
- 构造初始DFP方向。取 d 0 = − ∇ f ( x 0 ) d_0=-\nabla f(x_0) d0=−∇f(x0),令 k = 0 k=0 k=0
- 一维搜索确定步长 λ k \lambda_k λk,使得 f ( x k + λ k d k ) = min λ ≥ 0 f ( x k + λ d k ) (53) f(x_k+\lambda_k d_k)=\min_{\lambda\geq0}f(x_k+\lambda d_k) \tag{53} f(xk+λkdk)=λ≥0minf(xk+λdk)(53)并令 x k + 1 = x k + λ k d k (54) x_{k+1}=x_k+\lambda_k d_k \tag{54} xk+1=xk+λkdk(54)
- 检查终止条件。如果 ∣ ∣ ∇ f ( x k + 1 ) ∣ ∣ < ϵ ||\nabla f(x_{k+1})||<\epsilon ∣∣∇f(xk+1)∣∣<ϵ,终止迭代并输出 x k + 1 x_{k+1} xk+1;否则转步骤6。
- 检查迭代次数。如果 k + 1 = n k+1=n k+1=n,令 x 0 ← x n x_0 \leftarrow x_n x0←xn并转步骤(3);否则转步骤7。
- 构造DFP方向。令 G k + 1 = G k + Δ x k Δ x k T Δ x k T Δ g k − G k Δ g k Δ g k T G k Δ g k T G k Δ g k d k + 1 = − G k + 1 ∇ f ( x k ) (55) \begin{aligned} G_{k+1}&=G_k+\frac{\Delta x_k \Delta x_k^T}{\Delta x_k^T \Delta g_k}-\frac{G_k\Delta g_k \Delta g_k^T G_k}{\Delta g_k^T G_k \Delta g_k} \\ d_{k+1}&=-G_{k+1} \nabla f(x_k) \end{aligned} \tag{55} Gk+1dk+1=Gk+ΔxkTΔgkΔxkΔxkT−ΔgkTGkΔgkGkΔgkΔgkTGk=−Gk+1∇f(xk)(55)令 k ← k + 1 k \leftarrow k+1 k←k+1,转步骤4。
这里每 n n n次迭代就重新初始化可以减少计算误差的积累,有利于快速收敛。
共轭梯度法与牛顿法类似,具有二阶收敛速度。
2.5最小二乘法
若干的目标函数的平方和作为目标函数的优化问题就称为最小二乘问题。这是一类特殊的无约束最优化问题,如下式:
min
s
(
x
)
=
∑
i
=
1
m
f
i
2
(
x
)
(56)
\min s(x)=\sum_{i=1}^{m}f_i^2(x) \tag{56}
mins(x)=i=1∑mfi2(x)(56)
其中
x
∈
R
n
x\in R^n
x∈Rn且
m
≥
n
m\geq n
m≥n。
首先考虑线性最小二乘问题。当
f
i
(
x
)
f_i(x)
fi(x)均为线性函数时,最小二乘问题可以写成矩阵形式:
min
s
(
x
)
=
∣
∣
A
x
−
b
∣
∣
2
(57)
\min s(x)= ||Ax-b||^2 \tag{57}
mins(x)=∣∣Ax−b∣∣2(57)
其中
A
∈
R
m
×
n
A\in R^{m\times n}
A∈Rm×n,
b
∈
R
n
b\in R^n
b∈Rn。一般要求
A
A
A为列满秩矩阵(否则一定有冗余的列,去掉之)。令
0
=
∇
s
(
x
)
=
2
A
T
(
A
x
−
b
)
=
2
(
A
T
A
x
−
A
T
b
)
(58)
\begin{aligned} 0&=\nabla s(x) \\ &=2 A^T(Ax-b) \\ &=2(A^TAx-A^Tb) \end{aligned} \tag{58}
0=∇s(x)=2AT(Ax−b)=2(ATAx−ATb)(58)
注意到当
A
A
A列满秩时,
A
T
A
A^TA
ATA正定且可逆,故可推出最优解
x
ˉ
=
(
A
T
A
)
−
1
A
T
b
(59)
\bar{x}=(A^TA)^{-1}A^Tb \tag{59}
xˉ=(ATA)−1ATb(59)
这就是线性最小二乘问题的最优解。其中
(
A
T
A
)
−
1
A
T
(A^TA)^{-1}A^T
(ATA)−1AT称为矩阵
A
A
A的伪逆(广义逆)或者摩尔彭若斯伪逆(Moore–Penrose pseudoinverse),记为
A
+
A^+
A+。故最优解也可以记为
x
ˉ
=
A
+
b
(60)
\bar{x}=A^+b \tag{60}
xˉ=A+b(60)
接下来考虑非线性最小二乘问题
min
s
(
x
)
=
∑
i
=
1
m
f
i
2
(
x
)
(61)
\min s(x)=\sum_{i=1}^{m}f_i^2(x) \tag{61}
mins(x)=i=1∑mfi2(x)(61)
解决问题的思想与牛顿法相似,就是用目标函数的泰勒展开的最优解去逼近目标函数的最优解。由此产生了下面要介绍的高斯-牛顿法。将
f
i
(
x
)
f_i(x)
fi(x)在点
x
k
x_k
xk处做一阶泰勒展开,有
f
i
(
x
)
≈
f
i
(
x
k
)
+
∇
f
i
(
x
k
)
T
(
x
−
x
k
)
(62)
f_i(x)\approx f_i(x_k)+\nabla f_i(x_k)^T(x-x_k) \tag{62}
fi(x)≈fi(xk)+∇fi(xk)T(x−xk)(62)
则
s
(
x
)
≈
∑
i
=
1
m
[
f
i
(
x
k
)
+
∇
f
i
(
x
k
)
T
(
x
−
x
k
)
]
2
(63)
s(x) \approx \sum_{i=1}^{m}[f_i(x_k)+\nabla f_i(x_k)^T(x-x_k)]^2 \tag{63}
s(x)≈i=1∑m[fi(xk)+∇fi(xk)T(x−xk)]2(63)
这是一个线性最小二乘问题。令
A
=
(
∇
f
1
(
x
k
)
,
.
.
.
,
∇
f
m
(
x
k
)
)
b
=
−
(
f
1
(
x
k
)
,
.
.
.
,
f
m
(
x
k
)
)
T
(64)
A=(\nabla f_1(x_k), ..., \nabla f_m(x_k) ) \\ b=-(f_1(x_k),...,f_m(x_k))^T \tag{64}
A=(∇f1(xk),...,∇fm(xk))b=−(f1(xk),...,fm(xk))T(64)
如果
A
A
A为列满秩矩阵,则
x
−
x
k
=
(
A
T
A
)
−
1
A
T
b
(65)
x-x_k=(A^TA)^{-1}A^Tb \tag{65}
x−xk=(ATA)−1ATb(65)
如果记函数的向量
f
(
x
)
=
(
f
1
(
x
)
,
.
.
.
,
f
m
(
x
)
)
T
(66)
\mathbf{f}(x)=(f_1(x),...,f_m(x))^T \tag{66}
f(x)=(f1(x),...,fm(x))T(66)
则
f
(
x
)
\mathbf{f}(x)
f(x)的雅克比矩阵为
∇
f
(
x
)
=
(
∇
f
1
(
x
)
,
.
.
.
,
∇
f
m
(
x
)
)
(67)
\nabla \mathbf{f}(x)=(\nabla f_1(x), ... , \nabla f_m(x)) \tag{67}
∇f(x)=(∇f1(x),...,∇fm(x))(67)
则上面推导的公式变为
x
−
x
k
=
−
(
∇
f
(
x
k
)
T
∇
f
(
x
k
)
)
−
1
∇
f
(
x
k
)
T
f
(
x
k
)
(68)
x-x_k=-(\nabla \mathbf{f}(x_k)^T\nabla \mathbf{f}(x_k))^{-1}\nabla \mathbf{f}(x_k)^T \mathbf{f}(x_k) \tag{68}
x−xk=−(∇f(xk)T∇f(xk))−1∇f(xk)Tf(xk)(68)
将
x
x
x作为下一次迭代点,有
x
k
+
1
=
x
k
−
(
∇
f
(
x
k
)
T
∇
f
(
x
k
)
)
−
1
∇
f
(
x
k
)
T
f
(
x
k
)
(69)
x_{k+1}=x_k-(\nabla \mathbf{f}(x_k)^T\nabla \mathbf{f}(x_k))^{-1}\nabla \mathbf{f}(x_k)^T \mathbf{f}(x_k) \tag{69}
xk+1=xk−(∇f(xk)T∇f(xk))−1∇f(xk)Tf(xk)(69)
上述公式称为高斯-牛顿迭代公式。可以证明,当迭代初始点接近最优点时,高斯-牛顿法是二阶收敛的。为了解决收敛性问题,模仿阻尼牛顿法,我们可以用一维搜索来确定迭代的步长,从而得到阻尼高斯-牛顿法。它的步骤总结如下:
- 选取初始数据。 x 0 x_0 x0,允许误差 ϵ > 0 \epsilon>0 ϵ>0,令 k = 0 k=0 k=0
- 检查终止条件。如果 ∣ ∣ ∇ f ( x k ) T f ( x k ) ∣ ∣ < ϵ ||\nabla \mathbf{f}(x_k)^T \mathbf{f}(x_k)||<\epsilon ∣∣∇f(xk)Tf(xk)∣∣<ϵ,迭代终止并输出 x k x_k xk;否则转步骤3。
- 计算高斯-牛顿方向 d k = − ( ∇ f ( x k ) T ∇ f ( x k ) ) − 1 ∇ f ( x k ) T f ( x k ) (70) d_k=-(\nabla \mathbf{f}(x_k)^T\nabla \mathbf{f}(x_k))^{-1}\nabla \mathbf{f}(x_k)^T \mathbf{f}(x_k) \tag{70} dk=−(∇f(xk)T∇f(xk))−1∇f(xk)Tf(xk)(70)
- 进行一维搜索求 λ k \lambda_k λk,使得 s ( x k + λ k d k ) = min λ ≥ 0 s ( x k + λ d k ) (71) s(x_k+\lambda_k d_k)=\min_{\lambda\geq0}s(x_k+\lambda d_k) \tag{71} s(xk+λkdk)=λ≥0mins(xk+λdk)(71)并令 x k + 1 = x k + λ k d k k ← k + 1 (72) x_{k+1}=x_k+\lambda_k d_k \\ k \leftarrow k+1 \tag{72} xk+1=xk+λkdkk←k+1(72)转步骤2。
不要忘记高斯-牛顿法的推导过程中要求雅克比矩阵
∇
f
(
x
k
)
\nabla \mathbf{f}(x_k)
∇f(xk)为列满秩矩阵,这在实际计算中可能无法保证。为此,Levenberg和Marquardt对高斯-牛顿方向做了改进,即令
d
k
=
−
(
∇
f
(
x
k
)
T
∇
f
(
x
k
)
+
α
k
I
n
)
−
1
∇
f
(
x
k
)
T
f
(
x
k
)
(73)
d_k=-(\nabla \mathbf{f}(x_k)^T\nabla \mathbf{f}(x_k) + \alpha_k I_n)^{-1}\nabla \mathbf{f}(x_k)^T \mathbf{f}(x_k) \tag{73}
dk=−(∇f(xk)T∇f(xk)+αkIn)−1∇f(xk)Tf(xk)(73)
其中
α
k
>
0
\alpha_k>0
αk>0为参数,
I
n
I_n
In为
n
n
n阶单位矩阵。迭代步长仍为
1
1
1不变。这种方法被称为Levenberg-Marquart法或最小二乘LM法。
3.直接法
直接法是一类不需要函数导数信息的方法,通过试算来寻找优化方向。它的优势在于对函数没有可微假设。劣势是收敛速度往往很慢。
3.1坐标轮换法
坐标轮换法(univariate search technique)的思想极为朴素:从某个初始点 x 0 x_0 x0出发,沿着各个坐标轴依次做一维搜索得 x k x_k xk,直到 x k x_k xk不再显著变化。
坐标轮换法的步骤描述如下:
- 选取初始点。设定 x 0 x_0 x0,初始误差 ϵ > 0 \epsilon>0 ϵ>0,令 k ← 0 k \leftarrow 0 k←0。
- 进行一维搜索。从 x k − 1 x_{k-1} xk−1出发,沿坐标轴方向 e k e_k ek进行一维搜索,求 λ k − 1 \lambda_{k-1} λk−1和 x k x_k xk使得 f ( x k − 1 + λ k − 1 e k ) = min λ f ( x k − 1 + λ e k ) x k = x k − 1 + λ k − 1 e k (74) \begin{aligned}f(x_{k-1}+\lambda_{k-1}e_k)&=\min_{\lambda}f(x_{k-1}+\lambda e_k) \\ x_k&=x_{k-1}+\lambda_{k-1} e_k\end{aligned} \tag{74} f(xk−1+λk−1ek)xk=λminf(xk−1+λek)=xk−1+λk−1ek(74)
- 检查迭代次数。若 k = n k=n k=n,转步骤4;否则令 k ← k + 1 k\leftarrow k+1 k←k+1,返回步骤2。
- 检查终止准则。若 ∣ ∣ x n − x 0 ∣ ∣ < ϵ ||x_n-x_0||<\epsilon ∣∣xn−x0∣∣<ϵ,终止迭代输出 x n x_n xn;否则令 x 0 = x n , k ← 1 x_0=x_n,k\leftarrow 1 x0=xn,k←1,返回步骤2。
可以证明,如果每次一维搜索都有唯一最优解,则坐标轮换法是收敛的。
3.2 其他方法
直接法的其他方法还有:步长加速法(模式搜索法,pattern search method),Rosenbrock旋转方向法,Powell法,单纯形调优法等等,这里不介绍了。
4.结语
数学上有一个原则:把复杂问题转化为简单问题。例如:约束优化通过拉格朗日形式转化为无约束优化。将非线性问题通过近似和残差约束转化为线性问题。
无约束最优化问题是优化领域研究的基础性问题。基础性方法,如最速下降法、牛顿法、共轭梯度法,是其他复杂方法的基础。近年来学术和工业界更加常用的交替方向乘子法(ADMM)和原始对偶法(Primal Dual)也是在基础问题和方法上推演而来。因而掌握经典方法的思想对于理解新方法大有裨益。