最优化方法 课程设计报告 运用DFP算法解决无约束最优化问题

实验二、 无约束最优化【实验目的】1．了解无约束最优化方法的一些基本概念。

2．熟悉掌握用相关的命令来求解无约束最优化问题。

【实验内容】把题目和相应的完整命令写在实验报告上。

1：无约束最优化问题实际上是什么问题？求这类问题的最优解的基本思路是什么？2：求()5x f x e x =-在区间[1,2]内的极小值点和极小值。

3：已知22212312123(,,)3sin f x x x x x x x x =+-。

(1) 求123(,,)f x x x 在点(1,1,0)-附近的极小值；(2) 求123(,,)f x x x 在点(1,1,0)-附近的极小值点和极小值，要求优化算法用高斯－牛顿法，搜索方向用拟牛顿法的DFP 公式。

【相关知识说明】无约束最优化是指在没有约束条件下，求多变量实值函数极值。

无约束最优化问题的数学表达式为12min (),(,,,)n n f x x x x x R =∈。

一般f 为非线性函数，x 是n 维实变量，实际上这是一个多元函数无条件极值问题。

由于求极大值问题，可以用添加负号的方式转化为求极小值问题，因此通常只讨论求极小值问题。

应该注意的是，极值问题的解，即极值点，都是局部最优解，全局最优解只能从局部最优解的比较中得到。

如何求解无约束最优化问题的最优解呢？一般是采用迭代法，即先选择一个初始点，再寻找该点处的下降方向（我们称为搜索方向），在该方向上求极小点，得到一个新的点，然后在新点处再寻找下降方向和在该方向上的求极小点，……，如此下去，最终得到最优解。

我们先来看求一元函数y=f(x)在[x1,x2]内的极小值的命令：说明：其中'fun'是函数f(x)的表达式，当然也可以是关于f(x)的函数M-文件名。

返回值x 是极小值点。

现在我们来回答问题1。

问题1：求()2sin x f x e x -=在区间[0,6]内的极小值点和极小值.命令如下f='2*exp(-x)*sin(x)';x=fminbnd(f,0,6) %极小值点fval=2*exp(-x)*sin(x) %对应x 的极小值大家得到的结果是什么呢？这些是一元函数求极值，那么怎么求多元函数的极值呢？可以用下面的最简形式的命令：如果还必须满足更苛刻的要求，可以用下面的命令说明：(1) 返回值中，x 是极小值点。

问题2运用两种不同的方法求解非线性最小二乘问题 42min ||()||x RF x ，其中问题分析和解决框架：上述问题属于无约束最优化问题，具体又为其中应用很广的一类问题，即非线性最小二乘问题，可看作无约束极小化的特殊情形。

我们既可以采用一般的无约束最优化问题的典型算法，如牛顿法，拟牛顿法，共轭梯度法求解，也可采用根据目标函数的特殊结构，对一般的无约束最优化问题进行改造得到的一些针对非线性最小二乘问题的特殊方法，如高斯-牛顿法，Levenberg-Marquardt 法，最小非线性二乘的拟牛顿法等。

本次试验我采用了一般无约束最优化问题的共轭梯度法，拟牛顿法（DFP ）和针对非线性最小二乘的高斯牛顿法。

由于刚开始用的共轭梯度法和拟牛顿法（DFP ），没有采用解决非线性最小二乘的典型算法，觉得有些缺憾，故加上高斯-牛顿法，共采用三种方法，以求对具体算法有更进一步的认识。

由于要求的停机准则和求解一般无约束最优化问题的共轭梯度法，拟牛顿法（DFP ）有些差异，我对停机准则作了些修改，采用规定的表达式来确定停机与否。

对一维搜索求步长，采用的进退法(确定搜索步长区间)和黄金分割法（0.618法）。

具体程序和运行结果：共轭梯度法：共轭梯度法（Conjugate Gradient method ）是介于最速下降法与牛顿法之间的一种方法，是典型的共轭方向法，它的每一个搜索方向是互相共轭的，而这些搜索方向d 仅仅是负梯度方向与上一次迭代搜索方向的组合，因此，存储量少，计算方便。

它仅需利用一阶导数信息，克服了最速下降法收敛慢（“之”字形搜索）的缺点，又避免了牛顿法需要存储和计算Hesse 矩阵并求逆的缺点。

而且共轭梯度法不需要矩阵存储，且有较快的收敛速度和二次终止性等优点。

程序：%% conjugate gradient methodclc;clear;format long;error=1e-5; %停机门限syms x1 x2 x3 x4 r;P=(x1+10*x2)^2+(sqrt(5)*(x3-x4))^2+((x2-2*x3)^2)^2+(sqrt(10)*(x1-x4)^2)^2;f = [ x1+10*x2; 5^(1/2)*(x3-x4); (x2-2*x3)^2; 10^(1/2)*(x1-x4)^2]; v=[x1 x2 x3 x4];x=[3;-1;0;1];j=jacobian(f,v); %求jacobian行列式g=jacobian(P,v); %求目标函数梯度向量g=g';J=subs(j,v,x); %初值带入表达式F=subs(f,v,x);G=subs(g,v,x);k=0;ticwhile (sum((J'*F).^2))^(1/2)>error %判断停机与否if k==0d=-G; %初始搜索方向为最速下降方向，即负梯度方向elsebeta=G_new'*G_new/(G_old'*G_old);d=-G_new+beta*d; %从x（k）出发搜索方向介于从x（k-1）出发搜索方向和点x（k）负梯度方向之间endy=subs(P,v,x+r*d);interval=jintuifa(y,r);step=gold(y,r,interval); %一维搜索确定步长G_old=subs(g,v,x); %前一点x（k-1）梯度x=x+step*d;G_new=subs(g,v,x); %新一点x（k）梯度J=subs(j,v,x); %新的迭代点数值带入表达式F=subs(f,v,x);k=k+1; %迭代次数加1end;tocdisp('Conjugate gradient method');kxsigma=(sum((J'*F).^2))^(1/2)F=(sum(F.^2))^(1/2) %显示迭代次数，变量取值，停机表达式值，目标函数值%%运行结果：由运行结果可知，对非线性最小二乘问题，共轭梯度法由于没有包含二阶导数信息，所需迭代次数多，收敛慢，这是一个主要缺点。

北方民族大学课程设计报告系（部、中心）信息与计算科学学院专业信息与计算科学班级 09信计（3）班小组成员课程名称最优化方法设计题目名称运用DFP算法解决无约束最优化问题提交时间2012年6月26日成绩指导教师变尺度法是在牛顿法的基础上发展起来的，它和梯度法亦有密切关系.变尺度法避免了Newton法在每次迭代都要计算目标函数的Hesse矩阵和它的逆矩阵而导致随问题的维数增加计算量迅速增加.DFP算法是变尺度法中一个非常好的算法.DFP算法首先是1959年由Davidon提出的后经Fletcher和Powell改进，故名之为DFP算法，它也是求解无约束优化问题最有效的算法之一.DFP变尺度法综合了梯度法、牛顿法的优点而又避弃它们各自的缺点，只需计算一阶偏导数，无需计算二阶偏导数及其逆矩阵，对目标函数的初始点选择均无严格要求，收敛速度快.本文主要分析DFP算法原理及运用Matalb软件编程解决实际数学问题.最后运算结果符合计算精度且只用了一次迭代，由此可见收敛速度快.关键词：Newton法变尺度法Hesse矩阵Matlab软件一、课程设计目的 (4)二、课程设计要求 (4)三、课程设计原理 (4)（1）变尺度法基本原理 (4)（2）DFP算法 (5)四、实验内容 (6)五、数学建模及求解 (7)1.DFP算法迭代步骤 (7)2.DFP算法的流程图 (7)六、程序实现 (8)七、数值实验的结果与分析 (11)八、实验总结与体会 (12)1.DFP公式恒有确切解 (12)2.DFP算法的稳定性 (12)参考文献 (13)一、 课程设计目的：1、掌握无约束优化问题DFP 算法的数值求解思路；2、训练分析DFP 算法的运算存储量及收敛速度的能力，了解算法的优缺点；3、通过运用DFP 算法求解实际无约束优化问题的意义；4、熟悉应用Matlab 求解无约束最优化问题的编程方法.二、 课程设计要求熟悉了解DFP 算法原理及求解无约束优化问题的步骤，并运用Matlab 件编程实现求解问题.三、 课程设计原理（1）变尺度法基本原理在Newton 法中，基本迭代公式k k k k P t X X +=+1，其中，1=k t ，)()]([12k k k X f X f P ∇∇-=-，于是有2,1,0,11=-=-+k g G X X k k k k (1)其中0X 是初始点，k g 和k G 分别是目标函数)(X f 在点k X 的梯度和Hesse 矩阵.为了消除这个迭代公式中的Hesse 逆矩阵1-k G ，可用某种近似矩阵)(k k X H H =来替换它，即构造一个矩阵序列}{k H 去逼近Hesse 逆矩阵序列}{1-k G此时式（1）变为k k k k g H X X -=+1事实上，式中k k k g H P -=无非是确定了第k 次迭代的搜索方向，为了取得更大的灵活性，我们考虑更一般的的迭代公式k k k k k g H t X X -=+1 (2)其中步长因子k t 通过从k X 出发沿k k k g H P -=作直线搜索来确定.式（2）是代表很长的一类迭代公式.例如，当I H k ≡（单位矩阵）时，它变为最速下降法的迭代公式.为使kH确实与1-k G 近似并且有容易计算的特点，必须对k H 附加某些条件：第一，为保证迭代公式具有下降性质，要求}{k H 中的每一个矩阵都是对称正定的.理由是，为使搜索方向k k k g H P -=是下降方向，只要0<-=-k k T k k T k g H g P g成立即可，即0>k k T k g H g成立.当k H 对称正定时，此公式必然成立，从而保证式（2）具有下降性质.第二，要求k H 之间的迭代具有简单形式.显然，k k k E H H +=+1 (3)是最简单的形式了.其中k E 称为校正矩阵，式（3）称为校正公式.第三，必须满足拟Newton 条件.即：)()(111k k k k k X X g g H -=-+++ (4)为了书写方便也记k k k g g y -=+1k k k X X S -=+1于是拟Newton 条件可写为k k k S y H =+1 (5)有式（3）和（5）知，k E 必须满足k k k k S y E H =+)(或k k k k k y H S y E == (6)（2）DFP 算法DFP 校正是第一个拟牛顿校正是1959年由Davidon 提出的后经Fletcher 和Powell 改进故名之为DFP 算法它也是求解无约束优化问题最有效的算法之一.DFP 算法基本原理考虑如下形式的校正公式T k k k T k k k k k V V U U H H βα++=+1 (7)其中k k V U ,是特定n 维向量，k k βα,是待定常数.这时，校正矩阵是T k k k T k k k k V V U U E βα+=.现在来确定k E .根据拟Newton 条件，k E 必须满足（6），于是有k k k k T k k k T k k k y H S y V V U U -=+)(βα或k k k k T k k k k T k k k y H S y V V y U U -=+βα.满足这个方程的待定向量k U 和k V 有无穷多种取法，下面是其中的一种：k k T k k k S y U U =α，k k k T k k k y H y V V -=β注意到k T ky U 和k T k y V 都是数量，不妨取 k k S U =，k k k y H V =，同时定出k T k k y S 1=α，kk T k k y H y 1-=β. 将这两式代回（5.32）得 kk T k k T k k k k T k T k k k k y H y H y y H y S S S H H -+=+1. （8） 这就是DFP 校正公式.四、 实验内容采用课本P102页例5.3和P107页例5.4进行数值计算；1，求22212125),(m in x x x x f +=，取初始点T X ]2,2[0=.2，求2221214),(m in x x x x f +=，取初始点T X ]1,1[0=.五、 数学建模及求解1.DFP 算法迭代步骤在拟Newton 算法中，只要把第五步改为计算式（8）而其他不变，该算法就是DFP 算法了.但是由于计算中舍去误差的影响，特别是直线搜索不精确的影响，可能要破坏迭代矩阵k H 的正定性，从而导致算法失效.为保证k H 的正定性，采取以下重置措施：迭代1+n 次后，重置初始点和迭代矩阵，即I H X X n ==+010,以后重新迭代.已知目标函数)(X f 及其梯度)(X g ，问题的维数n ，终止限ε.（1） 选定初始点.计算)(00X f f =，)(00X g g =.（2） 置I H =0，00g P -=，0=k .（3） 作直线搜索),(1k k k P X ls X =+；计算)(11++=k k X f f ，)(11++=k k X g g .（4） 判别终止准则是否满足：若满足，则打印),(11++k k f X ，结束；否转（5）.（5） 若n k =，则置10+=k X X ，10+=k f f ，10+=k g g ，转（2）；否则转（6）.（6） 计算k k k X X S -=+1，k k k g g y -=+1，kk T k k T k k k k T k T k k k k y H y H y y H y S S S H H -+=+1， 111+++-=k k k g H P ，置1+=k k ，转（3）.2.DFP 算法的流程图六、程序实现七、 数值实验的结果与分析由上述运行结果可得出：第一题迭代一次就解得：]0664.0,2220.0[*0150.1---=e X 与精确解]0,0[=X 误差远小于610-=ε，符合要求.第二题同样迭代一次就解得：]0555.0,1110.0[*0150.1--=e X 与精确解]0,0[=X误差远小于610-=ε，符合要求.且所计算的k H 矩阵和梯度与精确计算所得一样，这也表明DFP 算法的matalb 程序完全符合要求.八、 实验总结与体会DFP 变尺度法综合了梯度法、牛顿法的优点而又避弃它们各自的缺点，只需计算一阶偏导数，无需计算二阶偏导数及其逆矩阵，对目标函数的初始点选择均无严格要求，收敛速度快，这些良好的性能已作阐述。

第七章 无约束最优化的解析法本章主要内容：最速下降法及其收敛性与收敛速度 Newton 切线法及其收敛性与收敛速度 阻尼Newton 法 共轭梯度法及其收敛性 变度量法、最小二乘法教学目的及要求：掌握最速下降法并理解其收敛性与收敛速度，掌握Newton 切线法并理解其收敛性与收敛速度，了解阻尼Newton 法；掌握共轭梯度法并理解其收敛性；了解变度量法、最小二乘法。

教学重点：最速下降法．教学难点：变度量法．教学方法：启发式．教学手段：多媒体演示、演讲与板书相结合．教学时间：6学时．教学内容：§7.1 最速下降法考虑无约束最优化问题m i n ()f x ， （7.1.1）其中:n f R R →具有一阶连续偏导数．算法7-1（最速下降法）Step1 选取初始数据．选取初始点0x ，给定允许误差0ε>，令0k =． Step2 检查是否满足终止准则．计算()k f x ∇，若()k f x ε∇<，迭代终止，k x 为问题（7.1.1）的近似最优解；否则，转Step3．Step3 进行一维搜索．取()k k d f x =-∇，求k λ和1k x +，使得()min ()k k k k k f x d f x d λλλ≥+=+， 1k k k k x x d λ+=+．令:1k k =+，返回Step2．特别地，考虑1m i n ()2T T f x x Qx b x c =++， （7.1.2） 其中,n n n x R Q R ⨯∈∈为正定矩阵，,n b R c R ∈∈．设第k 次迭代点为k x ，从点k x 出发沿()k f x -∇作一维搜索，得1()k k k k x x f x λ+=-∇，其中k λ为最优步长．根据定理 6.1.1，有1()()0T k k f x f x +∇∇=．而(),n f x Q x b x R∇=+∀∈， 所以1()()()k k k k f x f x Q f x λ+∇=∇-∇，从而(()())()0T k k k k f x Q f x f x λ∇-∇∇=，而Q 正定，即()()0T k k f x Q f x ∇∇>，故由上式解出()()()()T k k k T k k f x f x f x Q f x λ∇∇=∇∇， （7.1.3） 于是1()()()()()T k k k k k T k k f x f x x x f x f x Q f x +∇∇=-∇∇∇， （7.1.4） 这是最速下降法用于问题（7.1.2）的迭代公式．例1 用最速下降法求解问题2212min ()4f x x x =+， （7.1.5）其中12(,)T x x x =．取初始点(0)(1,1)T x =，允许误差0.1ε=．解 问题（7.1.5）中的f 是正定二次函数，且800,,0020Q b c ⎛⎫⎛⎫=== ⎪ ⎪⎝⎭⎝⎭．f 在点12(,)T x x x =处的梯度12()(8,2)T f x x x ∇=．第一次迭代：令搜索方向(0)(0)()(8,2)T d f x =-∇=--，(0)d ε==>，从点(0)x 出发沿(0)d 作一维搜索，由（7.1.3）式和（7.1.4）式有0680.130769520λ==， (1)(1,1)0.130769(8,2)(0.046152,0.738462)T T T x =+--=-．第二次迭代：令(1)(1)()(0.369216, 1.476924)T d f x =-∇=-，(1) 1.522375d ε==>，从点(1)x 出发沿(1)d 作一维搜索，按（7.1.4）式得(2)(0.101537,0.147682)T x =．第三次迭代：令(2)(2)()(0.812296,0.295364)T d f x =-∇=--，(2)0.864329d ε==>，按（7.1.4）式求得(3)(0.009747,0.107217)T x =-．第四次迭代：令(3)(3)()(0.077976,0.214434)T d f x =-∇=-，(3)0.228171d ε==>，按（7.1.4）式求得(4)(0.019126,0.027816)T x =．第五次迭代：令(4)(4)()(0.153008,0.055632)T d f x =-∇=--，(4)0.162807d ε==>，按（7.1.4）式求得(5)(0.001835,0.020195)T x =-．此时，(5)()f x ε∇=<，已满足精度要求，故得问题（7.1.5）的近似最优解(5)(0.001835,0.020195)T x =-．实际上问题（7.1.5）的最优解为(0,0)T x =．定理7.1.1 设:n f R R →具有一阶连续偏导数，0n x R ∈，记0()f x α=，假定水平集(,)S f α有界，令{}k x 是由最速下降法求解问题（7.1.1）产生的点列，则（1）当{}k x 是有穷点列时，其最后一个点是f 的平稳点；（2）当{}k x 是无穷点列时，它必有极限点，并且任一极限点都是f 的平稳点．定理7.1.2 设:n f R R →具有二阶连续偏导数，由最速下降法解问题（7.1.1）产生的点列{}k x 收敛于x ．若存在0ε>和0M m >>，使得当x x ε-<时，有222(),T n m y y f x y M y y R ≤∇≤∀∈， （7.1.7） 则{}k x 线性收敛于x ．§7.2 Newton 法21()()()()()()()()2T T k k k k k k f x x f x f x x x x x f x x x ϕ≈=+∇-+-∇-． 令 ()0x ϕ∇=，即2()()()0k k k f x f x x x ∇+∇-=，解之得211[()]()k k k k x x f x f x -+=-∇∇．算法7-2（Newton 法）Step1 选取初始数据．选取初始点0x ，给定允许误差0ε>，令0k =． Step2 检查是否满足终止准则．计算()k f x ∇，若()k f x ε∇<，迭代终止，k x 为问题（7.1.1）的近似最优解；否则，转Step3．Step3 构造Newton 方向．计算21[()]k f x -∇，取21[()]()k k k d f x f x -=-∇∇． Step4 求下一个迭代点．令1k k k x x d +=+，:1k k =+，返回Step2．例2 用Newton 法求解问题（7.1.5），仍取初始点(0)(1,1)T x =，允许误差0.1ε=．解 (0)()(8,2)T f x ∇=，2(0)80()02f x ⎛⎫∇= ⎪⎝⎭，故2(0)11/80[()]01/2f x -⎛⎫∇= ⎪⎝⎭，(0)2(0)1(0)1[()]()1d f x f x -⎛⎫=-∇∇=- ⎪⎝⎭，(1)(0)(0)(1,1)(1,1)(0,0)T T T x x d =+=-=．由于 (1)()00.1f x ∇=<，迭代结束，得(1)x 为问题（7.1.5）的最优解．定理7.2.1 设:n f R R →具有三阶连续偏导数，,()0n x R f x ∈∇=，若存在0ε>和0m >，使得当x x ε-≤时，有22(),T n m y y f x y y R ≤∇∀∈， （7.2.2） 则当初始点0x 充分接近x 时，由Newton 法解问题（7.1.1）产生的点列{}k x 收敛于x ，并有二阶收敛速度．算法7-3（阻尼Newton 法）Step1 选取初始数据．选取初始点0x ，给定允许误差0ε>，令0k =． Step2 检查是否满足终止准则．计算()k f x ∇，若()k f x ε∇<，迭代终止，k x 为问题（7.1.1）的近似最优解；否则，转Step3．Step3 构造Newton 方向．计算21[()]k f x -∇，取21[()]()k k k d f x f x -=-∇∇． Step4 进行一维搜索．求k λ和1k x +，使得()min ()k k k k k f x d f x d λλλ≥+=+， 1k k k k x x d λ+=+．令:1k k =+，返回Step2．例3 用阻尼Newton 法求解下面问题：222121min ()(1)2()f x x x x =-+-，（7.2.6） 其中12(,)T x x x =．取初始点(0)(0,0)T x =，允许误差0.1ε=．解 第一次迭代：212112212(1)8()()4()x x x x f x x x ⎛⎫----∇= ⎪-⎝⎭， 22212111168()28()84x x x x f x x ⎛⎫--+-∇= ⎪-⎝⎭， 故 (0)(0)()(2,0),()2T f x f x ε∇=-∇=>，2(0)2(0)1201/20(),[()]0401/4f x f x -⎛⎫⎛⎫∇=∇= ⎪ ⎪⎝⎭⎝⎭．于是，Newton 方向 (0)2(0)1(0)[()]()(1,0)T d f x f x -=-∇∇=，从(0)x 出发沿(0)d 作一维搜索，即求(0)(0)2400min ()min(1)2f x d λλλλλ≥≥+=-+的最优解，得到012λ=．令 (1)(0)(0)0(1/2,0)T x x d λ=+=．(1)(1)()(0,1),()1T f x f x ε∇=-∇=>．第二次迭代：2(1)2(1)184111(),[()]48124f x f x --⎛⎫⎛⎫∇=∇= ⎪ ⎪-⎝⎭⎝⎭， (1)2(1)1(1)[()]()(1/4,1/2)T d f x f x -=-∇∇=．从(1)x 出发沿(1)d 作一维搜索，即求(1)(1)24001min ()min [8(2)(2)]128f x d λλλλλ≥≥+=-+- 的最优解，得到12λ=．令(2)(1)(1)1(1,1)T x x d λ=+=．(1)(1)()(0,1),()1T f x f x ε∇=-∇=>．此时，(2)(2)()(0,0),()0T f x f x ε∇=∇=<，得问题（7.2.6）的最优解为(2)(1,1)T x =，这是惟一的最优解．定理7.2.2 设:n f R R →具有二阶连续偏导数，0n x R ∈，记0()f x α=，假定水平集(,)S f α有界，并且对一切2,()n x R f x ∈∇正定．若{}k x 是由阻尼Newton法求解问题（7.1.1）产生的点列，则（1）当{}k x 是有穷点列时，其最后一个点是f 的唯一极小点；（2）当{}k x 是无穷点列时，它必收敛于f 的唯一极小点．§7.3 共轭梯度法最速下降法和Newton 法是最基本的无约束最优化方法，它们的特性各异：前者计算量较小而收敛速度慢；后者虽然收敛速度快，但需要计算目标函数的Hesse 矩阵及其逆矩阵，故计算量大．本节介绍一类无需计算二阶导数并且收敛速度快的方法．定义 设n n Q R ⨯∈为正定矩阵．若n R 中的向量组011,,,m d d d - 满足0,,0,1,,1,T i j d Qd i j m i j =∀=-≠ ，则称011,,,m d d d - 是Q 共轭的．定理7.3.1 设n n Q R ⨯∈是正定矩阵，n R 中非零向量组011,,,m d d d - 是Q 共轭的，则这m 个向量线性无关．定理7.3.2 设n p R ∈，011,,,n d d d - 是n R 中线性无关的向量组，若p 与每个i d 都正交，则0p =．考虑正定二次函数的无约束最优化问题1min ()2T T f x x Qx b x c =++， （7.3.3） 其中n n Q R ⨯∈为正定矩阵，,n b R c R ∈∈．定理7.3.3 设n n Q R ⨯∈为正定矩阵，011,,,n d d d - 是n R 中一组Q 共轭的非零向量．对于问题（7.3.3），若从任意点0n x R ∈出发依次沿011,,,n d d d - 进行一维搜索，则至多经过n 次迭代可得问题（7.3.3）的最优解．算法7-4（共轭方向法）给定一个正定矩阵n n Q R ⨯∈．Step1 选取初始数据．选取初始点0x ，给定允许误差0ε>． Step2 选取初始搜索方向．计算0()f x ∇，求出0d ，使00()0T f x d ∇<，令0k =． Step3 检查是否满足终止准则．若()k f x ε∇<，迭代终止；否则，转Step4． Step4 进行一维搜索．求k λ和1k x +，使得()min ()k k k k k f x d f x d λλλ≥+=+， 1k k k k x x d λ+=+．Step5 选取搜索方向．求1k d +使10,0,1,,T k j d Qd j k +== ，令:1k k =+，返回Step3．如果用共轭方向法求解正定二次函数的无约束最优化问题1min ()2T T f x x Qx b x c =++， （7.3.3） 其中n n Q R ⨯∈为正定矩阵，,n b R c R ∈∈（此时算法中的正定矩阵应与二次函数的正定矩阵一致），那么容易推出迭代公式为1()T k k k k k T k kf x d x x d d Qd +∇=-． （7.3.10） 对于求解问题（7.1.1），我们还有如下一些方法．Fletcher-Reeves （FR ）公式：212()()k k k f x f x α+∇=∇；Dixon-Myers （DM ）公式：11()()()T k k k T k k f x f x d f x α++∇∇=-∇； Polak-Ribiere-Polyak （PRP ）公式：11()[()()]()()T k k k k T k k f x f x f x f x f x α++∇∇-∇=∇∇． 算法7-5（FR 共轭梯度法）Step1 选取初始数据．选取初始点0x ，给定允许误差0ε>． Step2 检查是否满足终止准则．计算0()f x ∇，若0()f x ε∇<，迭代终止，0x为（7.1.1）的近似最优解；否则，转Step3．Step3 构造初始搜索方向．计算00(),0d f x k =-∇=．Step4 进行一维搜索．求k λ和1k x +，使得()min ()k k k k k f x d f x d λλλ≥+=+， 1k k k k x x d λ+=+．Step5 检查是否满足终止准则．计算1()k f x +∇，若1()k f x ε+∇<，迭代终止，1k x +为（7.1.1）的近似最优解；否则，转Step6．Step6 检查迭代次数．若1k n +=，令0:n x x =，返回Step3；否则，转Step7．Step7 构造共轭方向．用FR 公式取11()k k k k d f x d α++=-∇+，212()()k k k f x f x α+∇=∇，令:1k k =+，返回Step4． 注意，如果算法7-4的Step7中k α的形式改为DM 公式或PRP 公式，则分别得到DM 共轭梯度法和PRP 共轭梯度法．定理7.3.5 设:n f R R →具有一阶连续偏导数，0n x R ∈，记0()f x α=，并假设水平集(,)S f α有界．若{}k x 是由共轭梯度法（包括任何一种仅仅与算法7-5中k α的形式不同的共轭梯度法）解问题（7.1.1）所产生的点列，则（1）当{}k x 是有穷点列时，其最后一个点是f 的平稳点；（2）当{}k x 是无穷点列时，它必有极限点，并且任一极限点都是f 的平稳点．可以证明：共轭梯度法产生的点列{}k x 是n 步二阶收敛的，即2lim sup k n k k x x q x x +→∞-≤<∞-．例4 用FR 共轭梯度法求解问题（7.2.6）222121min ()(1)2()f x x x x =-+-，仍取初始点(0)(0,0)T x =，允许误差0.1ε=．222121min ()(1)2()f x x x x =-+-解 因为212112212(1)8()()4()x x x x f x x x ⎛⎫----∇= ⎪-⎝⎭， 所以(0)(0)()(2,0),()2T f x f x ε∇=-∇=>．令(0)(0)()(2,0)T d f x =-∇=，从(0)x 出发，沿(0)d 进行一维搜索，得(1)(0)(0)001/4,(1/2,0)T x x d λλ==+=．从而(1)(1)()(0,1),()1T f x f x ε∇=-∇=>．由FR 公式有2(1)02(0)()14()f x f x α∇==∇， 因此，新的搜索方向为 (1)(1)(0)0()(1/2,1)T d f x d α=-∇+=．从(1)x 出发，沿(1)d 进行一维搜索，得(2)(1)(1)111,(1,1)T x x d λλ==+=．此时(2)(2)()(0,0),()0T f x f x ε∇=∇=<．得问题（7.2.6）的最优解为(2)(1,1)T x =．§7.4 变度量法前面介绍的最速下降法和阻尼Newton 法，它们的迭代公式可以统一为1,(),k k k k k k k x x d d G f x λ+=+⎧⎨=-∇⎩ （7.4.1）其中,n n k k G R λ⨯∈是从点k x 出发沿k d 进行一维搜索的最优步长．当k n G I =（n 阶单位矩阵）时，（7.4.1）式即为最速下降法的迭代公式；当21[()]k k G f x -=∇时，（7.4.1）式就是阻尼Newton 法的迭代公式．因此，如果能够使k G 的选取既不需要计算Hesse 矩阵及其逆矩阵，又能很好地近似于21[()]k f x -∇，则由（7.4.1）式确定的迭代算法将会保持Newton 法收敛速度快的优点，同时又具有计算简单的特性．设问题（7.1.1）中目标函数f 具有二阶连续偏导数，且2()k f x ∇正定．为使由（7.4.1）中第2式确定的搜索方向是f 在点k x 处的下降方向，根据定理1.2.3，应当要求()0T k k f x d ∇<，或即()()0T k k k f x G f x ∇∇>，所以我们应要求（7.4.1）式中的k G 是正定矩阵．设在第1k +次迭代后得到1k x +，将f 在点1k x +处作Taylor 展开，取二阶近似，得21111111()()()()()()()2TT k k k k k k f x f x f x x x x x f x x x ++++++≈+∇-+-∇-， 对上式两边求梯度，有2111()()()()k k k f x f x f x x x +++∇≈∇+∇-，令k x x =，得到2111()()()()k k k k k f x f x f x x x +++∇-∇≈∇-，即21111[()][()()]k k k k k f x f x f x x x -+++∇∇-∇≈-．易知，当f 为正定二次函数时，上式成为等式，即 21111[()][()()]k k k k k f x f x f x x x -+++∇∇-∇=-． （7.4.2）因为具有正定Hesse 矩阵的函数在极小点附近可用二次函数很好地近似，所以如果我们迫使1k G +满足类似于（7.4.2）式的关系式，即令111[()()]k k k k k G f x f x x x +++∇-∇=-， （7.4.3）则1k G +就可以很好地近似于211[()]k f x -+∇．因此称（7.4.3）式为拟Newton 条件．为方便起见，记1k k k x x x +∆=-，1()()k k k g f x f x +∆=∇-∇，1k k k G G G +∆=-，并称k G ∆为校正矩阵，则拟Newton 条件可以写成k k k k k G g x G g ∆∆=∆-∆． （7.4.4）综上所述，我们得到如下的一类算法．算法7-6（拟Newton 法）Step1 选取初始数据．选取初始点0n x R ∈和初始矩阵0G ，要求0G 为正定矩阵（可取0n G I =），给定允许误差0ε>，令0k =．Step2 检查是否满足终止准则．计算()k f x ∇，若()k f x ε∇<，迭代终止，k x 为问题（7.1.1）的近似最优解；否则，转Step3．Step3 构造搜索方向．令()k k k d G f x =-∇．Step4 进行一维搜索．求k λ和1k x +，使得()min ()k k k k k f x d f x d λλλ≥+=+， 1k k k k x x d λ+=+．Step5 产生校正矩阵．求出满足（7.4.4）的校正矩阵k G ∆，令1,:1k k k G G G k k +=+∆=+，返回Step2．算法7-7（DFP 法）Step1 选取初始数据．选取初始点0x 和初始矩阵0n G I =，给定允许误差0ε>．Step2 检查是否满足终止准则．计算0()f x ∇，若0()f x ε∇<，迭代终止，0x 为（7.1.1）的近似最优解；否则，转Step3．Step3 构造初始DFP 方向．取00()d f x =-∇，令0k =．Step4 进行一维搜索．求出k λ和1k x +，使得01()min (),.k k k k k k k k k f x d f x d x x d λλλλ≥++=+⎧⎪⎨=+⎪⎩ Step5 检查是否满足终止准则．计算1()k f x +∇，若1()k f x ε+∇<，迭代终止，1k x +为（7.1.1）的近似最优解；否则，转Step6．Step6 检查迭代次数．若1k n +=，令0:n x x =，返回Step3；否则，转Step7．Step7 构造DFP 方向．用DFP 公式1T T k k k k k k k k T T k k k k kx x G g g G G G x g g G g +∆∆∆∆=+-∆∆∆∆，算出1k G +，取111()k k k d G f x +++=-∇，令:1k k =+，返回Step4．定理7.4.4 设:n f R R →具有一阶连续偏导数，0n x R ∈，记0()f x α=，并假设水平集(,)S f α有界．若{}k x 是由DFP 法求解问题（7.1.1）产生的点列，则（1）当{}k x 是有穷点列时，其最后一个点是f 的平稳点；（2）当{}k x 是无穷点列时，它必有极限点，并且其任一极限点都是f 的平稳点．可以证明DFP 法具有超线性收敛性．算法7-8（BFGS 法）Step1 选取初始数据．选取初始点0x 和初始矩阵0n G I =，给定允许误差0ε>．Step2 检查是否满足终止准则．计算0()f x ∇，若0()f x ε∇<，迭代终止，0x 为（7.1.1）的近似最优解；否则，转Step3．Step3 构造初始BFGS 方向．取00()d f x =-∇，令0k =．Step4 进行一维搜索．求出k λ和1k x +，使得01()min (),.k k k k k k k k k f x d f x d x x d λλλλ≥++=+⎧⎪⎨=+⎪⎩ Step5 检查是否满足终止准则．计算1()k f x +∇，若1()k f x ε+∇<，迭代终止，1k x +为（7.1.1）的近似最优解；否则，转Step6．Step6 检查迭代次数．若1k n +=，令0:n x x =，返回Step3；否则，转Step7． Step7 构造BFGS 方向．用BFGS 公式111()T T T T k k k k k k k k k k k k k T T T k k k k k kx x g G g G G x g G G g x x g x g x g +⎛⎫∆∆∆∆=++-∆∆+∆∆ ⎪∆∆∆∆∆∆⎝⎭， 算出1k G +，取111()k k k d G f x +++=-∇，令:1k k =+，返回Step4．可以证明BFGS 法具有超线性收敛性．§7.5 最小二乘法在实际应用中，我们经常遇到目标函数为若干个函数的平方和的最优化问题： 21min ()()mi i s x f x ==∑， （7.5.1）其中n x R ∈，一般假设m n ≥，这类问题称为最小二乘问题．当每个()i f x 都是线性函数时，问题（7.5.1）称为线性最小二乘问题，否则，称为非线性最小二乘问题．由于最小二乘问题相对于一般无约束最优化问题而言具有特殊形式，因此除能运用本章前面介绍的一般求解方法外，还应有更为简便有效的方法．一、线性最小二乘法当()i f x 为线性函数时，即1(),1,2,,ni ij j i j f x a x b i m ==-=∑ ，问题（7.5.1）就成为线性最小二乘问题．如令11111,n m mn m a a b A b a a b ⎛⎫⎛⎫ ⎪ ⎪== ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭，12()((),(),,())T m f x f x f x f x = ，则 2()()()()()T T s x f x f x Ax b Ax b Ax b ==--=-， 从而问题（7.5.1）可表示为2min ()s x Ax b =-． （7.5.2） 因为()()()2T T T T T s x Ax b Ax b x A Ax b Ax b b =--=-+，所以()s x 为二次函数，且2()22,()2T T T s x A Ax A b s x A A ∇=-∇=．显然，对一切n y R ∈，均有20T T y A Ay Ay =≥，即知，T A A 是半正定矩阵，从而由定理2.3.7知，()s x 是n R 上的凸函数．定理7.5.1 n x R ∈是问题（7.5.2）的最优解的充要条件是x 满足方程组T T A Ax A b =． （7.5.3） 容易知道，矩阵T A A 正定的充要条件是对于一切非零向量n y R ∈，有20T T y A Ay Ay =>．记1212(,,,),(,,,)T n n y y y y A p p p == ，则上式等价于10nj j j Ay p y ==≠∑． （7.5.4）而（7.5.4）式成立的充要条件是向量组12,,,n p p p 线性无关，这又等价于rank A n =．从而T A A 为正定矩阵的充要条件是rank A n =，换句话说，()s x 为正定二次函数的充要条件是rank A n =．定理7.5.2 若rank A n =，则1()T T x A A A b -= （7.5.5）为问题（7.5.2）的唯一最优解．显然，方程组Ax b =有解的充分必要条件是问题（7.5.2）的最优值min ()0s x =．例5 求解线性最小二乘问题 2min Ax b -， （7.5.6）其中31223,3141A b ⎛⎫⎛⎫ ⎪ ⎪=-=- ⎪ ⎪ ⎪ ⎪--⎝⎭⎝⎭．解 因为11472671,()726714350T T A A A A --⎛⎫⎛⎫== ⎪ ⎪-⎝⎭⎝⎭，所以，由公式（7.5.5），求得问题（7.5.6）的最优解22673213/14137141343/103501x ⎛⎫-⎛⎫⎛⎫⎛⎫ ⎪=-= ⎪⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭ ⎪-⎝⎭． 由于问题（7.5.6）的最优值为()()11.454285710T Ax b Ax b --=≠．因此，方程组12121232,233,41x x x x x x +=⎧⎪-=-⎨⎪-+=-⎩无解．二、Gauss-Newton 法现在讨论非线性最小二乘问题2m i n ()()s x f x =， （7.5.7） 其中12,()((),(),,()),()n T m i x R f x f x f x f x f x ∈= 不全是线性函数，且假定()i f x 具有一阶连续偏导数，1,2,,i m = ．算法7-9（阻尼Gauss-Newton 法）Step1 选取初始数据．选取初始点0x ，给定允许误差0ε>，令0k =． Step2 检查是否满足终止准则．计算()k f x 及Jacobi 矩阵()k f x ∇，若()()T k k f x f x ε∇<，迭代终止，k x 为问题（7.5.7）的近似极小点；否则，转Step3．Step3 构造Gauss-Newton 方向．计算1[()()]T k k f x f x -∇∇，取1[()()]()()T T k k k k k d f x f x f x f x -=-∇∇∇．Step4 进行一维搜索．求k λ和1k x +，使得()min ()k k k k k s x d s x d λλλ≥+=+， 1k k k k x x d λ+=+．令:1k k =+，返回Step2．同阻尼Newton 法一样，阻尼Gauss-Newton 法具有全局收敛性．三、Levenberg-Marquardt 法算法7-10（LM 法）Step1 选取初始数据．选取初始点0x ，给定初始参数00α>，放大因子1β>及允许误差0ε>，令0k =．Step2 求目标函数值．计算()k f x 及()k s x ． Step3 求Jacobi 矩阵．计算()k f x ∇． Step4 检查是否满足终止准则．若()()T k k f x f x ε∇<，迭代终止，k x 为问题（7.5.7）的近似最优解；否则，转Step5．Step5 构造LM 方向．计算1[()()]T k k k n f x f x I α-∇∇+，令1[()()]()()T T k k k k n k k d f x f x I f x f x α-=-∇∇+∇．Step6 检查目标函数是否下降．计算()k k f x d +和()k k s x d +， 若()()k k k s x d s x +<，转Step8；否则，转Step7．Step7 放大参数．令:k k αβα=，返回Step5．Step8 缩小参数．令11,/,:1k k k k k x x d k k ααβ++=+==+，返回Step2． 初始参数0α和放大因子β应取适当数值，例如，根据经验可取00.01,10αβ==．可以证明，LM 法具有全局收敛性．。

最优化DFP算法姓名：施政学号：1010010125 班级：1 班专业：通信与信息系统目录1 算法流程图 (1)1.1DFP算法的流程图 (1)1.2 黄金分割法流程图 (1)1.3 回退法计算初始区间的算法 (2)2 测试函数 (3)2.1 二维、二次函数 (3)2.2 二维、高次函数 (3)2.3 高维、二次函数 (4)2.4 高维、高次函数 (4)3 运行结果及分析 (5)4 Matlab源程序 (6)4.1 主函数 (6)4.2 DFP算法函数 (8)4.3 黄金分割法函数 (9)4.4 回退法求解初始区间 (10)4.5 计算测试函数的值 (11)4.6 计算测试函数的梯度 (12)5 参考文献 (12)1 算法流程图对于DFP算法主要涉及到3个主要的算法，分别是：利用回退法计算初始区间、利用黄金分割法进行一维搜索、然后利用DFP算法计算最小点对应的自变量的值。

下面分别画出了这三个算法流程图。

1.1DFP算法的流程图设定控制误差为ε；输入的初始点坐标是0x；0E是与0x同维的单位阵。

图 11.2 黄金分割法流程图给定精确度ε>0;当区间长度小于等于ε时，即停止运行，同时取x=(a+b)/2作为最小点坐标。

在给定初始区间[a,b]内，求最小点时对应的α值，要保证α是大于等于零的，否则函数值就不是朝下降方向递降的了。

在本算法中保证初始区间的端点是大于等于零的，就可以满足这一条件了。

算法如下图所示：图 21.3 回退法计算初始区间的算法针对这个算法,参考文献[1]上面利用的回退法不能保证a ，b 两个端点的值大于零，因为利用黄金分割法求α时，α肯定是大于等于零的，所以可以对书上的算法适当的改进。

初始的点是0x ；步长是x Δ；算法如下：图 3注意：上面的算法是针对一维的情况，所以在计算0x ；1x ；2x 时，应该注意使用0000*x a t p =+；0a 代表的是对应DFP 算法上次迭代的自变量的坐标值，0p 代表的是0a 点的梯度，0t 开始的值是0；同样有1010*x a t p =+；2020*x a t p =+；10t t t =+Δ；21t t t =+Δ。

无约束优化方法1.坐标轮换法2.鲍威尔法3.梯度法4.牛顿法5.变尺度法1.坐标轮换法坐标轮换发是一种不计算函数梯度,而是通过函数值本身,即可求出寻优方向,因而也称为直接寻优法.在以后提到的鲍威尔法(Powell)法也属于直接寻优法。

对于坐标轮换法，我们做个比喻：如果我们在北京的老城区找一个地方，我们可以沿着经纬线去找。

这个比喻为我们提供了一种思路，既可以取坐标的方向为寻优的方向，这就是坐标轮换法。

它在每次搜索中，只允许一个变量的变化，其余量保持不变，即沿着坐标方向轮流进行搜索的方法。

该方法把多变量的优化问题轮流转化成一系列单变量的优化问题。

对应于n 个变量的寻优函数，若在第轮沿第k 个坐标第i 个坐标方向ki i S e =进行搜索，则迭代公式为1(0,1,...,1,2,...,)k k k k i i i i X X S k i n α-=+==其中搜索方向取坐标方向，即k i i S e =（1,2,...,i n =）。

若0k k n X X -‖‖<ε，则*kn X X ←，否则10k kn X X +←，进行下一轮的搜索，一直到满足精度要求为止。

其搜索路径如图所示这种方法的收敛效果与目标函数等值线形有很大关系。

如果目标函数为二元二次函数，其等值线为圆或长轴平行于坐标轴的椭圆时，此方法很有效，经过两次搜索即可以达到最优点，如图所示。

如果等值线为长轴不平行于坐标轴的椭圆，则需多次迭代才能达到最优点，但因坐标轮换法是坐标方向搜索而不是沿脊线搜索，所以就终止到脊线上而不能找到最优解。

从上述分析可以看出，采用坐标轮换法只能轮流沿着坐标的方向搜索，尽管也能使函数值步步下降，但经过曲折迂回的路径才能达到极值点；尤其极值点附近步长很小，收敛很慢，所以坐标轮换法不是一种很好的搜索方法。

但是可以构造很好的搜索策略，下面讨论的鲍威尔法就是这种情况。

例题：已知22121212()10460f X x x x x x x =+---+，设初始点：(0)[0,0]T X=，精度0.1=ε，用最优步长法的坐标轮换法求目标函数的最优解。