第四部分
数学基础
291
第二十五章 矩阵论
矩阵是高等代数学中的一个常见工具，也常见于统计学、物理学、计算机科学等学科。矩阵运
算是数值分析领域的一个重要问题，将矩阵分解为简单矩阵的组合可以在理论和实际应用中简化矩
阵运算。一些应用广泛而形式特殊的矩阵，如稀疏矩阵和准对角矩阵，都有特定的快速运算方法。
假设矩阵A是m × n维的实数矩阵，则A ∈ R
m×n
，并记A
T
为A的转置矩阵（Transpose Ma-
trix）。 如 果 m = n，则称A为方块矩阵，简称方阵，并记|A| 或det(A)为其行列式（Determinant）。
如果方阵A ∈ R
n×n
对称，则有A
T
= A，比如非对角元都是零元的对角阵（Diagonal Ma-
trix）。 如 果 对 角 阵 对 角 元 都 是 1，则称其为单位矩阵（Identity Matrix）， 并 记 作 I。如果A可逆
（Invertible）， 记 A
−1
为其逆矩阵（Inverse Matrix）， 则 有 AA
−1
= A
−1
A = I，此时A也被称作非
奇异阵（Nonsingular Matrix）， 并 且 满 足 |A| ̸=0。方阵A 所有对角元素之和称作矩阵A 的迹
（Trace），记作tr(A)=
(
i
a
ii
。
定义25.1 (分块矩阵). 如果矩阵A ∈ R
m×n
，利用矩阵行或列向量，可分块表示成如下两种形式：
⎧
⎨
⎩
A =[A
•1
,A
•2
,...,A
•n
],
A =[A
1•
,A
2•
,...,A
m•
]
T
,
其中，A
•i
=[a
1i
,...,a
mi
]
T
∈ R
m
为矩阵A的第i列，A
i•
=[a
i1
,...,a
in
] ∈ R
1×n
是A的第i行。
25.1 特征系统
25.1.1. 特征值与特征向量
定义25.2. 对于矩阵A ∈ C
n×n
，如果存在λ ∈ C与非零向量ω ∈ C
n
，满足Aω = λω，则称λ是矩阵A
的一个特特特征征征值值值（Eigenvalue）， ω是其对应特特特征征征向向向量量量（Eigenvector）。 矩 阵 A对应于复空间C
n×n
上的
一个线性变换。
定义25.3. 如果矩阵A ∈ C
n×n
存在特征值λ
1
,...,λ
n
，则矩阵A的谱谱谱半半半径径径（Spectral Radius）有定
义：
ρ(A) = max
1≤i≤n
|λ
i
|. (25.1)
293
搜索与排名 Searching and Ranking
!"#
如果特征值满足关系|λ
1
| > |λ
2
| ≥ ··· ≥ |λ
n
|，则称λ
1
为A的主主主特特特征征征值值值，其对应特征向量ω
1
为A的主主主
特特特征征征向向向量量量，矩阵A的谱半径ρ(A)=|λ
1
|。
命题25.1. 如果矩阵A存在特征值λ
1
,...,λ
n
，则特征值满足：
"
i
λ
i
= tr(A), (25.2)
!
i
λ
i
= |A|. (25.3)
如果特征值彼此不相同，则对应特征向量ω
1
,...,ω
n
也线性无关。
定义25.4 (Gershgorin圆盘). 对于矩阵A ∈ C
n×n
，如果记r
i
=
(
j̸=i
|a
ij
|，1 ≤ i ≤ n，我们称
D(a
ii
; r
i
)=
?
z ∈ C
6
6
6
6
|z − a
ii
| ≤ r
i
@
(25.4)
是矩阵A的Gershgorin圆圆圆盘盘盘（Gershgorin Disc）。
定理25.1 (Gershgorin圆盘定理). 复数域上矩阵的任意特征值都至少落在一个Gershgorin圆盘内。
定义25.5 (矩阵范数). 如果矩阵A ∈ R
m×n
，则它存在如下几种形式的范数∥A∥：
• ∥A∥
1
= max
1≤j≤n
m
(
i=1
|a
ij
|表示矩阵的最大列和（绝对值），称作列列列和和和范范范数数数。
• ∥A∥
2
= σ
1
(A)=
)
ρ(A
T
A)，表示矩阵最大的奇异值，称作谱谱谱范范范数数数。
• ∥A∥
∞
= max
1≤i≤m
n
(
j=1
|a
ij
|，表示矩阵的最大行和（绝对值），称作行行行和和和范范范数数数。
• ∥A∥
F
=
)
tr(A
T
A)=
Q
m
(
i=1
n
(
j=1
|a
ij
|
2
，称作Frobenius范范范数数数。
• ∥A∥
⋆
=
(
i
σ
i
(A)，表示矩阵奇异值之和，称作核核核范范范数数数。
定理25.2. 任意复数域上的矩阵，其谱半径不大于其任意一种诱导范数。
证明： 假设矩阵A ∈ C
n×n
的特征值为λ
i
，1 ≤ i ≤ n。任选一个特征值λ ∈ {λ
1
,...,λ
n
}，设ω ̸=0是
对应特征向量。我们构造矩阵X =[ω, 0,...,0] ∈ C
n×n
，则∥X∦=0且AX = λX。根据矩阵范数的
一致性（也称次可乘性）
|λ|∥X∥ = ∥AX∥≤∥A∥∥X∥
可知|λ∥≤∥A∥。由λ的任意性证得ρ(A) ≤∥A∥，即A的谱半径是其任意一种诱导范数的下界。
25.1.2. 二次型
定义25.6. 多元变量x
1
,...,x
n
的二二二次次次型型型（Quadratic Form）是形如
Q(x)=a
11
x
2
1
+ ···+ a
ij
x
i
x
j
+ ···+ a
nn
x
2
n
=
"
1≤i≤n
"
1≤j≤n
a
ij
x
i
x
j
= x
T
Ax (25.5)
的n元二次函数，其中A =(a
ij
)
n×n
是一个二次型矩阵。
$%"&'#(
294
)!*+",$
25.1. 特征系统
!"#
二次型矩阵A存在多种表示形式，比如Q(x)=2x
2
1
+ x
1
x
2
+5x
2
x
1
− 4x
2
2
等价于
Q(x)=2x
2
1
+2x
1
x
2
+4x
2
x
1
− 4x
2
2
,
Q(x)=2x
2
1
+8x
1
x
2
− 2x
2
x
1
− 4x
2
2
,
.
.
.
.
.
.
.
.
.
Q(x)=2x
2
1
+3x
1
x
2
+3x
2
x
1
− 4x
2
2
相应地二次型矩阵也会发生变化。假设原始二次型矩阵是A，由于Q(x
1
,...,x
n
)=x
T
(A + A
T
)x/2，
则(A + A
T
)/2是对称矩阵，即每个二次型都存在唯一的对称型二次型矩阵，并作为基本二次型矩
阵，如前例中的基本二次型矩阵为
A =
⎡
⎣
23
3 −4
⎤
⎦
.
定义25.7. 一个对称矩阵A ∈ R
n×n
是正正正（（（负负负）））定定定矩矩矩阵阵阵（Positive or Negative Definite Matrix）， 当 且
仅当对所有非零向量x ∈ R
n
都有x
T
Ax > 0(< 0)。 A 是半半半正正正（（（负负负）））定定定矩矩矩阵阵阵（Positive or Negative
Semi-definite Matrix）当且仅当对所有向量x ∈ R
n
都有x
T
Ax ≥ 0(≤ 0)。
命题25.2. 假设矩阵A的主特征向量是ω
1
，最小特征值对应的特征向量是ω
n
，则二次型Q(x)=
x
T
Ax在ω
1
方向上变化速率最大，在ω
n
方向上变化速率最小。
25.1.3. 严格对角占优矩阵
定义25.8 (严格对角占优矩阵). 对于矩阵A ∈ C
n×n
，如果对任意的1 ≤ i ≤ n都有
|a
ii
| >
"
j̸=i
|a
ij
|, (25.6)
则称A是严严严格格格对对对角角角占占占优优优矩矩矩阵阵阵（Strictly Diagonally Dominant Matrix）。
定理25.3. 严格对角占优矩阵是非奇异矩阵。
25.1.4. 随机矩阵
随机矩阵（Stochastic Matrix）， 也 称 概率矩阵（Probability Matrix）、 转移矩阵（Transition
Matrix）或者马尔科夫矩阵 （Markov Matrix）是一种用以述马尔科夫链状态转移的一种矩阵，
主要应用于概率论、统计学、线性代数、金融数学及计算机领域。
定义25.9 (随机矩阵). 对于非负矩阵A ∈ R
n×n
，如果A的行和等于1，则称它为行随机矩阵或右随机
矩阵，矩阵的每一行对应一个概率分布；如果矩阵A的列和等于1，则称其列随机矩阵或左随机矩
阵，矩阵的每一列对应一个概率分布；如果A的行和与列和都等于1，则称其为双随机矩阵，矩阵
的每一行、每一列都对应一个概率分布。
定理25.4. 随机矩阵的主特征值是1。
$%"&'#(
295
)!*+",$
搜索与排名 Searching and Ranking
!"#
证明： 假设矩阵A ∈ C
n×n
是一个行随机矩阵，根据定理25.2可知
ρ(A) ≤∥A∥
1
=1.
由于A1 =(1,...,1)
T
= 1，表明矩阵A存在一个特征向量λ =1，对应特征向量是1。显然，特征
值λ =1是矩阵A的主特征值。列随机矩阵与双随机矩阵类似可证。
25.1.5. 特征系统的解
Algorithm 19 幂法
Input: 给定迭代次数T 或误差阈值ϵ，t =1，任意非零向量x
0
∈ R
n
while t<T 或∥x
t
− x
t−1
∥ > ϵ do
y
t
= Ax
t−1
x
t
= y
t
/∥y
t
∥
t = t +1
end while
Output: 主特征向量ω
1
= x
T
定理25.5. 如果A ∈ R
n×n
是一个正定矩阵，则幂法收敛于A的主特征向量ω
1
。
证明： 如果A是正定矩阵，则n个特征向量一定线性无关，并且可以使用它们线性表示任意非零向
量x
0
，即存在非零向量a =(a
1
,...,a
n
)
T
，满足
x
0
=
"
i
a
i
ω
i
.
由幂法迭代步骤可知
A
t
x
0
=
"
i
a
i
A
t
ω
i
=
"
i
a
i
λ
t
i
ω
i
.
对等式左右两端同除以λ
t
1
，则有
(A
t
x
0
)/λ
t
1
= a
1
ω
1
+
"
2≤i≤n
#
a
i
(λ
i
/λ
1
)
t
ω
i
$
. (25.7)
由于λ
1
是主特征值，对任意i ≥ 2，当t →∞时都有(λ
i
/λ
1
)
t
→ 0，则(A
t
x
0
)/λ
t
1
→ a
1
ω
1
，表明幂法
收敛且收敛于A的主特征向量。
注解25.1. 幂法是Google核心搜索算法PageRank使用的一种迭代方法，由于Google网页链接矩阵是
一个随机矩阵，其收敛性可以保证。幂法迭代过程矩阵-向量乘积运算是主要的计算开销，其计算
复杂度为O(n
2
)。在实际应用时，由于海量的网页数目，计算复杂度不容小觑。如果A是稀疏矩阵，
矩阵-向量的乘积运算执行效率非常高，幂法仍然是一个不错的选择。此外，随机矩阵主特征值与
第二大特征值之差，也称谱谱谱隙隙隙（Eigengap or Spectral Gap），会影响到幂法的收敛速度。一般地，谱
隙越大幂法的收敛速度越快。
$%"&'#(
296
)!*+",$
25.2. 矩阵微分
!"#
25.2 矩阵微分
矩阵微分（Matrix Differentiation）是一个实用的数学工具，人们在机器学习理论证明、优化
算法推导时都可看到其身影。本节对各种形式的矩阵微分定义进行总结并供详细推导 [298]。
25.2.1. 函数y = Ax的偏导
定义25.10. 假设函数ψ : R
n
-→ R
m
，x ∈ R
n
，则y = ψ(x) ∈ R
m
关于x的一阶偏导为
∂ψ(x)
∂x
= J(x, ψ)=
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
∂y
1
/∂x
1
∂y
1
/∂x
2
··· ∂y
1
/∂x
n
∂y
2
/∂x
1
∂y
2
/∂x
2
··· ∂y
2
/∂x
n
.
.
.
.
.
.
.
.
.
.
.
.
∂y
m
/∂x
1
∂y
m
/∂x
2
··· ∂y
m
/∂x
n
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
(25.8)
称作是变换ψ(·)的雅雅雅克克克比比比矩矩矩阵阵阵（Jacobian Matrix）
J
ij
=
T
∂y
i
∂x
j
U
, ∀1 ≤ i ≤ m, 1 ≤ j ≤ n. (25.9)
如果y是一个标量，即m =1，则雅克比矩阵J(x, ψ)=
#
∂y/∂x
1
, ∂y/∂x
2
, ··· , ∂y/∂x
n
$
∈ R
1×n
；如
果x是一个标量，即n =1，则雅克比矩阵J(x, ψ)=
#
∂y
1
/∂x, ∂y
2
/∂x, ··· , ∂y
m
/∂x
$
T
∈ R
m
。
命题25.3. 假设y = Ax，其中A ∈ R
m×n
，x ∈ R
n
，并且A与x无关，则有
∂y/∂x = ∂(Ax)/∂x = A. (25.10)
证明： 由于y
i
=
(
k
a
ik
x
k
，则有
∂y
i
/∂x
j
= a
ij
,
对任意的1 ≤ i ≤ m和1 ≤ j ≤ n都成立，显然∂y/∂x = A。
命题25.4. 假设y = Ax，其中A ∈ R
m×n
，x ∈ R
n
，并且A与x无关。如果x还是向量z的函数，
且A与z无关，则有
∂y/∂z = ∂(Ax)/∂z = A(∂x/∂z). (25.11)
证明： 由于y
i
=
(
k
a
ik
x
k
，则有
∂y
i
/∂z
j
=
"
k
#
a
ik
(∂x
k
/∂z
j
)
$
,
等式右侧是A(∂x/∂z)的第(i, j)个元素，则有∂y/∂z =(∂y/∂x)(∂x/∂z)=A(∂x/∂z)。
$%"&'#(
297
)!*+",$
搜索与排名 Searching and Ranking
!"#
25.2.2. 函数c = y
T
Ax的偏导
命题25.5. 假设A ∈ R
m×n
，x ∈ R
n
，y ∈ R
m
，并且A和y都与x无关，则有
∂(y
T
Ax)/∂x = y
T
A, (25.12)
∂(y
T
Ax)/∂y =(Ax)
T
. (25.13)
特别地，y
T
Ax对y
T
的一阶偏导是y
T
Ax关于y一阶偏导的转置，即有
∂(y
T
Ax)/∂y
T
=
#
∂(y
T
Ax)/∂y
$
T
= Ax. (25.14)
证明： 记z
T
= y
T
A，则∂(y
T
Ax)/∂x = ∂(z
T
x)/∂x = z
T
，即第一个等式得证。由于y
T
Ax =
(y
T
Ax)
T
=(Ax)
T
y，利用第一个等式可以证明∂(y
T
Ax)/∂y = ∂[(Ax)
T
y]/∂y =(Ax)
T
。
推论25.1. 假设A ∈ R
n×n
，x ∈ R
n
，并且A与x无关，则有
∂(x
T
Ax)/∂x = x
T
(A + A
T
), (25.15)
∂
2
(x
T
Ax)/∂x∂x
T
= A + A
T
. (25.16)
如果A是对称矩阵A = A
T
，则有
∂(x
T
Ax)/∂x =2x
T
A, (25.17)
∂
2
(x
T
Ax)/∂x∂x
T
=2A. (25.18)
证明： 根据二次型的定义x
T
Ax =
(
i
(
j
a
ij
x
i
x
j
，可知
∂(x
T
Ax)/∂x
k
=
"
j
a
kj
x
j
+
"
i
a
ik
x
i
= x
T
(A
T
k•
+ A
•k
)
对任意1 ≤ k ≤ n都成立，显然∂(x
T
Ax)/∂x = x
T
(A + A
T
)。此外，
∂
2
(x
T
Ax)/(∂x∂x
T
)=∂[x
T
(A + A
T
)]/∂x
T
= ∂[(A + A
T
)x]/∂x = A + A
T
.
命题25.6. 假设x ∈ R
n
，y ∈ R
n
，并且x与y都是关于向量z的函数，则有
∂(x
T
y)/∂z = x
T
(∂y/∂z)+y
T
(∂x/∂z). (25.19)
证明： 根据x
T
y =
(
i
(x
i
y
i
)可知
∂(x
T
y)/∂z
k
=
"
i
∂(x
i
y
i
)/∂z
k
=
"
i
[x
i
(∂y
i
/∂z
k
)+y
i
(∂x
i
/∂z
k
)] = x
T
(∂y/∂z
k
)+y
T
(∂x/∂z
k
),
对任意k ≥ 1都成立，从而可得
∂(x
T
y)/∂z =[∂(x
T
y)/∂x](∂x/∂z)+[∂(x
T
y)/∂y](∂y/∂z)=y
T
(∂x/∂z)+x
T
(∂y/∂z). (25.20)
如果y = x，则有∂(x
T
x)/∂z = ∂∥x∥
2
2
/∂z =2x
T
(∂x/∂z)。
$%"&'#(
298
)!*+",$
25.2. 矩阵微分
!"#
命题25.7. 假设x ∈ R
n
，y ∈ R
m
，A ∈ R
m×n
，且x与y都是关于向量z的函数，但A与z无关，则有
∂(y
T
Ax)/∂z = y
T
A(∂x/∂z)+(Ax)
T
(∂y/∂z). (25.21)
当y = x时，则有
∂(x
T
Ax)/∂z = x
T
(A + A
T
)(∂x/∂z). (25.22)
如果A是对称矩阵，从而可以得到
∂(x
T
Ax)/∂z =2x
T
A(∂x/∂z). (25.23)
证明： 由链式法则可得
∂(y
T
Ax)/∂z =[∂(y
T
Ax)/∂x](∂x/∂z)+[∂(y
T
Ax)/∂y](∂y/∂z )=y
T
A(∂x/∂z)+(Ax)
T
(∂y/∂z),
证毕。
定义25.11. 假设ψ : R
n
-→ R，x ∈ R
n
，则c = ψ(x)关于x的二阶导数
∂
2
c
∂x
2
= H(x, ψ)=
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
∂
2
c/(∂x
1
∂x
1
) ∂
2
c/(∂x
1
∂x
2
) ··· ∂
2
c/(∂x
1
∂x
n
)
∂
2
c/(∂x
2
∂x
1
) ∂
2
c/(∂x
2
∂x
2
) ··· ∂
2
c/(∂x
2
∂x
n
)
.
.
.
.
.
.
.
.
.
.
.
.
∂
2
c/(∂x
n
∂x
1
) ∂
2
c/(∂x
n
∂x
2
) ··· ∂
2
c/(∂x
n
∂x
n
)
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
(25.24)
称作是变换ψ(·)的黑黑黑塞塞塞矩矩矩阵阵阵（Hessian Matrix）：
H
ij
=
T
∂
2
c
∂x
i
∂x
j
U
, ∀1 ≤ i ≤ m, 1 ≤ j ≤ n. (25.25)
定义25.12. 假设ψ : R
m×n
-→ R，X ∈ R
m×n
，则c = ψ(X)关于X的一阶偏导
∂c
∂X
=
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
∂c/∂x
11
∂c/∂x
12
··· ∂c/∂x
1n
∂c/∂x
21
∂c/∂x
22
··· ∂c/∂x
2n
.
.
.
.
.
.
.
.
.
.
.
.
∂c/∂x
m1
∂c/∂x
m2
··· ∂c/∂x
mn
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
(25.26)
是c对矩阵逐元素偏导：
∂c/∂X =
%
∂c/∂x
ij
&
m×n
. (25.27)
命题25.8. 假设α ∈ R
m
，β ∈ R
n
，则有
∂(α
T
Xβ)/∂X = αβ
T
, (25.28)
∂(α
T
Xα)/∂X = αα
T
. (25.29)
$%"&'#(
299
)!*+",$
搜索与排名 Searching and Ranking
!"#
25.2.3. 矩阵迹的偏导
如果矩阵A，B，X都是R
n×n
上的方阵，则有
∂ tr(X)/∂X = I, (25.30)
∂ tr(AX)/∂X = A
T
, (25.31)
∂ tr(AX
T
)/∂X = A, (25.32)
∂ tr(AX
T
B)/∂X = BA, (25.33)
∂ tr(AXB)/∂X =(BA)
T
, (25.34)
∂ tr(AX
−1
)/∂X = −(X
−1
X
−1
)
T
, (25.35)
∂ tr(AX
−1
B)/∂X = −(X
−1
BAX
−1
)
T
, (25.36)
∂ tr(X
T
AX)/∂X =(A + A
T
)X, (25.37)
∂ tr(XAX
T
)/∂X = X(A + A
T
), (25.38)
∂ tr(AXBX)/∂X = A
T
X
T
B
T
+ B
T
X
T
A
T
, (25.39)
∂ tr(AXBX
T
)/∂X = AXB + A
T
XB
T
, (25.40)
∂ tr(AXX
T
B)/∂X =(A
T
B
T
+ BA)X. (25.41)
25.2.4. 行列式的偏导
如果矩阵A，B，X都是R
n×n
上的方阵，则有
∂|X|/∂X = |X|(X
T
)
−1
, (25.42)
∂|AXB|/∂X = |AXB||X|(X
T
)
−1
, |A| ̸=0, |X| ̸=0, | B| ̸=0, (25.43)
∂ log |X|/∂X =(X
T
)
−1
, |X| > 0, (25.44)
∂ log |X
T
AX|/∂X =(A + A
T
)X(X
T
AX)
−1
, |X
T
AX| > 0, (25.45)
∂ log |X
T
X|/∂X =2X(X
T
X)
−1
, |X
T
X| > 0. (25.46)
25.2.5. 逆矩阵的偏导
定义25.13. 假设A ∈ R
m×n
的每个元素都是一个标量c的函数，则A关于c的导数为
∂A
∂c
=
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
∂a
11
/∂c ∂a
12
/∂c ··· ∂a
1n
/∂c
∂a
21
/∂c ∂a
22
/∂c ··· ∂a
2n
/∂c
.
.
.
.
.
.
.
.
.
.
.
.
∂a
m1
/∂c ∂a
m2
/∂c ··· ∂a
mn
/∂c
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
=
T
∂a
ij
∂c
U
m×n
. (25.47)
$%"&'#(
300
)!*+",$
25.2. 矩阵微分
!"#
命题25.9. 假设A ∈ R
n×n
是一个非奇异矩阵，其各个元素都是标量c的函数，则有
∂A
−1
/∂c = −A
−1
(∂A/∂c)A
−1
. (25.48)
证明： 由于A是非奇异矩阵，则存在逆矩阵A
−1
，满足AA
−1
= I，两边同时取关于c的导数，则有
(∂A/∂c)A
−1
+ A(∂A
−1
/∂c)=0,
整理可得∂A
−1
/∂c = −A
−1
(∂A/∂c)A
−1
。
$%"&'#(
301
)!*+",$
第二十六章 优化理论
优化理论与算法是数学的一个重要分支，是一门应用广泛、实用性强的学科，它所研究的问题
是讨论在众多的方案中什么样的方案最优以及怎样找出最优方案。
最优化在航空航天、生命科学、水利科学、地球科学、工程技术等自然科学领域和经济金融等
社会科学领域有着广泛和重要的应用，其研究和发展一直受到广泛关注。最优化研究包含理论、方
法和应用：最优化理论主要研究解的最优性条件、灵敏度分析、解的存在性和一般复杂性等，而最
优化方法研究包括构造新算法、证明解的收敛性、算法的比较和复杂性等，最优化的应用研究则包
括算法的实现、算法的程序、软件包及商业化、在实际问题的应用。
现代优化理论中最重要的未解难题是发现通用的全局最优化条件。如果没有全局最优化条件，
我们无从知道何时找到最优解，现有解是否最优，从而无法有效地组织优化过程或及时中断搜索。
任何全局最优化条件都具有理论意义和实用价值。
)!*+",$%"&'#(
26.1 线性规划
线性规划（Linear Programming）是优化理论的重要组成部分，在各个领域都有广泛应
用，比 如 在 生 产 领 域，通 过 合 理 分 配 有 限 的 资 源， 以 最 小 的 成 本 获 得 最 高 的 利 润。 一 般 认
为，线 性 规划 肇 始 于 斯坦 福 大 学 教授George Dantzig，实际上早在1826年左右，法国大数学
家Jean Fourier已开始钻研求解线性不等式组，并出Fourier-Motzkin消元法。1939年，荷兰数
学家Tjalling Koopmans与苏联数学家Leonid Kantorovich联合发展出一种新的数学理论，称之
为Linear Programming，两人由于在资源最优配置理论方面的贡献荣膺诺贝尔经济学奖（1975年）。
同一年，George Stigler[299]在研究美军士兵食品的最佳营养搭配时，出一种启发式算法求解线
性规划问题，并荣获1982 年诺贝尔经济学奖。1947年，George Dantzig受到Stigler的启发，出了
单纯形法（Simplex Method）[300]，线性规划真正步入公众的视野。1979年，苏联数学家Leonid
Khachiyan证明了线性规划问题可以在多项式时间内求解，并出了基于收缩球体的非线性几何
理论的椭球法（Ellipsoid Algorithm）[301]，一时轰动全球，并掀起了新一轮线性规划的研究高
303
搜索与排名 Searching and Ranking
!"#
潮。遗憾的是广泛的数值试验表明，椭球算法的计算性能远逊于单纯形法。在1984年，印度裔科学
家Narendra Karmarkar引入非线性规划中的牛顿投影发明了另外一种多项式时间的求解方法，被
称作Karmarkar方法[302]。 Karmarkar 方法无论是理论还是数值计算上都优于椭球法，并极大地
激发了人们的研究热情，自从Karmarkar方法诞生，关于内点法的研究成果不断推陈出新，至今已
经有2000多份研究结果、上百个改进算法面世。实际上，无论是Ellipsoid 算法还是Karmarkar方
法都属于内点法（Interior Point Method），与单纯形法截然不同。
单纯刑法与内点法孰优孰劣，目前仍然难以分晓。单纯形法沿着边界由一个顶点移动到相邻的
顶点，内点法却是穿过可行域内部逼近最优解。由于每一个顶点的相邻顶点数目有限，因此前一步
移动时的计算量不大。反观一个内点却有无穷多的相邻内点，所以内点法需要更加周详地考虑每
一步移动的方向与步长。综合而言，对于中小型问题，由于可行域顶点不多，单纯形法的表现大
体上优于内点法。对于大型或超大型问题而言，内点法的优势比较明显。线性规划理论日趋成熟，
已经成为现代科学的一个重要工具。1961年Charnes和Cooper[303]继Dantzig之后又出目标规划
（Goal Programming）的概念，推动优化理论研究进一步发展。
26.1.1. 基本概念
线性规划是一类特殊的规划问题，目标函数和约束条件都是线性的，约束为等式或者不等式。
满足不等式线性约束的点集称为半平面（half-plane），规划的约束集合是所有半平面的交集。简单
的线性规划问题可以使用作图法求解。通常，最值点可以在约束集的边（entire edge）或者角点
（corner point）取得。
给定一个m维列向量b =(b
1
,b
2
,...,b
m
)
T
，n维列向量c =(c
1
,c
2
,...,c
n
)
T
，一个m × n的矩
阵A。标准的最大化问题即寻找一个n维向量x =(x
1
,x
2
,...,x
n
)
T
以使目标函数c
T
x最大，约束条
件可以表示为Ax ≤ b，x ≥ 0。这个问题可以表示成如下形式的规划问题：
max z = c
T
x
s.t. Ax ≤ b
x ≥ 0
(26.1)
类似地，可以定义标准的最小化问题：
min z = y
T
b
s.t. y
T
A ≥ c
T
y ≥ 0
(26.2)
在标准的最大化问题(26.1)中的x和标准的最小化问题(26.2)中的y，如果分别满足两个问题的约
束，则称其为可行解（feasible），由所有可行解构成的集合称为约束集（constrain set）。如果一个
规划问题的约束集非空，则称该问题是可行的；否则，则不可行（infeasible）。
对于可行的最值问题，如果目标函数可以在约束集上取得任意大（小）的值，则称其为无界的
（unbounded）， 否 则 就 是 有 界 的 （ bounded）。 因 此 ， 对 于 任 意 一 个 线 性 规 划 问 题 ， 存 在 三 种 可 能
$%"&'#(
304
)!*+",$
26.1. 线性规划
!"#
情形：有界可行，无界可行，不可行（约束集为空）。 使得目标函数取得最值的可行解又称为最优
解（optimal）。
26.1.2. 标准型
所有的线性规划问题都可以通过一定的变换，转化为标准形式。常见的有如下几种情形：
（1） 目标函数最小化和最大化相互转化；
（2） 消去等式约束，将可以显式表达的变量代入到目标函数和其他约束条件，可以消去至少一个
变量和一个约束；
（3） 变量没有非负约束，可以使用两个非负约束的差替换。
26.1.3. 对偶问题
任意线性规划问题都存在一个对偶问题（Dual Problem）与之对应。对偶问题可以做如下定
义，假设原始线性规划问题是形如(26.1)的最大化问题，那么其对偶问题为形如(26.2)的最小化问
题，标准最大化问题都有一个标准最小化的对偶问题与之对应，二者互为对偶问题。
根据标准最大化问题(26.1)，我们可以改造乘法形式的CCR模型(28.3)如下：
max (0,...,0,y
k1
,...,y
km
)
V WX Y
c
(ν,µ)
T
s.t.
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
−x
11
−x
12
··· −x
1s
y
11
y
12
··· y
1m
−x
21
−x
22
··· −x
2s
y
21
y
22
··· y
2m
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
−x
n1
−x
n2
··· −x
ns
y
n1
y
n2
··· y
nm
x
k1
x
k2
··· x
ks
00··· 0
−x
k1
−x
k2
··· −x
ks
00··· 0
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
V WX Y
A
(ν,µ)
T
≤
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
0
0
.
.
.
0
1
−1
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
VWXY
b
(ν,µ) ≥ 0
(26.3)
根据标准最大化模型的对偶形(26.2)，可得下式成立：
min (λ
1
,...,λ
n
, θ
1
, θ
2
)b
s.t. (λ
1
,...,λ
n
, θ
1
, θ
2
)A ≥ c
λ
i
≥ 0,i=1,...,n,θ
1
≥ 0, θ
2
≥ 0
(26.4)
$%"&'#(
305
)!*+",$
搜索与排名 Searching and Ranking
!"#
将其展开，化简可得：
min (θ
1
− θ
2
)
s.t. (θ
1
− θ
2
)x
k
−
n
(
i=1
λ
i
x
i
≥ 0
−y
k
+
n
(
i=1
λ
i
y
i
≥ 0
λ
i
≥ 0,i=1,...,n,θ
2
≥ 0, θ
2
≥ 0
(26.5)
可重新写作下面形式的模型：
min θ
s.t. θx
k
−
n
(
i=1
λ
i
x
i
≥ 0
−y
k
+
n
(
i=1
λ
i
y
i
≥ 0
λ
i
≥ 0,i=1,...,n
(26.6)
标准最值问题与其对偶问题的关系可由以下几个定理及推论予以说明。
定理26.1. 如果x是标准最大化问题(26.1)的可行解，y是其对偶问题(26.2)的可行解，则c
T
x ≤ y
T
b
（c
T
≤ y
T
A，Ax ≤ b，则c
T
x ≤ y
T
，Ax ≤ y
T
b）。
定理26.2. 如果标准最值问题与其对偶问题都是可行的，那么一定都是有界可行的。
证明： 假设y是最小化问题(26.2)的某个可行解，x是其对偶问题的任意一个可行解，那么根据定
理26.1可知，对于其对偶问题，y
T
b 都是c
T
x的上界，说明c
T
x是有界的；反之，可以证明假设x是
其对偶问题的某个可行解，y 是最小化问题的任意可行解，那么对于最小化问题，c
T
x 都是y
T
b的
下界，说明y
T
b 也是有界的。
定理26.3. 如果问题(26.1)和(26.2)分别存在可行解x
∗
、y
∗
，满足c
T
x
∗
= y
∗T
b，则称x
∗
、y
∗
分别为问
题(26.1)和(26.2)的最优解。
定理26.4 (对偶定理). 如果一个标准线性规划问题是有界可行的，则其对偶问题也是有界可行的，
二者的最优目标函数值是相同的，并且两个问题都存在最优解。
定理26.5 (均衡定理). 假设x
∗
、y
∗
分别是问题(26.1)、(26.2)的可行解，则x
∗
、y
∗
是最优解的充要条
件为：
1. 对于任意i，若满足
n
(
j=1
a
ij
x
∗
j
<b
i
，都有y
i
=0；
2. 对于任意j，若满足
m
(
i=1
y
∗
i
a
ij
>c
j
，都有x
j
=0。
$%"&'#(
306
)!*+",$
26.1. 线性规划
!"#
有时，这两个条件称作互互互补补补松松松弛弛弛条条条件件件（（（Complementary Slackness Condition））） 。
证明： 充分性：已知对于i，只有当
n
(
j=1
a
ij
x
∗
j
= b
i
，才有y
i
̸=0，那么
m
"
i=1
y
∗
i
b
i
=
m
"
i=1
y
∗
i
n
"
j=1
a
ij
x
∗
j
=
m
"
i=1
n
"
j=1
a
ij
x
∗
j
y
∗
i
(26.7)
同理，对于j，只有当
m
(
i=1
y
∗
i
a
ij
= c
j
，才有x
j
̸=0，那么
n
"
j=1
x
∗
j
c
j
=
n
"
j=1
x
∗
j
m
"
i=1
a
ij
y
∗
i
=
m
"
i=1
n
"
j=1
a
ij
x
∗
j
y
∗
i
(26.8)
也就是说
m
"
i=1
y
∗
i
b
i
=
n
"
j=1
x
∗
j
c
j
(26.9)
根据定理26.3，x
∗
、y
∗
分别为问题(26.1)和(26.2)的最优解。
必要性：若x
∗
、y
∗
分别为问题(26.1)和(26.2)的最优解，则根据定理26.1，
n
"
j=1
x
∗
j
c
j
≤
m
"
i=1
n
"
j=1
a
ij
x
∗
j
y
∗
i
≤
m
"
i=1
y
∗
i
b
i
由于
n
(
j=1
x
∗
j
c
j
与
m
(
i=1
y
∗
i
b
i
都是有界的，根据定理26.4，上式左侧等于右侧，从而三项相等，对于前两
项有
n
"
j=1
(c
j
−
m
"
i=1
y
∗
i
a
ij
)x
∗
j
=0
由于x
∗
j
≥ 0，c
j
≤
m
(
i=1
y
∗
i
a
ij
，从而当c
j
<
m
(
i=1
y
∗
i
a
ij
时，x
j
=0；当
n
(
j=1
a
ij
x
∗
j
<b
i
时，y
i
=0。
26.1.4. 基
考虑线性规划问题
max z =
n
(
i=1
c
i
x
i
s.t.
n
(
k=1
a
ik
x
k
= b
i
,i=1,...,m
x
i
≥ 0,i=1,...,n
(26.10)
假设A ∈ R
m×n
是约束条件的系数矩阵，秩为m。假设A经过初等变换，分解为两个部分
A =(B | N)
其中，B ∈ R
m×m
是A的一个非奇异子阵，称为线性规划的一个“基”（Basis），列向量线性无关。
$%"&'#(
307
)!*+",$
搜索与排名 Searching and Ranking
!"#
不失一般性，设矩阵A的前m列就是一个“基”：
B =
⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎣
a
11
a
12
··· a
1m
a
21
a
22
··· a
2m
.
.
.
.
.
.
.
.
.
.
.
.
a
m1
a
m2
··· a
mm
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎦
=(P
1
,P
2
, ··· ,P
m
) (26.11)
矩阵的列P
i
,i =1,...,m称作是基向量，则剩余列构成非基向量N =(P
m+1
,...,P
n
)。基向量对应
的变量x
i
,i=1,...,m称为基变量，否则，称为非基变量。
假设规划问题(26.10)的前m列是线性无关的，构成基向量，则约束条件可以写作：
m
"
i=1
P
i
x
i
= b −
n
"
i=m+1
P
i
x
i
(26.12)
矩阵形式为BX
B
= b − NX
N
，当X
N
=0时，解X =(x
1
,...,x
m
, 0,...,0)称为基解，其中非零元
数目不大于方程个数。如果解X ≥ 0，称它是基可行解。基可行解空间的基，称为可行基。
26.1.5. 拉格朗日乘子与影子价格
在进行项生产活动的过程中，如果投入的生产要素为x
1
,...,x
n
，产出为u = f(x
1
,...,x
n
)，
则在资源总量为a，且满足限制ϕ(x
1
,...,x
n
)=a时，要求最大的产出，可以运用条件极值的方法，
通过构造拉格朗日函数
F (x
1
,...,x
n
; λ)=f(x
1
,...,x
n
)+λ[a − ϕ(x
1
,...,x
n
)] (26.13)
来求解，这里资源总量a是一个常数。
现在考虑：如果资源总量a是一个变量，那么a的变化将会对产出u = f (x
1
,...,x
n
)产生什么样
的影响呢？
为讨论简单起见，不妨设产出为二元函数u = f( x, y)，其中x, y为两种生产要素，约束条件
为ϕ(x, y)=a，则构造拉格朗日函数：
F (x, y; λ)=f(x, y)+λ[a − ϕ(x, y)] (26.14)
产出最大化的必要条件为：
F
′
x
= f
′
x
− λϕ
′
x
=0
F
′
y
= f
′
y
− λϕ
′
y
=0
ϕ(x, y)=a
(26.15)
假设该问题存在最优解是x
0
= x(a),y
0
= y(a), λ
0
= λ(a)，它们都是a的函数，且满足：
f
′
x
(x
0
,y
0
)=λ
0
ϕ
′
x
(x
0
,y
0
)
f
′
y
(x
0
,y
0
)=λ
0
ϕ
′
y
(x
0
,y
0
)
ϕ(x
0
,y
0
)=a
(26.16)
$%"&'#(
308
)!*+",$
26.1. 线性规划
!"#
则最优产出u
0
= f(x
0
,y
0
)也是a的函数，为了考察a的变化对u
0
的影响，只需要求u
0
相对于a的边际
函数：
du
0
da
= f
′
x
(x
0
,y
0
)
dx
0
da
+ f
′
x
2
(x
0
,y
0
)
dy
0
da
= λ
0
ϕ
′
x
(x
0
,y
0
)
dx
0
da
+ λ
0
ϕ
′
x
(x
0
,y
0
)
dy
0
da
= λ
0
[ϕ
′
x
(x
0
,y
0
)
dx
0
da
+ ϕ
′
y
(x
0
,y
0
)
dy
0
da
]
(26.17)
由于ϕ(x, y)=a，则有：
ϕ
′
x
(x
0
,y
0
)
dx
0
da
+ ϕ
′
y
(x
0
,y
0
)
dy
0
da
=1 (26.18)
即
du
0
da
= λ
0
(26.19)
这个结果表明，产出最大化时的拉格朗日乘子λ
0
，正是资源总量a对最优目标函数值u
0
的边
际贡献。如果资源总量a再增加一个单位，产出将随之增加λ
0
个单位。换句话说，此时的资源投
入如果再增加一个单位，将能够带来λ
0
个单位的追加效益。不难看出，拉格朗日乘子（Lagrange
Multipliers）是有着非常明确的经济意义的。在经济学上，我们把λ
0
称为产出最大化时资源的影子
价格（Shadow Price）。
影子价格又称会计价格、最优计划价格。假设种资源的市场价格为p，如果我们将一个单位
的这种资源投入到项生产活动中可以产生p单位的效益，则数量p 就反映了这种资源在该项生产
活动中的“价值”。 在经济学上，我们就把数量p称为这种资源在该项生产活动中的影子价格。显
然，影子价格不同于市场价格，且对于同一种资源来说，在不同的企业、不同的时期，其影子价格
也是不同的。从影子价格的经济学意义可以看出，影子价格实际上是资源投入项生产活动的潜在
边际效益，它反映了产品的供求状况和资源的稀缺程度。而且资源的数量、产品的价格都影响着影
子价格的大小。一般说来，资源越丰富，其影子价格就越低，反之亦然。正因为如此，企业的管理
者在进行科学决策的时候，影子价格是必须要参考的重要依据之一。
➊
注解26.1. 在Lingo结果报告中，有关于资源对对对偶偶偶价价价格格格（（（Dual Price））） 的 列 ， 实 际 上 正 是 上 文 所 述
的影影影子子子价价价格格格。对于最大值问题，影子价格表示增加单位的资源总量（或右手边（Right Hand Side，
RHS））引起目标函数的增加量；对于最小值问题，影子价格表示减少单位的资源总量引起目标函
数下降的量。 如果在最优生产计划下，某种资源的使用小于相应的资源总量限制（资源有剩余），
则该资源的影子价格必然为零。Lingo报告中，还有一列产品差差差额额额成成成本本本（（（Reduced Cost））），有时又称
作是机机机会会会成成成本本本（（（Opportunity Cost））），实际上产品的“差额成本”等于生产该产品的“机会成本”
同销售产品的“市场价格”之差，对于标准最大化问题，其对偶问题的负的松弛变量即表示差额成
本。
在 利 润 最 大 化 的 生 产 计 划 问 题 中（见 图26.1）， 生 产 一 个 单 位 的 第 2种产品消耗的资源
为a
12
,...,a
m2
，各种资源的影子价格由对偶问题26.2 所反映，即y
1
,y
2
,...,y
m
，则生产第2种产品
➊ http://blog.sciencenet.cn/blog-383627-307859.html
$%"&'#(
309
)!*+",$
搜索与排名 Searching and Ranking
!"#
图 26.1: 原始问题：利润最大化的生产计划问题 图 26.2: 对偶问题：成本最小化的资源定价问题
的机会成本（或隐含成本）可以计算：
a
12
y
1
+ ...+ a
m2
y
m
(26.20)
如图26.2所示，根据市场定价，一个单位第2种产品的价格是c
2
，则生产一个单位产品的隐含成
本1单位该产品的价格之间的差即是差额成本：
y
m+2
= a
12
y
1
+ ...+ a
m2
y
m
− c
2
= −σ
2
(26.21)
只有当单位产品的市场价格大于隐含成本时，有利可图，（单纯形法中的）检验数（σ
2
）为正，
始安排生产（进基）；否则不安排。如果价格等隐含成本，检验数为0，表明通过生产本产品增加总
利润的可能性已不存在（变量出基后不再进基）。
26.2 单纯形法
1947年， 美 国 数 学 家George Dantzig[300]出一种解线性规划的迭代算法– 单纯形法
（Simplex Method），它被评为20世纪十大算法之一。
➊
单纯形法在初始表上重复进行初等行变换（Elementary Row Operation，ERO）、 高 斯 消 元 ，
由于线性规划方程数目小于变量个数时出现不定解，单纯形法从线性方程组中寻找单纯形，
➋
从每
一个单纯形中确定一组解，根据单纯形对应的解对目标函数的改变方向，选择下一个单纯形，逐步
趋于最值点。
➊ 20世纪十大算法包括：Metropolis Algorithm for Monte Carlo[304], Simplex Method for Linear Programming[305], Krylov Sub-
space Iteration Methods[306], The Decompositional Approach to Matrix Computations[307], The Fortran Optimizing Compiler[308],
QR Algorithm for Computing Eigenvalues[309], Quicksort Algorithm for Sorting[310], Fast Fourier Transform[311], Integer Relation
Detection[312], Fast Multipole Method[313].
➋ 一个d维空间R
d
中的单纯形S定义为由d +1个顶点x
1
,...,x
d
∈ R
d
所构成的凸包（Convex Hull）。单纯形就是一种简单的几何形状，
通过将d 维空间中的d +1个顶点相连形成。比如，一维单纯形是一个边，二维单纯形是一个三角形，三维单纯形就是一个四面体。
$%"&'#(
310
)!*+",$
26.2. 单纯形法
!"#
首先考虑如下包含n个变量，m个约束的标准形式的线性规划问题：
max z = c
T
x
s.t. Ax = b
x ≥ 0
(26.22)
其中，A ∈ R
m×n
，一般地，m ≤ n。对于约束条件中包含不等式的情况，可以通过添加松弛变量
（Slacks）或剩余变量（Surplus）将其转化为等式形式。
为方便计算，向模型中添加m个“人工变量（Artificial Variables）”：
x
n+i
= b
i
−
n
"
j=1
a
ij
x
j
则人工变量非负。我们将x
1
,...,x
n
称作非基变量，它们构成的集合记做N，x
n+1
,...,x
n+m
称作基
变量，它们构成的集合记做B。由人工变量的构造可知，非基变量可由基变量线性表出。
在此基础上，可以构建如下形式的规划问题：
max z =
(
i∈N
c
i
x
i
s.t. x
i
= b
i
−
(
n
j=1
a
ij
x
j
,i∈ B
x ≥ 0
(26.23)
使用单纯形法求解线性规划问题的基本步骤：
（1） 确定初始基可行解；
（2） 做最优性检验，通过，即是最优解，否则，转入下一步；
（3） 转换可行解，得到相邻的基可行解，转入上一步。
在单纯形法中，原问题的最优解满足：
（1） 是基本解
（2） 可行（X
B
= B
−1
b ≥ 0）
（3） 检验数C − C
B
B
−1
A ≤ 0,YA ≤ C，即对偶解可行
根据线性规划的基本性质可知，当线性规划方程的数目小于变量个数时，规划问题存在不定
解。使用单纯形法解线性规划时，如果线性规划变量个数小于约束条件的个数，则将其转换为对偶
问题再求解。
单纯形法易于程序实现，而且容易使用，但是它只适用于可以转化为标准型的规划问题，此外
算法迭代次数将随着约束和变量数目的增加而迅速上升。一般地，单纯形法需要2n − 3n步迭代，
其中n表示原问题中变量的个数。Karmarkar方法通常仅仅需要迭代次数在10 − 100次即可求解模
型，显然在求解大规模问题时更具优势。
$%"&'#(
311
)!*+",$
搜索与排名 Searching and Ranking
!"#
26.3 内点法
内点法（Interior Point Method）， 又 称 牛顿障碍法（Newton Barrier Method）， 由 John von
Neumann出，后经由Narendra Karmarkar发展成为解线性规划的一个重要方法。
例26.1. 考虑如下线性规划问题：
max x
1
+ x
2
s.t. 2px
1
+ x
2
≤ p
2
+1,p= {0, 0.1, 0.2,...,0.9, 1}
(26.24)
使用内点法求解它的过程可见图26.3： 图中，蓝色的线表示约束条件，红色线是内点法求解路径。
图 26.3: 内点法求解示例
26.4 非线性规划
在实际应用问题中，无论是问题规模还是变量个数都是庞大的，非线性程度越来越高，难以使
用简单的线性模型予以刻画。通常，可以使用目标函数或约束条件至少有一个是非线性函数的数学
规划表示，即非线性规划（Nonlinear Programming, NLP），一般可以表示成下面形式：
min f(x)
s.t. g
i
(x) ≥ 0,i =1,...,m
h
j
(x)=0,j =1,...,l
(26.25)
一般认为，非线性规划诞生的重要标志是1951年Harold W. Kuhn与Albert W. Tucker发表的最
优化条件（KT条件）[314, 315]。此后的50 年代主要是对梯度法和牛顿法的研究。在50 年代还发展
$%"&'#(
312
)!*+",$
26.4. 非线性规划
!"#
了可分离规划和二次规划的多种解法，它们大都是以单纯形法为基础。50年代末到60年代末出现了
许多解非线性规划问题的有效的算法，70年代又得到进一步的发展。非线性规划在工程、管理、经
济、科研、军事等方面都有广泛的应用，为最优设计供了有力的工具。20世纪80年代以来，计算
机的飞速发展使非线性规划的研究如虎添翼，陆续诞生了信赖域法、稀疏拟牛顿法、内点法和有限
存储法，在大规模非线性问题求解与并行计算方面也取得长足的发展。
非线性规划数值算法大致可分为直接搜索算法和 间接搜索算法。直接搜索算法仅仅利
用目标函数值信息直接建立搜索求解的方法，无须借助目标函数的梯度信息，因此又称作
无梯度（Non-gradient）算法或零阶（Zeroth Order）算法。直接搜索算法包括随机搜索方
法（Random Search）、 网 格 搜 索 方 法 （ Grid Search）、 坐 标 轮 换 法 （ Univariate Search）[316]、
旋转坐标法（Rotating Coordinates）[317]、模式方向法（Pattern Direction）[318]、 Powell方
法[319]和Nelder-Mead法[320]。间接搜索算法需要利用一阶导数或黑塞矩阵（Hessian Matrix）信
息，故又称为基于梯度的方法，典型的包括最速下降法、共轭梯度法（Fletcher-Reeves）、牛顿法、
Marquardt方法、拟牛顿法、序列二次规划方法（Sequential Quadratic Programming, SQP）、 增
广拉格朗日方法（Augmented Lagrangian Method）和非线性内点法。非线性规划目前还没有适
用于各种问题的一般算法，各个方法都有自己特定的适用范围。
26.4.1. 基本概念
由线性规划的性质我们知道，如果线性规划的最优解存在，其最优解只能在其可行域的边界上
（特别是可行域的顶点上）达到。非线性规划的最优解（如果存在）则可能在其可行域的任意一点
达到。
定义26.1 (极小值点). 给定集合S上的点ˆx，如果存在ε ≥ 0，对于任意的x ∈S，当∥x − ˆx∥≤ε时，
都有f(x) ≥ f(ˆx)，则称ˆx是函数f 在S上的局局局部部部极极极小小小值值值点点点。如果对任意的x ∈ S，都有f(x) ≥ f(ˆx)，
则ˆx是函数f在S上的全全全局局局极极极小小小值值值点点点。
定义26.2 (驻点、临界点、拐点、鞍点). 如果函数f在某点ˆx处不可微或f
′
(ˆx)=0，则称ˆx是临临临界界界点点点
（Critical Point）， 对 应 函 数 值 称 作 是 临 界 值 。 如 果 函 数 是 一 阶 可 导 的 ， 且 f
′
(ˆx)=0，则称点ˆx 是驻驻驻
点点点 （ Stationary Point）。 如 果 函 数 在 ˆx处由凹变凸或由凸变凹，则称ˆx是函数的一个拐拐拐点点点 （Inflection
Point），同时如果f(ˆx)=0，则称ˆx 是鞍鞍鞍点点点（Saddle Point）。
定理26.6 (费马定理). 给定函数f :(a, b) -→ R，如果ˆx ∈ (a, b)是f 的一个局部极值点，若f在ˆx处是可
微的，则f
′
(ˆx)=0。
定义26.3 (有效约束和无效约束). 假设x是非线性规划的一个可行解，它自然满足所有的约束。考
虑某不等式约束条件g
i
(x) ≥ 0。当g
i
(x)=0时，点x处于此约束条件形成的可行域边界上，约束条
件对点x的微小摄动构成限制，则称此约束条件是点x的有有有效效效约约约束束束（active）； 反 之 ， 当 g
i
(x) > 0时，
$%"&'#(
313
)!*+",$
搜索与排名 Searching and Ranking
!"#
约束条件对点x的微小摄动（small variation）不会构成限制，则称此约束条件是点x 的无无无效效效约约约束束束
（inactive）。所有的等式约束条件都是点x的有效约束。
定义26.4 (可行方向). 假设x
0
是非线性规划的一个可行解，若存在δ > 0，对任意的λ ∈ (0, δ)，
x
0
+ λd仍处于可行域上，则称d是x
0
的一个可可可行行行方方方向向向。
将g
i
(x)在点x
0
处作一阶泰勒展开，由
g
i
(x
0
+ λd)=g
i
(x
0
)+λd
T
∇g
i
(x
0
)+O(λ)
可知满足d
T
∇g
i
(x
0
) ≥ 0,i =1,...,m 的方向d是点x
0
的一个可行方向。
定义26.5 (下降方向). 假设x
0
是非线性规划的一个可行解，若存在δ > 0，对任意的λ ∈ (0, δ)都
有f(x
0
+ λd) <f(x
0
)，则称d是x
0
的一个下下下降降降方方方向向向。
将f(x)在点x
0
处作一阶泰勒展开，由
f(x
0
+ λd)=f(x
0
)+λd
T
∇f(x
0
)+O(λ)
可知满足d
T
∇f(x
0
) < 0的方向d是点x
0
的一个下降方向。
定义26.6 (正则点). 给定某可行解x，如果所有有效的不等式约束函数与等式约束函数在x处的梯度
线性无关，则称点x是正正正则则则点点点（Regular Point）。
根据定义可知，从点出发，沿任意可行方向（距离很小）都不能使目标函数值减少，则该点
即局部极小值点。
定理26.7 (最优化的必要性条件). 假设ˆx是目标函数f : R
n
-→ R的无约束最小值点，f是含ˆx的开
集S上的一阶连续可微函数，则满足一一一阶阶阶必必必要要要性性性条条条件件件，即∇f(ˆx)=0。如果f又是二阶连续可微的，
则满足二二二阶阶阶必必必要要要性性性条条条件件件，即∇
2
f(ˆx)半正定。
26.5 拉格朗日乘子法
拉格朗日乘子法（Lagrange Multiplier Method）主要用于解决等式约束的优化问题
min f(x)
s.t. h
i
(x)=0,i=1,...,m
(26.26)
其中，x ∈ R
n
，目标函数f : R
n
-→ R是一阶连续可微函数。
我们知道，可行域上目标函数f的极小值点应当是一个驻点
➊
，对于可行域上驻点存在一个有
用的性质：
➊ Stationary Point:函数在该点处的梯度为零
$%"&'#(
314
)!*+",$
26.5. 拉格朗日乘子法
!"#
定理26.8. 点x是可行域上驻点的充要条件是：沿任意可行方向移动点x，目标函数值都不会发生改
变。
给定点x，设V 是保持目标函数值不变的方向集，则从x出发不改变目标函数值的任意移动方向
皆垂直于目标函数在该点处的梯度。对任意的v ∈ V ⊂ R
n
，都有等式
v
T
∇f(x)=0 (26.27)
成立。
规划问题中的约束条件限制了点x可以移动的方向：同时满足所有的等式约束，不会改变函
数g
i
的值。设V
′
是x的可行方向集合，那么从x出发不“违背”约束条件的任何移动方向都垂直于约
束函数的梯度。对任意的v ∈ V
′
⊂ R
n
，都有等式
v
T
∇h
i
(x)=0,i=1,...,m (26.28)
成立。我们由此得出结论：
定理26.9. 点x是可行域上驻点的充要条件是：沿任意方向移动点x，如果它改变了目标函数的值，
则它至少会违背一个约束条件。
假设存在个从x 出发的移动方向，如果它是可行方向，即满足所有的约束条件，而且在移动
后，目标函数值发生改变，无论是变大还是变小，皆表明x 并非可行域上的驻点，矛盾产生。换句
话说，目标函数在可行域上驻点x处的梯度垂直于由V
′
生成的子集，而后者又与约束函数在x处的
梯度垂直。
定理26.10 (拉格朗日乘子定理). 假设x
∗
是目标函数f的局部极小值点，若约束函数的梯度向
量组h
1
(x
∗
),...,h
m
(x
∗
)线线线 性性性 无无无 关关关 （ x
∗
是一 个 正则 点），则 存 在唯 一的 拉 格朗 日乘 子 向量Λ
∗
=
(λ
∗
1
,...,λ
∗
m
)使得等式
∇f(x
∗
)+
m
"
i=1
λ
∗
i
∇h
i
(x
∗
)=0 (26.29)
成立。
如果f与h
i
,i=1,...,m都是二阶连续可微的，则对任意v ∈ V (x
∗
)
v
T
(∇
2
f(x
∗
)+
m
"
i=1
λ
∗
i
∇
2
h
i
(x
∗
))v ≥ 0 (26.30)
其中，V (x
∗
)=
,
v | v
T
∇h
i
(x
∗
)=0,i=1,...,m
-
。
证明： 我们下面利用罚函数法（Penalty Approach）给出证明 [321]。对于含有等式约束的非线性
规划问题，我们通过对违反约束条件的点施加惩罚，将有约束最优化问题转化为无约束最优化问
题。为了方便表示，令h =(h
1
,h
2
,...,h
m
)
T
。
$%"&'#(
315
)!*+",$
搜索与排名 Searching and Ranking
!"#
对于k =1, 2,...，定义函数：
F
k
(x)=f(x)+
k
2
∥h(x)∥
2
+
α
2
∥x − x
∗
∥
2
(26.31)
其中，x
∗
表示约束问题的局部极小值点，α > 0。第二项是对违反约束条件的点h(x)=0施加的惩
罚项，第三项则确保x
∗
是
min f(x)+
α
2
∥x − x
∗
∥
2
s.t. h(x)=0
严格局部极小值点。
由于x
∗
是一个局部极小值点，我们可以选择半径ε > 0，使得封闭球
S = {x | ∥x − x
∗
∥≤ε} (26.32)
内所有可行解x，都有f(x
∗
) ≤ f(x )。
对于约束优化问题
min F
k
(x)
s.t. x ∈ S
(26.33)
假设其最优解是x
k
，对于所有的k，则有
F
k
(x
k
)=f(x
k
)+
k
2
∥h(x
k
)∥
2
+
α
2
∥x
k
− x
∗
∥
2
≤ F
k
(x
∗
)=f(x
∗
) (26.34)
由f(x
k
)在S上有界，可知 lim
k→∞
h(x
k
)=0（否则，左端无界，矛盾）。对于每个极限值¯x =lim
k→∞
x
k
，
都有h(¯x)=0，从而对不等式 26.34 两边同时取极限可得：
f(¯x)+
α
2
∥¯x − x
∗
∥
2
≤ f(x
∗
) (26.35)
由于¯x ∈ S上的可行解，所以f(x
∗
) ≤ f(¯x)，那么¯x = x
∗
，表明当k充分大时，x
k
是F
k
(x)的无约束
局部极小值点。
下面我们利用无约束优化的相关性质证明拉格朗日乘子定理：根据一阶必要条件，当k充分大
时，有
0=∇F
k
(x
k
)=∇f(x
k
)+k∇h(x
k
)h(x
k
)+α(x
k
− x
∗
) (26.36)
成立。
由于∇h(x
∗
)的秩为m，当k充分大时∇h(x
k
)有相同的结论，所以∇h(x
k
)
T
∇h(x
k
) 可逆，那么
对于上式，整理可得：
kh(x
k
)=−(∇h(x
k
)
T
∇h(x
k
))
−1
∇h(x
k
)
T
#
∇f(x
k
)+α(x
k
− x
∗
)
$
(26.37)
两边同时取极限可知
lim
k→∞
kh(x
k
)=Λ
∗
= −(∇h(x
∗
)
T
∇h(x
∗
))
−1
∇h(x
∗
)
T
∇f(x
∗
) (26.38)
$%"&'#(
316
)!*+",$
26.5. 拉格朗日乘子法
!"#
根据等式 (26.36)可得
∇f(x
∗
)+Λ
∗
∇h(x
∗
)=0 (26.39)
由此证明了一阶优化条件。
对于F
k
(x)，根据其二阶优化条件可知，当k充分大时，对任意的α > 0，矩阵
∇
2
F
k
(x
k
)=∇
2
f(x
k
)+k∇h(x
k
)(∇h(x
k
))
T
+ k
m
"
i=1
h
i
(x
k
)∇
2
h
i
(x
k
)+αI (26.40)
是半正定的。给定v ∈ V (x
∗
)，则v
T
∇h(x
∗
)=0，取v
k
表示v在(∇h(x
k
))
T
零子空间上的投影，则
v
k
= v −∇h (x
k
)
#
(∇h(x
k
))
T
∇h(x
k
)
$
−1
v
T
∇h(x
k
) (26.41)
已知(v
k
)
T
∇h(x
k
)=0，∇
2
F
k
(x
k
)半正定，可得
0 ≤ (v
k
)
T
∇
2
F
k
(x
k
)v
k
=(v
k
)
T
#
∇
2
f(x
k
)+k
m
"
i=1
h
i
(x
k
)∇
2
h
i
(x
k
)
$
v
k
+ α∥v
k
∥
2
(26.42)
由于x
k
→ x
∗
，kh
i
(x
k
) → λ
∗
i
，并且v
T
∇h(x
∗
)=0，则v
k
→ v，对任意的v ∈ V (x
∗
)，都有
0 ≤ v
T
#
∇
2
f(x
∗
)+
m
"
i=1
λ
∗
i
∇
2
h
i
(x
∗
)
$
v + α∥v∥
2
(26.43)
参数α可以取任意小的正数，则对任意的v ∈ V (x
∗
)有
0 ≤ v
T
#
∇
2
f(x
∗
)+
m
"
i=1
λ
∗
i
∇
2
h
i
(x
∗
)
$
v (26.44)
证毕。
根据一阶必要条件和等式约束条件，由m + n个方程构建的包含m + n 个未知数的方程组，并
确定出唯一解。对于一阶必要条件，要求约束函数的梯度向量组线性无关，如果是线性相关的，那
么会产生什么影响呢？我们下面给出一个例子说明这个问题：
min f(x)=x
1
+ x
2
s.t. h
1
(x)=(x
1
− 1)
2
+ x
2
2
− 1=0
h
2
(x)=(x
1
− 1)
2
+ x
2
2
− 4=0
(26.45)
它只包含一个可行点(0, 0)，并且是极小值点，两个约束函数的梯度是线性相关的，因此对于任意
的λ
1
, λ
2
，等式
∇f(x
∗
)+λ
1
∇h
1
(x
∗
)+λ
2
∇h
2
(x
∗
)=0
始终不成立。
定理26.11 (充分条件). 假设f和g是二阶连续可微函数，定义拉格朗日函数
L(x, Λ)=f(x)+Λ
T
h(x) (26.46)
$%"&'#(
317
)!*+",$
搜索与排名 Searching and Ranking
!"#
设x
∗
∈ R
n
，Λ
∗
∈ R
m
满足
∇
x
L(x
∗
, Λ
∗
)=0
∇
Λ
L(x
∗
, Λ
∗
)=0
(26.47)
并且对任意的v ̸=0，当v
T
∇h(x
∗
)=0时，
v
T
∇
2
xx
L(x
∗
, Λ
∗
)v>0 (26.48)
则x
∗
是非线性等式约束优化问题的严格局部最小值点。
充分条件可以使用罚函数法或者可行方向法予以证明 [321]。
26.5.1. 对偶问题与KKT条件
KKT（Karush-Kuhn-Tucker）条件，又称KT（Kuhn-Tucker）条件，是非线性规划的一个重
要结论。早在1939 年，William Karush[314] 在其博士论文中就得到这一结论，然而直到Harold
W. Kuhn与Albert W. Tucker在1951 年首次公开发表[315]，人们才意识到这一工作的重要性。
KKT条件是对拉格朗日乘子法（Lagrange Multipliers Method）的推广，对于包含不等式约
束条件的非线性规划问题：
min f(x)
s.t. g
i
(x) ≤ 0,i =1,...,m
h
j
(x)=0,j =1,...,l
(26.49)
其中，
x
∈
R
n
，f : R
n
-→ R是目标函数，g
i
: R
n
-→ R, (i =1,...,m)是不等式约束函数，
h
j
: R
n
-→ R, (j =1,...,l)是等式约束函数，都是一阶连续可微函数。
构造拉格朗日函数
L(x, α, β)=f(x)+
m
"
i=1
α
i
g
i
(x)+
l
"
j=1
β
j
h
j
(x) (26.50)
其中，α ≥ 0, β为拉格朗日乘子向量。
记Θ
p
(x) = max
α≥0,β
L(x, α, β)，由于
max
α≥0,β
L(x, α, β)=f(x) + max
α≥0,β
7
m
"
i=1
α
i
g
i
(x)+
l
"
j=1
β
j
h
j
(x)
8
(26.51)
• 当x是可行解时，满足所有约束条件，则h
j
(x)=0,j =1,...,l，函数Θ
p
(x)可简化为
Θ
p
(x)=f(x) + max
α≥0,β
m
"
i=1
α
i
g
i
(x)
由于g
i
(x) ≤ 0,i =1, . . . , m, α ≥ 0，则
m
(
i=1
α
i
g
i
(x) ≤ 0，当且仅当α =0时等式成立，从
而Θ
p
(x)=f(x)。
$%"&'#(
318
)!*+",$
26.5. 拉格朗日乘子法
!"#
• 当x不可行时，违反至少一个约束条件，则通过灵活选择α, β可得Θ
p
(x)=∞。
综上可得
Θ
p
(x) = max
α≥0,β
L(x, α, β)=
⎧
⎨
⎩
f(x) x是可行解
∞ 否则
(26.52)
利用拉格朗日函数可以建立如下无约束最优化问题
ˆp =min
x
Θ
p
(x)=min
x
max
α≥0,β
L(x, α, β) (26.53)
它与原始问题(26.49)等价（解相同），我们称其为原问题。如果记Θ
d
(α, β)=min
x
L(x, α, β)，容易
证明Θ
d
(α, β) ≤ ˆp。假设x
0
是问题(26.53) 的任意可行解，即g
i
(x
0
) ≤ 0,h
j
(x
0
)=0, α ≥ 0，则有
m
"
i=1
α
i
g
i
(x
0
)+
l
"
j=1
βh
j
(x
0
) ≤ 0 (26.54)
从而
L(x
0
, α, β)=f(x
0
)+
m
"
i=1
α
i
g
i
(x
0
)+
l
"
j=1
βh
j
(x
0
) ≤ f(x
0
) (26.55)
也就是说
Θ
d
(α, β)=min
x
L(x, α, β) ≤ L (x
0
, α, β) ≤ f (x
0
) (26.56)
由x
0
的任意性可知，Θ
d
(α, β) ≤ ˆp成立。
对于下面形式的问题
ˆ
d = max
α≥0,β
Θ
d
(α, β) = max
α≥0,β
min
x
L(x, α, β) (26.57)
构成原问题的对偶问题，必然有
ˆ
d ≤ ˆp（弱对偶性，Weak Duality） (26.58)
差值ˆp −
ˆ
d称为对偶间隙（Duality Gap），当对偶间隙为零时
ˆ
d =ˆp（强对偶性，Strong Duality） (26.59)
当强对偶性成立时，原问题的最优解是ˆx，对偶问题的最优解是ˆα,
ˆ
β，则拉格朗日函数在可行
解(ˆx, ˆα,
ˆ
β)的函数值同原始问题目标函数值f(ˆx)相等，即
L(ˆx, ˆα,
ˆ
β)=f(ˆx)+
m
"
i=1
ˆα
i
g
i
(ˆx)+
l
"
j=1
ˆ
β
j
h
j
(ˆx)=f(ˆx) (26.60)
• ˆx是Θ
p
(x)的极值点，由拉格朗日函数与Θ
p
(x)的定义，极值点的必要条件，可得
∇f(ˆx)+
m
"
i=1
ˆα
i
∇g
i
(ˆx)+
l
"
j=1
ˆ
β
j
∇h
j
(ˆx)=0 (26.61)
我们称之为平稳性（Stationarity）条件。
$%"&'#(
319
)!*+",$
搜索与排名 Searching and Ranking
!"#
• ˆx是原始问题的可行解，满足h
j
(ˆx)=0,j =1,...,l，所以
m
"
i=1
ˆα
i
g
i
(ˆx)=0 (26.62)
既然ˆα
i
g
i
(ˆx) ≤ 0, ∀i，则有
ˆα
i
g
i
(ˆx)=0,i=1,...,m (26.63)
我们称之为互补松弛（Complementary Slackness）条件，当ˆα
i
> 0时，必然有g
i
(ˆx)=0，
表明此不等式约束条件是有效的。
I. 必要条件
如果ˆx是极小值点，并且满足一些约束规范条件（Constraint Qualification，Regularity Con-
ditions）， 则 存 在 拉 格 朗 日 乘 子 向 量 ˆα ∈ R
m
, β ∈ R
l
，满足下面几条性质（统称为KKT优化条
件）：
• 主可行性条件（Primal Feasibility Conditions）
g
i
(ˆx) ≤ 0
h
j
(ˆx)=0
(26.64)
• 对偶可行性条件（Dual Feasibility Condition）
ˆα
i
≥ 0 (26.65)
• 互补松弛条件
ˆα
i
g
i
(ˆx)=0 (26.66)
• 平稳性条件
∇f(ˆx)+
m
"
i=1
ˆα
i
∇g
i
(ˆx)+
l
"
j=1
ˆ
β
j
∇h
j
(ˆx)=0 (26.67)
其中，i =1, . . . , m, j =1,...,l。特别地，当m =0时，非线性规划仅含等式约束，则KTT条件退
化为拉格朗日条件。
II. 约束规范条件
对于极值点ˆx，若满足KKT必要性条件，必需满足一些约束规范条件。常用的约束规范条件有
• 线性无关约束规范条件（Linear Independence Constraint Qualification, LICQ）： 有 效 不 等
式约束函数与等式约束函数在ˆx处的梯度线性无关（正则点）
• Slater约束规范条件：对于凸优化问题，所有的有效不等式约束在点ˆx处是严格不等式
$%"&'#(
320
)!*+",$
26.6. 梯度下降法
!"#
III. 充分条件
假设目标函数f与不等式约束函数g
i
(i =1,...,m)都是凸函数，等式约束函数h
j
(j =1,...,l)是
仿射函数，ˆx是一个可行解，如果存在0 ≤ ˆα ∈ R
m
, β ∈ R
l
，满足互补松弛与平稳性条件，则ˆx是全
局极小值点。
IV. Fritz John必要条件
对于非线性优化问题
min f(x)
s.t. g
i
(x) ≤ 0,i =1,...,m
(26.68)
根据Farkas引理，设ˆx是一个局部极小值点，则存在不全为零的拉格朗日乘子
ˆ
λ
i
,i=0, 1,...,m，满
足下面几个性质（统称Fritz John必要条件）
ˆ
λ
0
∇f(ˆx)+
m
(
i=1
ˆ
λ
i
∇g
i
(ˆx)=0
ˆ
λ
i
≥ 0,i=0, 1,...,m
ˆ
λ
i
g
i
(ˆx)=0,i=1, 2,...,m
(26.69)
当原问题比较复杂，而对偶问题具有良好的结构时，利用原问题的强对偶性“通过求解对偶问
题实现对原问题的求解”是一种常用的优化方法，比如当目标函数与不等式约束函数都是凸函数，
等式约束函数是仿射函数，而不等式约束条件是严格的，则相应的优化问题是强对偶的。一般而
言，KKT条件是强对偶性成立的必要条件，特别地，当原问题是凸优化问题时，KKT条件是强对
偶性成立的充分必要条件。
26.6 梯度下降法
梯度下降法（Gradient Descent Method）， 又 称 最速下降法（Steepest Descent Method）， 由
大数学家Cauchy于1847年最先使用。它是最古老的一种解析方法，而其他解析方法大多承其衣钵，
并构成最优化方法的基础。梯度下降法的优点是工作量少，存储变量少，对初始点不敏感。当然，
它也有缺点。比如，对于多峰形式（Multi-Modal）的参数分布，它的收敛速度很慢，容易陷入局
部最优点。
如果我们希望利用迭代算法寻找二阶可微函数f的最小点，假设当前最优点是x
k
∈ R
n
，在选
择下一个点x
k+1
= x
k
+ αh 时，遵循的原则或目标是能够最小化函数值f(x
k+1
)，即
min
h∈R
n
f(x
k+1
)=f(x
k
+ αh)
(26.70)
其中，h表示搜索的方向，一般地∥h∥ =1，α > 0表示步长或者学习率（Learning Rate）。
根据泰勒展开公式有：
f(x
k
+ αh)=f(x
k
)+αg(x
k
)
T
h + O(α)
$%"&'#(
321
)!*+",$
搜索与排名 Searching and Ranking
!"#
其中，g(x
k
)=∇f(x
k
)。
不考虑高阶部分，等式(26.70)等价于最小化g(x
k
)
T
h。根据内积公式可知，当h = −g(x
k
)=
−∇f(x
k
) （负梯度方向）时，f (x
k
)下降的速度最快（最速下降法由此而来）。
梯度下降法在机器学习中有广泛的应用，比如在监督型学习算法中，根据经验风险最小化
原则训练模型就可能用到梯度下降法。假设预测模型一阶可微f = f (ω,x)，则预测模型在数据
集S = {x
i
,y
i
}
n
i=1
上的经验损失
L(ω,S)=
1
n
n
"
i=1
ℓ(y
i
,f(ω,x
i
)) (26.71)
其中，ℓ(y, f(ω,x))是预测模型在单个样本上的损失函数，假设它是一阶可微的。
根据梯度下降法有下面的参数更新法则有：
ω
t+1
= ω
t
+ α
1
n
n
"
i=1
∇L(y
i
,f(ω,x
i
)) (26.72)
由于每次更新模型参数时，都要用到所有的训练样本，梯度下降法也称批量梯度下降法。批量梯度
下降法在更新参数时，用到了数据集中全部的样本，时间复杂度为O(n)。
在大数据处理时，对于些类型的损失函数，计算所有样本的平均梯度值是不可行的。为此，
研究人员议使用随机抽取的单个样本计算近似的平均梯度值，利用它来更新模型参数，能够在大
幅减少单次迭代计算量的同时，还能取得同批量梯度下降法相似的效果，这种方法被称作随机梯度
下降法，也称增量随机梯度下降法。
假设随机选取的样本是x
0
，则使用下面的公式更新参数：
ω
t+1
= ω
t
+ α
t
∇L(y
0
,f(ω,x
0
)) (26.73)
为保证收敛，通常要求学习率α
t
序列满足
T
"
t=1
α
t
= ∞,
T
"
t=1
α
2
t
< ∞ (26.74)
一般选择α
t
=1/t或α
t
= α
0
(1 + λα
0
)
−1
。
随机梯度下降法由于单次迭代计算量小，效率很高，远优于经典的优化算法，如牛顿法，特别
适合大数据机器学习[322]。它也存在一些缺陷，比如在迭代过程中产生许多噪声、每次迭代的搜
索方向未必是整体梯度下降的方向等。随机梯度下降法已经被用于训练线性的支持向量机[85]、条
件随机场模型[323]。
1992年，Polyak与Juditsky[324]在随机梯度下降法的基础上，每次更新参数时，增加一个计算
参数平均值的过程
¯ω
t+1
=
t
t +1
¯ω
t
+
1
t +1
ω
t+1
(26.75)
并使用收敛的参数平均值作为最优参数，因此称为平均随机梯度下降法（Averaged Stochastic
Gradient, ASGD）。 平 均 随 机 梯 度 下 降 法 的 优 点 是 实 现 简 单 ， 无 需 计 算 目 标 函 数 Hessian 矩阵逆的
$%"&'#(
322
)!*+",$
26.6. 梯度下降法
!"#
情况下，就可以达到与二阶随机梯度下降法相似的渐进收敛速度，比随机梯度下降法还快。然而，
在实际应用中，平均随机梯度下降法却遭受冷遇[325, 326]，部分原因在于其渐进收敛速度要求训
练集足够大，在中小型训练集上的表现逊色于随机梯度下降法，部分原因是对于参数调优的要求较
高。 2010年，Wei Xu[325]设计了一套新的学习率调整规则（Learning Rate Schedule）
α
t
= α
0
(1 + γα
0
)
−c
(26.76)
扩大了平均随机梯度下降法的适用范围，使平均随机梯度下降法逐渐被接受和使用，比如Lin等
人[326]设计了一种并行随机梯度下降法，求解线性支持向量分类模型，做大规模图片分类。通常，
α
0
选择的数值很小（比如α
0
=0.01）。参数γ、c 的选取则需要具体问题具体分析。对于γ，[325]分
析认为最好选择使用目标函数Hessian矩阵的最小特征值；对于c，常见的选择有{1, 3/4, 2/3}。
26.6.1. 共轭梯度下降法
共轭梯度下降法（Conjugate Gradient Descent Method）是一种共轭方向法，我们首先熟悉
一下共轭与共轭方向的含义。
定义26.7 (共轭). 假设A ∈ R
n×n
是一个对称正定矩阵，如果存在两个向量x, y ∈ R
n
，使得
x
T
Ay =0, (26.77)
我们称x和y关于矩阵A共共共 轭轭轭 。 如 果 A是单位矩阵，则x
T
y =0，即x, y正交，可知正交是共轭的一个
特例。
定义26.8 (共轭方向). 如果存在一组非零向量{x
k
}
n
k=1
，对任意两个向量x
i
,x
j
都有
x
T
i
Ax
j
=0, (26.78)
我们称向量族{x
k
}
n
k=1
关于矩阵A共轭，每个向量对应一个共共共轭轭轭方方方向向向。
共轭方向在寻找二次型函数f(x)=c + b
T
x +
1
2
x
T
Ax 极小值点时，存在一个重要的性质：从任
意初始点出发，依次沿x
1
,x
2
,...,x
n
方向做一维最优搜索，至多n步即可收敛到极小值点。共轭方
向法有限步收敛的性质，因此可以称作是二次收敛算法。
理论与实践证明，将二次收敛算法用于非二次型目标函数，亦有很好的效果，但不能保证经过
有限次迭代达到极小点。此时，可以对目标函数在极小点附近进行泰勒展开，略去二次项之后的高
阶项对目标函数做二次逼近。
26.6.2. 坐标下降法
坐标下降法（Coordinate Descent Method）是一种非梯度优化算法，用于寻找一个多元函数
的极值点。坐标下降法在每次迭代过程，基于当前点沿个坐标方向进行一维搜索，并在后续优化
过程循环使用不同的坐标方向。
$%"&'#(
323
)!*+",$
搜索与排名 Searching and Ranking
!"#
给定多元函数
f : R
n
-→ R
搜索空间的坐标基是e
1
,...,e
n
，当前所在的点是x
t
=(x
t
1
,...,x
t
n
)。假设下一次搜索选择的方向
是e
i
，坐标下降法则固定其他坐标方向上的数值，在方向e
i
做一维搜索，以优化目标函数：
x
t+1
i
= arg min
x
F (x)=f(x
t
1
,...,x
t
i−1
,x,x
t
i+1
,...,x
t
n
) (26.79)
坐标下降法可以归结为下面的简单流程：
Algorithm 20 坐标下降法
Input: 初始点x
0
，偏差阈值ϵ
0
> 0
for k =1, 2,...do
1. 选择坐标方向e
1
进行一维搜索
x
k
1
= arg min
x
1
f(x
1
,x
k−1
2
,...,x
k−1
n
)
2. 将x
k
1
代入函数，选择坐标方向e
2
进行一维搜索
x
k
2
= arg min
x
2
f(x
k
1
,x
2
,x
k−1
3
,...,x
k−1
n
)
3. 将x
k
2
代入函数，选择坐标方向e
2
进行一维搜索
x
k
3
= arg min
x
3
f(x
k
1
,x
k
2
,x
3
,x
k−1
4
,...,x
k−1
n
)
4. 依照上步方式依次做下一个坐标方向上的搜索，直至在所有坐标方向上执行一遍，得到
新的搜索点
x
k
=(x
k
1
,x
k
2
,...,x
k
n
)
5. 计算偏差ϵ = |f( x
k−1
) − f(x
k
)|，如果ϵ < ϵ
0
，则停止迭代，否则重复上述步骤。
end for
26.6.3. 次梯度法
二十世纪八十年代，Naum Shor[327, 328]研究发明了次梯度法（Subgradient Method），它是
解凸优化问题的一种迭代算法，也适用于目标函数不可微的情形。当目标函数可微时，其搜索方向
同最速下降法相同。从速度上来看，次梯度法比牛顿法略逊一筹，但是由于它内存消耗少，适合处
理大型问题。
$%"&'#(
324
)!*+",$
26.7. 牛顿法
!"#
26.7 牛顿法
梯度法在些情况下收敛速度会很慢，无法有效地确定最优搜索方向。如果在梯度法所
使用的一阶偏导信息外，我们利用其他更多信息，比如目标函数的二阶偏导，则可以有效改
善梯度法，高效地确定最优搜索方向，此法称作牛顿法（Newton’s Method）或牛顿-拉弗森法
（Newton-Raphson Method）
➊
。
假设目标函数f(x)在点x = x
k
有定义并且连续，则我们可以在此处做二阶泰勒展开，从而有
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
).
为了搜索到目标函数的最小值点，我们使用迭代算法思想，从当前点x = x
k
出发，搜索下一个
最优点x = x
k+1
，使得目标函数值不断下降。如果迭代算法先后两次优化结果数值上比较接近，则
可以忽略高阶项O((x − x
k
)
2
)，并得到目标函数的近似表达式
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
). (26.80)
根据连续函数极值条件可得∇f(x)=∇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+1
) −∇f(x
k
)). (26.81)
如果黑塞矩阵是奇异矩阵，需要对算法进行修正。由于牛顿法利用目标函数的二阶偏导信息修正搜
索方向，收敛速度很快。同时，牛顿法在处理黑塞矩阵，特别是计算黑塞矩阵的逆矩阵时，会产生
大量的计算开销。牛顿法利用二阶泰勒展开式逼近目标函数，如果f(x)对应一个二次函数，可以直
接确定最优解析解。多数情况下，目标函数不属于二次函数，由于泰勒展开存在逼近误差，算法不
可能直接搜索到最优值点。当目标函数属于高次函数时，优化速度会受到很大影响。
假设目标函数f(x)二阶可微，则函数的一阶偏导，又称梯度向量，形式如下：
g(x)=∇f =(
∂f
∂x
1
,
∂f
∂x
2
,...,
∂f
∂x
n
)
T
. (26.82)
函数的二阶偏导矩阵，又称黑塞矩阵（Hessian Matrix），形式如下：
H(x)=∇
2
f =
⎡
⎢
⎢
⎢
⎣
∂
2
f/∂x
2
1
··· ∂
2
f/∂x
1
∂x
n
.
.
.
.
.
.
.
.
.
∂
2
f/∂x
n
∂x
1
··· ∂
2
f/∂x
2
n
⎤
⎥
⎥
⎥
⎦
= H
T
(x) (26.83)
黑塞矩阵H(x)亦称为梯度向量g( x)的雅克比矩阵（Jacobian Matrix）。
➊ 牛顿法最初由艾萨克·牛顿发明，后于1690年由约瑟夫·拉弗森再次出。
$%"&'#(
325
)!*+",$
搜索与排名 Searching and Ranking
!"#
26.8 拟牛顿法
1959年，物理学家William C. Davidon[329]设计出拟牛顿法（Quasi-Newton Methods），
也称变尺度法（Variable Metric Method），以解决非线性优化问题。 根据牛顿法最优搜索路径：
x
k+1
− x
k
=(∇
2
f(x
k
))
−1
(∇f(x
k+1
) −∇f(x
k
))，黑塞矩阵反映了目标函数的二阶信息，有利于
加快算法的迭代速度，但计算黑塞矩阵逆矩阵的运算却十分复杂。William Davidon[329]设计出
一种新的算法，既能够利用目标函数的二阶信息，还能避免直接进行逆矩阵运算：通过构造新
的矩阵，逼近黑塞矩阵的逆矩阵。矩阵近似逼近的策略是拟牛顿算法的核心思想，并以此衍生
出多种形式的优化算法，如DFP更新方程[330]（DFP Updating Formula）、 BHHH方法SR1方程
（Symmetric Rank One）、 BFGS法与L-BFGS（Limited-memory BFGS）[331, 332, 333, 334]等高效
的优化方法。
拟牛顿法使用
¯
H
k+1
表示黑塞矩阵逆矩阵[∇
2
f(x
k
)]
−1
的近似矩阵，每一步迭代都充分利用当前
信息确定下一个搜索方向，确保目标函数值随迭代次数增加而稳步下降，近似矩阵最终将收敛于最
优解对应的黑塞矩阵的逆。对于目标函数f(x)，拟牛顿法根据下面形式的搜索规则
x
k+1
− x
k
=
¯
H
k+1
(∇f(x
k+1
) −∇f(x
k
)) (26.84)
实现高效搜索。为方便表示，我们记
⎧
⎨
⎩
∆x
k
= x
k+1
− x
k
∆G
k
= ∇f(x
k+1
) −∇f(x
k
)
(26.85)
则拟牛顿法的搜索路径可以简单写作：
∆x
k
=
¯
H
k+1
∆G
k
. (26.86)
假设近似矩阵
¯
H
k+1
（尺度矩阵）有下面形式的结构
¯
H
k+1
=
¯
H
k
+ ∆
¯
H
k
, (26.87)
在上次迭代的近似矩阵
¯
H
k
基础上进行调整，∆
¯
H
k
即校正矩阵。一般地，算法要求近似矩阵为对称
正定矩阵，则校正矩阵也是对称矩阵。如果将拆分形式的近似矩阵
¯
H
k+1
代入(26.84)可得：
∆x
k
=(
¯
H
k
+ ∆
¯
H
k
)∆G
k
=
¯
H
k
∆G
k
+ ∆
¯
H
k
∆G
k
. (26.88)
此时，我们的目标就是构造合适的校正矩阵，以满足：
∆
¯
H
k
∆G
k
= ∆x
k
−
¯
H
k
∆G
k
. (26.89)
由于校正矩阵是对称矩阵，我们利用两个列向量∆x
k
与∆G
k
，构造如下形式的校正矩阵可以设
想，∆
¯
H
k
存在一种简单的形式：
∆
¯
H
k
= ∆x
k
(η
k
∆x
k
)
T
−
¯
H
k
∆G
k
(ξ
k
¯
H
k
∆G
k
)
T
. (26.90)
$%"&'#(
326
)!*+",$
26.8. 拟牛顿法
!"#
等式左右两端同时乘以∆G
k
，则有：
∆
¯
H
k
∆G
k
= ∆x
k
(η
k
∆x
k
)
T
∆G
k
−
¯
H
k
∆G
k
(ξ
k
¯
H
k
∆G
k
)
T
∆G
k
= ∆x
k
−
¯
H
k
∆G
k
(26.91)
根据对应关系，我们可以取
⎧
⎨
⎩
(η
k
∆x
k
)
T
∆G
k
=1
(ξ
k
¯
H
k
∆G
k
)
T
∆G
k
=1
(26.92)
如果(∆x
k
)
T
∆G
k
，(
¯
H
k
∆G
k
)
T
∆G
k
都不等于零，我们可以立即确定两个待定参数
⎧
⎨
⎩
η
k
=1/(∆x
k
)
T
∆G
k
ξ
k
=1/(
¯
H
k
∆G
k
)
T
∆G
k
(26.93)
从而得到校正矩阵
∆
¯
H
k
=
∆x
k
(∆x
k
)
T
(∆x
k
)
T
∆G
k
−
(
¯
H
k
∆G
k
)
T
(
¯
H
k
∆G
k
)
T
(
¯
H
k
∆G
k
)
T
∆G
k
(26.94)
及尺度矩阵
¯
H
k+1
=
¯
H
k
+
∆x
k
(∆x
k
)
T
(∆x
k
)
T
∆G
k
−
(
¯
H
k
∆G
k
)
T
(
¯
H
k
∆G
k
)
T
(
¯
H
k
∆G
k
)
T
∆G
k
. (26.95)
迭代时常取第一个尺度矩阵
¯
H
0
为单位矩阵。
定理26.12. 如果x
k
不是目标函数f(x)的极小值点，且
¯
H
k
是正定矩阵，则有
(∆x
k
)
T
∆G
k
̸=0, (
¯
H
k
∆G
k
)
T
∆G
k
̸=0,
根据式子(26.95)构造尺度矩阵。如果
¯
H
k
是对称正定矩阵，则根据式子(26.95)确定的尺度矩
阵
¯
H
k+1
也是对称正定矩阵。拟牛顿法确定的搜索方向是梯梯梯度度度下下下降降降方方方向向向。
26.8.1. DFP变尺度法
拟牛顿法以较小的代价构造黑塞矩阵的近似矩阵，不要求黑塞矩阵的近似矩阵是非奇异矩
阵（可逆），能够产生超线性收敛性，性能优于牛顿法与梯度法。 根据拟牛顿法思想，Roger
Fletcher与Michael Powell[330]发展出DFP变尺度法（Davidon-Fletcher-Powell Variable Metric
Method）。
26.8.2. BFGS
BFGS（Broyden-Fletcher-Goldfarb-Shanno）算法是CG Broyden[331]、Roger Fletcher[332]、
Donald Goldfarb[333]、David Shanno[334] 各自独立出的一种拟牛顿法，它是目前最流行的也
是最有效的一种校正方法。BFGS 是DFP 法的对偶形式，在处理非平滑优化问题时，BFGS也能够
取得很好的性能。它还存在一个变体，有限内存的BFGS 算法（Limited-memory BFGS, L-BFGS），
实现占用的内存远小于BFGS 算法。BFGS 在计算过程中，需要存储一个n × n的矩阵H。对于高维
矩阵，如n = 100, 000，若使用双精度类型存储矩阵H，占用内存100, 000 ×100, 000 ×8B≈ 74.5GB，
因此L-BFGS在实际应用中的意义显而易见。
$%"&'#(
327
)!*+",$
搜索与排名 Searching and Ranking
!"#
Algorithm 21 DFP变尺度法
Input: 初始点x
0
，梯度偏差阈值ϵ > 0，初始尺度矩阵
¯
H
0
= I（单位矩阵）
for k =1, 2,...do
1. 如果梯度偏差∥∇f(x
k
) −∇f(x
k−1
)∥≤ϵ，则认定x
k−1
是近似极小值点，停止迭代。否则，
转向下一步。
2. 计算尺度矩阵
¯
H
k
=
¯
H
k−1
+
∆x
k−1
(∆x
k−1
)
T
(∆x
k−1
)
T
∆G
k−1
−
(
¯
H
k−1
∆G
k−1
)
T
(
¯
H
k−1
∆G
k−1
)
T
(
¯
H
k−1
∆G
k−1
)
T
∆G
k−1
3. 取搜索方向d
k
= −
¯
H
k
(∇f(x
k
) −∇f(x
k−1
))，并在d
k
方向上进行一维搜索，确定最佳步
长λ
k
：
λ
k
= arg min
λ>0
f(x
k
+ λd
k
)
4. 确定下一个搜索点x
k+1
= x
k
+ λ
k
d
k
end for
图 26.4: Army of Four: Broyden-Fletcher-Goldfard-Shanno
$%"&'#(
328
)!*+",$
26.9. 分枝定界法
!"#
26.9 分枝定界法
1960年，Ailsa Land与Alison Doig[335]创立分枝定界法（Branch and Bound Method），它是
解决整数规划问题最常用的一种算法，也可求解混合整数线性规划问题（MILP）。 以 最 小 化 目 标
函数的线性规划问题为例，假设S是问题的可行解集合，分枝定界法主要包括三个基本环节：分枝
（Branching）、 定界 （Bounding）与剪枝（Pruning）。 在 分 枝 阶 段 ， 将 S 分割成多个相互不交的
子集S
1
,S
2
, ···，满足S =
P
i
S
i
，并计算目标函数在每个子集上的最小值点；定界阶段计算目标函
数在每个子集的上界与下界；分枝定界法法最关键的环节在于对可行解子集的筛选，即剪枝环节。
在剪枝环节，根据目标函数在定界阶段的计算结果可以判定，如果目标函数值在子集甲上的下界比
在子集乙的上界都要大，可以直接舍弃子集甲。分枝定界法通过多次调用分枝、定界与剪枝三个环
节，逐渐逼近最优解。
在实际使用时，可以设计一定的启发式规则，当目标函数值上界与下界的差值比较接近时，则
前终止分枝步骤，从而大幅减少搜索时间。
26.10 割平面法
20世纪50年代，Ralph Gomory出割平面法（Cutting-Plane Method）。
26.11 随机搜索算法
随机搜索算法（Random Search）对于目标函数的连续性、可微性无任何限制，是最简单、最
直观的一种搜索算法。随机搜索一说最早源自[336, 337, 338]。随机算法迭代地从一个点移动到另
一个点，以期后一个点相对于目标函数优于前者。随机搜索算法可概括为如下基本流程：
Algorithm 22 随机搜索算法
Input: 初始点x，最大迭代次数T ，偏差阈值ϵ
repeat
• 从超球面S = {x|x ∈ R
n
, ∥x − x
t−1
∥ = δ
t
} 中抽取数据点y
• 如果f(y) <f(x)，则设置临时最优点x = y
until t ≤ T 或者f(x) − f (y) > ϵ
Output: x是最优点
随机搜索算法在机器学习领域中应用广泛，经常使用的随机搜索方法有：模拟退火算法
（Simulated Annealing）、 遗 传 算 法 （ Genetic Algorithms）、 进 化 优 化 （ Evolutionary Program-
ming）、 粒 子 群 优 化 算 法 （ Particle Swarm Optimization）、 蚁 群 优 化 （ Ant Colony Optimiza-
$%"&'#(
329
)!*+",$
搜索与排名 Searching and Ranking
!"#
tion）、 交 叉 熵 （ Cross Entropy）、 随 机 近 似 算 法 （ Stochastic Approximation）、 多 点 搜 索 算 法
（Multi Start）、禁忌搜索算法（Tabu Search, Taboo Search）[339]等。
随机搜索算法因其简单易于实现的特性，得到广泛应用。对于大型优化问题，相比确定性算法
（如分支界定法）随机搜索算法则略胜一筹，遗憾之处在于它无法保证渐进收敛到最优解，只能通
过大量的试错搜索，依概率收敛到个未必全局最优的点。
26.11.1. 模拟退火算法
在热力学中，“退火”（Annealing）是指物体缓慢降温的物理现象，温度越低，则物体的能量
就越低。当温度下降到一定程度后，液体开始冷凝与结晶，在结晶状态，系统能量达到最低。如果
物体快速降温，亦称“淬火”（Quenching）， 则 可 能 产 生 能 量 不 是 最 低 的 非 结 晶 态 。 退 火 过 程 使 得
物体分子在任何温度下，都能有充足的时间、以较大的概率找到使其内能更低的位置，最终系统也
最稳定。
1953年，Metropolis等人[340]出使用蒙特卡洛模拟（Monte Carlo Simulation）方法确定热
槽中分子的稳态分布。 他们在模拟中使用的规则是：给定一个能量为E
1
的状态S
1
，如果此时移
动一个分子的位置，系统状态发生转移，从S
1
转移到S
2
，能量从E
1
变成E
2
。如果E
2
− E
1
≤ 0，
则接 受此状 态。 反 之，如 果能 量升高（恶 化状 态），系统 不会 直 接拒 绝新状 态S
2
，而是以概
率exp(
E
2
−E
1
κT
)有条件地接受它，其中κ > 0是Boltzmann常数，T 表示热槽的温度。通过不断重复这
一流程，他们推断：在给定的温度T 下，分子能量满足标准的Boltzmann分布：
P (E = E(S)) =
1
Z(T )
exp
*
−
E(S)
κT
+
(26.96)
其中，E是代表分子能量的随机变量，Z(T )是标准化因子：
Z( T )=
"
S∈S
exp
*
−
E(S)
κT
+
这个公式就是著名的Metropolis准则。
根据推断，在温度T 时，分子从状态S
1
到S
2
，能量从E
1
变成E
2
。根据Boltzmann分布可以计算
两种状态下的概率差值：
∆P = P (E = E
2
) − P (E = E
1
)=
1
Z( T )
exp
*
−
E
2
κT
+#
1 − exp
*
E
2
− E
1
κT
+$
如果E
1
<E
2
，由于∆E = E
2
− E
1
> 0，则exp
*
E
2
−E
1
κT
+
> 1，从而∆P<0。由此可以说明，温度
不变时，分子停留在低能量状态S
1
的概率比停留在高能量状态S
2
的概率要大。根据Boltzmann分布
可知，温度越高，不同能量状态对应的概率相差就越小，当温度足够高时，不同能量状态下的概率
基本相同。随着温度下降，低能量状态对应的概率越来越大，当温度趋于0时，其概率达到1。
1983年，S cott Kirkpatrick等人[341]首次将模拟退火（Simulated Annealing）同组合优化
（Combinatorial Optimization）联系起来。他们将能量对应到成本函数，将系统状态对应到组合
$%"&'#(
330
)!*+",$
26.11. 随机搜索算法
!"#
优化问题的解，系统中粒子的一次移动等价地视为组合优化问题中的一次试验（Trial）。 系 统 温 度
在组合最小化问题中就是一个控制参数，首先“冷却”处于高温条件下的解空间，然后，逐渐降温
直至系统在一个稳定解“结晶”。 1985年，Vlado
˘
Cern
´
y[342]独立地发现可将模拟退火原理应用于
货郎担问题（Traveling Salesman Problem, TSP），并取得巨大成功。
模拟退火算法（Simulated Annealing Algorithm）能够实现全局最优或近似全局最优，奥秘
就在于其接受状态的策略：除了接受优化解，还使用一个随机接受准则（Metropolis准则）有限度
地接受恶化解（Deterioration）。 模 拟 退 火 算 法 能 够 逃 脱 局 部 最 优 点 的 园 囿 ， 很 大 程 度 上 归 功 于 那
些看似不利状态的功劳，实现全局最优或接近全局最优。
Algorithm 23 模拟退火算法
Input: 系统温度T ，允许最低温度T
0
，迭代最大次数N，初始解x
1
，Boltzmann常数设为0 < κ < 1
repeat
1. 从当前解x
1
的邻域中随机选择点x
2
，计算∆E ! f (x
2
) − f(x
1
)
2. 如果∆E<0，则接受x
2
为新解；否则，以概率
P =exp(
−∆E
T
)
接受x
2
。具体地，判定P>rand(0, 1)是否成立。若成立，则接受x
2
；否则，拒绝x
2
。
3. 以速率κ降温：T ← κT
4. i −−
until i<0或T<T
0
模拟退火算法至少有两个重要的控制参数：初始温度T 和Boltzman常数κ > 0，需谨慎选择与
调整。它们主要从下面几点影响算法搜索的进程：
（1） 初始温度T的设置问题：初始温度的设置是影响模拟退火算法全局搜索性能的重要因素之
一，初始温度高，则搜索到全局最优解的可能性大，但需要消耗大量的计算时间；反之，效
率增高但搜索性能可能受到影响。在实际应用中，初始温度一般需根据多次实验结果进行选
择，比如David Fogel[343]在处理TSP问题时，根据100条随机产生的路线，取出最短与最长
距离d
0
,d
1
，并使用下面公式设定初始温度：
T = −(d
1
− d
0
)/ log 0.9
（2） 退火管理问题（Annealing Schedule）： 模 拟 退 火 算 法 全 局 搜 索 性 能 也 和 退 火 管 理 方 式 密 切
相关。一般而言，在同一温度下“充分”搜索（退火）是相当必要的，但需要计算开销。在
实际应用中，需要针对具体问题的不同性质和特征设置合理的退火平衡条件。
$%"&'#(
331
)!*+",$
搜索与排名 Searching and Ranking
!"#
在模拟退火算法构造新解、接受新解过程，存在几个细节性的问题需要注意：（1）新解的产生
一般是对当前解的一种简单变换，比如改变当前解的部分元素（逆序、交换、插入等），相比重新
完全生成时间开销要小。（2）计算当前解与新解的目标函数之差时，如果能够分离出变换部分对目
标函数差的独立影响，则目标函数差的计算就可以采用增量方式。（3）当确定接受新解时，需要以
新解代替当前解。若已知产生新解的变换流程，可只就变换部分进行替换，从而节省完全替换的时
间。
模拟退火算法确定的最优解与初始解无关，具有渐近收敛性，从理论上被证明是一种依概
率1收敛于全局最优解的全局优化算法，此外模拟退火算法具有并行性。
26.11.2. 遗传算法
26.12 网格搜索算法
网格搜索算法（Grid Search Algorithm）将搜索空间划分成网格状，根据各个格点处的目标
函数值，选择出最优的解。在二维搜索空间（两个参数）中，网格搜索算法可以迅速确定最佳格
点，但代价不菲，极易带来维数灾难问题：对于含有多个参数（比如d个）的优化问题，若沿着搜
索空间每个坐标选取1, 000的网格点，则产生的格点总数达到m =1.0 × 10
3d
，相应地需要对目标函
数评估m次。
如果待估的目标函数仅含少许几个参数，并且评估目标函数的代价低廉，那么网格搜索就是最
佳的选择。反之，则可考虑使用网格游走（Grid Walk）或Nelder-Mead 算法。
26.12.1. 网格游走算法
网格游走算法（Grid Walk）同网格搜索算法类似，但它不会搜索整个网格空间，而是从当前
位置出发，根据其周围有限个邻居点上的目标函数值，选择“走下山”（Walk Down Hill）的下一
站。
假设当前位置是x =(x
1
,x
2
,...,x
d
)，每次移动的步长为s，网格游走算法保持其他坐标值不
变，在每个坐标当前位置往前或往后一步的地方评估目标函数：
⎧
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩
f
−
1
= f(x
1
+ s, x
2
,...,x
d
)
f
−
1
= f(x
1
− s, x
2
,...,x
d
)
f
+
i
= f(x
1
,...,x
i−1
,x
i
+ s, x
i+1
,...,x
d
), 2 ≤ i ≤ d − 1
f
−
i
= f(x
1
,...,x
i−1
,x
i
− s, x
i+1
,...,x
d
), 2 ≤ i ≤ d − 1
f
+
d
= f(x
1
,x
2
,...,x
d
+ s)
f
−
d
= f(x
1
,x
2
,...,x
d
− s)
(26.97)
如果当前位置是最佳的，则减小步长（比如减半），再重复上述步骤。如果个邻居格点是最佳的，
$%"&'#(
332
)!*+",$
26.13. 坐标轮换法
!"#
则游走到此位置，继续依此策略进行搜索
➊
。网格游走算法每次仅需评估2d个邻居格点，维数灾难
问题自然可以避免，优化的效率远高于网格搜索算法，在高维参数空间，更是如此。
网格游走算法通过不断缩小搜索范围，逐渐逼近局部最优点。在具体实现时，需要设置一个步
长下限s
0
> 0，当步长缩减至s<s
0
，则可以停止游走，结束搜索过程。为有效反映不同坐标（对
应一个参数）的特征，对于每个坐标可以灵活地选择不同步长。
美中不足的是，网格游走算法易陷入局部最优点，为此，可以与随机搜索算法融合，通过引入
随机因素弥补此一缺憾。
26.13 坐标轮换法
坐标轮换法（Univariate Search）是由D’Esopo[316]在1959年发明的。坐标轮换法每次搜索只
允许一个变量变化，其他变量保持不变，沿坐标方向轮流进行搜索的寻优方法。使用坐标轮换法搜
索目标函数f(x)的最小值点，基本流程可概括如下：
Algorithm 24 坐标轮换法
Input: 初始点x
1
，迭代次数T ，允许误差ϵ > 0
for t =1,...,T do
1. 取y
1
= x
t
，从y
1
出发，依次沿第k(1 ≤ k ≤ n)个坐标轴方向做最优一维搜索：
f(y
k
+ λ
k
e
k
)=min
λ>0
f(y
k
+ λe
k
)
其中，y
k+1
= y
k
+ λ
k
e
k
。
2. 如果∥y
n+1
− y
1
∥≥ϵ，取x
t+1
= y
n+1
，结束本次迭代。
3. 如果∥y
n+1
− y
1
∥ < ϵ，则迭代终止，y
n+1
是f(x)的近似最小值点。
end for
26.14 共轭方向法
1964年，Michael Powell[319]对坐标轮换法进行改进，出了共轭方向法（Conjugate Direc-
tion Method）， 也 称 Powell方法。Powell方法主要包括基本搜索、加速搜索和搜索方向调整三部
分，具体步骤如下：
➊ 著名数学家George P
´
olya于1921年证明了一个有趣的定理：喝醉的酒鬼总能找到回家的路，喝醉的小鸟则可能永远也回不了家。游走
空间的维度越高，回到出发点的概率则越低。事实上，酒鬼在网格状分布的街道中随机游走，始终能够回到出发点，然而小鸟在三维网格
中随机游走，就不那么幸运了，它最终回到出发点的概率只有大约34%。在四维网格中随机游走，最终能回到出发点的概率是19.3%，在
八维空间概率仅为7.3%。
$%"&'#(
333
)!*+",$
搜索与排名 Searching and Ranking
!"#
Algorithm 25 Powell方法
Input: 初始点x
0
，由n个相互线性无关的搜索方向构成的初始方向组{h
0
,h
1
,...,h
n−1
}，终止阈
值ϵ > 0，取k =0
for k =1, 2,...do
1. 基本搜索：令y
0
= x
k
，依次沿着h
0
,h
1
,...,h
n−1
的方向进行一维搜索，得到相应的辅助
迭代点{y
1
,...,y
n
}：
y
i
= y
i−1
+ λ
i−1
h
i−1
其中，λ
i−1
= arg min
λ≥0
f(y
i−1
+ λ
i−1
h
i−1
),i=1,...,n
2. 加速搜索：令h
n
= y
n
−y
0
，若∥h
n
∥≤ϵ，停止迭代，取下个迭代点x
k+1
= y
n
，否则转到
下一步。
3. 搜索方向调整：取∆
m−1
! f(y
m−1
) − f(y
m
) = max
i∈{1,...,n}
f(y
i−1
) − f(y
i
)，如果
f(2y
n
− y
0
) − [2f(y
n
) − f(y
0
)] ≥ 2∆
m−1
成立，则搜索方向无需调整，令x
k+1
= y
n
；否则，转入下一步。
4. 搜索方向组调整：令x
k+1
= y
n
+ λ
n
h
n
，其中
λ
n
= arg min
λ≥0
(y
n
+ λh
n
)
丢弃第m个搜索方向向量：
{h
0
,...,h
m−1
, h
m
,h
m+1
,...,h
n−1
} ← {h
0
,...,h
m−1
,h
m+1
,h
m+2
,...,h
n
}
end for
$%"&'#(
334
)!*+",$
26.15. Nelder-Mead法
!"#
Powell方法是一种高效的零阶优化方法，但事无绝对，在些条件下其表现不尽如人意。当一
个搜索方向不能供进一步的改善，则接下来的搜索方向就不再是共轭的。有时，在迭代几步以
后，搜索方向可能平行。
26.15 Nelder-Mead法
1965年，John Nelder和Roger Mead[320]出一种启发式局部搜索算法– Nelder-Mead法，有
时也称“下山单纯形法”（Downhill Simplex Method）。 一 般 地 ， Nelder-Mead法只需计算一个或
两个点处的函数值，因此广泛应用于函数值计算开销较大的情形。
Nelder-Mead法使用单纯形（见26.2节：Dantzig的单纯形法）的概念优化含有多个参数的目
标函数。在搜索过程中，始终维护一个大小为d +1的顶点集，利用四种基本的变形规则（反射、
扩展、收缩和压缩），灵活调整单纯形搜索方向。Nelder-Mead单纯形算法主要步骤如下：
1. 初始化单纯形S = S
0
：任意选取点x
1
∈ R
n
，根据给定步长s
i
，选取x
i
= x
1
+ s
i
e
i−1
,i =
2,...,n+1，其中e
i
表示第i维坐标的单位坐标基。
2. 排序（Order）：对单纯形各顶点按函数值排序，不失一般性，假设
f(x
1
) ≤ f(x
2
) ≤ ···≤ f(x
n
) ≤ f(x
n+1
)
3. 中心点（Centroid）：计算除f (x
n+1
)以外的其他顶点的中心点
c =
1
n
n
"
i=1
x
i
4. 变形（Transformation）：根据各顶点函数值，改变当前单纯形形状。 首先尝试利用“反射
（Reflection）”、“扩展（Expansion）”、“收缩（Contraction）” 搜 索 一 个 比 x
n
更好的顶点替换
它，搜索的点皆位于由x
n+1
和c确定的直线上。如果搜索成功，将使用新点替换x
n+1
，结束本
次迭代；否则，向顶点x
1
方向“压缩（Shrink, Reduce）”单纯形S。具体做法如下：
• 反射：选择反射点x
r
= c + α(c − x
n+1
)，其中α > 0。如果f(x
1
) ≤ f (x
r
) <f(x
n
)，则将
使用x
r
从S中替换出x
n+1
，并结束本次迭代。
• 扩展：如果f(r) <f(x
1
)，表明此方向利于搜索，于是选择扩展点x
e
= c + γ(c − x
n+1
)，
其中，γ > α。如果f(x
e
) <f(x
r
)，则使用x
e
替换x
n+1
；否则，使用x
r
替换x
n+1
。结束本
次迭代。
• 收缩：如果f (r) ≥ f(x
n
)，表明此方向不利于搜索，转向x
n+1
方向收缩。选择收缩
点x
c
= c + ρ(c − x
n+1
)，其中ρ < 0。如果f(x
c
) <f(x
n+1
)，则使用x
c
替换x
n+1
，并结束
本次迭代；否则，“压缩”单纯形。
$%"&'#(
335
)!*+",$
搜索与排名 Searching and Ranking
!"#
• 压缩：如果利用“反射”、“扩展”、“收缩”搜索失败，表明最佳点x
1
可能临近局部最优
点，可以缩小搜索范围，选择向顶点x
1
方向“压缩”单纯形，更新除x
1
以外的所有顶点
（σ > 0）：
x
i
← σ(x
i
+ x
1
),i=2,...,n+1
更新单纯形S后，结束本次迭代。
Nelder-Mead单纯形算法中使用的参数，一般分别取α =1, γ =2, ρ = −0.5, σ =0.5。算法的
终止条件，可选择在单纯形最短的边低于个阈值。如果搜索空间有界，阻止单纯形法越界搜索最
简单的法子，可以令越界点的函数值最不可能最优，比如最小化问题则越界点处函数值越大越好，
最大化问题则函数值越小越好。
Nelder-Mead算法同网格搜索算法类似，都属于局部搜索方法，但它从多个角度对改进了网格
搜索方法。比如，适应性地调整步长，使用更少的目标函数估计。单纯形算法与网格游走相似，对
于初始点（单纯形）的选择十分敏感，选择的初始点不同，则可能得到不同的局部最优点。对于单
纯形法，如果初始单纯形太小，很容易陷入局部最优。随机搜索算法、网格搜索算法比较擅长寻找
起始点，但收敛速度不容乐观。加州伯克利大学教授Jonathan Shewchuk
➊
出融合两种类型的搜
索算法，实现优势互补：（1）使用网格搜索（粗糙网格）寻找初始点，然后运行网格游走算法，并
将网格搜索算法使用的步长减半。（2）使用网格搜索寻找多个（比如20个）表现好的起始点，然后
将这些点分散到网格中，并选择其中一个点作为起始点，开始单纯形搜索，重复执行多次，优选出
表现最好的点。
26.16 Levenberg-Marquardt算法
Levenberg-Marquardt优化算法，也称阻滞的最小二乘法，由Levenberg [344]和Marquardt
[345]两人分别在1944 年和1963年出。它是使用最广泛的非线性最小二乘算法，主要处理非线性
最优化问题。它介于牛顿法和梯度下降法之间，比牛顿法更健壮但收敛速度要慢，在非线性曲线拟
合、计算机视觉等领域得到广泛应用。
给定一组数据{x
i
,y
i
}
n
i=1
，x
i
∈ R
m
，y
i
∈ R，曲线拟合的目的是使用一个含参函数
f
θ
: R
m
-→ R (26.98)
拟合指定数据。θ ∈ R
k
表示参数向量。函数对数据的拟合效果可以使用下面形式的误差平方和衡
量：
L(θ)=
1
2
n
"
i=1
(f
θ
(x
i
) − y
i
)
2
(26.99)
其中，系数1/2只是方便计算，不影响结果。模型参数的估计是通过最小化误差平方和实现的，此
类方法因此被称作最小二乘法，当f是非线性函数时，称为非线性最小二乘法。
➊ Jonathan Shewchuk：http://www.cs.berkeley.edu/simjrs
$%"&'#(
336
)!*+",$
26.17. 信赖域方法
!"#
在实际应用中，人们有时会使用加权误差平方和（卡方）
χ
2
(θ)=
1
2
n
"
i=1
(
f
θ
(x
i
) − y
i
ω
i
)
2
(26.100)
估计模型参数
θ = arg mi n
θ∈Θ
χ
2
(θ) (26.101)
其中，ω表示误差权值向量，Θ表示参数空间。
由于目标函数是可微的，我们计算它关于参数θ的导数：
∂χ
2
(θ)
∂θ
=
n
"
i=1
(
f
θ
(x
i
) − y
i
ω
2
i
)
∂f
θ
(x
i
)
∂θ
=(ˆy
θ
− y)
T
WJ (26.102)
其中，ˆy
θ
∈ R
n
表示函数估计结果，y ∈ R
n
表示真实数据向量，W 表示权值矩阵，对角线部分
是1/ω
2
i
，J表示雅克比矩阵。
根据梯度下降法，可行的搜索方向为
d
gd
= J
T
W (y − ˆy
θ
) (26.103)
根据Gauss-Newton法，搜索方向d
gn
有下面等式成立：
J
T
WJd
gn
= J
T
W (y − ˆy
θ
) (26.104)
Levenberg-Marquardt算法则介于两者之间，选择的搜索方向d
lm
满足：
(J
T
WJ + λI)d
lm
= J
T
W (y − ˆy
θ
) (26.105)
其中的调节参数λ越大，则越接近于梯度下降法，λ越小，则越接近于Newton法。
Marquardt [345]使用下面的形式更新模型参数：
(J
T
WJ + λ diag(J
T
WJ))d
lm
= J
T
W (y − ˆy
θ
) (26.106)
26.17 信赖域方法
26.18 罚函数法
利用罚函数法（Penalty Function Method）， 可 将 非 线 性 规 划 问 题 的 求 解 ， 转 化 为 求 解 一 系
列无约束极值问题，因而也称序列无约束最小化技术（Sequential Unconstrained Minimization
Technique, SUMT）。罚函数法求解非线性规划问题的思想：利用问题中的约束函数作出适当的罚
函数，由此构造出带参数的增广目标函数，把问题转化为无约束非线性规划问题。罚函数法有两种
形式，一种是外罚函数法，另一种是内罚函数法。
$%"&'#(
337
)!*+",$
搜索与排名 Searching and Ranking
!"#
26.19 障碍函数
障碍函数（Barrier Function）是一个连续函数，当数据点趋近于可行域（Feasible Region）
边界时，函数值增加到无穷。在有约束优化（Constrained Optimization）问题中，障碍函数可以
用作违反约束条件的惩罚项，以防迭代过程搜索路径逃离出可行域。最常见的两种障碍函数是倒数
障碍函数（Inverse Barrier Functions）和对数障碍函数（Logarithmic Barrier Functions）
➊
。
我们首先给出障碍函数的基本数学定义：
定义26.9. 集合K ⊂ R
n
称为锥锥锥，如果对于每一个x ∈ K，都有tx ∈ K，∀t ≥ 0。如果锥K是凸的，
则称它为凸锥。
定义26.10. 假设K是一个凸锥，如果连续函数f : K -→ R
P
{+∞}满足条件
⎧
⎨
⎩
f(x) < +∞,x∈ int(K),
f(x)=+∞,x∈ bdry(K).
我们称它是障障障碍碍碍函函函数数数，其中int(K)表示K的内部（Interior）， bdry(K)表示K的边界（Boundary）。
如果障碍函数形如f(x)=−log φ(x)，并且φ : K -→ R
+
也是一个连续函数，则我们称它对对对数数数障障障碍碍碍
函函函数数数[346]。
对于下面的原始优化问题：
min c
T
x
s.t. Ax ≤ b
(26.107)
其中，A =(a
1
,...,a
m
)
T
，则可以定义如下形式的障碍函数：
f(x)=
⎧
⎪
⎨
⎪
⎩
−
m
(
i=1
log(b
i
− a
T
i
x), Ax < b,
+∞, Ax ≥ b.
(26.108)
将有约束优化问题转化成下面的无约束优化问题：
min c
T
x + λf(x)
26.20 范数逼近
范数逼近（Norm Approximation）在信号处理、压缩感知（Compressed Sensing）
➋
等领域
中有很多应用[347]，[348]对解决范数最优化问题有过概述，主题是将无约束的范数最小化问题转
化为线性规划问题，本节我们做出简要介绍。
➊ 由于对数障碍函数与主对偶内点法（Primal-dual Interior Point Method）存在联系，人们开始加深对数障碍函数的研究；倒数障碍函
数相比而言，计算开销较大，因此，对数障碍函数更为实用
➋ 压缩感知也称压缩采样或稀疏采样，属于一种新的采样理论，2004年由斯坦福大学教授Emmanuel Candes、美国科学院院士David
Donoho、菲尔兹奖获得者陶泽轩等学者出。它能够利用信号的稀疏性特征，随机采样获取信号的离散样本，通过非线性重建算法重建
原始信号，在信息论、图像处理、光学成像、模式识别等领域受到高度关注，并被美国科技评论评委2007年度十大科技进展。
$%"&'#(
338
)!*+",$
