目标函数的几种极值求解方法

目标函数极值求解的几种方法

题目：()()

2

22

1

122min -+-x x

，取初始点()()

T

x 3,11

=，分别用最速下降法，

牛顿法，共轭梯度法编程实现。

一维搜索法：

迭代下降算法大都具有一个共同点，这就是得到点()k x 后需要按某种规则确定一个方向()k d ，再从()k x 出发，沿方向()k d 在直线（或射线）上求目标函数的极小点，从而得到()k x 的后继点()1+k x ，重复以上做法，直至求得问题的解，这里所谓求目标函数在直线上的极小点，称为一维搜索。

一维搜索的方法很多，归纳起来大体可以分为两类，一类是试探法：采用这类方法，需要按某种方式找试探点，通过一系列的试探点来确定极小点。另一类是函数逼近法或插值法：这类方法是用某种较简单的曲线逼近本来的函数曲线，通过求逼近函数的极小点来估计目标函数的极小点。本文采用的是第一类试探法中的黄金分割法。原理书上有详细叙述，在这里介绍一下实现过程：

⑴ 置初始区间[11,b a ]及精度要求L>0,计算试探点1λ和1μ，计算函数值

()1λf 和()1μf ，计算公式是：()1111382.0a b a -+=λ，()1111618.0a b a -+=μ。令

k=1。

⑵ 若L a b k k <-则停止计算。否则，当()K f λ>()k f μ时,转步骤⑶;当

()K f λ≤()k f μ时,转步骤⑷ 。

⑶ 置k k a λ=+1,k k b b =+1,k k μλ=+1,()1111618.0++++-+=k k k k a b a μ,计算函数值

()1+k f μ,转⑸。

⑷ 置k k a a =+1,k k b μ=+1,k k μμ=+1,()1111382.0++++-+=k k k k a b a λ,计算函数值()1+k f λ,转⑸。

⑸ 置k=k+1返回步骤 ⑵。

1. 最速下降法

实现原理描述：在求目标函数极小值问题时,总希望从一点出发,选择一个目

标函数值下降最快的方向,以利于尽快达到极小点,正是基于这样一种愿望提出的最速下降法,并且经过一系列理论推导研究可知,负梯度方向为最速下降方向。

最速下降法的迭代公式是()()()k k k k d x x λ+=+1，其中()k d 是从()k x 出发的搜索方

向，这里取在点()k x 处最速下降方向，即()()k k x f d -?=。k λ是从()k x 出发沿方向()k d 进行的一维搜索步长，满足()()()()()()k k k k k d x f d x f λλλ+=+≥0

min 。

实现步骤如下：

⑴ 给定初点()n R x ∈1 ，允许误差0>ε，置k=1。 ⑵ 计算搜索方向()()k k x f d -?=。

⑶ 若()ε≤k d ，则停止计算；否则，从()k x 出发，沿方向()k d 进行的一维搜索，求k λ，使()()()()()()k k k k k d x f d x f λλλ+=+≥0

min 。

⑷ ()()()k k k k d x x λ+=+1，置k=k+1返回步骤 ⑵。

2. 拟牛顿法

基本思想是用不包括二阶导数的矩阵近似牛顿法中的Hesse 矩阵的逆矩阵，

因构造近似矩阵的方法不同，因而出现了不同的拟牛顿法。

牛顿法迭代公式：()()()k k k k d x x λ+=+1，()k d 是在点()k x 处的牛顿方向，

()()

()

()()

k k k x f x f d ?-?=-1

2，k λ是从()k x 出发沿牛顿方向()k d 进行搜索的最优步长。

用不包括二阶导数的矩阵k H 近似取代牛顿法中的Hesse 矩阵的逆矩阵

()

()

1

2-?k x f ，1+k H 需满足拟牛顿条件。

实现步骤：

⑴ 给定初点()1x ，允许误差0>ε。

⑵ 置n I H =1（单位矩阵），计算出在()1x 处的梯度()()11x f g ?=，置k=1。 ⑶ 令()k k k g H d -=。

⑷ 从()k x 出发沿方向()k d 搜索，求步长k λ，使它满足

()()()()()()

k k k k k d x f d x f λλλ+=+≥0

min ，令()()()k k k k d x x λ+=+1。

⑸ 检验是否满足收敛标准，若()

(

)ε≤+1k y f ，

则停止迭代，得到点()

1+-

=k x x ，

否则进行步骤⑹。

⑹ 若k=n ，令()()11+=k x x ，返回⑵；否则进行步骤⑺。 ⑺

令

()

()

11++?=k k x f g ，

()()()

k k k x x p -=+1，

()k

k k g g q -=+1，

()()()()()()()()k k T

k k

T k k k k T k T k k k k q H q H q q H q p p p H H -+=+1，置k=k+1 。返回⑶。

3. 共轭梯度法

若()()()k d d d ,,,21 是n R 中k 个方向，它们两两关于A 共轭，即满足

()()k j i j i Ad d j T i ,,1,;,0 =≠=，称这组方向为A 的k 个共轭方向。共轭梯度法的

基本思想是把共轭性与最速下降法相结合，利用已知点处的梯度构造一组共轭方向，并沿这组方向进行搜索，求出目标函数的极小点，根据共轭方向的基本性质这种方法具有二次终止性。

实现步骤如下：

⑴ 给定初点()1x ，允许误差0>ε，置 ()()11x y =，()()()11y f d -?=，k=j=1。

⑵ 若()()

ε≤j y f ，则停止计算；否则，作一维搜索，求j λ，满足 ()()()()()()j j j j j d y f d y f λλλ+=+≥0

min ，令()()()j j j j d y y λ+=+1。

⑶ 若n j <，则进行步骤⑷，否则进行步骤⑸

⑷ 令()()()

()j j j j d y f d β+-?=++11，其中()

()

()()

2

2

1j

j j y f y f ??=

+β，置j=j+1，转⑵。

⑸ 令()()11++=n k y x ，()()11+=k x y ，()()()11y f d -?=，置j=1，k=k+1，转⑵ 。

4. 实验结果

1=。迭用以上三种方法通过Matlab编程得到实验数据。初始值()()T

x3,1

代精度sum(abs(x1-x).^2)<1e-4。

实验结果分析：

由上表格可以看到最速下降法需要四次迭代实现所要求的精度，拟牛顿法和共轭梯度法需要三次。

程序：

%精确一维搜索法的子函数,0.618（黄金分割）法,gold.m

%输入的变量x为初始迭代点是二维的向量,d为初始迭代方向是二维的向量

%输出变量是在[0,10]区间上使函数取得极小值点的步长因子

function alfa=gold(x,d)

a=0;b=10;tao=0.618;

lanmda=a+(1-tao)*(b-a);

mu=a+tao*(b-a);alfa=lanmda;%初始化

f=((x(1)+alfa*d(1))-2)^2+2*(x(2)+alfa*d(2)-1)^2;%目标函数

m=f;alfa=mu;n=f;

while 1

if m>n

if abs(lanmda-b)<1e-4

alfa=mu; return

else

a=lanmda; lanmda=mu; m=n;

mu=a+tao*(b-a); alfa=mu;

n=((x(1)+alfa*d(1))-2)^2+2*(x(2)+alfa*d(2)-1)^2;

end

else

if abs(mu-a)<1e-4

alfa=lanmda; return

else

b=mu; mu=lanmda; n=m;

lanmda=a+(1-tao)*(b-a); alfa=lanmda;

m=((x(1)+alfa*d(1))-2)^2+2*(x(2)+alfa*d(2)-1)^2;

end

end

end

%梯度子函数，tidu.m,输入的变量为二维的向量，返回梯度在x处的数值向量function g=tidu(x)

%待求解的函数

f=(x(1)-2)^2+2*(x(2)-1)^2;

%求函数的梯度表达式

g=[2*(x(1)-2) 4*(x(2)-1)];

x1=x(1); x2=x(2);

%最速下降法极小化函数的通用子函数zuisu.m

%输入变量为初始的迭代点，输出变量为极小值点

function x0=zuisu(x)

%判断梯度范数是否满足计算精度1e-4的要求.是,标志变量设为1，输出结果; %否,标志变量设为0

if sum(abs(tidu(x)).^2)<1e-4

flag=1; x0=x;

else

flag=0;

end

%循环求解函数的极小点

while flag==0

d=-tidu(x); a=gold(x,d); x=x+a*d

%判断梯度范数是否满足计算精度的要求.是,标志变量设为1，输出结果;

%否,标志变量设为0，继续迭代

if sum(abs(tidu(x)).^2)<1e-4

flag=1; x0=x;

else

flag=0;

end

End

%拟牛顿法极小化函数的通用子函数，gonge.m

%输入变量为初始的迭代点，输出变量为极小值点

function x0=ninewton(x)

%判断梯度范数是否满足计算精度的要求.是,标志变量设为1，输出结果;

%否,标志变量设为0，继续迭代

if sum(abs(tidu(x)).^2)<1e-4

flag=1； x0=x;

else

flag=0;

end

%初始的H矩阵为单位矩阵

h0=eye(2);

%循环求解函数的极小点

while flag==0

%计算新的迭代方向

d=-h0*tidu(x)'; a=gold(x,d);

x1=(x'+a*h0*d)'; s=x1-x;

y=tidu(x1)-tidu(x); v=s*y';

%校正H矩阵

h0=(eye(2)-s'*y./v)*h0*(eye(2)-y'*s./v)+s'*s./v;

%判断下一次和上一次迭代点之差是否满足计算精度的要求.是,标志变量设为1，输出结果;否,标志变量设为0，继续迭代

if sum(abs(x-x1).^2)<1e-4

flag=1; x0=x;

else

flag=0;

end

x=x1

end

%共轭剃度法极小化函数的通用子函数，gonge.m

%输入变量为初始的迭代点，输出变量为极小值点

function x0=gonge(x)

%判断梯度范数是否满足计算精度的要求.是,标志变量设为1，输出结果;

%否,标志变量设为0，继续迭代

if sum(abs(tidu(x)).^2)<1e-4

flag=1; x0=x;

else

flag=0;

end

%第一次的迭代方法为负梯度方向

d1=-tidu(x);a=gold(x,d1);x1=x+a*d1;

%循环求解函数的极小点

while flag==0

g1=tidu(x); g2=tidu(x1);

%利用FR公式求解系数bata

bata=(g2*g2')/(g1*g1');

d2=-g2+bata*d1;

a=gold(x1,d2);

x=x1;

x1=x+a*d2

%判断下一次和上一次迭代点之差是否满足计算精度的要求.是,标志变量设为1，输出结果;否,标志变量设为0，继续迭代

if sum(abs(x1-x).^2)<1e-4

flag=1; x0=x1;

else

flag=0;

end

end