系统工程最优化

最优化
1 最优化问题概论
1.1 最优化问题的类型
许多实际问题都可以归结为最优化问题, 其种类繁多, 一般来说分类如下:
最优化问题⎪
⎪⎪⎩
⎪⎪⎪⎨⎧⎪⎪⎪⎩⎪⎪⎪
⎨⎧⎪⎩⎪⎨⎧⎩⎨⎧⎩⎨⎧动态问题非凸规划凸规划非线性规划线性规划约束问题多维问题一维问题
无约束问题静态问题 图1.1 最优化问题的分类
下面是两个有名的无约束优化问题. 例1.1.1 求Rosentrock 函数的极小点, 即
minimize 100(x 2 - x 12)2 + (1 - x 1)2.
例1.1.2 求Dixon 函数的极小点, 即
minimize (1 - x 1)2
+ (1 - x 10)2
+
∑=+-9
1
212)(i i i x x
.
这两个问题常用于检验各种计算方法的优劣.
1.2 非线性规划的例子
例1.2.1 边长为a 正方形铁皮，在四个角处剪去长相等的正方形制成无盖容器. 问如何剪法使容器容积最大?
解 设剪去正方形的边长为x , 与此相应的容器容积
f (x ) = (a - 2x )2x .
令
f ’(x ) = 2(a - 2x )(-2)x + (a - 2x )2 = (a - 2x )(a - 6x ) = 0.
由此解出两个驻点
x =
21a , x =6
1a . 第一个驻点不合实际意义,现在判断第二个驻点是否是极大点. 因为
f ’’(x ) = 24x - 8a , f ’’(
6
1
a ) = -4a < 0, 说明x =61a 是极大点,即将铁皮的四个角剪去长为6
1
a 的正方形. 制成容器的容积为
f (61a ) = (a - 2⨯61a )2 ⨯61a = 27
2a 3. 例1.2.2 把半径为1的实心金属球熔化,铸成一个实心圆柱体.问圆柱体的半径和高各为
多少才能使它的表面积最小?
解 设圆柱体的半径为r , 高为h , 本题的数学模型为
minimize f (r , h ) = 2πrh +2πr 2
subject to
πr 2h =
34π or
r 2h =
3
4. 现在利用Lagrange 乘子法解这个问题. Lagrange 函数为
L (r , h , λ) = 2πrh +2πr 2
-λ(r 2
h -
3
4). 分别对r , h , λ求偏导数,并令它们等于零, 有
r L
∂∂= 2πh +4πr -λ2rh = 0, h L
∂∂= 2πr - λr 2 = 0, λ∂∂L = -r 2h +3
4= 0. 由第一、二两个方程解得h = 2r , 再与第三个方程联立解出
r =
3
32, h = 233
2
. 此圆柱体的表面积是6π3
232
⎪⎭
⎫ ⎝⎛.
例1.2.3 已知n 个年份某经济指标的数值y t , t = 1, 2, … , n . 试用逻辑增长曲线模型
y t =
bt
ae
k
-+1 (1.2.1) 拟合该组样本数据, 其中k , a , b 是待定参数.
将n 个样本数据代入(1.2.1)得
y 1 =b ae k -+1, y 2 =b
ae
k 21-+, … , y n = nb ae k -+1. 按最小二乘法确定参数k , a , b , 此问题需要求解
minimize f (k , a , b ) =
∑
=-+-
n
t bt
t ae
k
y 1
2)1(. 例1.2.4 在数据分类的支持向量机方法中有个最小球问题, 即已知n 维空间m 个点x (1), x (2), …, x (m ), 要用一个半径最小的球将这些点包含在内(此问题作了简化). 求球的中心和半径.
设球的中心为a 半径为R, 问题的模型为
minimize f (a , R ) = R 2
subject to (x (i )
- a )T
(x (i )
- a ) ≤ R 2
, i = 1, 2, … , m .
例1.2.5 Markowitz 资产组合选择问题. 有n 种资产可供选择，这n 种资产收益率的均值为r 1, r 2, … ,r n , 协方差矩阵为H = (σij )n ×n , 其中σij 是第i , j 两种资产收益率的协方差, 投资者的预期收益率为r p . 问如何确定每种资产的比例, 使得资产组合的预期收益率达到r p 并且方差风险最小?
设每种资产占总资产的比例分别为x1, x2, …, x n. 根据概率论知识, 资产组合收益率的均值为r1x1 + r2x2 + … + r n x n, 方差风险为x T Hx, 所以此问题需求解数学模型
minimize f(x) = x T Hx
subject to r1x1 + r2x2 + … + r n x n = r p,
x1 + x2 + … + x n = 1,
x1, x2, … , x n ≥0.
其中x = (x1, x2, … , x n)T.
上面的例 1.2.1和 1.2.3是无约束极值问题, 其他三个问题是约束极值问题, 其中的例1.2.3也称为最小二乘问题.
一般来说, 无约束优化问题的形式为
min f(x),
约束优化问题的形式为
min f(x)
s.t. x∈D.
其中x∈ R n, D⊂ R n, min是minimize的缩写, s.t.是subject to的缩写.
1.3 二维约束极值问题的图解法
例1.3.1用图解法求解
min f(x) = (x1- 2.5)2 + (x2- 3)2
s.t. g1(x) = x1-x22≥ 0,
g2(x) = -2x1-x2 + 6 ≥ 0,
g3(x) = x2≥ 0.
解首先画出问题的可行域, 如图1.1的阴影区域所示. 然后以(2.5, 3)为圆心画一些同心圆, 与可行域相交于x* = (2.25, 1.5)T的点即是最优解.
图1.1 例1.3.1的图示图1.2 例1.3.1的图示
例1.3.2用图解法求解
min f(x) = (x1 + 2)2 + (x2- 2)2
s.t. g1(x) = 1.6x13- 2x2 + 0.875 = 0,
g2(x) = 8x1 + 7x2- 24 ≤ 0,
g3(x) = x1≥ 0,
g4(x) = x2≥ 0.
解首先画出问题的可行域, 见图1.2的曲线段. 然后以(-2, 2)为圆心画一些同心圆, 与曲线段相交于x* = (0, 0.4375)T的点即是最优解.
1.4 基础知识
1.4.1 梯度与Hesse矩阵
n 元实值函数f (x ) = f (x 1,x 2,…,x n )在x =(x 1, x 2,…, x n )T
处的n 个1阶偏导数构成的向量称为f (x )在x 处的梯度(gradient)，记为∇x f (x )或∇f (x )，即
∇f (x ) = T
n x x f x x f x x f ⎪⎪⎭
⎫
⎝⎛∂∂∂∂∂∂)(,,)(,)(21 . f (x )在x 处的Hesse 矩阵(Hessian matrix)∇2f (x )是指其2阶偏导数构成的下述对称矩阵：
∇2
f (x ) =⎪⎪
⎪⎪⎪⎪
⎪
⎪
⎪⎭
⎫
⎝⎛∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂22
221
22222
21
22
1221221
2)()()()()()()()()
(n n n n n x x f x x x f x x x f x x x f x x f x x x f x x x f x x x f x x f . 若f (x ) =
2
1x T
Ax + c T x , 其中A 是n 阶实对称矩阵, c 是n 维列向量，x 是n 维未知列向量.容易证明, ∇f (x ) = Ax + c, ∇2f (x ) = A .
设g (x ) = (g 1(x ), g 2(x ), … , g m (x ))T
, 其中g i (x )是n 元实值函数, i = 1, 2, … , m . g (x )的Jacobi 矩阵是指下述m ⨯ n 矩阵.
Dg (x ) = ⎪⎪⎪⎪⎪⎪⎪⎪
⎭
⎫
⎝
⎛∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂n m
m m n n x x g x x g x
x g x x g x x g x x g x x g x x g x x g )()()()()()
()()()(212221
212111 =⎪
⎪⎪⎪⎪⎭
⎫
⎝⎛∇∇∇T T
2T 1)()()(x g x g x g m . 例1.4.1 求二元函数
f (x ) = f (x 1, x 2) = -x 13 + 4x 12 - 2x 1x 2 + x 22
在x = (-6, -3)T
处的梯度和Hesse 矩阵.
解 因为
1)(x x f ∂∂= -3x 12 + 8x 1 - 2x 2, 2
)
(x x f ∂∂= -2x 1 + 2x 2, 所以
1)(x x f ∂∂= -150, 2
)
(x x f ∂∂= 6, 于是∇f (x ) = (-150, 6)T
. 另外
212)(x x f ∂∂= -6x 1 + 8, 212)(x x x f ∂∂∂= -2, 2
2
2)(x x f ∂∂= 2, 所以
212)(x x f ∂∂= 42, 212)(x x x f ∂∂∂= -2, 2
2
2)
(x x f ∂∂= 2,
于是∇2
f (x ) = ⎪
⎪⎭
⎫
⎝⎛--22242. 1.4.2 T aylor 展开式及有关定理
常用的Taylor 展开式有以下几种.
设f (x )是n 元实值函数, 具有1阶连续偏导数, s ∈ R n
. 那么
f (x + s ) = f (x ) + ∇f (x )T
s + o (|| s ||) (1.4.1a)
或者
f (x + s ) = f (x ) + ∇f (x + θs )T s , θ ∈ (0, 1). (1.4.1b)
若f (x )具有直到2阶的连续偏导数, 那么
f (x + s ) = f (x ) + ∇f (x )T s +
2
1s T ∇2
f (x )s + o (|| s ||2) (1.4.2a) 或者
f (x + s ) = f (x ) + ∇f (x )T s +
2
1s T ∇2
f (x + θs )s , θ ∈ (0, 1). (1.4.2b) 定理1.4.1设
g (x ) = (g 1(x ), g 2(x ), … , g m (x ))T , 其中g i (x )是n 元实值函数, i = 1, 2, … , m , x , s ∈ R n
, t ∈ R. 那么
g (x + s ) - g (x ) =
⎰+1
)(sdt ts x Dg . (1.4.3)
证 令ϕi (t ) = g i (x + ts ), 则ϕi (0) = g i (x ), ϕi (1) = g i (x + s ), ϕi ’(t ) = ∇g i (x + ts )T
s . 将以上式子代入
ϕi (1) - ϕi (0) =
⎰1
)('dt t i
ϕ
得到
g i (x + s ) - g i (x ) =
⎰+∇1
T
)(sdt ts x g i , i = 1, 2, … , m . 以上m 个等式用矩阵符号表示即为(1.4.3).
设A (x ) = (a ij (x ))是m ⨯ n 函数矩阵, 其中a ij (x )是定义在开集D ⊂ R n 上的实值函数. 如果存在正数L , 使得对于任何x (1), x (2)
∈ D , 均有 || A (x (2)) - A (x (1)) || ≤ L || x (2) - x (1) ||,
则称A (x )在D 上Lipschitz 连续, L 称为Lipschitz 常数.
定理1.4.2 设g (x ) = (g 1(x ), g 2(x ), … , g m (x ))T , 其中g i (x )是定义在开凸集D ⊂ R n 上的实值函数, 具有直到2阶的连续偏导数; D g (x )在D 上Lipschitz 连续, L 为Lipschitz 常数. 那么, 对于任何x , x + s ∈ D , 均有
|| g (x + s ) - g (x ) - D g (x )s || ≤
2
L
|| s ||2. (1.4.4) 证 利用定理1.4.1可得
g (x + s ) - g (x ) - D g (x )s = ⎰+1
)(sdt ts x Dg - D g (x )s
= ⎰-+1
])()([dt s x Dg s ts x Dg .
于是
|| g (x + s ) - g (x ) - D g (x )s || ≤
⎰
⋅-+1
)()(dt s x Dg ts x Dg
≤
⎰
⋅1
dt s ts L =
2
L || s ||2
. 设f (x )是定义在开凸集D ⊂ R n 上的实值函数, 具有直到2阶的连续偏导数, 在定理1.4.1
中令g (x ) = ∇f (x ), Dg (x ) = ∇2
f (x ), 便有
∇f (x + s ) = ∇f (x ) +⎰
+∇1
2)(sdt ts x f . (1.4.5)
若f (x )具有直到3阶的连续偏导数, ∇2f (x )在D 上Lipschitz 连续, 则根据定理1.4.2,有
|| ∇f (x + s ) - ∇f (x ) - ∇2
f (x )s || ≤
2
L || s ||2
. (1.4.6) 1.4.3梯度的性质
函数f (x )的梯度有以下重要性质:
(i) 若∇f (x ) ≠ 0, 则∇f (x )与过x 的等值面(2维的情形为等值线) f (x ) = f (x )垂直;
(ii) ∇f (x )是f (x )在x 处变化率最大的方向, 或者说-∇f (x )是f (x )在x 处的最陡下降方向(the steepest descent direction).
设x = (x 1, x 2, … , x n )T . 根据微分学知识, f (x )在x =(1x ,2x ,…,n x )T 处的n 个偏导数不全为0时, 曲面f (x ) = f (x )上过x 的切平面方程为
1)(x x f ∂∂(x 1 -1x ) +2)(x x f ∂∂(x 2 -2x )+ … +n x x f ∂∂)
((x n -n x ) = 0. 用向量符号表示为
∇f (x )T
(x -x ) = 0,
所以∇f (x )是曲面f (x ) = f (x )上过x 的切平面方程的法向量.
例如, f (x ) = f (x 1,x 2) = x 12
+ x 22的梯度函数为
∇f (x ) = (2x 1, 2x 2)T ,
在x = (2, 1)T 处的梯度
∇f (x ) = (4, 2)T .
由于f (x ) = 5, 过x 的等值线为 x 12 + x 22 = 5.
将x 2看作x 1的函数, 两边对x 1求导数得2x 1 + 2x 2
12dx dx = 0. 于是12dx dx = -2
1x x . 在x = (2, 1)T
处,
1
2
dx dx = -2. 所以等值线x 12 + x 22 = 5在x = (2, 1)T 处的切线方程为 x 2 - 1= -2(x 1 - 2),
即2x 1 + x 2 = 5. 其法向量(2, 1)T 与∇f (x ) = (4, 2)T 同方向, 见图1.7.
图1.3
x = (2, 1)T 处的梯度 图1.4 x = (1, 1)T 处的梯度
又如 f (x ) = f (x 1,x 2) = x 12
- x 2的梯度函数为
∇f (x ) = (2x 1, -1)T
,
在x = (1, 1)T 处的梯度
∇f (x ) = (2, -1)T .
由于f (x ) = 0, 过x 的等值线为 x 12 - x 2 = 0.
由此得x 2 = x 12
,
12dx dx = 2x 1. 在x = (1, 1)T 处, 1
2dx dx = 2. 所以等值线x 12 - x 2 = 0上过x = (1, 1)T
的切线方程为x 2 - 1= 2(x 1 - 1), 或-2x 1 + x 2 = -1. 其法向量(-2, 1)T 与∇f (x ) = (2, -1)T 方向相反, 见图1.4.
在3维的情形, ∇f (x )与等值面f (x ) = f (x )的关系如图1.5所示, 其中的虚线表示x 处的切平面.
图1.5梯度与等值面的关系
再考虑梯度的第二个性质. 设0 ≠ s ∈R n
, 那么{x + αs | α ≥ 0}从x 出发沿s 方向的射线,
见图1.6.
图1.6 射线的图示
f (x )在x 处沿s 方向的方向导数(变化率)定义为
s x f ∂∂)(=ααα)
()(lim 0
x f s x f -++→. 记ϕ(α) = f (x + αs ), 则ϕ’(α) = ∇f (x + αs )T s , ϕ’(0) = ∇f (x )T s . 所以
s x f ∂∂)(=αϕαϕα)0()(lim 0
-+→=ϕ’(0) = ∇f (x )T
s . 设|| s || = 1, s 与∇f (x )的夹角为θ, 那么 ∇f (x )T s = ||∇f (x )||cos θ.
当θ = 0, 即s 与∇f (x )同方向时, cos θ = 1, 这时∇f (x )T
s 取得最大值. 所以s =∇f (x )是f (x )在x 处变化率最大的方向, 或者说-∇f (x )是f (x )在x 处的最陡下降方向.
1.4.4凸函数及有关定理
设f (x )是定义在n 维欧氏空间R n 的某个凸集D 上的函数. 如果对于D 中的任意两点x (1)
和x(2)，以及任何实数α (0 < α < 1), 有
f(αx(1) + (1 -α)x(2))≤αf(x(1)) + (1 -α)f(x(2)), (1.4.7)
就称f(x)是定义在D上的凸函数(convex function). 如果x(1) ≠x(2) 时总有
f(αx(1) + (1 -α)x(2)) <αf(x(1)) + (1 -α)f(x(2)) , (1.4.8)
则称f(x)是定义在D上的严格凸函数.
将(1.4.7)和(1.4.8)两式中的不等号反向，则得到凹函数(concave function)和严格凹函数的定义. 显然，若f(x)是凸函数(严格凸函数)则-f(x)是凹函数(严格凹函数).
在1维的情形, 凸函数和凹函数的图形分别如图1.7的(a)和(b)所示,(c)则为非凸非凹的情形.
(a) (b) (c)
图1.7 几种形式的函数
关于凸（凹）函数有以下定理.
定理1.4.3若干个凸函数的非负线性组合是凸函数.
定理1.4.4若f(x)是凹函数，则f(x)≥0的解集是凸集.
定理1.4.5设f(x)在R n的某个开凸集D上具有一阶连续偏导数，则f(x)为D上凸函数的必要充分条件是：对于D中的任何两个不同点x(1)和x(2)有
f(x(2)) -f(x(1))≥∇f(x(1))T(x(2) -x(1)).
定理1.4.6设f(x)在R n的某个开凸集D上具有二阶连续偏导数，则f(x)为D上凸函数的必要充分条件是，f(x)的海色矩阵∇2f(x)在D上处处正半定.
定理1.4.7若A是n阶（对称）正半定矩阵，c是n维列向量，x是n维未知列向量，则二次函数f(x) = x T Ax + c T x是凸函数.
2 无约束极值问题
2.1 局部极小点的概念及有关定理
设f(x)为定义在n维欧氏空间R n的某个子集D上的n元实值函数，x*∈D.若存在正数δ，使得D中所有满足|| x -x*||<δ的x均满足f(x)≥f(x*)，则称x*为f(x)在D上的局部极小点，f(x*)为局部极小值. 若以上等号仅当x = x*时成立，则称x*为f(x)在D上的严格局部极小点，f(x*)为严格局部极小值. 若对于所有x∈D都有f(x)≥f(x*)，则称x*为f(x)在D上的全局极小点，f(x*)为全局极小值. 若以上等号仅当x = x*时成立，则称x*为f(x)在D上的严格全局极小点，f(x*)为严格全局极小值.
将上述的不等号反向，即可得到相应的极大点和极大值的定义.函数的极小点和极大点统称为极值点. 显然，凸函数的极值点一定是极小点.由于函数f(x)的极大点（值）相当-f(x)的极小点（值），下面仅讨论极小问题. 全局极小点又称为最优解，全局极小值称为最优目标函数值.
可微函数f(x)的驻点(stationary point)是满足方程组∇f(x) = 0的x. 许多算法都是先求驻点,然后判断其是否是局部极小点或全局极小点.
定理2.1.1 (局部极小点的1阶必要条件) 设f(x)在x*∈R n的一个邻域内1阶连续可微. 如
果x *是f (x )的局部极小点, 那么∇f (x *) = 0.
证 用反证法,设∇f (x *
) ≠ 0. 令s = -∇f (x *
), 则∇f (x *
+ ts )T
s 是t 的一元函数并且t = 0时为-|| ∇f (x *) ||2 < 0. 所以存在正数T 使得
∇f (x *
+ ts )T s < 0, t ∈ [0, T ].
对于任何t ∈ (0, T ], 根据Taylor 公式有
f (x *
+t s ) = f (x *
) + ∇f (x *
+ θt s )T
s . (2.1.1)
由于θ ∈ (0,1), 0<θt <t , 所以当t 是充分小的正数时, ∇f (x * + θt s )T s <0, 因而f (x * +t s ) < f (x *).这与x *是f (x )局部极小点的前提相矛盾, 所以∇f (x *) = 0.
定理2.1.2 (局部极小点的2阶必要条件) 设f (x )在x *∈R n 的一个邻域内2阶连续可微. 如果x *是f (x )局部极小点, 那么∇f (x *) = 0, ∇2f (x *)正半定.
证 定理2.1.1已证明∇f (x *) = 0, 下面仅证明∇2f (x *)正半定. 假如∇2f (x *
)不是正半定的, 那么存在s ∈R n
使得s T
∇2
f (x *)s <0. 考虑
f (x * + ts ) = f (x *) + ∇f (x *)T ts +
2
1t 2s T ∇2f (x *
+ θts )s = f (x *) +
2
1t 2s T ∇2f (x *
+ θts )s (2.1.2) 其中t ∈ R, θ ∈ (0,1). 由于s T
∇2
f (x *
+ θts )s 是t 的连续函数并且t = 0时为负数, 所以| t |充分小时, (2.1.2) 式右边的第二项是负数,使得f (x * + ts ) < f (x *). 这与x *是f (x )局部极小点的前提相矛盾, 所以∇2f (x *)正半定.
定理2.1.3 (局部极小点的2阶充分条件) 设f (x )在x *
∈R n
的一个邻域内2阶连续可微. 如果∇f (x *
) = 0, ∇2
f (x *
)正定, 那么x *
是f (x )的严格局部极小点.
证 考虑
f (x * + ts ) = f (x *) + ∇f (x *)T ts +
2
1t 2s T ∇2f (x *
+ θts )s = f (x *
) +
2
1t 2s T ∇2f (x *
+ θts )s , (2.1.3) 其中s 是R n 中任意一个非零的向量, t ∈ R, θ ∈ (0,1). 由于s T ∇2f (x * + θts )s 是t 的连续函数并且t = 0时为正数, 所以当 t 充分接近0时, (2.1.3) 式右边的第二项是正数, 因而f (x * + ts ) > f (x *), 即x *
是f (x )的严格局部极小点.
例2.1.1 考虑二元函数
f (x ) = x 14
+ x 24
- x 13
- x 23
.
由于
∇f (x ) =⎪⎪⎭⎫ ⎝⎛--22322
2313434x x x x , ∇2f (x ) =⎪⎪⎭
⎫
⎝
⎛--22222
1126612x x x x , x * = (0,0)T 时, ∇f (x *) =⎪
⎪⎭⎫ ⎝⎛00, ∇2f (x *
) =⎪⎪⎭
⎫ ⎝⎛0000正半定, 所以不能利用定理2.1.3断定x *
= 0是局部极小点.
将x = (t , t )T
代入f (x )可得到
f (t , t ) = 2t 4
- 2t 3
= 2t 3
(t - 1),
可见, 当t 是充分小的正数时, f (t , t )<0 = f (0, 0). 所以x * = (0, 0)不是局部极小点.
下面介绍有关凸函数的几个定理.
定理2.1.4 设f (x )是凸函数. 那么
(i) f (x )的任何局部极小点x *
是其全局极小点.
(ii) 如果f (x )可微, 那么任何驻点x *
是全局极小点.
证 (i) 假设x *是局部极小点但不是全局极小点. 那么可找到一个点x 使得f (x ) < f (x *). 考虑x *
和x 之间的直线段
x = t x + (1 - t )x *, 0 < t ≤ 1. 由于f (x )是凸函数,
f (x ) ≤ tf (x ) + (1 - t )f (x *) = t [f (x ) - f (x *)] + f (x *) < f (x *).
当t 充分接近0时, x 充分接近x *
, 所以上式与x *
是局部极小点相矛盾. 故x *
是全局极小点.
(ii) 假若∇f (x *) = 0, 但x *不是全局极小点. 那么可找到一个点x 使得f (x ) < f (x *). 令
ϕ(t ) = f (t x + (1 - t )x *) = f (x *+ t (x - x *)), 那么
ϕ’(t ) = ∇f (x *+ t (x - x *))T (x - x *).
于是
∇f (x *)T (x - x *) = ϕ’(0) = t
t t )
0()(lim 0
ϕϕ-+
→
=t x f x x t x f t )
())((lim ***0--++→≤ =t
x f x f t x tf t )()()1()(lim **0--++→ = f (x ) - f (x *
) < 0.
这样一来∇f (x *) ≠ 0. 所以f (x )的驻点x *是全局极小点.
定理2.1.5设f (x )2阶连续可微, 水平集Ω = {x ∈ R n | f (x ) ≤ f (x (0))是凸集, 其中x (0) ∈ R n , ∇f (x (0)) ≠ 0. 如果∇2f (x )在Ω上处处正定, 那么, Ω内部有个点x *使得∇f (x *) = 0.
证 Ω应当是个有界集. 要不然存在x ∈ Ω, 0 ≠ s ∈ R n , 对于任意正数t , 都有f (x + ts ) ∈ Ω. 但是
f (x + ts ) = f (x ) + ∇f (x )T ts +
2
1t 2s T ∇2
f (x + θts )s , 其中θ ∈ (0, 1). 由于∇2f (x + θts ) 在Ω上处处正定, 当t 充分大时, f (x + ts )将会超过f (x (0)), 这与f (x + ts ) ∈ Ω相矛盾. 所以Ω是个有界集, 于是f (x )在Ω内部某个点x *
处取得最小值. 根据定理2.1.1, ∇f (x *) = 0.
推论2.1.1 考虑二次函数
f (x ) =
2
1x T Ax + c T
x , 其中A 是n 阶实对称矩阵, c 是n 维列向量, x 是n 维未知列向量.
(i) 如果A 正半定并且方程组∇f (x ) = Ax + c = 0有解, 那么其任何解是f (x )的全局极小点.
(ii) 如果A 正定, 那么x * = -A -1
c 是f (x )的严格全局极小点.
证 (i) 若存在x *∈R n 使得∇f (x *) = Ax * + c = 0, 由于f (x )是凸函数, 根据定理2.1.4(ii), x *是f (x )的全局极小点.
(ii) 令∇f (x ) = Ax + c = 0, 得到x *
= -A -1
c . 由于∇2
f (x *
) = A 正定, 定理2.1.3和2.1.4, x *
= -A -1c 是f (x )的严格全局极小点.
2.2 函数极小点的计算原理
根据定理2.1.1, 可微非线性函数f(x)在x*∈R n处取得局部极小的必要条件是∇f(x*) = 0. 因而可考虑解方程组
∇f(x) = 0 (2.2.1)
求出f(x)的驻点, 然后再判断其是否是局部或全局极小点. 但(2.2.1)是一个具有n个等式n个变量的非线性方程组, 一般难以求解. 而且, 即使求出一个解, 很可能也不是局部极小点. 实际中常采用所谓的下降(迭代)算法, 就是说从n维空间某个点x(0)出发, 在其附近找一个点x(1), 使得f(x(1)) < f(x(0)). 然后在x(1)附近找一个点x(2), 使得f(x(2)) < f(x(1)). 如此等等, 产生一个序列
x(0), x(1), …, x(k), x(k+1), …,
它们对应的函数值满足
f(x(0)) > f(x(1)) > … > f(x(k)) > f(x(k+1)) > …,
一直到趋近于一点x*使得∇f(x*)充分接近零向量, 即x*是驻点. 由于计算过程中函数值一直在减少, x*极有可能是一个局部极小点.
假设现在有了一点x(k), 根据Taylor展开式, x(k)与附近的一点x(k+1)有如下关系式:
f(x(k+1)) = f(x(k)) + ∇f(x(k))T(x(k+1)-x(k)) + o(|| x(k+1)- x(k) ||).
令s(k) = x(k+1)-x(k), 上式成为
f(x(k+1)) = f(x(k)) + ∇f(x(k))T s(k)+ o(|| s(k) ||).
如果|| s(k) ||充分小, ∇f(x(k))T s(k) <0, 即s(k)与∇f(x(k))的夹角大于90︒, 那么令x(k+1) = x(k) + s(k)就有
f(x(k+1)) < f(x(k)).
但|| s(k) ||较大时, 上式可能不成立. 不过, 只要α是充分小的正数, 令
x(k+1) = x(k) + αs(k), (2.2.2)
仍会有f(x(k+1)) < f(x(k)). 使得∇f(x(k))T s(k) < 0的s(k)称为f(x)在x(k)处的下降方向(descent direction).
(2.2.2)式中的α称为步长(step length). 在x(k)和s(k)已知的条件下, 集合
{x(k) + αs(k) | α > 0}
是一条以x(k)为端点的射线. 在这条射线上找一点x=x(k)+ αs(k)使得f(x) < f(x(k))称为线搜索(line search)或一维搜索.
算法2.2.1计算n元函数f(x)驻点的下降迭代算法
(i) 给定x(0)∈ R n, 置k = 0;
(ii) 如果|| ∇f(x(k)) || < ε (如ε = 10-6~10-4)则停止; 否则
(iii) 寻找一个s(k)∈ R n, 使得∇f(x(k))T s(k) < 0;
(iv) 进行线搜索, 找一个αk并令x(k+1) = x(k) + αs(k),使得f(x(k+1)) < f(x(k));
(v) 令k := k + 1, 返回(ii).
2.3 线搜索的方法
2.3.1 几种下降方向
前面讲到, 下降算法是从某一点x(k)出发沿着一个下降方向s(k)移动到下一点x(k+1)使得f(x(k+1)) < f(x(k)). 这种方法的关键有两点, 一是要确定下降方向s(k), 即找一个非零向量s(k)使得∇f(x(k))T s(k) <0; 二是要确定合适的步长αk使得f(x(k) + αk s(k)) < f(x(k)).
设∇f(x(k)) ≠0, 典型的下降方向有
s(k) = -∇f(x(k)),
即s(k)是f(x)在x(k)处的负梯度方向或最陡下降方向. ∇2f(x(k))正定时,
s(k) = -∇2f(x(k))-1∇f(x(k))
也是下降方向, 称为牛顿方向. 事实上,只要A是正定矩阵,
s(k) = -A∇f(x(k))
就是下降方向, 因为∇f(x(k))T s(k) = -∇f(x(k))T A∇f(x(k)) < 0.
图2.1显示s(k)与∇f(x(k))以及x(k) + αk s(k)与x(k)和s(k)之间的关系. s(k)和∇f(x(k))这两个向量的起点都应该是原点, 但把它们看作”自由”向量, 画在x(k)处.
O
图2.1 x(k) + αk s(k)的几何意义
以负梯度作为搜索方向的算法称为梯度法, 以牛顿方向作为搜索方向的算法称为牛顿法, 还有其他许多算法, 在第3章讨论. 本章假定搜索方向已经确定, 仅介绍如何确定步长.
2.3.2 线搜索的步长条件
第k次迭代的起始点x(k)和搜索方向s(k)已知的条件下, f(x(k) +αs(k))是α的一元函数. 记
ϕ(α) = f(x(k) +αs(k)), α > 0.
确定步长α一方面希望函数值有足够的减少,另一方面则希望计算量不要太大,因此需要做出一种折衷的选择. 找出一个α使得ϕ(α)取得局部或全局极小的方法称为精确线搜索, 否则称为非精确线搜索. 根据复合函数求导法, 精确线搜索的必要条件是
ϕ’(α) = ∇f(x(k) +αs(k))T s(k) = 0. (2.2.3)
精确线搜索计算量大, 下面仅讨论非精确线搜索的方法.
一种非精确线搜索的方法是寻找一个步长α, 满足
f(x(k) +αs(k)) ≤f(x(k)) + c1α∇f(x(k))T s(k), (2.2.4a)
∇f(x(k) +αs(k))T s(k)≥c2∇f(x(k))T s(k), (2.2.4b)
其中0<c1<c2<1, 如c1 = 10-4, c2 = 0.9. (2.2.4a)称为充分减少条件(sufficient decrease condition)或Armijo条件, (2.2.4b)称为曲率条件(curvature condition), 两者统称为Wolfe条件.
图2.2反映条件(2.2.4a), 其中的曲线是ϕ(α) = f(x(k) +αs(k))的图形, 虚线是线性函数l(α) = f(x(k)) + c1α∇f(x(k))T s(k)的图形. 虚线下方的曲线对应的区间, 即箭线表示的范围内的α值满足(2.2.4a)式.
图2.2 充分减少条件图示
条件(2.2.4b)也可写为
ϕ’(α) ≥c2ϕ’(0).
图2.3反映此条件, 其中虚线的斜率为c2ϕ’(0), 切点对应的两个区间, 即箭线表示的范围内的α值满足(2.2.4b). 由(2.2.4a)式限定的α值实际上可以很小, (2.2.4b)则迫使其进入一个α值较大的区间.
图2.3 曲率条件图示
但是(2.2.4b)也可能使α远离局部极小点. 为使其不至于离局部极小点太远, 可使用强Wolfe 条件. 此条件要求α满足
f (x (k ) + αs (k )) ≤ f (x (k )) + c 1α∇f (x (k ))T s (k ), (2.2.5a)
|∇f (x (k )
+ αs (k ))T s (k )
| ≤ c 2|∇f (x (k ))T s (k )
|, (2.2.5b)
其中0<c 1<c 2<1. 与Wolfe 条件相比, 强Wolfe 条件的第二个式子除了要求α满足(2.2.4b)外, 还要求α满足
∇f (x (k ) + αs (k ))T s (k ) ≤ -c 2∇f (x (k ))T s (k ).
下述定理表明, 线搜索可以找到一个步长满足Wolfe 条件和强Wolfe 条件.
定理2.3.1 设f (x )是连续可微的n 元实值函数. 令s (k )是x (k )处的一个下降方向, f (x )沿着射线{x (k ) + αs (k ) | α > 0}有下界, 0<c 1<c 2<1. 那么存在一个步长区间满足Wolfe 条件和强Wolfe 条件.
证 由于ϕ(α) = f (x (k ) + αs (k ))在α >0时有下界, 0 < c 1 < 1, 直线l (α) = f (x (k )) + c 1α∇f (x (k ))T s (k )与ϕ(α)的图形至少相交于一点. 令α是交点对应的最小α值, 这时
f (x (k ) + αs (k )) = f (x (k )) + c 1α∇f (x (k ))T s (k ),
所有小于α的步长都满足充分减少条件. 根据一元函数微分学的Lagrange 定理, 存在
α
~∈(0,α)使得 f (x (k ) +αs (k )) - f (x (k )) =α∇f (x (k ) +α
~s (k ))T s (k ). 以上两式合在一起, 考虑到∇f (x (k ))T s (k )< 0和0< c 1< c 2, 得
∇f (x (k ) +α~
s (k ))T s (k ) = c 1∇f (x (k ))T s (k ) > c 2∇f (x (k ))T s (k ).
所以α~满足Wolfe 条件并且使严格的不等式成立. 由于f (x )是光滑函数, 存在一个包含α~的
区间使Wolfe 条件成立, 由于∇f (x (k ) +α~s (k ))T s (k ) <0, 而-c 2∇f (x (k ))T s (k ) >0, 所以强Wolfe 条件也
成立.
与Wolfe 条件类似, 满足Goldstein 条件的步长也可使f (x )减少并防止步长αk 不至于太小. 该条件为
f (x (k )) + (1 - c )αk ∇f (x (k ))T s (k ) ≤ f (x (k ) + αk s (k )) ≤ f (x (k )) + c 1αk ∇f (x (k ))T s (k ), (2.2.6)
其中0< c < 1/2. 第二个不等式就是(2.2.4a), 第一个不等式则从ϕ(α)图形的下方控制步长, 见图2.4所示.
图2.4 Goldstein 条件图示
2.3.3 步长的确定
为使步长满足Wolfe 条件或Goldstein 条件并非易事, 每次迭代都需从一个初始步长开始, 反复计算以找到一个满足Wolfe 条件或Goldstein 条件的步长, 我们将这一过程称为步长试验. 用αk 表示第k 次迭代实际使用的步长, 为了获取αk 而设定的初始步长为αk 0. 如果只要求αk 满足充分减少条件(2.2.4a), 可使用以下程序进行步长试验.
算法2.3.1 计算满足充分减少条件的步长αk . (i) 选择αk 0 >0, ρ, c ∈(0, 1), α = αk 0;
(ii) 若f (x (k ) + αs (k )) ≤ f (x (k )) + c α∇f (x (k ))T s (k ), 则αk =α, 停止, 否则 (iii) 令α := ρα, 转(ii).
以上算法的每次迭代需要计算f (x )和∇f (x )的数值, 迭代次数较多时计算量较大而且产生的步长αk 可能很小. 作为这一算法的替代, 可使用插值法作步长试验. 用αk 和αk -1表示第k 和第k -1次迭代实际使用的步长, 用αi , αi -1和αj 等表示为了获取αk 和αk -1而在试验阶段产生的步长. 仍用αk 0表示第k 次迭代最初设定的步长. 用ϕ(α)表示的充分减少条件为
ϕ(αk ) ≤ ϕ(0) + c 1αk ϕ’(0).
如果
ϕ(αk 0) ≤ ϕ(0) + c 1αk 0ϕ’(0), (2.2.7)
就令αk := αk 0, 停止. 否则区间[0, αk 0]包含满足充分减少条件的步长. 利用ϕ(0), ϕ’(0)和ϕ(αk 0)
构造二次插值函数
ϕq (α) =))0()0()((2
'00k k k αϕαϕαϕ--α2
+ϕ’(0)α + ϕ(0). 它满足ϕq (0) = ϕ(0), ϕq ’(0) = ϕ’(0), ϕq (αk 0) = ϕ(αk 0). ϕq (α)的极小点为
α1 = ]
)0()0()([2)0(0'020
'k k k αϕϕαϕαϕ---. (2.2.8)
如果α1满足充分减少条件则停止, 否则利用ϕ(0), ϕ’(0), ϕ(αk 0)和ϕ(α1)构造三次插值函数
ϕc (α) = a α3 + b α2 + αϕ’(0) + ϕ(0),
其中
⎪⎪⎭⎫ ⎝⎛b a =)
(101
2120k k αααα-⎪⎪⎭
⎫ ⎝⎛--31302120ααααk k
⎪⎪⎭⎫ ⎝
⎛----0'01'1)0()0()()0()0()(k k αϕϕαϕαϕϕαϕ
该三次插值函数的极小点为
α2 = a
a b b 3)
0(3'2ϕ-+-. (2.2.9)
如果α2满足充分减少条件则停止, 否则利用ϕ(0), ϕ’(0)和最新的两个ϕ(α∙)构造三次插值函数以获取新的α值, 直到其满足充分减少条件(2.2.4a). 如果αi 太接近αi -1或者比αi -1小很多, 就令αi -1 = αi -1/2, 以保证最终产生的αk 不至于太小.
步长试验中的初始步长αk 0可以人为给定, 至少α00必需人为给定. 进行了一次迭代后, 以后的初始步长可按下式计算:
αk 0 = αk -1)
()()
1()1()()(k T k k T k s x f s x f ∇∇--.
或者利用f (x
(k -1)
), f (x (k ))和ϕ’(0) = ∇f (x (k ))T s (k )
构造一个插值函数, 以其极小点
αk 0 = )
0()]()([2'
)
1()(ϕ--k k x f x f
作为步长试验的初始步长.
下面是一个获取满足强Wolfe 条件的步长试验程序, 其中c 1 = 10-4
, c 2 = 0.9. 算法2.3.2 计算满足强Wolfe 条件的步长αk . (i) 令αk 0 = 0, 选择正数α1和αmax , 令i = 1;
(ii) 若ϕ(αi ) > ϕ(0) + c 1αi ϕ’(0), 或者当i > 1时有ϕ(αi ) > ϕ(αi -1), 令
αk = zoom(αi -1, αi ),
停止; 否则
(iii) (a) 若|ϕ’(αi )| ≤ -c 2ϕ’(0), 令αk = αi , 停止; 否则 (b) 若ϕ’(αi ) ≥ 0, 令
αk = zoom(αi , αi -1),
停止; 否则
(c) 选择αi +1 ∈ (αi , αmax ), 令i := i +1, 转(ii).
以上算法的zoom(∙, ∙)是子程序zoom 返回的一个α值.
程序zoom 如下, 其中的αlo 是计算过程中产生的所有满足充分减少条件的步长对应的函数值的最小值, αhi 满足ϕ’(αlo )(αhi - αlo ) < 0, c 1 = 10-4, c 2 = 0.9.
算法2.3.3 获取步长αj 的子程序zoom
(i) 利用二次或三次插值法在αhi 和αlo 之间找一个试验步长αj ;
(ii) 若ϕ(αj ) > ϕ(0) + c 1αj ϕ’(0), 或者ϕ(αj ) ≥ ϕ(αlo ), 令αhi = αj , 转(i); 否则 (iii) (a) 若|ϕ’(αj )| ≤ -c 2ϕ’(0), 令zoom(∙, ∙) = αj , 停止; 否则 (b) 若ϕ’(αj )( αhi - αlo ) ≥ 0, 令αhi = αlo ; (c) 令αlo = αj , 转(i).
以上内容详见Nocedal&Wright(2006)的著述.
3 无约束最优化的计算方法
3.1 梯度法
设f (x )是定义在R n
上的光滑函数. 下面考虑如何使用下降算法寻找f (x )的一个驻点x *
. 设第k 次迭代的初始点为x (k ), ∇f (x (k )) ≠ 0. 梯度法以s (k ) = -∇f (x (k ))为搜索方向, k = 0, 1, 2,等等. 迭代公式为
x (k +1) = x (k ) - αk ∇f (x (k )). (3.1.1)
其中αk 由精确或非精确线搜索方法确定, 使得∇f (x (k +1))<∇f (x (k )). 如果αk 是
ϕ(α) = f (x (k ) + αs (k ))
在α >0时的一个局部极小点, 那么αk 满足
ϕ’(α) = ∇f (x (k ) + αs (k ))T s (k ) = 0. (3.1.2)
这时的梯度法称为最速下降法.
一般来说, 最速下降法的步长需要使用精确线搜索获取,计算量较大. 但对于某些特殊的函数, 如严格凸二次函数, 容易得到步长的解析表达式.
考虑
min f (x ) =
2
1x T Ax + c T
x , (3.1.3)
其中A 是n 阶正定(对称)矩阵, c 是n 维列向量, x 是n 维未知列向量.
由于
∇f (x ) = Ax + c ,
因此
∇f (x (k ) + αs (k )) = A (x (k ) - α∇f (x (k ))) + c .
将上式和s (k ) = -∇f (x (k )
)代入(3.1.2), 得到
0 = -∇f (x (k ))T
[A (x (k )
- α∇f (x (k )
)) + c ] = -∇f (x (k ))T [Ax (k ) + c - αA ∇f (x (k ))] = -∇f (x (k ))T [∇f (x (k )) - αA ∇f (x (k ))],
所以
αk = )
()()()()
(T )()(T )(k k k k x f A x f x f x f ∇∇∇∇. (3.1.4) 例3.1.1 用最速下降法解
min f (x ) =
2
1(x 12
- 4x 1x 2 + 5x 22) - x 1 + 2x 2, 初始点x (0)
= (-1, 1)T
.
解 目标函数的Hesse 矩阵为A =⎪
⎪⎭
⎫
⎝⎛--5221, 初始值f (x (0)
) = 8. 由于 ∇f (x ) =⎪⎪⎭
⎫ ⎝⎛++---252122121x x x x ,
所以∇f (x (0)) = (-4, 9)T . 利用(3.1.4)有
α0 = )
()()
()()
0(T )0()0(T )0(x f A x f x f x f ∇∇∇∇= 0.1716814. 于是
x (1)
= x (0)
- α0∇f (x (0)
) = (-0.3132743, -0.5451328)T
.
以后的计算过程类似, 共进行了6次迭代, 全部计算结果见表3.1.
表3.1 最速下降法的计算结果(x (0) = (-1, 1)T )
点列x (0), x (1), x (2)等等如图3.1所示, 其中的椭圆是令目标函数f (x )等于一些正数得到的等值线. 可以看到前后两个下降方向s (k ) = -∇f (x (k ))和s (k +1) = -∇f (x (k +1))总是直交的, 这种情况称为”锯齿”现象. 椭圆越扁,需要的迭代次数越多. 一般来说, 最速下降法的收敛速度不快.
图3.1 梯度法中的锯齿现象
3.2 牛顿法
3.2.1 一元非线性方程的牛顿法
牛顿法起源于一元非线性方程求根. 现在设f (x )是一元实值函数. 如果存在实数x *, 使得f (x *
) = 0, 称x *
是方程f (x ) = 0的根, 或者说x *
是函数f (x )的零点. f (x )是线性或二次函数时, 很容易求出其零点(若存在的话). 当f (x )是更复杂的非线性函数时, 通常从某一个实数x 0出发, 通过一系列迭代产生一个序列
x 0, x 1, … , x k , x k +1, …
使其趋近于f (x )的一个零点. 产生点列的方法很多, 牛顿法是十分有效的一种. 设f (x )一阶可微, 其图形如图3.2的曲线所示.
图3.2 牛顿法的几何原理
其中的x k 为当前解. 过曲线上的点(x k , f (x k ))作切线, 该切线与横轴的交点x k +1即是下一次迭代的起始点. 由图可见
f ’(x k ) =
1
)
(+-k k k x x x f ,
所以
x k +1 = x k -
)
()('
k k x f x f (3.2.1)
此式即为牛顿迭代公式. 迭代终止的条件为|f (x k +1)|或|x k +1 - x k |是充分小的正数. (3.2.1)式也可将f (x )在x k 处线性化得到, 即由
0 = f (x ) ≈ f (x k ) + f ’(x k )(x - x k )
得出
x = x k -
)
()('
k k x f x f .
3.2.2 非线性方程组的牛顿法
考虑方程组
g i (x ) = 0, i = 1, 2, … , n , (3.2.2)
其中g i (x )是可微的n 元实值函数.
对于某个x (k ) ∈ R n , 将每个g i (x )在x (k )处线性化得
0 = g i (x ) ≈ g i (x (k )) + ∇g i (x (k ))T (x - x (k )), i = 1, 2, … , n . (3.2.3)
记
g (x (k )) = ⎪⎪⎪⎪
⎪⎭
⎫
⎝⎛)()()()()(2)(1k n k k x g x g x g , D g (x (k )) = ⎪⎪
⎪⎪
⎪⎭
⎫
⎝⎛∇∇∇T )(T )(2T )(1)()()(k n k k x g x g x g 则(3.2.3)可写为
g (x (k )
) + D g (x (k )
) (x - x (k )
) = 0. (3.2.4)
若D g (x (k ))可逆, 则有
x = x
(k +1)
= x (k ) - D g (x (k ))-1g (x (k )
). (3.2.5)
3.2.3 多元函数极值的牛顿法
设f (x )是可微的n 元实值函数, 现在要求一个x *∈R n , 使得∇f (x *) = 0.
由于∇f (x ) = 0是一个具有n 个变量n 个方程的方程组, 即
i
x x f ∂∂)
( = 0, i = 1, 2, … , n , 而∇f (x )的Jacobi 矩阵D (∇f (x )) = ∇2f (x ), 与(3.2.4)对应有
∇f (x (k )
) + ∇2
f (x (k )
)(x - x (k )
) = 0.
若∇2
f (x (k )
)可逆, 则
x (k +1)
= x (k )
- ∇2
f (x (k ))-1
∇f (x (k )
). (3.2.6)
公式(3.2.6)常用于求f (x )的局部极小点, 这时常假定∇2f (x (k )
)正定. 向量 s (k ) = -∇2f (x (k ))-1∇f (x (k )) (3.2.7)
称为x (k )处的牛顿方向. 由于
∇f (x (k ))T s (k ) = -∇f (x (k ))T ∇2f (x (k ))-1∇f (x (k )) < 0,
牛顿方向确实是下降方向.
也可以换一个角度看待牛顿法. 将f (x )在x (k )处作Taylor 展开, 忽略高于二次项的部分得
f (x ) ≈ f (x (k )) + ∇f (x (k ))T (x - x (k )) +
2
1
(x - x (k ))T ∇2f (x (k )) (x - x (k )). 代之以min f (x )的是
min q (x ) = ∇f (x (k ))T (x - x (k )) +
2
1
(x - x (k ))T ∇2f (x (k )) (x - x (k )). (3.2.8) 令
∇q (x ) = ∇f (x (k )) + ∇2f (x (k )) (x - x (k )) = 0,
便得到
x (k +1) = x (k ) - ∇2f (x (k ))-1∇f (x (k )).
此式便是(3.2.6). 令s (k ) = x (k +1) - x (k ), 上式成为
s (k )
= -∇2
f (x (k ))-1
∇f (x (k )
).
此式便是(3.2.7).
由于q (x )仅是f (x )在x (k )
处的近似, 只有当x (k +1)
距离x (k )很近, 即s (k )
很短时, 近似程度才
较好, 所以常按下式确定下一个点.
x (k +1) = x (k ) + αs (k ), (3.2.9)
其中的α即为步长, 通常0<α≤ 1. 有的学者将按公式(3.2.9)进行迭代的方法称为阻尼牛顿法.
再考虑严格凸二次规划(3.1.3). 将f (x )配完全平方:
f (x ) =
21x T Ax + c T x = 21(x + A -1c )T A (x + A -1
c ) -2
1c T A -1c , 可见最优解为x = -A -1c . 现在设x (0)是R n 的任何一点, 令步长为1, 用牛顿迭代公式求x (1)
. 由于
∇f (x (0)) = Ax (0) + c , ∇2f (x (0)) = A ,
代入(3.2.6)得
x (1)
= x (0)
- A -1
(Ax (0)
+ c ) = -A -1
c ,
所以用牛顿法只需一次迭代便可得到一个严格凸二次函数的最优解
x * = -A -1c .
3.2.4 牛顿法的计算步骤及收敛性
算法3.2.1 用牛顿法求函数f (x )驻点的计算步骤
设∇2
f (x )处处非奇异.
(i) 给定x (0)∈R n , 置k = 0;
(ii) 令x (k +1) = x (k ) - ∇2f (x (k ))∇f (x (k ));
(iii) 若|| ∇f (x (k )) || < ε (如ε = 10-6~10-4), 停止, x (k +1)为所求; 否则令k := k +1, 返回步骤(ii). 定理 3.2.1 设f (x )是具有3阶连续偏导数的实值函数, ∇2f (x )在x *的一个邻域N 上 Lipschitz 连续, ∇f (x *) = 0, ∇2f (x *)正定. 那么,由牛顿迭代公式x (k +1) = x (k ) - ∇2f (x (k ))-1∇f (x (k ))所产生的点列如果收敛到x *, 则收敛的阶为2.
证 根据定理1.4.2, 当k 充分大时, 存在正数L 使得
|| ∇f (x *) - ∇f (x (k )) - ∇2f (x (k ))(x - x (k ))|| ≤
2
1
L || x - x (k ) ||2, 并且∇2f (x (k ))正定. 于是
|| x (k +1) - x * || = || x (k ) - x * - ∇2f (x (k ))-1∇f (x (k ))|| = || ∇2f (x (k ))-1[∇f (x *) - ∇f (x (k )) - ∇2f (x (k ))(x * - x (k ))]|| ≤
2
1
L || ∇2f (x (k ))-1 || || x * - x (k ) ||2 ≤ L || ∇2f (x *)-1 || || x (k ) - x * ||2 .
所以定理成立.
以上定理指出, 牛顿法具有非常理想的收敛速度, 实际计算效果也非常好. 例3.2.1 用牛顿法解Rosenbrock 问题
min f (x ) = 100(x 2 - x 12)2 + (1 - x 1)2,
初始点为x (0) = (0, 0)T .
解 目标函数的梯度和Hesse 矩阵为
∇f (x ) =⎪⎪⎭⎫ ⎝⎛-+---)(200)22)(40021212121x x x x x x , ∇2f (x ) =⎪⎪⎭
⎫
⎝⎛--+--2004004002)3(40011212x x x x . 由于x (0)
= (0, 0)T
, 可算得
∇f (x (0)
) =⎪⎪⎭⎫ ⎝⎛-02, ∇2f (x (0)) =⎪⎪⎭
⎫ ⎝⎛200002, ∇2f (x (0))-1 =⎪⎪⎭⎫
⎝⎛005.0005
.0. 由此可得到
x (1) = x (0) - ∇2f (x (0))-1∇f (x (0)) = (1, 0)T , || ∇f (x (0)) ||2 = 4.
∇f (x (1)
) =⎪⎪⎭
⎫
⎝⎛-200400, ∇2f (x (1)
) =⎪⎪⎭
⎫
⎝⎛--2004004001202, ∇2
f (x (1))-1
=⎪⎪⎭⎫ ⎝
⎛⨯±⨯⨯----233310495025.110975124.410975124.410487562.2. 于是
x (2) = x (1) - ∇2f (x (1))-1∇f (x (1)) = (1, 1)T ,
∇f (x (2)) = (4.768372⨯10-5, -2.384186⨯10-5)T ,
|| ∇f (x (2)) ||2 = 2.842171⨯10-9. 其中|| ∇f (x (2)) || < 10-4, 停止计算. 问题的解为x = (1, 1)T
.
如果初始点x (0)
= (10, 10)T
, 经过6次迭代得到 x = (1.000005, 1.000011)T
.
初始点x (0) = (100, 100)T 时, 经过10次迭代得到
x = (0.9999962, 0.9999925)T .
3.3 共轭梯度法
3.3.1 共轭方向及有关定理
设A 是n 阶实对称正定矩阵, s (0), s (1), …, s (m )是n 维非零向量. 如果
(s (i ))T
As (j )
= 0, i , j = 0, 1, … , m , i ≠ j ,
称s (0)
, s (2)
, …, s (m )
A -共轭, 简称共轭.
定理3.3.1 考虑二次规划
min f (x ) =
2
1x T
Ax + cx (3.3.1) 其中A 是n 阶正定矩阵, c 是n 维列向量, x 是n 维未知列向量. 设s (0), s (1), …, s (n -1)A -共轭. 从某一x (0)出发依次沿s (0), s (1), …, s (k )作精确线搜索到达x (k +1). 那么 ∇f (x (k +1))T s (i ) = 0, i = 0, 1 , … , k ,
即每个迭代点处的梯度与前面的搜索方向直交.
证 由于是精确线搜索, 有 ∇f (x (i +1))T s (i ) = 0, i = 0, 1, … , k ,
即定理对于i = k 成立.
现在设i 至少比k 少1. 由于 ∇f (x (j +1)) - ∇f (x (j )) = A (x (j +1) - x (j )) = αj As (j ),
其中αj 是第j 次精确线搜索的步长, 所以
∇f (x
(k +1)
) =
∑+=+∇-∇k
i j j j x f x f 1
)()1())()((+ ∇f (x (i +1))
=
∑+=k
i j j j
As 1
)(α
+ ∇f (x (i +1)).
上式两边的转置由乘以s (i ), 考虑到s (i )与s (j )共轭(i ≠ j ), 有
∇f (x
(k +1))T s (i )
=
∑+=k
i j i j j
As s 1
)(T )()(α
+ ∇f (x (i +1))T s (i ) = 0.。