无约束最优化方法学习笔记

这一段时间学习优化理论中的一些经典方法，写下一些要点和体会，帮助自己和有需要的人理解最优化方法。本文介绍了如下优化方法的原理和计算步骤：

• 最速下降法。
• 牛顿法。
• 共轭梯度法。
• 变尺度法（变度量法）。
• 线性最小二乘法。
• 非线性最小二乘问题的高斯牛顿法和Levenberg-Marquardt法。
• 坐标轮换法。

１.基础知识

首先来看无约束最优化问题：

minf(x)

其中函数 f:Rn→R.求解此问题的方法方法分为两大类：最优条件法和迭代法。

所谓的最优条件法，是指当函数存在解析形式，能够通过最优性条件求解出显式最优解。对于无约束最优化问题，如果f(x)在最优点x¯附近可微，那么x¯是局部极小点的必要条件为：

∇f(x¯)=0

我们常常就是通过这个必要条件去求取可能的极小值点，再验证这些点是否真的是极小值点。当上式方程可以求解的时候，无约束最优化问题基本就解决了。实际中，这个方程往往难以求解。这就引出了第二大类方法：迭代法。

迭代法，也称为“搜索”法，主要思想是通过简单的运算构造点列，逐渐逼近问题的最优解。这里说的“点”是多维空间中的点，也称为“向量”。还有少部分算法通过构造点的集合来逼近问题的最优解，如单纯形调优法。

用于求解无约束最优化问题的方法可以分为解析法和直接法两大类。解析法在构造迭代公式的过程中往往使用了泰勒展开来作近似或者推导，因此迭代步骤中含有梯度∇f(x)或黑塞(Hessian)矩阵∇2f(x)，在问题的解析形态较好的情况下使用往往能获得比较快的收敛速度。而直接法则从物理角度思考如何递推，不会用到梯度或者黑塞矩阵，它对问题的解析形态几乎没有要求，只要能计算出函数值即可。当然，它的收敛速度往往难于保证。

2.解析法

求解无约束最优化问题的解析法主要有：最速下降法、牛顿法、共轭梯度法（DFP法）和变尺度法（变度量法）。对于特殊的最小二乘问题，有最小二乘法。这些方法各有千秋，除了最小二乘法，后面的方法都针对前面方法的某个问题做了改进。这些方法的核心就是研究如何确定每一步迭代的方向和步长。

2.1最速下降法

最速下降法(Gradient Descent Method)的基本思想是：因为连续函数沿着负梯度方向的下降速度是最快的（这一结论由梯度和方向导数的定义可以推出），所以每次迭代我们都从当前点出发，沿着负梯度方向前进一个最优步长，可以期望能较快逼近函数的极值。

最速下降法仅有三个步骤：

1. 设置初始值。设置迭代起点x0∈Rn,允许误差ϵ>0和迭代变量初值k←0。
2. 检查终止条件。如果||∇f(xk)||<ϵ，停止迭代输出xk作为近似最优解；否则转步骤3。
3. 迭代，也就是通过一维搜索求下一个迭代点。取搜索方向为负梯度方向dk=−∇f(x)，求λk使得
   f(xk+λkdk)=minλ≥0f(xk+λdk)
   再令
   xk+1=xk+λdk
   转步骤２。

步骤中隐含了条件：函数f(x)必须可微，也就是说函数f(x)的梯度必须存在。这个步骤描述是最速下降发最抽象的形式。其中最关键的步骤是求解问题

minλ≥0f(xk+λdk)

我们写出它的最优性必要条件

df(xk+λdk)dλ=∇f(xk+λdk)Tdk=0

当f(x)的形式确定，我们可以通过求解这个一元方程来获得迭代步长λ。当此方程形式复杂，解析解不存在，我们就需要使用“一维搜索”来求解λ了。一维搜索是一些数值方法，有0.618法、Fibonacci法、抛物线法等等，这里不详细解释了。

在实际使用中，为了简便，也可以使用一个预定义的常数而不用一维搜索来确定步长λ。这时步长的选择往往根据经验或者通过试算来确定。步长过小则收敛慢，步长过大可能震荡而不收敛。

最速下降法是最基本的迭代优化方法。它最简单，最基础，通常是收敛的。可以证明，最速下降法是一阶收敛的，往往需要多次迭代才能接近问题最优解。这是它的不足。

2.2牛顿法

牛顿法是从函数的二阶泰勒展开式推导而来，其思想就是利用目标函数二阶近似的解去逼近目标函数的解。将函数f在迭代点xk附近作二阶泰勒展开，有

f(x)≈ϕ(x)=f(xk)+∇f(xk)T(x−xk)+12(x−xk)T∇2f(xk)(x−xk)

为了计算ϕ(x)的极小值，令它的梯度为零，即

∇ϕ(x)=0

有

∇f(xk)+∇2f(xk)(x−xk)=0

从而推出

x=xk−[∇2f(xk)]−1∇f(xk)

我们将这里的x作为第k+1次迭代的估值，就有了迭代公式

xk+1=xk−[∇2(fxk)]−1∇f(xk)

牛顿法(Newton Method)利用了函数的二阶信息，即黑塞矩阵，来加速迭代收敛。具有二阶收敛速度是它的显著优势。牛顿法的步骤为：

1. 设置初始值。给定迭代初值x0∈Rn,ϵ>0,令k←0。
2. 检查终止条件。如果||∇f(x)||<ϵ,迭代终止，xk为近似最优解；否则，转步骤３。
3. 迭代计算。取迭代方向
   dk=−[∇2f(xk)]−1∇f(xk)
   令
   xk+1=xk+dkk←k+1
   转步骤２。

牛顿法要求初始点在最优点附近(泰勒展开的前提就是在邻域内)，否则可能不收敛。为了使得牛顿法能够全局收敛，提出了阻尼牛顿法（Damped Newton Method）。阻尼牛顿法的改进在于每次的搜索步长不固定为１，而是通过一维搜索来确定步长。其步骤如下：
1. 设置初始值。给定迭代初值x0∈Rn,ϵ>0,令k←0。
2. 检查终止条件。如果||∇f(x)||<ϵ,迭代终止，xk为近似最优解；否则，转步骤３。
3. 取迭代方向

dk=−[∇2f(xk)]−1∇f(xk)

4. 进行一维搜索确定步长λk使得

f(xk+λkdk)=minλ≥0f(xk+λdk)

令

xk+1=xk+λdkk←k+1

转步骤２。

可以证明，当f(x)具有二阶连续偏导数且∇2f(x)正定，阻尼牛顿法是全局收敛的。

2.3共轭梯度法

最速下降法收敛速度慢，牛顿法虽然收敛速度快但是需要计算黑塞矩阵的逆矩阵，计算复杂度比较高。有的时候我们希望有一种收敛快且不用计算二阶导数的方法，共轭梯度法（Conjugate Gradient Method）应运而生了。当然，计算量减小的代价是理论推导更加繁琐。

要了解共轭梯度法，首先要了解共轭的概念。设有两个n维向量d1、d2和一个n阶正定矩阵Q，如果它们满足

dT1Qd2=0

则称向量d1、d2关于Q共轭。当Q为单位阵时，共轭就成了正交。可将共轭看成是正交概念的推广。共轭的概念还可以推广到含有m个向量的向量组，即

dTiQdj=0,∀i,j=1,...m,i≠j

我们用正定二次函数

f(x)=12xTQx+bTx+c

来推导，最后推广到一般问题。注意这里的Q∈Rn×n为正定矩阵，b∈Rn，c∈R。

我们从x0出发，沿着某个下降方向d0作一维搜索得到下一个迭代点x1，那么有

f(x1)=f(x0+λ0d0)=minλ≥0f(x0+λd0)

它的最优性条件是对λ求导，于是有

∇f(x0+λ0d0)Td0=0

即

∇f(x1)Td0=0

新的迭代方向d1如果取负梯度方向，就成了最速下降法。负梯度方向虽然是当前函数值下降最快的方向，却未必是直指函数最优点的方向。设函数最优点为x¯，我们希望搜索方向d1能直接指向最优值，即

x¯=x1+λ1d1

而作为最优点，x¯应满足

0=∇f(x¯)=Qx¯+b=Q(x1+λ1d1)+b=(Qx1+b)+λ1Qd1=∇f(x1)+λ1Qd1

两边左乘dT0并注意到dT0∇f(x1)=0就有

dT0Qd1=0

上式表明，搜索方向d1与d0是关于Q共轭的。

可以证明，如果能构造出与Q共轭的向量组，至多通过n次迭代就可以求得正定二次函数最优化问题的最优解。由此产生一类称为共轭方向法的方法，当我们利用梯度信息来构造与Q共轭的向量组，就产生了共轭梯度法。

下面推导如何利用梯度构造与Q共轭的向量组。首先初始搜索方向取为负梯度方向，即

d0=−∇f(x0)

新的搜索方向由负梯度方向和上一次搜索方向的线型组合产生，即

dk+1=−∇f(xk+1)+αkdk,k=0,1,...,n−2

其中αk待定。我们利用dk+1与dk关于Q共轭作为约束条件，有

0=dTk+1Qdk=(−∇f(xk+1)T+αkdTk)Qdk=−∇f(xk+1)TQdk+αkdTkQdk

从而解得

αk=∇f(xk+1)TQdkdTkQdk

总结一下，我们得到了n个搜索方向的递推公式如下：

⎧⎩⎨⎪⎪⎪⎪d0dk+1αk=−∇f(x0)=−∇f(xk+1)+αkdk=−∇f(xk+1)TQdkdTkQdk

可以证明，上式得到的dk是一组关于Ｑ共轭的方向向量。

将此公式向一般无约束最优化问题推广，必须消去Q。这里直接给出三个消去Q后αk的表达式：
(1)FR公式（Fletcher-Reeves）

αk=||∇f(xk+1)||2||∇f(xk)||2

(2)DM公式（Dixon-Myers）

αk=−||∇f(xk+1)||2dTk∇f(xk)

(3)PRP公式（Polak-Ribiere-Polyak）

αk=∇f(xk+1)T[∇f(xk+1)−∇f(xk)]||∇f(xk)||2=||∇f(xk+1)||2−∇f(xk+1)T∇f(xk)||∇f(xk)||2

当问题为正定二次函数优化问题时，这三个公式是等价的。对非二次函数最优化问题，它们将产生不同的搜索方向。

共轭梯度法的步骤总结如下：

1. 设置初始值。迭代初始点x0和允许误差ϵ。
2. 检查终止条件。如果||∇f(x0)||<ϵ，迭代终止并输出x0；否则转步骤3。
3. 构造初始搜索方向。
   d0=−∇f(x0)
   　并令k=0
4. 进行一维搜索。求λk使得
   f(xk+λkdk)=minλ≥0f(xk+λdk)
   并令
   xk+1=xk+λkdk
5. 检查终止条件。如果||∇f(xk+1)||<ϵ，迭代终止并输出xk+1；否则转步骤6。
6. 检查迭代次数。如果k+1=n，令
   x0=xn
   转步骤3；否则转步骤7。
7. 构造共轭方向。令
   dk+1αk=−∇f(xk+1)+αkdk=||∇f(xk+1)||2||∇f(xk)||2
   
   这里αk的计算取了FR公式，也可以选取另外两个公式。再令k=k+1，转步骤4。

共轭梯度法在n次迭代后将无法产生新的共轭方向（因为Rn中的共轭方向至多只能有n个），故将第n次迭代点xn作为新的起点，重新进行一轮共轭梯度迭代。共轭梯度法比最速下降法的收敛条件更弱。当f具有一阶连续偏导数且有极值点时，共轭梯度法就是收敛的。对于二次函数最优化，理论上共轭梯度法n次迭代即可收敛。在一定条件下，共轭梯度法具有二次收敛速度。

共轭梯度法的Matlab实现请参考我的另一篇文章：Matlab实现FR共轭梯度法。

2.4变尺度法

最速下降法和阻尼牛顿法的迭代公式可以写成

{xk+1=xk+λkdkdk=−Gk∇f(xk)

当Gk取单位阵时是最速下降法，当Gk=[∇2f(xk)]−1时是阻尼牛顿法。事实上，我们把上式确定的迭代法统称为拟牛顿法(Quasi-Newton Method)。因此改进的思路可以让Gk近似[∇2f(xk)]−1，从而避免计算二阶导数和矩阵求逆，既简化计算又保持快速收敛性。

我们同样利用二阶泰勒展开来考虑如何近似。将f在点xk+1处作二阶泰勒展开，有

f(x)≈f(xk+1)+∇f(xk+1)T(x−xk+1)+12(x−xk+1)T∇2f(xk+1)(x−xk+1)

对上式两边求梯度，有

∇f(x)≈∇f(xk+1)+∇2f(xk+1)(x−xk+1)

令x=xk，有

∇f(xk)≈∇f(xk+1)+∇2f(xk+1)(xk−xk+1)

可得如下近似关系

[∇2f(xk+1)]−1[∇f(xk+1)−∇f(xk)]=xk+1−xk

我们可以令Gk+1满足此关系式，即

Gk+1[∇f(xk+1)−∇f(xk)]=xk+1−xk

从而近似阻尼牛顿法。为了简化后续推导，记

ΔxkΔgkΔGk=xk+1−xk=∇f(xk+1)−∇f(xk)=Gk+1−Gk

则上式可以改写为

ΔGkΔgk=Δxk−GkΔgk

为了求出ΔGk的表达式，我们采用待定系数法，设

ΔGk=ΔxkpTk−GkΔgkqTk

其中pk,qk∈Rn为待定向量。两边右乘Δgk有

ΔGkΔgk=(pTkΔgk)Δxk+(qTkΔgk)GkΔgk

与前面推导的近似式子对比，则

pTkΔgk=1qTkΔgk=1

考虑到Gk近似于黑塞矩阵的逆（为对称矩阵），也应该为对阵矩阵，可构造性地设

pkqk=αkΔxk=βkGkΔgk

解得

αkβk=1ΔxTkΔgk=1ΔgkGTkΔgk

回代可得

ΔGk=ΔxkΔxTkΔxTkΔgk−GkΔgkΔgTkGkΔgTkGkΔgk

上述公式由Davidon提出，被Fletcher和Powell改进得到，称为DFP公式。应用DFP公式的拟牛顿法称为DFP法或者变尺度法（变度量法）。

最后总结一下变尺度法的步骤：

1. 选取初始值。初始点x0∈Rn，允许误差 ϵ>0和G0=In为单位阵。
2. 检查终止条件。如果||∇f(x0)||<ϵ，终止迭代并输出x0；否则转步骤6。
3. 构造初始DFP方向。取d0=−∇f(x0)，令k=0
4. 一维搜索确定步长λk，使得
   f(xk+λkdk)=minλ≥0f(xk+λdk)
   并令
   xk+1=xk+λkdk
5. 检查终止条件。如果||∇f(xk+1)||<ϵ，终止迭代并输出xk+1；否则转步骤6。
6. 检查迭代次数。如果k+1=n，令x0←xn并转步骤(3)；否则转步骤7。
7. 构造DFP方向。令
   Gk+1dk+1=Gk+ΔxkΔxTkΔxTkΔgk−GkΔgkΔgTkGkΔgTkGkΔgk=−Gk+1∇f(xk)
   令k←k+1，转步骤4。

这里每n次迭代就重新初始化可以减少计算误差的积累，有利于快速收敛。

共轭梯度法与牛顿法类似，具有二阶收敛速度。

2.5最小二乘法

若干的目标函数的平方和作为目标函数的优化问题就称为最小二乘问题。这是一类特殊的无约束最优化问题，如下式：

mins(x)=∑i=1mf2i(x)

其中x∈Rn且m≥n。

首先考虑线性最小二乘问题。当fi(x)均为线性函数时，最小二乘问题可以写成矩阵形式：

mins(x)=||Ax−b||2

其中A∈Rm×n，b∈Rn。一般要求A为列满秩矩阵(否则一定有冗余的列，去掉之)。令

0=∇s(x)=2AT(Ax−b)=2(ATAx−ATb)

注意到当A列满秩时，ATA正定且可逆，故可推出最优解

x¯=(ATA)−1ATb

这就是线性最小二乘问题的最优解。其中(ATA)−1AT称为矩阵A的伪逆（广义逆）或者摩尔彭若斯伪逆(Moore–Penrose pseudoinverse)，记为A+。故最优解也可以记为

x¯=A+b

接下来考虑非线性最小二乘问题

mins(x)=∑i=1mf2i(x)

解决问题的思想与牛顿法相似，就是用目标函数的泰勒展开的最优解去逼近目标函数的最优解。由此产生了下面要介绍的高斯－牛顿法。将fi(x)在点xk处做一阶泰勒展开，有

fi(x)≈fi(xk)+∇fi(xk)T(x−xk)

则

s(x)≈∑i=1m[fi(xk)+∇fi(xk)T(x−xk)]2

这是一个线性最小二乘问题。令

A=(∇f1(xk),...,∇fm(xk))b=−(f1(xk),...,fm(xk))T

如果A为列满秩矩阵，则

x−xk=(ATA)−1ATb

如果记函数的向量

f(x)=(f1(x),...,fm(x))T

则f(x)的雅克比矩阵为

∇f(x)=(∇f1(x),...,∇fm(x))

则上面推导的公式变为

x−xk=−(∇f(xk)T∇f(xk))−1∇f(xk)Tf(xk)

将x作为下一次迭代点，有

xk+1=xk−(∇f(xk)T∇f(xk))−1∇f(xk)Tf(xk)

上述公式称为高斯－牛顿迭代公式。可以证明，当迭代初始点接近最优点时，高斯－牛顿法是二阶收敛的。为了解决收敛性问题，模仿阻尼牛顿法，我们可以用以为搜索来确定迭代的步长，从而得到阻尼高斯－牛顿法。它的步骤总结如下：

1. 选取初始数据。x0，允许误差ϵ>0，令k=0
2. 检查终止条件。如果||∇f(xk)Tf(xk)||<ϵ，迭代终止并输出xk；否则转步骤3。
3. 计算高斯－牛顿方向
   dk=−(∇f(xk)T∇f(xk))−1∇f(xk)Tf(xk)
4. 进行一维搜索求λk，使得
   s(xk+λkdk)=minλ≥0s(xk+λdk)
   并令
   xk+1=xk+λkdkk←k+1
   转步骤2。

不要忘记高斯－牛顿法的推导过程中要求雅克比矩阵∇f(xk)为列满秩矩阵，这在实际计算中可能无法保证。为此，Levenberg和Marquardt对高斯－牛顿方向做了改进，即令

dk=−(∇f(xk)T∇f(xk)+αkIn)−1∇f(xk)Tf(xk)

其中αk>0为参数，In为n阶单位矩阵。迭代步长仍为1不变。这种方法被称为Levenberg-Marquart法或最小二乘LM法。

3.直接法

直接法是一类不需要函数导数信息的方法，通过试算来寻找优化方向。它的优势在于对函数没有可微假设。劣势是收敛速度往往很慢。

3.1坐标轮换法

坐标轮换法(univariate search technique)的思想极为朴素：从某个初始点x0出发，沿着各个坐标轴依次做一维搜索得xk，直到xk不再显著变化。

坐标轮换法的步骤描述如下：

1. 选取初始点。设定x0，初始误差ϵ>0，令k←0。
2. 进行一维搜索。从xk−1出发，沿坐标轴方向ek进行一维搜索，求λk−1和xk使得
   f(xk−1+λk−1ek)xk=minλf(xk−1+λek)=xk−1+λk−1ek
3. 检查迭代次数。若k=n，转步骤4；否则令k←k+1，返回步骤2。
4. 检查终止准则。若||xn−x0||<ϵ，终止迭代输出xn；否则令x0=xn,k←1，返回步骤2。

可以证明，如果每次一维搜索都有唯一最优解，则坐标轮换法是收敛的。

3.2 其他方法

直接法的其他方法还有：步长加速法(模式搜索法，pattern search method)，Rosenbrock旋转方向法，Powell法，单纯形调优法等等，这里不介绍了。

4.结语

数学上有一个原则：把复杂问题转化为简单问题。例如：约束优化通过拉格朗日形式转化为无约束优化。将非线性问题通过近似和残差约束转化为线性问题。

无约束最优化问题是优化领域研究的基础性问题。基础性方法，如最速下降法、牛顿法、共轭梯度法，是其他复杂方法的基础。近年来学术和工业界更加常用的交替方向乘子法(ADMM)和原始对偶法(Primal Dual)也是在基础问题和方法上推演而来。因而掌握经典方法的思想对于理解新方法大有裨益。