http://www.ma-xy.com
目录
第一章 无约束非线性规划 1
1.1 引言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 问题的引入与分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.3 模型规范化及基本理论 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3.1 点的性质 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3.2 目标函数的性质 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3.3
最优性条件
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.4 算法框架 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.5 搜索步长的确定：一维搜索 (线搜索) . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.5.1 黄金分割法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5.2 Fibonacci 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.5.3 二次插值法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.5.4 三次插值法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.5.5 Armijo-Goldstein 准则 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.6 MATLAB 应用实例 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.7 搜索方向的确定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.7.1 最速下降法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.7.2 牛顿法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.7.3 修正牛顿法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.7.4 信赖域方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.7.5 共轭梯度法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.7.6 拟牛顿法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
1.7.7 模式搜索法：不使用导数的最优化方法 . . . . . . . . . . . . . . . . . . . . 26
1.8 MATLAB 应用实例 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
1
http://www.ma-xy.com
第一章 无约束非线性规划
1.1 引言
许多模型 (如支持向量机，线性回归等) 最终都会转化为一个优化问题。这是因为我们在处理
问题时，总是希望在众多可选择的情况中选择最优的。回想一下，在数学分析中的优化问题：我
们会求一个函数 f (x) 的极大值、极大值点、最大值、最大值点等等。又或者，我们在“Lagrange
乘子法”中讨论了在等式约束“h(x) = 0”下 f (x) 最优化。现在，我们将这些所谓的函数“放
宽”：我们希望在“一类事物”中，找到某个指标下的最优个体 (或者最优群体)。明显，这一类
事物应该是一个集合。不妨记集合为 ϕ，集合中的个体 x 的特点 (特征) 记为 I(x)。
当然，我们不能一眼看出来集合中哪个个体最好，所以，我们需要东奔西跑的去找。但是，
我们一般不会盲目的去寻找。我们将寻找最优解 (个体) 的方法分为有指导 (有方向) 的寻找和随
机性寻找以及二者相结合的 3 大类方法。这种寻找过程 (优化算法) 必然是一步一步反复迭代的：
这一步找到了张三，下一步找到了李四，最终找到最优结果。有指导寻找：A → B → ···；随机
寻找：A, C, F, ···。一个显著的特点是：有指导的寻找应该是现在的结果比上一次要好，随机则
不一定，而二者的结合方法亦不确定。
我们不可能设计一个算法 (寻找方法)，说这个算法对所有优化问题都适用，毕竟各种优化问
题会有它自身的特点，因此要具体问题具体分析。后面，我们将在分析具体的某一类优化问题时，
给出其适应的算法。
1.2 问题的引入与分析
考虑如下不含参数的无约束非线性优化问题
min
x∈R
2
f(x) = x
1
exp(−(x
2
1
+ x
2
2
)) + (x
2
1
+ x
2
2
)/20 (1.1)
在上述优化问题中外加参数，形成如下含参数的无约束非线性规划问题
min
x∈R
2
f(x; a, b, c) = (x
1
− a) exp(−(x
1
− a)
2
+ (x
2
− b)
2
) + ((x
1
− a)
2
+ (x
2
− b)
2
)/c
2
(1.2)
其中：x = (x
1
, x
2
) ∈ R
2
，a, b, c 为参数。我们的目标是寻找 x = (x
1
, x
2
) ∈ R
2
，使得 f(x) 最小。
在 MATLAB 中，Optimization Toolbox 使用下面 3 种方法求解上述无约束非线性最小化问
题：
1
http://www.ma-xy.com
1.3 模型规范化及基本理论 第一章 无约束非线性规划
1. 拟牛顿法：使用二次或三次线搜索技术并使用 BFGS 公式来计算海赛矩阵；
2. Nelder-Mead 算法：只使用目标函数值来直接寻找最优点，可用于处理非平滑目标函数；
3. 信赖域法：特别适用于有稀疏矩阵或其他结构的大规模问题。
用 MATLAB 解上述问题 (1.2)：
1 y = @(x , a , b , c )
2 (x ( 1 )−a) .* exp(−((x ( 1 )−a) .^2+(x (2)−b) . ^ 2 ) ) +((x (1)−a)^2+(x ( 2 )−b) ^2) / c ;
3 a = 0;
4 b = 0;
5 c = 20;
6 obj = @(x ) y ( x , a , b , c ) ;
7 e 2s u rf c ( obj , [ −2 ,2]) ;
8 x0 = [ −0 . 5 ; 0 ] ;
9 opti o ns = optimset ( ’ fminunc ’ , ’ Algorithm ’ , ’ quaci – newton ’ )
10 opti o ns . Display = ’ i t e r ’ ; %点 索 引
11 [ x , fv a l , e x i t f l a g , output ] = fminunc ( obj , x0 , optio n s )
12 %% 绘 制 等 高 线 图 以 及 迭 代 动 图
13
1.3 模型规范化及基本理论
将前面的例子 (1.1) 和 (1.2) 忽略参数，写成一般形式
min
x∈R
n
f(x)
上式即为无约束非线性优化的一般形式，我们的目标是在 x ∈ R
n
中求最优 x
∗
，使得 f 最
小。我们在数学分析中已经学过当 n = 1, 2, 3 时，某些特定的 f 的极值点的求解。要注意的是，
我们是求某一类型的 f(如 f 连续可微) 的极值点而并非任何一个 f 。下面，我们将在 R
n
中讨论
f(x) 极值问题。
1.3.1 点的性质
首先，我们给出 n 维空间 R
n
中点 x 的一些定义：
定义 (半范数) 定义映射 ∥ · ∥ : R
n
→ R 为 R
n
上的半范数，当且仅当它具有如下性质
(i) ∥x∥ ≥ 0, ∀x ∈ R
n
；
(ii) ∥αx∥ = |α|||x||, ∀α ∈ R, ∀x ∈ R
n
；
(iii) ∥x + y∥ ≤ ∥x∥ + ∥y∥, ∀x, y ∈ R
n
。
定义 (范数) 映射 ∥ · ∥ : R
n
→ R 为 R
n
上的范数，当且仅当它是半范数，且
∥x∥ = 0 ⇔ x = θ
其中：θ 为 R
n
中的 0 点。
http://www.ma-xy.com 2 http://www.ma-xy.com
http://www.ma-xy.com
第一章 无约束非线性规划 1.3 模型规范化及基本理论
设 x = (x
1
, x
2
, . . . , x
n
)
T
∈ R
n
，常用的向量范数有
∥x∥
∞
= max
i
|x
i
| (l
∞
范数)
∥x∥
1
=
n
i=1
|x
i
| (l
1
范数)
∥x∥
2
=
n
i=1
x
2
i
(l
2
范数)
∥x∥
p
=
n
i=1
|x
i
|
p
1
p
(l
p
范数)
类似于向量范数的定义，我们可以定义矩阵范数。设 A ∈ R
m×n
，其诱导矩阵范数为
∥A∥ = max
x=0
∥Ax∥
∥x∥
其中：∥x∥ 是某一向量范数。特别地，有
∥A∥
1
= max
j
{∥a
·j
∥
1
} = max
j
n
i=1
|a
ij
|
∥A∥
∞
= max
i
{∥a
i·
∥
1
} = max
i
n
j=1
|a
ij
|
∥A∥
2
= ( λ
A
T
A
)
1
2
其中：λ
A
T
A
表示 A
T
A 的最大特征值，a
·j
表示 A 的第 j 列，a
i·
表示 A 的第 i 行。显然
∥A
−1
∥ =
1
∥A∥
∥I∥ = 1
1.3.2 目标函数的性质
下面来讨论某些特定的 f : R
n
→ R。首先，我们给出函数导数的定义。
定义 (一阶导数) 如果
∂f
∂x
i
(x), i = 1, 2, . . . , n 存在且连续，则函数 f : R
n
→ R 在 x ∈ R
n
连续可微。如果 f 在开集 D ⊂ R
n
中的每一点连续可微，则称 f 在 D 中连续可微，记为
f ∈ C
1
(D)。定义 f 在 x 处的梯度为
g = ∇f(x) =
∂f
∂x
1
(x), ··· ,
∂f
∂x
n
(x)
T
定义 (二阶导数) 如果
∂
2
f
∂x
i
∂x
j
(x), i = 1, 2, . . . , n 存在且连续，则函数 f : R
n
→ R 在
x ∈ R
n
二次连续可微。如果 f 在开集 D ⊂ R
n
中的每一点二次连续可微，则称 f 在 D 中二次
连续可微，记为 f ∈ C
2
(D)。定义 f 在 x 处的 Hessian 矩阵为
G = [∇
2
f(x)]
ij
=
∂f(x)
∂x
i
∂x
j
1 ≤ i, j ≤ n
http://www.ma-xy.com 3 http://www.ma-xy.com
http://www.ma-xy.com
1.3 模型规范化及基本理论 第一章 无约束非线性规划
设 f : R
n
→ R 在开集 D ⊂ R
2
上连续可微。对于 x ∈ R
n
, d ∈ R
n
，f 在 x 点关于方向 d 的
方向导数定义为
∂f
∂d
(x) = lim
θ→0
f(x + θd) − f(x)
θ
上述定义的方向导数等于 ∇f (x)
T
d。其中，∇f(x) 为梯度，它是 f 的导数 f
′
(x) 的转置，是 n ×1
向量。
设 f : R
n
→ R 在开集 D ⊂ R
2
上连续可微，对于 x ∈ R
n
, d ∈ R
n
，f 在 x 点关于方向 d 的
方向导数定义为
∂
2
f
∂d
2
(x) = lim
θ→0
∂f
∂d
(x + θd) −
∂f
∂d
(x)
θ
上述定义的二阶方向导数等于 d
T
∇
2
f(x)d。其中，∇
2
f(x) 是 f 在 x 点的 Hessian 矩阵。
1.3.3 最优性条件
上面给出了向量 (矩阵) 范数、f 连续可微、二次连续可微的定义。下面，我们将在这些定义
的基础上给出无约束非线性优化问题的最优性条件。
回到我们前面的目标：求 x ∈ R
n
，使 f (x) 最小。我们自然会问：什么是最小，什么情况
下最小 (极小值存在性)。首先，我们给出极小值的两种类型的定义：局部极小点和总体极小点。
f : x ∈ R → R 的极小点类型示意图如图 (1.1) 所示
图 1.1: 极小点类型示意图
定义 (局部极小点) 如果 ∃f > 0, ∀x ∈ N(x
∗
, δ) = {x ∈ R
n
|∥x − x
∗
∥ < δ)}，有
f(x
∗
) ⩽ f (x)
则称 x
∗
为 f 的一个局部极小点。若 f (x
∗
) < f (x) 且 x = x
∗
，则称 x
∗
为 f 的一个严格局部极
小点。
定义 (总体极小点) 如果 ∀x ∈ R
n
都有
f(x
∗
) ⩽ f (x)
则称 x
∗
为 f 的一个总体极小点。若 f (x
∗
) < f (x) 且 x = x
∗
，则称 x
∗
为 f 的一个严格总体极
小点。
http://www.ma-xy.com 4 http://www.ma-xy.com
http://www.ma-xy.com
第一章 无约束非线性规划 1.4 算法框架
应该指出，实践中，我们只是求一个局部极小点而非全局 (总体) 极小点，因为总体极小点
往往是难以求解的。后面大部分章节讨论局部极小点，在全局优化及智能优化中讨论的是全局极
小点。并且，应当注意的是：当问题具有某种凸性时，局部极小点就是总体极小点。
下面给出 (局部) 极小点存在的充分必要条件 (即最优性条件/极小点解的存在性)。为方便书
写，记
g(x) = ∇f(x), g
k
= ∇f(x
k
)
G(x) = ∇
2
f(x), G
k
= ∇
2
f(x
k
)
一阶必要条件：设 f : D ⊂ R
n
→ R 在开集 D 上连续可微。若 x
∗
∈ D 是局部极小点，则
g(x
∗
) = 0
二阶必要条件：设 f : D ⊂ R
n
→ R 在开集 D 上二次连续可微。若 x
∗
∈ D 是局部极小点，
则
g(x
∗
) = 0, G(x
∗
) ≥ 0
二阶充分条件：设 f : D ⊂ R
n
→ R 在开集 D 上二次连续可微 (f ∈ C
2
(D))，若 g(x
∗
) = 0
并且 G(x
∗
) 是正定矩阵，则 x
∗
∈ D 是 f 的一个严格极小点。
一般的，目标函数的稳定点 (g(x) = 0) 不一定是极小点。但当 f 为凸函数时，其稳定点、局
部极小点和总体极小点是等价的。
设 f : R
n
→ R 是凸函数，且 f ∈ C
1
，则 x
∗
是总体极小点的充分必要条件是 g(x
∗
) = 0。
1.4 算法框架
最优化方法通常采用迭代方法进行求解。我们给定一个初始点 x
0
，然后按照某一方向以某
个大小的步子去靠近 x
∗
。当然，我们会迭代许多次以靠近 x
∗
。在这个过程中会产生一系列的点，
我们将其放在一起，记为 {x
k
}
n
k=1
(设迭代 n 次，x
0
为初始点)。后面我们会研究序列 {x
k
} 的性
质。
设 x
k
是第 k 次迭代的搜索点，d
k
是第 k 次迭代的搜索方向，α
k
是第 k 次迭代的步长，则
第 k + 1 次的搜索点可表示为
x
k+1
= x
k
+ α
k
d
k
从这个表达式中可以看出，不同的 α
k
, d
k
决定了 x
k+1
，并由此构成了不同的方法 (这个留在后面
讨论)。在最优化方法中，搜索方向 d
k
是 f 在 x
k
处的下降方向，即
∇f(x
k
)d
k
< 0
或者说
f(x
k
+ α
k
d
k
) < f (x
k
)
http://www.ma-xy.com 5 http://www.ma-xy.com
http://www.ma-xy.com
1.4 算法框架 第一章 无约束非线性规划
按照上面的设计思路，我们给出如下的最优化算法的基本结构：
step1. 确定初始点 x
k
，设置算法的终止准则 A；
step2. 确定搜索方向 d
k
；
step3. 确定搜索步长 α
k
；
step4. 更新搜索点
x
k+1
= x
k
+ α
k
d
k
step5. 更新 x
k+1
是否满足终止准则 A，不满足就返回 step2，k ← k + 1。
按照上面的算法框架，我们可以得到一系列搜索点 {x
k
}。下面，我们来分析一下 {x
k
}。我们
希望通过最优化算法求解的最优解能够足够接近真实极小点 x
∗
(虽然更多时候极小点是未知的，
但我们设其为 x
∗
)。下面的问题是如何衡量一个序列 x
k
与 x
∗
的接近程度？
如果
lim
k→∞
∥x
k
− x
∗
∥
p
= 0
我们称 {x
k
} 在 p 阶范数下收敛于 x
∗
。并且更多情况，我们讨论 p = 2 时的 L
2
范数。上面给出
了确定序列的收敛性 (当然，有随机下有概率收敛的概念)。下面考虑如果我们有多个算法，那么
哪个算法收敛的精度高且收敛速度快呢？
精度的问题即为上述的收敛性。而对于收敛速度，我们需要讨论序列 {x
k
} 收敛到 x
∗
的速
度。考虑 x ∈ R 的情况。定义 {x
k
} 到 x
∗
的收敛速度：设 {x
k
} 收敛于 x
∗
，且 ∃α > 0 和常数
C，使得
lim
k→∞
∥x
k+1
− x
∗
∥
∥x
k
− x
∗
∥
α
= C
则称 {x
k
} 具有 Q − α 阶收敛速度。特别地
1. 当 α = 1, C > 0 时，叫做线性收敛速度；
2. 当 1 < α < 2, C > 0 或者 α = 1, q = 0 时，叫做超线性收敛速度；
3. 当 α = 2 时，叫做二阶收敛速度。
一般认为具有二阶收敛速度或超线性收敛速度的算法是比较好的。当然，许多时候，我们需
要知道 x
∗
，因此，有必要再讨论收敛性和收敛速度。在前面的算法框架中，我们提到了算法的终
止准则，下面，我们来具体看一些终止准则：
¬当目标值的差量足够小时，终止。
|f(x
k+1
) − f(x
k
)| ≤ ε
或者
|f(x
k+1
) − f(x
k
)|
f(x
k
)
≤ ε
http://www.ma-xy.com 6 http://www.ma-xy.com
http://www.ma-xy.com
第一章 无约束非线性规划 1.5 搜索步长的确定：一维搜索 (线搜索)
当搜索值的差量足够小时，终止。
∥x
k+1
− x
k
∥ ≤ ε
或者
|(x
k+1
)
i
− (x
k
)
i
| ≤ ε
i
∀i ∈ n
或者
∥x
k+1
− x
k
∥
∥x
k
∥
≤ ε
®n 足够大时，终止。
n ≥ T
max
¯f, x 的差量都足够小时，终止。
∥x
k+1
− x
k
∥ ≤ ε
1
& |f(x
k+1
) − f(x
k
)| < ε
1
其中：ε, ε
1
为允许误差。
上面，我们讨论了算法的基本框架、算法收敛性、算法收敛速度和一些停机准则，但在算法的
基本框架中，我们仅给出了数值优化的整体框架，并没有涉及到细节：如何设计 x
k+1
= x
k
+α
k
d
k
。
其中，d
k
是方向向量，α
k
是步长向量。下面，我们将来讨论 α
k
和 d
k
的确定。首先是 α
k
，这
部分内容称为线搜索或一维搜索。因为我们要在每步 k 中寻找最优的 α
k
，使得目标值最小。然
后，我们讨论 d
k
，d
k
是使目标函数下降的方向。
1.5 搜索步长的确定：一维搜索 (线搜索)
考虑最优算法框架中搜索点的更新公式
x
k+1
= x
k
+ α
k
d
k
现在，我们来确定 α
k
。这里，我们假设 d
k
与 α
k
无关，并且 d
k
已知，即我们已经知道 d
k
的具
体方向 (在后面部分将介绍 d
k
的求解)。那么，我们的目标就变为求 α
k
，使目标函数 f 最小。
min
α>0
f(x
k
+ αd
k
)
为了书写简单，由于 x
k
, d
k
为已知量，故将 f(x
k
+ αd
k
) 记为 φ(α)
α
k
= arg min
α>0
φ(α)
如果 α
k
满足 min φ(α)，则称 α
k
为最优步长。这样的一维搜索为最优一维搜索或精确一维
搜索。但很明显的是：在实际计算中，精确 α
k
是难以求解的，或者需要很大的计算量。所以我
们需要一个非精确的 α
k
：我们自然希望 α
k
只要能使目标函数值下降就好了，即
φα
k
< φ (0)
http://www.ma-xy.com 7 http://www.ma-xy.com
http://www.ma-xy.com
1.5 搜索步长的确定：一维搜索 (线搜索) 第一章 无约束非线性规划
或者
f(x
k
+ α
k
d
k
) < f (x
k
)
如果 α
k
满足上面要求，则称 α
k
为粗步长。这样一维搜索为近似一维搜索或不精确一维搜索又
或称可接受一维搜索。
毫无疑问，在实际求解过程中 φ(α) 会有许多不同的形式：例如 φ(α) 光滑、φ(α) 可以写出
具体的表达式、可以关于 α 求导、φ(α) 单调或者 φ(α) 存在唯一极值 (凸性)，再或 φ(α ) 有许多
极值等等。
上面给出了求 α
k
的精确线性搜索和非精确线性搜索。下面，我们先来讨论精确线搜索，然
后讨论非精确线性搜索。精确线性搜索可以分为两类：一类是使用导数 φ
′
(α) 的搜索，如牛顿法，
插值法等等；另一类是不使用导数，仅依靠函数值 φ(α) 的搜索，如黄金分割法和 Fibonacci 法。
而非精确线性搜索有依靠 Wolfe 准则和 Armijo 准则的方法。
1.5.1 黄金分割法
黄金分割法也称 0.618 法，是一种基于区间收缩技术的搜索方法。所谓的区间收缩技术是指：
先确定包含最小值的搜索区间，然后不断分割这个区间，使最小值所在的区间越来越小。当区间
长度缩小至一定程度时，区间上各点的函数值均接近极小值，可作为最小值的近似。这种方法尤
其适合于非光滑 φ(α) 和导数 φ
′
(α) 复杂甚至写不出的情形。黄金分割法如图 (1.2) 所示
图 1.2: 黄金分割法示意图
设 φ(α) 是初始搜索区间 [a
1
, b
1
] 上的凸函数，如图 (1.2) 所示。我们要不断收缩 [a
1
, b
1
] 以
接近 α
∗
(问：[a
1
, b
1
] 如何确定)。设第 k 次分割的区间为 [a
k
, b
k
](α
∗
∈ [a
k
, b
k
])。取两个试操点
λ
k
, µ
k
∈ [a
k
, b
k
]，且 λ
k
< µ
k
，计算 φ(λ
k
) 和 φ(µ
k
)
(1) 若 φ(λ
k
) ≤ φ(µ
k
) 则令 a
k+1
= α
k
, b
k+1
= µ
k
，以形成 [α
k
, µ
k
] = a
k+1
, b
k+1
的新搜索空
间；
(2) 若 φ(λ
k
) > φ(µ
k
) 则令 a
k+1
= λ
k
, b
k+1
= b
k
，以形成 [λ
k
, b
k
] = a
k+1
, b
k+1
的新搜索空
间。
我们要求 λ
k
, µ
k
满足以下条件：
(1) b
k
− λ
k
= µ
k
− a
k
；
(2) b
k+1
− a
k+1
= τ(b
k
− a
k
)，其中，τ 为收缩率。
http://www.ma-xy.com 8 http://www.ma-xy.com
http://www.ma-xy.com
第一章 无约束非线性规划 1.5 搜索步长的确定：一维搜索 (线搜索)
由 (1)(2) 可知
λ
k
= a
k
+ (1 − τ)(b
k
− a
k
)
µ
k
= a
k
+ τ(b
k
− a
k
)
0.618 法即 τ = 0 .618 的分割方法。下面给出黄金分割法的步骤：
step1. 初始化。[a
1
, b
1
], k = 1，精度要求 ϵ > 0，τ = 0.618。
step2. 计算 λ
1
, µ
1
。
λ
1
= a
1
+ (1 − τ)(b
1
− a
1
)
µ
1
= a
1
+ τ(b
1
− a
1
)
step3. 确定新的搜索空间。
3.1). 计算 φ(λ
k
), φ(µ
k
)。
3.2). 若 φ(λ
k
) > φ(µ
k
)：
若 b
k
− λ
k
≤ ϵ，则停止，得到 µ
k
；
否则，令 a
k+1
:= λ
k
, b
k+1
:= b
k
, λ
k+1
:= µ
k
φ(λ
k+1
) = φ(µ
k
), µ
k+1
:= a
k+1
+ 0.618(b
k+1
− a
k+1
)
计算 φ(µ
k+1
)；
令 k = k + 1，返回 step3.2)。
若 φ(λ
k
) ⩽ φ(µ
k
)：
若 µ
k
− a
k
≤ ϵ，则停止，得到 λ
k
；
否则，令 a
k+1
:= a
k
, b
k+1
:= µ
k
, µ
k+1
:= λ
k
φ(µ
k+1
) = φ(λ
k
), λ
k+1
:= a
k+1
+ 0.382(b
k+1
− a
k+1
)
计算 φ(λ
k+1
)；
令 k := k = 1，返回 step3.2)。
前面，我们假设 φ(α) 在 [a
,
b
1
] 上是某种凸函数，如果不是凸函数，则需要一定的改进，参
考《最优化理论与方法》P72. 袁亚湘，孙文瑜。
1.5.2 Fibonacci 法
Fibonacci 法与 0.618 法相近，区别在于它的收缩率 τ 不是黄金分割数 0.618，而是使用
Fibonacci 数。Fibonacci 数列满足
F
0
= F
1
= 1
F
k+1
= F
k
+ F
k−1
k = 1, 2, . . .
http://www.ma-xy.com 9 http://www.ma-xy.com
http://www.ma-xy.com
1.5 搜索步长的确定：一维搜索 (线搜索) 第一章 无约束非线性规划
令 τ =
F
n−k
F
n−k+1
，故
λ
k
= a
k
+
1 −
F
n−k
F
n−k+1
(b
k
− a
k
)
= a
k
+
F
n−k−1
F
n−k+1
(b
k
− a
k
) k = 1, 2, . . . , n − 1
µ
k
= a
k
+
F
n−k
F
n−k+1
(b
k
− a
k
) k = 1, 2, . . . , n − 1
1.5.3 二次插值法
插值简介
我们只知道 φ(α) 的点列，如果我们假设 φ(α) 是一个多项式函数 (如二次函数)，即用多项
式来逼近 φ(α)，则称该方法为 φ(α) 的插值方法。
设 φ(α) 的形式是一个二次多项式的形式，并将其设为
q(α) = aα
2
+ bα + c
其中：a, b, c 为待定参数 (待求参数)。3 个参数要 3 个方程，而这 3 个方程的信息来自于已知
φ(α) 点列。
(1) 如果我们仅给出 1 个点 α
1
，并且知道了 φ(α
1
)，φ
′
(α
1
)，φ
′′
(α
1
) 的值，则要求 q(α) 在
点 α
1
满足这 3 个值的情况。并称该方法为 1 点二次插值法，又称牛顿法。
(2) 如果我们给出 2 个点 α
1
, α
2
，并且知道了 φ(α
1
)，φ(α
2
)，φ
′
(α
1
) 或者 φ(α
1
), φ
′
(α
1
), φ
′
(α
2
)
的值，则要求 q(α) 在 2 点处满足相应的数值条件。称该方法为 2 点二次插值法又称割线法。
(3) 如果我们给出 2 个点 α
1
, α
2
, α
3
，并且知道了 φ(α
1
), φ(α
2
), φ(α
3
), . . . 称该方法为 3 点二
次插值法又称割线法。
1 点二次插值法 (牛顿法)
设 q(α) = aα
2
+ bα + c，由
q(α
1
) = φ(α
1
)
q
′
(α
1
) = φ
′
(α
1
)
q
′′
(α
1
) = φ
′′
(α
1
)
可以求解 a, b, c
a = φ
′′
(α
1
)/2
b = φ
′
(α
1
) − φ
′′
(α
1
)α
1
则 q(α) 的最小值为
α
2
:= −
b
2a
= α
1
−
φ
′
(α
1
)
φ
′′
(α
1
)
http://www.ma-xy.com 10 http://www.ma-xy.com
http://www.ma-xy.com
第一章 无约束非线性规划 1.5 搜索步长的确定：一维搜索 (线搜索)
由此即得牛顿法的迭代公式：
α
k+1
= α
k
−
φ
′
(α
1
)
φ
′′
(α
1
)
注：这里的 k 并非 x
k+1
= x
k
+ α
k
d
k
，而是 α 的搜索序列。
下面给出牛顿法的局部二阶收敛速度。假定 φ : R → R , φ ∈ C
2
, φ
′
(α
∗
) = 0，则当初始点 α
0
充分接近 α
∗
时，由牛顿法迭代
α
k+1
= α
k
−
φ
′
(α
k
)
φ
′′
(α
k
)
k = 0, 1, . . .
产生的点列 {α
k
} 收敛，即 α
k
→ α
∗
。若 φ ∈ C
3
，则
lim
k→∞
|α
k+1
− α
∗
|
|α
k
− α
∗
|
=
1
2
φ
′′′
(α
∗
)
φ
′′
(α
∗
)
这表明
|α
k+1
− α
∗
| = O(|α
k
− α
∗
|
2
)
2 点二次插值法
¬设 q(α) = aα
2
+ bα + c，已知 α
1
, α
2
处的函数值 φ(α
1
), φ(α
2
), φ
′
(α
1
)，构建三元方程组，
解 a, b, c 有
a
=
φ
1
− φ
2
− φ
′
1
(α
1
− α
2
)
−(α
1
− α
2
)
2
b = φ
′
1
+ 2
φ
1
− φ
2
− φ
′
1
(α
1
− α
2
)
(α
1
− α
2
)
2
α
1
并且有
α
3
:= α
1
−
(α
1
− α
2
)φ
′
1
2[φ
′
1
−
φ
1
−φ
2
α
1
−α
2
]
写为迭代格式，有
α
k+1
= α
k
−
(α
k
− α
k−1
)φ
′
k
2[φ
′
k
−
φ
k
−φ
k−1
α
k
−α
k−1
]
设 q(α) = aα
2
+ bα + c，已知 α
1
, α
2
处的函数值 φ(α
1
)，φ
′
(α
1
)，φ
′
(α
2
) 的，构建三元方
程组，解 a, b, c，有
α
3
:= −
b
2a
= α
1
−
α
1
− α
2
φ
′
(α
1
) − φ
′
(α
2
)
φ
′
(α
1
)
写为迭代格式，有
α
k+1
= α
k
−
α
k
− α
k−1
φ
′
(α
k
) − φ
′
(α
k−1
)
φ
′
(α
k
)
下面，我们给出二点二次插值法的收敛速度。设 φ(α) 存在三阶连续导数，φ(α) ∈ C
3
，α
∗
满足 φ
′
(α
∗
) = 0, φ
′′
(α
∗
) = 0，则割线法迭代产生的序列 {α
k
} 收敛到 {α
∗
}，且其收敛速度的阶
为
1+
√
5
2
≈ 1.618。
http://www.ma-xy.com 11 http://www.ma-xy.com
http://www.ma-xy.com
1.5 搜索步长的确定：一维搜索 (线搜索) 第一章 无约束非线性规划
3 点二次插值法
设 q(α) = aα
2
+ bα + c，已知 α
1
, α
2
, α
3
处的函数值 φ(α
1
), φ(α
2
), φ(α
3
)，构建三元方程组，
解 a, b, c 有 α
4
α
4
=
1
2
(φ
1
+ φ
2
) +
1
2
(φ
1
− φ
2
)(φ
2
− φ
3
)(φ
3
− φ
1
)
(α
2
− α
3
)φ
1
+ (α
3
− α
2
)φ
2
+ (α
1
− α
2
)φ
3
上面的公式可以直接由拉格朗日插值公式 得到
L(α) =
(α − α
1
)(α − α
3
)
(α
1
− α
2
)(α
1
− α
3
)
φ
1
+
(α − α
1
)(α − α
3
)
(α
2
− α
1
)(α
2
− α
3
)
φ
2
+
(α − α
1
)(α − α
2
)
(α
3
− α
1
)(α
3
− α
1
)
φ
3
令 L
′
(α) = 0 即可得到。
下面讨论三点二次插值法的收敛速度。设 φ(α) 存在四阶连续导数，φ(α) ∈ C
4
，φ
′
(α
∗
) =
0, φ
′′
(α
∗
) = 0，则三点二次插值法产生的序列 {α
k
} 收敛速度约为 1.32。
1.5.4 三次插值法
三次插值法使用一个三次四项式 q(α) = c
1
(α −a)
3
+ c
2
(α −a)
2
+ c
3
(α −a) + c
4
来逼近 φ(α)。
这个 q(α) 中的待定系数 φ 要有 4 个条件来确定，我们可以利用四点的函数值 φ
1
, φ
2
, φ
3
, , φ
4
，也
可以利用三点函数值加一点的导数值 φ
1
, φ
2
, φ
3
, φ
′
4
等等。三次插值法比二次插值法有较好的收
敛效果。但通常要求计算导数组，且工作量比二次插值法大，所以，当导数容易求的时候，用三
次插值法较好。
上面讨论了基于区间收缩技术的黄金分割法和 Fibonacci 法，以及基于插值法的二次插值和
三次插值法，并讨论了算法的收敛性。但上面讨论的都是精确一维搜索方法，求 α
∗
往往需要很
大的计算量。下面，我们介绍非精确的一维搜索方法：Armijo-Goldstein 准则和 Wolfe-Powell 准
则。
1.5.5 Armijo-Goldstein 准则
Armijo(1966) 和 Goldstein(1965) 分别提出不精确一维搜索过程，我们现在不要求 α
∗
使
φ(α) 最小。我们要求 α 使得 φ(α) < φ(0) 即可 (φ(α) = f (x
k
+ αd
k
))。设
J = {α > 0
φ(α) < φ(0)}
J 是一个 α 的区间集合，则 J 中的点 α 是我们可以取值的。虽然 J 中的每一点皆是可取的，但
我们这里要选择一个“恰当的”。为此，我们需要设置一些准则来约束 J。如图 (1.3) 所示
http://www.ma-xy.com 12 http://www.ma-xy.com
http://www.ma-xy.com
第一章 无约束非线性规划 1.5 搜索步长的确定：一维搜索 (线搜索)
图 1.3: Armijo 准则示意图
¬如图 (1.3) 所示，J = [0, a]，为了保证目标函数单调下降，同时 φ 的下降不是太小 (如果太
小，可能导致序列 {f (x
k
)} 的极限值不是极小值)。必须避免所选择的 α 太靠近区间 J 的端点。
我们要求
φ(α
k
) ⩽ φ(0) + ρφ(α
k
)φ
′
(0) (1.3)
其中：ρ ∈ (0,
1
2
)。满足式 (1.3) 的 α
k
指代的区间是 J
1
= [0 , c]。
为避免 α 太小，我们加上另外一个要求
φ(α
k
) > φ(0) + (1 − ρ)φ(α
k
)φ
′
(0) (1.4)
满足式 (1.4) 的多项式 J
2
= [b, a]。综合 (1.3)(1.4)，则有 J
∗
= [b, c]，我们将 (1.3)(1.4) 称为
Armijo-Goldstein 不精确线性搜索准则。一旦 α ∈ J
∗
，则称 α 为可接受步长，J
∗
为可接受区间。
®如图1.3所示，A-G 准则可能将 α
∗
排除在 J
∗
之外。为此，Wolfe-Powell 准则给出了一个
更简单的条件来代替 (1.4)
φ
′
(α
k
) = σφ
′
(0) > φ
′
(0) (1.5)
解释为：可接受点 α
k
处的斜率大于等于初始点斜率的 σ 倍。满足 (1.5) 的 α
k
构成区间为
J
3
= [ e, a]。
注：φ
′
(0) = g
T
k
d
k
。
¯式 (1.5) 的不足在于，即使 σ → ∞ 时，也不能带来精确的线性搜索。但是，如果要求
|φ
′
(α
k
)| ⩽ σ|φ
′
(0)| (1.6)
则可以满足。一般的，σ 值越小，线性搜索越精确，满足 (1.6) 的 α
k
的区间为 J
4
。
综合 (1.3)(1.5)，则有 J
∗
= [e, c]。我们称 (1.3)(1.5) 为强 Wolfe-Powell 准则；称 (1.3)(1.6)
为 α 强 Wolfe-Powell 准则。
下面，我们给出 Armijo-Goldstein 不精确一维搜索方法的步骤：
step1. 初始化。在搜索区间 J
∗
= [0 , a] 中取定初始点 α
0
，计算 φ(0), φ
′
(0)。给出 ρ ∈ (0,
1
2
), t ⩾ 1，
令 a
0
:= 0 , b
0
:= a, k := 0。
step2. 检验准则 (1.3)。计算 φ(α
k
)，若
φ(α
k
) ⩽ φ(0) + ρα
k
φ
′
(0)
http://www.ma-xy.com 13 http://www.ma-xy.com
http://www.ma-xy.com
1.6 MATLAB 应用实例 1 第一章 无约束非线性规划
转到 step3，否则，令 α
k+1
:= α
k
, b
k+1
:= α
k
，转到 step4。
step3. 检验准则 (1.4)，若
φ(α
k
) > φ(0) + (1 − ρ)α
k
φ
′
(0)
停止迭代，输出 α
k
；否则，令 α
k+1
:= α
k
, b
k+1
:= α
k
，若 b
k
< a 转到 step4，否则，令 α
k+1
:=
tα
k
, k := k + 1，转到 step2。
step4. 取新的搜索点，取
α
k+1
=
α
k+1
+ b
k+1
2
令 k := k + 1，转到 step2。
关于不精确以为搜索的收敛性，可以参考《最优化理论与方法》P102。
总结：前面我们给出了求解 x
k+1
= x
k
+ α
k
d
k
中 α
k
的精确线搜索方法：黄金分割法、Fibonacci
法、插值法以及非精确线搜索方法：A-G 准则，W-P 准则。后面，我们要给出求解方向 d
k
的方
法。在此之前，我们给出 MATLAB 求一维无约束极值函数 fminbnd 的用法。
1.6 MATLAB 应用实例 1
MATLAB 中使用 fminbnd 函数命令来求解无约束一维极值优化，其调用格式为
[x,fval,exitag,output]=fminbnd(fun,
x
1
,
x
2
,options)
其中：fun 为目标函数，可以使句柄，可以是匿名函数也可以是函数文件；x
1
, x
2
为 x 的搜索区
间；options 为结构体；x 为极小点；fval 为极小值；exitag 返回函数 fminbnd 的求解状态：成
功/失败；output 返回 fminbnd 的求解信息：迭代次数、优化算法等。
fminbnd 函数的应用实例如下
1 a = 9/7 ;
2 fun = @(x ) s in (x−a ) ;
3 x1 = 1 ;
4 x2 = 2* pi ;
5 opti o ns = o p tion s e t ( ’ Display ’ , ’ i t e r ’ ) ;
6 option . PlotFcns = @optimplotfval ;
7 [ x , fva l , e x i t f l a g , output ] = fminbnd ( fun , x1 , x2 , opti o ns )
8
1.7 搜索方向的确定
下面，将介绍一些求解方向 d
k
的方法。
http://www.ma-xy.com 14 http://www.ma-xy.com
http://www.ma-xy.com
第一章 无约束非线性规划 1.7 搜索方向的确定
1.7.1 最速下降法
最速下降法是以负梯度方向为搜索方向向量。设 f (x) 在 x
k
附近连续可微，且 g
k
= ∇f(x
k
) =
0，由泰勒公式展开
f(x) = f (x
k
) + (x − x
k
)
T
∇f(x
k
) + o(∥x − x
k
∥)
如果将 x 换为 x
k
+ αd
k
，有
f(x
k
+ αd
k
) = f (x
k
) + αg
T
k
d
k
+ o(∥αd
k
∥)
则 f 在 x
k
处沿方向 d
k
的变化率为
lim
α→0
f(x
k
+ αd
k
) − f(x
k
)
α
= g
T
k
d
k
由 Gauchy-Schwartz 不等式
|d
T
k
g
k
| ⩽ ∥d
k
∥ ∥g
k
∥
所以，当且仅当 d
k
= −g
k
时，d
T
k
g
k
最小，从而有 −g
k
是最速下降方向。迭代格式为
x
k+1
= x
k
− α
k
g
k
上面给出了 d
k
= −g
k
，若对 α
k
采用精确线搜索，则有
φ(α) ≡ f (x
k
+ α
k
d
k
) = min
α⩾0
f(x
k
+ αd
k
)
α
k
应满足 φ
′
(α) =
d
dα
f(x
k
+ αd
k
)
α=α
k
= ∇f(x
k
+ α
k
d
k
)
T
d
k
= 0 。由 d
k
= −g
k
= ∇f(x
k
) 有
g
T
k+1
g
k
= d
T
k+1
d
k
= 0
这表明，在相邻的两个迭代点 k, k + 1，函数 f(x) 的两个梯度方向是相互正交的，也就是迭
代点列所走的路线是锯齿形的。当接近极小点时步长愈小，前进愈慢，至多有线性收敛速度。关
于最速下降法的收敛性和收敛速度我们不再讨论。
1.7.2 牛顿法
牛顿法的基本思想是利用目标函数 f (x) 在收敛处的二次 Taylor 展开来逼近 (代替)f (x)。求
x
k+1
使二次泰勒展开函数最小，以做为更新公式。
假设 f (x) 是二次连续可微函数，x
k
∈ R
n
，且 Hesse 矩阵 H ≜ G(x) = ∇
2
f(x) 是正定的。
我们在 x
k
附近用二次泰勒展开近似 f
f(x
k
+ s) ≈ q
(k)
(s) = f (x
k
) + ∇f(x
k
)
T
s +
1
2
s
T
∇
2
f(x
k
)s
= f(x
k
) + g
T
k
s +
1
2
s
T
G
k
s
http://www.ma-xy.com 15 http://www.ma-xy.com
http://www.ma-xy.com
1.7 搜索方向的确定 第一章 无约束非线性规划
我们的目标是求 s，使 q
(k)
(s) 最小，求上式右边二次函数 q
(k)
(s) 的稳定点
∇q
(k)
(s) = g
k
+ G
k
s = 0
有 s = −
g
k
G
k
¬
。其中：s = α
k
d
k
。如果我们令 α
k
= 1 ，就可得到更新公式 (迭代公式)
x
k+1
= x
k
−
g
k
G
k
= x
k
− G
−1
k
g
k
下面，给出牛顿法的局部收敛型和二阶收敛速度。设 f ∈ C
2
(R)，x
k
充分靠近 x
∗
: ∇f(x
∗
) = 0。
如果 ∇
2
f(x) 正定，且 Hesse 矩阵 G(x) 满足 Lipsditz 条件
，即 ∃L > 0，使得对 ∀i, j，∀x, y ∈ R
n
，
有
|G
ij
(x) − G
ij
(y)| ≤ L∥x − y∥
其中：G
ij
(x) 是 Hesse 矩阵 G(x) 的 (i, j) 元素。则对一切 k，牛顿迭代法所得到的的序列 {x
k
}
收敛于 x
∗
，并且有二阶收敛速度。
在上述局部收敛性中，我们设 f 在 x
∗
的 Hessan 矩阵是正定的，但这样并不是牛顿法收敛
的必要条件。
值得一提的是，当 x
0
远离 x
∗
时，G
k
不一定是正定的，因此牛顿方向不一定是下降方向，
其收敛性不能保证。这说明恒取 α
k
= 1 的牛顿法是不合适的，所以我们用一维搜索方法来确定
α
k
，并且注意，仅当 {α
k
} 收敛到 1 时，牛顿法才是二阶收敛的，这时牛顿法的迭代公式为
d
k
= −a
−1
k
g
k
x
k+1
= x
k
+ α
k
d
k
称 α
k
= 1 的方法为阻尼牛顿法。
这里，我们给出阻尼牛顿法的算法流程
step1. 初始化。初始点 x
0
∈ R
n
，终止误差 ε > 0，令 k := 0。
step2. 计算 g
k
。若 ∥g
k
∥ ⩽ ε，终止迭代，输出 x
k
；否则转到 step3。
step3. 求解 d
k
。解方程组构造牛顿方向，即解
a
k
d = g
k
求出 d
k
。
step4. 求 α
k
。进行一维线搜索，求 α
k
。
step5. 令 x
k+1
= x
k−1
+ α
k
d
k
, k := k + 1，转到 step2。
下面，我们给出阻尼牛顿法的总体收敛性：设 f : R
n
→ R ∈ C
2
(D)，D 为开凸集。如果
∀x
0
∈ D, ∃m > 0，使得 f (x) 在水平集 L(x
0
) = {x|f (x) ⩽ f(x
0
)} 上满足
u
T
∇
2
f(x)u ⩾ m∥u∥
2
, ∀u ∈ R
n
, x ∈ L(x
0
)
则在精确一维搜索下，带步长 α
k
的牛顿法产生的迭代序列 {x
k
} 满足：
(1) 当 {x
k
} 为有穷点列时，∃k, g(k) = 0；
(2) 当 {x
k
} 为无穷点列时，{x
k
} 收敛到 f 的唯一极小点 x
∗
。
¬
这里要求 G
k
正定且非奇异
注：记 Lipsditz 条件为 Lip
r
(D)，其中，r 为 Lip 常数，x ∈ D。
http://www.ma-xy.com 16 http://www.ma-xy.com
http://www.ma-xy.com
第一章 无约束非线性规划 1.7 搜索方向的确定
1.7.3 修正牛顿法
牛顿法面临的主要困难是 Hesse 矩阵 G
k
不一定是正定的，这时候二次搜索 q
k
(s) 不一
定有极小点 (q
′
= 0, q
′′
> 0)，甚至没有稳定点 (q
′
= 0)。为了克服这个困难，Goldstein 和
Price(1967) 提出了一种修正 G
k
的方法：当 G
k
非正定时，采用最速下降方向 −g
k
，将这种处理
方法与角度准则结合，给出
d
k
=
d
N
k
cos
d
N
k
, −g
k
⩾ η
− g
k
其中，η > 0 是正常数。这种方法能够保证 d
k
总满足
cos
d
N
k
, −g
k
⩾ η
从而算法的收敛性是可以保证的。
Goldfold(1966) 提出了一种修正方法，将 Hesse 矩阵 G
k
改为 G
k
+ V
k
I，其中，V
k
> 0，
G
k
+ V
k
I 正定。较理想的 V
k
取值是：V
k
不要远大于使 G
k
+ V I 正定的最小 V 。
记 E
k
= V
k
I，修正后的 G
k
为
¯
G
k
，下面给出修正牛顿法的算法程序：
step1. 初始化。初始搜索点 x
0
∈ R
n
，容许误差 ε > 0，令 k := 0。
step2. 计算 g
k
。若 ∥g
k
∥ ⩽ ε，终止迭代，输出 x
k
；否则转到 step3。
step3. 计算
¯
G
k
。计算 Hesse 矩阵 G
k
，如果 G
k
正定，V
k
= 0,
¯
G
k
= G
k
+ V
k
I；如果 G
k
非正
定，给出 V
k
，
¯
G
k
= G
k
+ V
k
I。
step4. 计算 d
k
。解
¯
G
k
d = −g
k
，得到 d
k
。
step5. 计算 α
k
。
step6. 计算 x
k+1
。令 x
k+1
= x
k
+ α
k
d
k
, k := k + 1，返回 step2。
上述算法的关键是如何解
¯
G
k
，也即如何解修正矩阵 E。下面给出基于 Cholesky 分解的一种
求
¯
G
k
的策略：先形成 G
k
的 Cholesky 分解 LDL
T
。然后令
¯
G
k
= L
¯
DL。其中，
¯
d
ii
= max{|d
jj
|, δ}，
d
jj
为 D 的对角元素，δ 为某个给定的小正数。但是，如果 G
k
为对称不定矩阵，则其 Cholesky
分解可能不存在。另外，即使这种分解存在，其一般也是不稳定的，因为其矩阵分解的元素可能
是无界的。进一步，当 G
k
仅微小波动时，
¯
G
k
也可能与 G
k
相差很大。
为了克服 Cholesky 分解方法的不稳定性，Gill 和 Murray(1974) 提出了一个数值稳定的处
理方法。对称正定矩阵的 Cholesky 分解可以描述为
d
jj
= g
jj
−
j−1
s=1
d
ss
− l
2
js
l
ij
=
1
d
jj
g
ij
−
j−1
s=1
d
ss
− l
js
l
is
, i ⩾ j + 1
其中：g
ij
表示 G
k
的元素，d
jj
表示 D 的对角元素。
现在，我们要求 Cholesky 分解因子 L 和 D 满足条件：¬D 的所有元素是严格正的；L, D
http://www.ma-xy.com 17 http://www.ma-xy.com
http://www.ma-xy.com
1.7 搜索方向的确定 第一章 无约束非线性规划
的所有元素满足一致有界，也就是说，对 k = 1, 2, . . . , n 和某一个正数 β，要求
d
kk
> δ
|r
ik
| ⩽ β i > k (1.7)
其中：r
ik
= l
ik
√
d
kk
，δ 为某个给定的小正数。
下面，我们来描述一下这个分解的第 j 步，假设修改 Cholesky 分解的前 j − 1 列已经计算
出来，对于 k = 1, 2, . . . , j − 1，式 (1.7) 成立。先计算
r
j
=
ξ
j
−
j−1
s=1
d
ss
l
2
js
其中：ξ
j
取 g
jj
，试验值
¯
d = max r
j
, δ。
为了断定
¯
d
是否可以接受作为
D
的第
j
个元素。我们检验
r
ij
=
l
ij
√
¯
d
是否满足式
(1.7)。
并且由 l
ij
= r
ik
/
d
jj
得到 L 的第 j 列，否则
d
jj
=
ξ
j
−
j−1
s=1
d
ss
l
2
js
其中：ξ
j
= g
jj
+ e
jj
。选取正数 e
jj
使得 max |r
ij
| = β，并且也产生 L 的第 j 列。
上述过程完成时，我们得到了正定矩阵
¯
G
k
的 Cholesky 分解
¯
G
k
= LDL
T
= G
k
+ E
其中：E 是非负对角矩阵，对角元素为 e
jj
。对于给定的 G
k
，这个非负对角矩阵 E 依赖于 β，
Gill 和 Murray 证明：如果 n > 1，则
∥E(β)∥
∞
< (
ξ
β
+ (n − 1)β)
2
+ 2(r + (n − 1)β
2
) + δ
其中：ξ 是 G
k
的非对角元素的最大模，r 是 G
k
的对角元素的最大模。
1.7.4 信赖域方法
在目标函数 f 的鞍点 处，我们有时会遇到 ∇f(x
k
) = 0，∇
2
f(x
k
) 非半定这种特殊情形。此
时，前面的修正牛顿法就不好用了。一种解决策略是用函数的负曲率方向作为搜索方向。保证目
标函数值仍是下降的。
前面介绍的牛顿法的基本思想是在 x
k
附近用二次函数
q
(k)
(s) = f (x
k
) + g
T
k
s +
1
2
s
T
G
k
s
来逼近 f(x)，并以 q
(k)
(s) 的极小点 s
k
修正 x
k
，反复迭代。其迭代公式为
x
k+1
= x
k
+ s
k
http://www.ma-xy.com 18 http://www.ma-xy.com
http://www.ma-xy.com
第一章 无约束非线性规划 1.7 搜索方向的确定
牛顿法显然具有二阶收敛速度，但其仅具有局部收敛性，即只有当 s 充分小时，q
(k)
(s) 才能逼近
f(x)。我们后面说的阻尼牛顿法具有总体收敛性。下面，我们将介绍一种新的方法 - 信赖域法。
这种方法不仅具有总体收敛性，并且也可以解决 Hesse 矩阵 G
k
非正定及 x
k
为鞍点等困难。
在给定一个点 x
k
后，我们先定义一个步长上界 h
k
，并且由此给出 x
k
的一个邻域 Ω
k
Ω
k
= {x| ∥x − x
k
∥ ⩽ h
k
}
Ω
k
称为信赖域。我们假设在 Ω
k
中用 q
k
(s) 来逼近 f (x) 是恰当的，然后求搜索方向 s
k
，使 q
k
(s)
最小，即
min
s
q
k
(s) = f (x
k
) + g
T
k
s +
1
2
s
T
G
k
s
s.t. ∥s∥ ⩽ h
k
其中：h
k
是步长上界，s = ∥x − x
k
∥，范数 ∥ · ∥ 可以用 L
2
范数 L
∞
范数。
值得一提的是，G
k
为 Hesse 矩阵，如果难于计算，可采用后面介绍的方法做近似。假设我们
给出了 h
k
，我们来尝试给出 ∥·∥ = ∥·∥
2
情况的解。如果取 ∥·∥ 为 L
2
范数，则相应的模型变为
min
s
q
k
(s) = f (x
k
) + g
T
k
s +
1
2
s
T
G
k
s
s.t. s
T
s ⩽ h
2
k
¬如果 s
∗
在 s
T
s ⩽ h
2
k
内，则约束不起作用，所以只需要求解
min
s
q
k
(s) = f (x
k
) + g
T
k
s +
1
2
s
T
G
k
s
即可。那么就变为牛顿法：当 G
k
为正定时，s
k
= −G
−1
k
g
k
为解。
如果 s
∗
在边界上有 s
T
s = h
2
k
，则引进 Lagrange 函数。
L(s, λ) = q
(k)
(s) +
1
2
λ(s
T
s − h
2
k
)
上面，我们假设假设我们已经给出了 h
k
，那么，我们该如何求 h
k
？一般地，当 q
(k)
(s) 与
f(x
b
+ s) 之间的一致性满足某种要求时，应选取尽可能大的 h
k
。设 ∆f
k
是 f 在第 k 步的实际
下降量。
∆f
k
= f
k
− f(x
k
+ s
k
)
对应的 q
(k)
(s) 下降量为
∆q
(k)
= f
k
− q
(k)
(s)
定义二者比值
r
k
= ∆ f
k
/∆q
(k)
r
k
衡量了 q
(k)
(s
k
) 与目标 f(x
k
+ s
k
) 的近似程度，r
k
越接近 1，表明近似程度越高。下面给出
信赖域方法的算法程序：
http://www.ma-xy.com 19 http://www.ma-xy.com
http://www.ma-xy.com
1.7 搜索方向的确定 第一章 无约束非线性规划
step1. 初始化。x
0
∈ R
n
, h
0
= ∥g
0
∥, µ =
1
4
, η
=
3
4
, ε >
0
, k
:= 1
。
step2. 给出 x
k
∈ R
n
, h
k
∈ R，计算 g
k
。若 ∥g
k
∥ ⩽ ε，终止迭代，输出 x
k
；否则，计算 G
k
。
step3. 解信赖域模型，求出 s
k
。
step4. 求 f(x
k
+ s
k
) 和 r
k
的值。
step5. 如果 r
k
⩽ µ，令 h
k+1
= ∥s
k
∥/4；如果 r
k
> η 并且 ∥s
k
∥ = h
k
，令 h
k+1
= 2 h
k
；否则，令
h
k+1
= h
k
。
step6. 如果 r
k
⩽ 0，令 x
k+1
= x
k
，否则 x
k+1
= x
k
+ s
k
。令 k := k + 1，返回 step2。
下面，给出信赖域的总体收敛性和收敛速度。
(1) 设 B ⊂ R
n
是有界集，x
k
∈ B, ∀N ，若 f ∈ C
2
在有界集 B 上 ∥G
k
∥
2
⩽ N, M ⩾ 0，则
信赖域算法产生一个满足一阶二阶必要条件的聚点 x
∞
(2) 若零点 x
∞
还满足 f 的 Hesse 矩阵 G
∞
是正定的，那么，对于主序列，有 r
k
→ 1, x
k
→
x
∞
, g(b(x
k
)) > 0，以及充分大的 k，约束 ∥s∥
2
< h
k
，此外收敛速度是二阶的。
1.7.5 共轭梯度法
共轭方向法是介于最速下降法和牛顿法之间的一个方法，它仅需要一阶导数信息，但克服了
最速下降法收敛速度慢的特点，又避免了牛顿法计算存贮二阶导数信息的麻烦。典型的共轭方向
法有共轭梯度法和拟牛顿法。我们先来介绍共轭方向法。
设 G ∈ R
n×n
是对称正定矩阵，d
1
, d
2
是 n 维非零向量。如果 d
T
1
Gd
2
= 0，则称向量 d
1
, d
2
是 G 共轭的。类似地，设 d
1
, d
2
, . . . , d
m
是 R
n
中一组非零向量。如果 ∀ij, i = j 有 d
T
i
Gd
j
= 0，
则称 d
1
, d
2
, . . . , d
m
是 G-共轭的。
显然，如果 d
1
, d
2
, . . . , d
m
是 G-共轭的，那么它们是线性无关的。下面给出共轭方向法的算
法流程 (共轭方向法即方向 d
k
共轭)。
step1. 初始化。x
0
∈ R
n
，计算 g
0
= g(x
0
)，给一个 d
0
使 d
T
0
g
0
< 0，令 k := 0。
step2. 计算 α
k
, x
k+1
，min
α∈R
f(x
k
+ αd
k
)。若 ∇f(x
k+1
) = 0 或 k = n − 1，则算法停止；否则转
到 step3。
step3. 计算 d
k+1
，使 d
T
k+1
Gd
j
= 0 , j = 0, 1, . . . , k(共轭方向法)，d
k
∈ R
n
。
step4. 令 k := k + 1 转到 step2。
引理 (扩张子空间定理) 给定严格凸的二次正定函数
f(x) =
1
2
x
T
Gx + b
T
x + c
其中：G = G
T
≻ 0。若 {d
0
, d
1
, . . . , d
n−1
} 是 G-共轭向量，则 ∀x
0
∈ R
n
，共轭方向法至多经过
n 步的精确线搜索后，可以找到 f(x) 的最小点。特别地，对 i = 0, 1, . . . . 迭代点 x
i+1
都是 f (x)
在 x
0
和方向 d
0
, d
1
, . . . , d
i
所张成的线性流形 (仿射子空间)M 中的最小点。
M(x
0
; s
i
) ≡ M(x
0
; {d
0
, d
1
, . . . , d
i
}) =
x
x = x
0
+
i
j=0
λ
j
d
j
, λ
j
∈ R
其中：用 s
i
= span{d
0
, d
1
, . . . , d
i
} 表示基向量 d
0
, d
1
, . . . , d
i
张成的线性子空间，简记作 M
i
。
http://www.ma-xy.com 20 http://www.ma-xy.com
http://www.ma-xy.com
第一章 无约束非线性规划 1.7 搜索方向的确定
共轭梯度的引出
1952，Hestenes 和 Stiefel 在求解线性方程组时，提出共轭梯度法 (线性方程组等价于极小化
一个正定二次函数)。1964 年，Fletcher 和 Reeves 提出了无约束极小问题的共轭梯度法，共轭梯
度法就是使得最速下降方向 g
k
具有共轭性。下面，我们以正定二次函数为例来推导共轭梯度法。
设：
f(x) =
1
2
x
T
Gx + b
T
x + c
f 的梯度为
g(x) = Gx + b
我们令
d
0
= −g
0
则 x
1
= x
0
+ α
0
d
0
。由精确线搜索性质
g
T
1
d
0
= 0
令
d
1
= −g
1
+ β
0
d
0
(1.8)
选择 β
0
，使得
d
T
1
Gd
0
= 0
对 (1.8) 两边同乘以 d
T
0
G，有
β
0
=
g
T
1
Gd
0
d
T
0
Gd
0
=
g
T
1
(g
1
− g
0
)
d
T
0
(g
1
− g
0
)
=
g
T
1
g
1
g
T
0
g
0
由扩张子空间引理，g
T
2
d
2
= 0 , i = 0, 1，利用 d
0
= −g
0
, d
1
= −g
1
+ β
0
d
0
，可知
g
T
2
g
0
= 0 g
T
2
g
1
= 0
又令
d
2
= −g
2
+ β
0
d
0
+ β
1
d
1
选择 β
0
和 β
1
，使得 d
T
2
Gd
2
= 0 , i = 0, 1, . . .，从而有
β
0
= 0
β
1
=
g
T
2
(g
2
− g
1
)
d
T
1
(g
2
− g
1
)
=
g
T
2
g
2
g
T
1
g
2
http://www.ma-xy.com 21 http://www.ma-xy.com
http://www.ma-xy.com
1.7 搜索方向的确定 第一章 无约束非线性规划
一般的，在第 k 次迭代中，令
d
k
= −g
k
+
k−1
i=0
β
i
d
i
(1.9)
选择 β
i
，使 d
T
k
Gd
i
= 0 , i = 0, 1, . . . , k − 1。已假定
g
T
k
d
i
= 0 , g
T
k
g
i
= 0 , i = 0, 1, . . . , k − 1 (1.10)
对 (1.9) 式和 d
T
j
G, j = 0, 1, . . . , k − 1，则
β
j
=
g
T
k
Gd
j
d
T
j
Gd
j
=
g
T
k
(g
j+1
− g
j
)
d
T
j
(g
j+1
− g
j
)
j = 0, 1, . . . , k − 1
由式 (1.10)，有
g
T
k
g
j+1
= 0 j = 0, 1, . . . , k − 2
g
T
k
g
j
= 0 j = 0, 1, . . . , k − 1
故得 β
j
= 0 , j = 0, 1, . . . , k − 2，和
β
k−1
=
g
T
k
(
g
k
−
g
k−1
)
d
T
k−1
(g
k
− g
k−1
)
=
g
T
k
g
k
g
T
k−1
g
k−1
因此，共轭梯度法的公式为
x
k+1
= x
k
+ α
k
d
k
d
k+1
= −g
k+1
+ β
k
d
k
k = 0, 1, . . .
其中：d
0
= −g
0
，α
k
由线搜索得到，β
k
有以下几种计算公式：
β
k
=
g
T
k+1
g
k+1
g
T
k
g
k
[Fletcher-Reeves 公式]
β
k
=
g
T
k+1
(g
k+1
− g
k
)
g
T
k
g
k
[Polak-Ribion-Polyak 公式]
β
k
=
g
T
k+1
(g
k+1
− g
k
)
d
T
k
(g
k+1
− g
k
)
[Growder-Wolfe 公式]
β
k
=
g
T
k+1
G
k+1
d
k
d
T
k
G
k+1
d
k
[Doniel 公式]
β
k
= −
g
T
k+1
g
k+1
d
T
k
g
k
[Dixon 公式]
β
k
=
g
T
k+1
g
k+1
d
T
k
(g
k+1
− g
k
)
[Dai-Yuan 公式]
对上面的 β
k
计算公式，如果 α
k
采用精确线搜索，那么 g
T
k+1
d
k
= 0 。特别地，当 g
k+1
= 0 时，有
g
T
k+1
d
k+1
= −∥g
k+1
∥
2
< 0
http://www.ma-xy.com 22 http://www.ma-xy.com
http://www.ma-xy.com
第一章 无约束非线性规划 1.7 搜索方向的确定
此时，搜索方向 d
k+1
一定是目标函数下降方向。在实际应用中，Fletcher-Reeves 公式较为普通，
对于一些大型的优化问题，Polak-Ribion-Polyak 公式的结果较好。下面讨论共轭梯度法的收敛
性。
由于共轭梯度法的计算公式较多，我们仅给出精确线搜索下 FR 共轭梯度的总体收敛性。设
f : R
n
→ R 在有界水平集 L = {x ∈ R
n
|f(x) ⩽ f(x
0
)} 上连续可微，即 f ∈ C
′
(L)。那么，精确
线搜索下的 FR 共轭梯度法，产生的序列 {x
k
} 至少有一个聚点是驻点，即
(1) 当 {x
R
} 是有穷点列时，最后一个点 x
∗
是 f 的驻点。
(2) 当 {x
R
} 是无穷点列时，它必有极限点，且其任意极限点为 f 的驻点。
共轭梯度法具有二次终止性，即对于二次函数，采用精确一维搜索的共轭梯度法在 n 次迭代
后终止。此外，共轭梯度法是至少线性收敛的，且在适当条件下，共轭梯度法具有 n 步二阶收敛
性。
1.7.6 拟牛顿法
牛顿法具有较快的收敛速度，关键是利用了 Hesse 矩阵提供的曲率信息，但是计算 Hesse 矩
阵困难且其存储量极大。拟牛顿法利用目标函数
f
和一阶导数
g
来构造
Hesse
矩阵的近似。由
此获得一个搜索方向，生成新的迭代点。采用不同的 Hesse 矩阵近似，对应着不同的拟牛顿法。
设 f : R
n
→ R 在开集 D ⊂ R
n
上的二次函数连续可微函数 f ∈ C
2
(D)，假设我们已经知道
了 x
k+1
∈ R
n
的具体值，f 在 x
k+1
附近的二次函数近似为
f(x) ≈ f (x
k+1
) + g
T
k+1
(x − x
k+1
) +
1
2
(x − x
k+1
)
T
G
k+1
(x − x
k+1
)
上式两边对 x 求导，有
g(x) ≈ g
k+1
+ G
k+1
(x − x
k+1
) (1.11)
¬如果 G
k+1
已知，且令 g(x) ≡ g(x
k+2
) ≈ 0，则由上式 (1.11) 可得到牛顿法的迭代公式
x
k+2
= x
k+1
− G
−1
k+1
g
k+1
如果 G
k+1
未知，如何由式 (1.11) 求 x
k+2
呢？
令 x = x
k
, s
k
= x
k+1
− x
k
, y
k
= g
k+1
− g
k
，得
G
−1
k+1
y
k
≈ s
k
显然，对于二次函数 f ，上式是精确成立的。如果我们构造 Hesse 矩阵的 G
k+1
的近似，我们要
求其近似能够满足上述等式，即
H
k+1
y
k
= s
k
其中：H
k+1
是 G
−1
k+1
的近似。或者设 B
k+1
是 G
k+1
的近似，则
B
k+1
s
k
= y
k
http://www.ma-xy.com 23 http://www.ma-xy.com
http://www.ma-xy.com
1.7 搜索方向的确定 第一章 无约束非线性规划
我们把 B
k+1
或者 H
k+1
满足的等式称作拟牛顿条件或者拟牛顿方程。
一般的拟牛顿算法流程如下：
step1. 初始化。x
0
∈ R
n
, H
0
∈ R
n×n
, 0 ⩽ ε < 1, k := 0。
step2. 如果 ∥g
k
∥ ⩽ ε，输出 x
k
；否则计算 d
k
= −H
k
g
k
。
step3. 求 α
k
。计算 arg min
α
f(x
k
+ d
k
α)，令 x
k+1
= x
k
+ α
k
d
k
。
step4. 校正 H
k
产生 H
k+1
。其中 H
k+1
要使得 H
k+1
y
k
= s
k
成立
y
k
= g
k+1
− g
k
s
k
= x
k+1
− x
k
step5.k := k + 1，转 step2.
在上述拟牛顿算法中，H
0
通常取单位矩阵 H
0
= I。由于在每一次迭代中不定矩阵 H
k
总是
不断变化的，故拟牛顿法也称为变尺度方法。
上面介绍了拟牛顿法的思想及拟牛顿法的算法流程。下面给出 H
k
的一些计算公式。
对称秩 1 校正公式
假设 H
k
已知，在构造满足 H
k+1
y
k
= s
k
的矩阵 H
k+1
时，可以令
H
k+1
= H
k
+ E
k
其中：E
k
是一个低秩的校正矩阵。秩 1 校正是指
E
k
= uv
T
即 H
k+1
= H
k
+ uv
T
。由拟牛顿条件，有
H
k+1
y
k
= ( H
k
+ uv
T
)y
k
= s
k
⇒ (v
T
y
k
)u = s
k
− H
k
y
k
故 u 必定是在方向 s
k
−H
k
y
k
上。假设 s
k
−H
k
y
k
= 0 (否则，H
k
已满足拟牛顿条件)。向量 v 满
足 v
T
y = 0，则
H
k
+1
y
k
= H
k
+
1
v
T
y
k
(s
k
− H
k
y
k
)v
T
(1.12)
由于 Hesse 矩阵是对称的，故要求 Hesse 逆近似也是对称的，从而取 v = s
k
− H
k
y
k
，得
H
k+1
y
k
= H
k
+
(s
k
− H
k
y
k
)(s
k
− H
k
y
k
)
T
(s
k
− H
k
y
k
)
T
y
k
(1.13)
式 (1.12) 称为 Broyden 秩一校正公式。特别地，当 v = y
k
时，称为 Broyden 秩一校正公式。
http://www.ma-xy.com 24 http://www.ma-xy.com
http://www.ma-xy.com
第一章 无约束非线性规划 1.7 搜索方向的确定
对称秩 2 校正公式
SR1 校正不能保证 H
k
的正定性，仅当 (s
k
−H
k
y
T
k
)y
k
> 0 时，SR1 校正才具有正定性，而
这个条件往往很难保证。设对称秩二校正为
H
k+1
= H
k
+ auu
T
+ bvv
T
令 H
k+1
满足拟牛顿条件，则
H
k
y
k
+ auu
T
y
k
+ bvv
T
y
k
= s
k
这里 u 和 v 并不唯一确定，但 u 和 v 的明显选择是
u = s
k
v = H
k
y
k
于是
au
T
y
k
= 1 bv
T
y
k
= −1
确定出
a = 1/u
T
y
k
= 1/ s
T
k
y
k
b = −1/v
T
y
k
= −1/y
T
k
H
k
y
k
因此
H
k+1
= H
k
+
s
k
s
T
k
s
T
k
y
k
−
H
k
y
k
y
T
k
H
k
y
T
k
H
k
y
k
(1.14)
(1.14) 式称为 DFP 公式，它是由 Davidon(1959) 提出，后来由 Fletder 和 Powell(1963) 发展的。
DFP 校正能保证 H
k
的正定性，每次迭代需要 3n
2
+ O(n) 次乘法运算，且方法具有超线性收敛
速度。但 DFP 方法具有数值不稳定性，有时产生数值上奇异的 Hessen 矩阵。
BFGS 校正公式
类似于关于 H
k
我们得到的 DFP 修正公式。那样，我们也可以从 B
k
得到 BFGS 修正公式
B
(BF GS)
k+1
= B
k
+
y
k
y
T
k
y
T
k
s
k
−
B
k
s
k
s
T
k
B
k
s
T
k
B
k
s
k
由于 B
k
s
k
= −α
k
g
k
, B
k
d
k
= −g
k
，故上式也可以写为
B
(BF GS)
k+1
= B
k
+
g
k
g
T
k
g
T
k
d
k
−
y
k
y
T
k
αy
T
k
d
k
事实上，只要通过对 DFP 校正公式作简单的变换，H ↔ B, s ↔ y，就可以得到 B
k
的 BFGS 校
正公式。对 B
k
的 BFGS 应用两次逆的秩一校正的 Sherman—Morrison 公式，就可以得到 H
k
的 BFGS 校正公式
H
(BF GS)
k+1
=
I −
s
k
y
T
k
s
T
k
y
k
H
k
I −
y
k
s
T
k
s
T
k
y
k
+
s
k
s
T
k
s
T
k
y
k
BFGS 是非常好的拟牛顿公式，它具有 DFP 校正所具有的性质，并且当采用不精确线搜索时，
BFGS 还具有总体收敛性质。
http://www.ma-xy.com 25 http://www.ma-xy.com
http://www.ma-xy.com
1.7 搜索方向的确定 第一章 无约束非线性规划
1.7.7 模式搜索法：不使用导数的最优化方法
前面介绍了基于导数的最优化方法：最速下降法、牛顿法、修正牛顿法、信赖域方法、共轭
梯度法、拟牛顿方法，这些方法都是当目标函数 f 的导数易求时，才能较好的使用，但如果函数
f 不规则，我们又该如何处理？下面，我们来介绍一些不使用导数的最优化方法：模式搜索法、
Rosenbrock 方法和 Powell 方法。
1961 年，Hooke 和 Jeeves 设计了模式搜索方法，算法从初始基点开始，包括两种类型的移
动。¬探测移动：依次沿着 n 个坐标轴进行，用以确定新的基点和利于函数值下降的方向。模
式移动：沿相邻两个基点连线方向进行，试图顺着“山谷”使函数值更快的减小。
设 x ∈ R
n
，f : R
n
→ R，坐标方向为 e
j
, j = 1, 2, . . . , n。
e
j
= (0 , 0, . . . , 0,
j
1, 0, . . . , 0)
T
给定初始步长 δ 和加速因子 α，使任取初始点 x
1
作为第 1 个基点，x
j
为第 j 个基点，在每轮探
测移动中，自变量 x 用 y
j
表示，即 y
j
是沿着 e
j
探测的出发点，y
1
是沿 e
1
探测的出发点，y
n+1
是沿 e
n
探测得到的点。
首先，从 y
1
= x 出发，进行探测移动。先沿 e
1
探测，如果 f (y
1
+ δe
1
) < f (y
1
) 则探测成功，
令 y
2
= y
1
+ δe
1
，并从 y
2
出发，沿 e
2
进行探测；否则，沿 e
1
方向探测失败，再沿 −e
1
方向探
测。如果 f(y
1
−δe
1
) < f (y
1
) 则探测成功，令 y
2
= y
1
−δe
1
，并从 y
2
出发，沿 e
2
进行探测。如果
f(y
1
−δe
1
) ⩾ f (y
1
) 则沿 −e
1
探测失败，令 y
2
= y
1
，再从 y
2
出发，沿 e
2
进行探测。方法同上，
得到的点记作 y
3
，直到 y
n+1
终止。如果 f(y
n+1
) < f (x
1
) 则 y
n+1
作为新的基点，x
2
:= y
n+1
。
这时，d = x
2
− x
1
是利用函数值减小的方向。
然后，沿方向 x
2
− x
1
进行模式移动，令新的 y
1
为
y
1
= x
2
+ α(x
2
− x
1
)
模式移动之后，以 y
1
为起点进行探测移动，沿坐标轴方向进行，探测完毕后，得到 y
n+1
，如果
f(y
n+1
) ⩾ f (x
2
)，则表示此次模式移动成功，于是取新的基点。
x
∗
3
= y
n+1
再沿方向 x
3
− x
2
进行模式移动。如果 f (y
n+1
) ⩾ f (x
2
) 则表明此次模式移动失败，于是返回到
基点 x
2
，减小步长 δ。再从 x
2
出发，沿坐标轴方向进行探测移动。如此下去，直到 δ < ε 为止。
模式搜索的算法流程如下：
step1. 初始化。x
1
∈ R
n
, e
i
e
j
, δ, α ⩾ 1，缩减率 β ∈ (0, 1)，允许误差 ε > 0，置 y
1
:= x
1
, k :=
1, j := 1。
step2. 如果 f(y
j
+ δe
j
) < f (y
j
)，则令
y
j+1
= y
j
+ δe
j
转到 step4，否则转到 step3.
step3. 如果 f(y
j
− δe
j
) < f (y
j
)，则令
y
j+1
= y
j
− δe
j
http://www.ma-xy.com 26 http://www.ma-xy.com
http://www.ma-xy.com
第一章 无约束非线性规划 1.8 MATLAB 应用实例 2
转到 step4；否则，令 y
j+1
= y
j
转到 step4。
step4. 如果 j < n 则置 j := j + 1，转到 step2，否则转到 step5。
step5. 如果 f(y
n+1
) < f (x
k
) 则置 x
k+1
= y
k+1
令
y
1
= x
k+1
+ α(x
k+1
− x
k
)
置 k := k +1, j = 1，转到 step2；否则，如果 δ ⩽ ε，则停止迭代，得到 x
k
；否则，置 δ := βδ, y
1
=
x
k
, x
k+1
= x
k
, k := k + 1, j = 1，转 step2。
模式移动方向可以看作是最速下降方向的近似，因此模式搜索方法也可以看作是最速下降法
的一种近似，但这种方法的收敛速度是比较慢的，不适合 n 较大的情况。
关于 Rosenbrock 方法、Powell 方法和单纯形等方法，可以直接参考《最优化理论与算法》
陈宝林。下面，我们给出 MATLAB 的求解多元无约束非线性规划问题的示例。
1.8 MATLAB 应用实例 2
MATLAB 中使用 fminunc 和 fminsearch 函数用来求解多维无约束非线性规划问题，fminunc
是利用导数搜索算法，fminsearch 是不使用导数的搜索算法。
(1)fminunc 在没有梯度矩阵输入时，采用拟牛顿法，在有梯度矩阵输入时，采用信赖域算法。
其调用格式为
[x,fval,exitag,output,grad,hession]=fminunc(fun,x0,options)
其中：fun 为目标函数句柄；x0 为初始点；options 为结构体参数/参数结构体；x 为极小点；fval
为极小值；exitag 为返回求解状态；output 为返回求解信息：迭代次数和所用算法等；gval 为
返回 fun 在极小点 x 处的梯度；hessien 为返回 fun 在极小点 x 处的 Hesse 矩阵。
我们用 fminunc 来求解如下优化问题
min f(x) = 100(x
2
− x
2
1
)
2
+ (1 − x
1
)
2
s.t. g =
− 400(x
2
− x
2
1
)x
1
− 2(1 − x
1
)
200(x
2
− x
2
1
)
> 0
1 f u n c t i o n [ f , g]=Afun(x)
2 f =100*(x (2)−x ( 1 ) ^2)^2+(1−x (1) ) ^2;
3 i f nargout>1
4 g=[−400*(x ( 2 )−x(1) ^2)*x (1)−2*(1−x ( 1 ) ) ;
5 200*(x (2)−x (1) ^2) ] ;
6 end
7 x0 = [ −1 ,2] ;
8 fun = @Afun ;
9 opti o ns = optimoptions (@fminunc , ’ Algorithm ’ , ’ quasi−newton ’ , ’ Sp e c i f yobje c t i v eGri d e n t ’ , true ) ;
10 opti o ns . Display = ’ i t e r ’ ;
11 [ x , fva l , e x i t f l a g , output ] = fminunc ( fun , x0 , opti o ns )
12
http://www.ma-xy.com 27 http://www.ma-xy.com
http://www.ma-xy.com
1.8 MATLAB 应用实例 2 第一章 无约束非线性规划
(2)fminsearch 采用单纯形搜索方法进行搜索，其调用格式为
[x,fval,exitag,output] = fminsearch(fun,x0,options)
下面，我们给出 fminsearch 的应用示例：
1 fun = @(x ) 100*( x (2)−x ( 1 ) ^2)^2+(1−x (1) ) ^2;
2 x0 = [ −1. 2 , 1 ];
3 opti o ns = optimset ( ’ PlotFcns ’ , @optimplotfval , ’ Display ’ , ’ i t e r ’ ) ;
4 [ x , fva l , e x i t f l a g , output ] = fminsearch ( fun , x0 , opti o ns )
5
http://www.ma-xy.com 28 http://www.ma-xy.com
