实验2 最速下降法和共轭梯度法的程序设计
(完整)实验2 最速下降法和共轭梯度法的程序设计
实验2 最速下降法和共轭梯度法的程序设计一、实验目的1、熟悉无约束优化问题的最速下降算法和共轭梯度法。
2、培养matlab 编程与上机调试能力。
二、实验课时：2个课时三、实验准备1、预习无约束优化问题的最速下降算法和共轭梯度法.2、熟悉matlab 软件的基本操作及程序编写。
四、实验内容课堂实验演示根据最速下降法编写程序，求函数21222121342)(m in x x x x x x x f -++-=的极小值，其中初始点为()01,1Tx =算法步骤如下：Step1：：给出初始点0x ,和精度1;0k ε<<=;Step2：计算()k f x ∇，如果()k f x ε∇≤，则停止迭代，输出结果;否则转step3；Step3：令下降方向()k k d f x =-∇，计算步长因子k λ使得0()min ()k k k k k f x d f x d λλλ≥+=+，令1,1k k k k x x d k k λ+=+=+,转step2.其程序如下:function ［x,iter,val,dval] = Steepest_Descent_Method(x ，eps)k = 1;dy = grad_obj(x);x_mat （:,1) = x ；%存储每一次迭代得到的点xwhile norm （dy)〉epsd = —dy ； % 搜索方向lambda = line_search(x ，d ）;%步长x = x + d ＊lambda ；k = k + 1；x_mat （:,k ） = x;dy = grad_obj （x)；enditer = k — 1;val = obj(x);％目标函数在极值点处的函数值dval = grad_obj(x)；％目标函数在极值点处的梯度%———----—-————-—-————-—-———-——---—----——-—---——---—-—-——--——--—x1 = linspace(—1。
最速下降法和共轭梯度法
v b Ax v, v t v, Av x x tv
output k+1， x enddo
算法中也可以用残差作为循环控制，我们可以利用以下的事后误差估计式
|| x x* || || r || ( A) || x || || b ||
证明 按照残差的定义，并且由于假设 x 是方程组的准确解，我们有
(0)
for k=0 to M-1 do
v ( k 1) b Ax ( k ) v( k ) , v( k ) t v ( k ) , Av ( k ) x ( k 1) x ( k ) tk v ( k )
output k+1， x enddo
( k 1)
为了节省空间，实际的算法为 input x ，A, b, M output 0, x for k=0 to M-1 do
定理 1 （A-标准正交系定理）设 A 对称正定， u
(1)
，u
(2)
，…， u
(n)
是一组关于内积
x, y A 的标准正交系。定义 x ( i ) x ( i 1) b Ax (i 1) , u ( i ) u ( i )
其中 x
(0)
是 R 中的任意点，则 Ax
n
( k 1) r (k ) , v(k ) (k ) (k ) x x v (k ) (k ) v Av , ( k 1) b Ax ( k 1) r ( k 1) (i ) k A (i ) v ( k 1) r ( k 1) r (i ) , v v (i ) A i 0 v , v
2
v, b Ax q( x)
matlab数学实验
《管理数学实验》实验报告班级姓名实验1：MATLAB的数值运算【实验目的】（1）掌握MATLAB变量的使用（2）掌握MATLAB数组的创建，（3）掌握MA TLAB数组和矩阵的运算。
（4）熟悉MATLAB多项式的运用【实验原理】矩阵运算和数组运算在MA TLAB中属于两种不同类型的运算，数组的运算是从数组元素出发，针对每个元素进行运算，矩阵的运算是从矩阵的整体出发，依照线性代数的运算规则进行。
【实验步骤】（1）使用冒号生成法和定数线性采样法生成一维数组。
（2）使用MA TLAB提供的库函数reshape，将一维数组转换为二维和三维数组。
（3）使用逐个元素输入法生成给定变量，并对变量进行指定的算术运算、关系运算、逻辑运算。
（4）使用MA TLAB绘制指定函数的曲线图，将所有输入的指令保存为M文件。
【实验内容】（1）在[0,2*pi]上产生50个等距采样数据的一维数组，用两种不同的指令实现。
0:(2*pi-0)/(50-1):2*pi 或linspace(0,2*pi,50)（2）将一维数组A=1:18，转换为2×9数组和2×3×3数组。
reshape(A,2,9)ans =Columns 1 through 71 3 5 7 9 11 132 4 6 8 10 12 14Columns 8 through 915 1716 18reshape(A,2,3,3)ans(:,:,1) =1 3 52 4 6ans(:,:,2) =7 9 118 10 12 ans(:,:,3) =13 15 17 14 16 18（3）A=[0 2 3 4 ;1 3 5 0],B=[1 0 5 3;1 5 0 5]，计算数组A 、B 乘积，计算A&B,A|B,~A,A= =B,A>B 。
A.*Bans=0 0 15 121 15 0 0 A&Bans =0 0 1 11 1 0 0 A|Bans =1 1 1 11 1 1 1~Aans =1 0 0 00 0 0 1A==Bans =0 0 0 01 0 0 0A>=Bans =0 1 0 11 0 1 0（4）绘制y= 0.53t e -t*t*sin(t),t=[0,pi]并标注峰值和峰值时间,添加标题y= 0.53t e -t*t*sint ，将所有输入的指令保存为M 文件。
用MATLAB实现最速下降法_牛顿法和共轭梯度法求解实例——张小强
机电产品优化设计课程设计报告姓名：张小强学号：201222080633学院：机械电子工程学院实验的题目和要求一．课程名称：最优化设计方法二．实验日期：2013年6月27日三．实验目的：掌握最速下降法，牛顿法和共轭梯度法的算法思想，并能上机编程实现相应的算法。
四．实验要求：用MATLAB 实现最速下降法，牛顿法和共轭梯度法求解实例。
五．实验原理：最速下降法是以负梯度方向最为下降方向的极小化算法，相邻两次的搜索方向是互相直交的。
牛顿法是利用目标函数)(x f 在迭代点k x 处的Taylor 展开式作为模型函数，并利用这个二次模型函数的极小点序列去逼近目标函数的极小点。
共轭梯度法它的每一个搜索方向是互相共轭的，而这些搜索方向k d 仅仅是负梯度方向k g -与上一次搜索方向1-k d 的组合。
五．运行结果如下: 题目：f=(x-2)^2+(y-4)^2①．最速下降法：M 文件：function [R,n]=steel(x0,y0,eps)syms x ;syms y ;f=(x-2)^2+(y-4)^2;v=[x,y];j=jacobian(f,v);T=[subs(j(1),x,x0),subs(j(2),y,y0)];temp=sqrt((T(1))^2+(T(2))^2);x1=x0;y1=y0;n=0;syms kk ;while (temp>eps)d=-T;f1=x1+kk*d(1);f2=y1+kk*d(2);fT=[subs(j(1),x,f1),subs(j(2),y,f2)];fun=sqrt((fT(1))^2+(fT(2))^2);Mini=Gold(fun,0,1,0.00001);x0=x1+Mini*d(1);y0=y1+Mini*d(2);T=[subs(j(1),x,x0),subs(j(2),y,y0)];temp=sqrt((T(1))^2+(T(2))^2);x1=x0;y1=y0;n=n+1;endR=[x0,y0];调用黄金分割法：M文件：function Mini=Gold(f,a0,b0,eps)syms x;format long;syms kk;u=a0+0.382*(b0-a0);v=a0+0.618*(b0-a0);k=0;a=a0;b=b0;array(k+1,1)=a;array(k+1,2)=b;while((b-a)/(b0-a0)>=eps)Fu=subs(f,kk,u);Fv=subs(f,kk,v);if(Fu<=Fv)b=v;v=u;u=a+0.382*(b-a);k=k+1;elseif(Fu>Fv)a=u;u=v;v=a+0.618*(b-a);k=k+1;endarray(k+1,1)=a;array(k+1,2)=b;endMini=(a+b)/2;输入：[R,n]=steel(0,1,0.0001)输出：R = 1.99999413667642 3.99999120501463n = 1②．牛顿法：M文件：syms x1x2;f=(x1-2)^2+(x2-4)^2;v=[x1,x2];df=jacobian(f,v);df=df.';G=jacobian(df,v);epson=1e-12;x0=[0,0]';g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});G1=subs(G,{x1,x2 },{x0(1,1),x0(2,1)});k=0;mul_count=0;sum_count=0;mul_count=mul_count+12;sum_count=sum_count+6; while(norm(g1)>epson)p=-G1\g1;x0=x0+p;g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});G1=subs(G,{x1,x2},{x0(1,1),x0(2,1)});k=k+1;mul_count=mul_count+16;sum_count=sum_count+11; end;kx0mul_countsum_count结果：k = 1x0 =24mul_count = 28sum_count = 17③．共轭梯度法：M文件：function f=conjugate_grad_2d(x0,t)x=x0;syms xi yi af=(xi-2)^2+(yi-4)^2;fx=diff(f,xi);fy=diff(f,yi);fx=subs(fx,{xi,yi},x0);fy=subs(fy,{xi,yi},x0);fi=[fx,fy];count=0;while double(sqrt(fx^2+fy^2))>ts=-fi;if count<=0s=-fi;elses=s1;endx=x+a*s;f=subs(f,{xi,yi},x);f1=diff(f);f1=solve(f1);if f1~=0ai=double(f1);elsebreakx,f=subs(f,{xi,yi},x),countendx=subs(x,a,ai);f=xi-xi^2+2*xi*yi+yi^2;fxi=diff(f,xi);fyi=diff(f,yi);fxi=subs(fxi,{xi,yi},x);fyi=subs(fyi,{xi,yi},x);fii=[fxi,fyi];d=(fxi^2+fyi^2)/(fx^2+fy^2);s1=-fii+d*s;count=count+1;fx=fxi;fy=fyi;endx,f=subs(f,{xi,yi},x),count输入：conjugate_grad_2d([0,0],0.0001)结果：x = 0.24998825499785 -0.24999998741273f = 0.12499999986176count = 10ans = 0.12499999986176六．结论如下:最速下降法越接近极小值，步长越小，前进越慢。
共轭梯度法程序
一、共轭梯度法共轭梯度法（Conjugate Gradient）是共轭方向法的一种，因为在该方向法中每一个共轭向量都是依靠赖于迭代点处的负梯度而构造出来的，所以称为共轭梯度法。
由于此法最先由Fletcher和Reeves （1964）提出了解非线性最优化问题的，因而又称为FR 共轭梯度法。
由于共轭梯度法不需要矩阵存储，且有较快的收敛速度和二次终止性等优点，现在共轭梯度法已经广泛地应用于实际问题中。
共轭梯度法是一个典型的共轭方向法，它的每一个搜索方向是互相共轭的，而这些搜索方向d仅仅是负梯度方向与上一次迭代的搜索方向的组合，因此，存储量少，计算方便，效果好。
二、共轭梯度法的原理设有目标函数f（X）=1/2X T HX+b T X+c 式1 式中，H作为f（X）的二阶导数矩阵，b为常数矢量，b=[b1,b2,b3,...b n]T 在第k次迭代计算中，从点X（k）出发，沿负梯度方向作一维搜索，得S(K)=-∆f（X（k））式2 X(k+1)=X(k)+ɑ(k)S(k) 式3在式中，ɑ（k）为最优步长。
设与S(k）共轭的下一个方向S(k+1)由点S(k）和点X(k+1)负梯度的线性组合构，即S (k+1)=-∆f （X （k+1））+β(k)S (k) 式4 根据共轭的条件有[S (k)]T ∆2f （X (k)）S (k+1)=0 式5 把式2和式4带入式5，得-[∆f(X (k))]T ∆2f （X (k)）[-∆f （X （k+1））+β(k)S (k) ]=0 式6 对于式1，则在点X (k)和点X (k+1)的梯度可写为∆f(X (k))=HX (k)+b 式7 ∆f （X （k+1）)=HX (k+1)+b 式8 把上面两式相减并将式3代入得ɑ(k)H S (k)=∆f （X （k+1）)-∆f(X (k)) 式9 将式4和式9两边分别相乘，并代入式5得-[∆f （X （k+1）)+β(k)∆f(X (k))]T [∆f （X （k+1）)-∆f(X (k)]=0 式10 将式10展开，并注意到相邻两点梯度间的正交关系，整理后得 β（k ）=22||))((||||))1((||k X f k X f ∆+∆ 式11把式11代入式4和式3，得 S （k+1）=-∆f （X (k)）+β（k ）S （k ）X (k+1)=X （k ）+ɑ（k ）S （k ）由上可见，只要利用相邻两点的梯度就可以构造一个共轭方向。
最优化问题的算法迭代格式
最优化问题的算法迭代格式最优化问题的算法迭代格式最优化问题是指在一定的条件下，寻找使某个目标函数取得极值（最大值或最小值）的变量取值。
解决最优化问题的方法有很多种，其中较为常见的是迭代法。
本文将介绍几种常用的最优化问题迭代算法及其格式。
一、梯度下降法梯度下降法是一种基于负梯度方向进行搜索的迭代算法，它通过不断地沿着目标函数的负梯度方向进行搜索，逐步接近极值点。
该方法具有收敛速度快、易于实现等优点，在许多应用领域中被广泛使用。
1. 算法描述对于目标函数 $f(x)$，初始点 $x_0$ 和学习率 $\alpha$，梯度下降算法可以描述为以下步骤：- 计算当前点 $x_k$ 的梯度 $\nabla f(x_k)$；- 更新当前点 $x_k$ 为 $x_{k+1}=x_k-\alpha\nabla f(x_k)$；- 如果满足停止条件，则输出结果；否则返回第 1 步。
2. 算法特点- 沿着负梯度方向进行搜索，能够快速收敛；- 学习率的选择对算法效果有重要影响；- 可能会陷入局部极小值。
二、共轭梯度法共轭梯度法是一种基于线性方程组求解的迭代算法，它通过不断地搜索与当前搜索方向共轭的新搜索方向，并在该方向上进行一维搜索，逐步接近极值点。
该方法具有收敛速度快、内存占用少等优点，在大规模问题中被广泛使用。
1. 算法描述对于目标函数 $f(x)$，初始点 $x_0$ 和初始搜索方向 $d_0$，共轭梯度算法可以描述为以下步骤：- 计算当前点 $x_k$ 的梯度 $\nabla f(x_k)$；- 如果满足停止条件，则输出结果；否则进行下一步；- 计算当前搜索方向 $d_k$；- 在当前搜索方向上进行一维搜索，得到最优步长 $\alpha_k$；- 更新当前点为 $x_{k+1}=x_k+\alpha_k d_k$；- 计算新的搜索方向 $d_{k+1}$；- 返回第 2 步。
2. 算法特点- 搜索方向与前面所有搜索方向都正交，能够快速收敛；- 需要存储和计算大量中间变量，内存占用较大；- 可以用于非线性问题的求解。
非线性方程组的牛顿迭代法-最速下降法
数学软件实验任务书实验一 非线性方程组的牛顿迭代法 1 实验原理对于非线性方程11221212(,,,)(,,,)(,,,)n n n n f x x x f x x x f f x x x ⎛⎫ ⎪ ⎪= ⎪ ⎪⎝⎭ 在x (k )处按照多元函数的泰勒展开，并取线性项得到 ()()()()(1)()1111()()()()(1)()()122()()()()(1)()1(,,,)(,,,)()0(,,,)k k k k k k n n k k k k k k k n n n k k k k k k n n n nn f x x x x x f x x x x x f x f x x x x x +++⎛⎫⎛⎫- ⎪ ⎪- ⎪ ⎪'+= ⎪ ⎪ ⎪ ⎪ ⎪ ⎪-⎝⎭⎝⎭ 其中1111()n n n n f f x x f x f f x x ∂∂⎛⎫ ⎪∂∂ ⎪' ⎪= ⎪∂∂ ⎪ ⎪∂∂⎝⎭(1)()()()()()1111(1)()()()()()()1221(1)()()()()()1(,,,)(,,,)[()](,,,)k k k k k k n n k k k k k k k n n n k k k k k k n n n n n x x f x x x x x f x x x f x x x f x x x ++-+⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪'=- ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭ 2 数据来源计算非线性方程组22220.50440x x y x y ⎧--+=⎨+-=⎩ 初值取11x y ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦3 实验步骤步骤一：编写牛顿迭代法的基本程序。
打开Editor 编辑器，输入以下语句：function [x,n,data]=new_ton(x0,tol)if nargin==1tol=1e-10;endx1=x0-f1(x0)/df1(x0);MATLAB 241 数值分析与应用n=1;%迭代过程while (norm(x1-x0)>tol)x0=x1;x1=x0-f1(x0)/df1(x0);n=n+1;%data 用来存放中间数据data(:,n)=x1;endx=x1;以文件名new_ton.m保存。
最速下降法与共轭梯度法实验
x1=a;x2=b;
Q=fsxhesse(f,x1,x2);
x0=[x1 x2]';
fx1=diff(f,'x1'); %对x1求偏导数
fx2=diff(f,'x2'); %对x2求偏导数
g=[fx1 fx2]'; %梯度
g1=subs(g); %把符号变量转为数值
a(i)=-(d(:,i)'*g1(:,i))/(d(:,i)'*G*d(:,i));
X(:,i+1)=X(:,i)+a(i)*d(:,i);
g1=[g1 subs(subs(g,x1,X(1,i+1)),x2,X(2,i+1))];
m=norm(g1(:,i+1));
i=i+1;
end
case 'WYL'
i=i+1;
end
case 'DY'
while m>=e
k(i-1)=g1(:,i)'*g1(:,i)/(d(:,i-1)'*(g1(:,i)-g1(:,i-1)));
d(:,i)=-g1(:,i)+k(i-1)*d(:,i-1);
a(i)=-(d(:,i)'*g1(:,i))/(d(:,i)'*G*d(:,i));
fxx=subs(fxx); %将符号变量转化为数值
fxy=subs(fxy);
fyx=subs(fyx);
fyy=subs(fyy);
x=[fxx,fxy;fyx,fyy]; %求hesse矩阵
运行函数
syms x1 x2;
最速下降法和共轭阶梯法
最速下降法和共轭阶梯法一、实验目的对矩阵A=[2 1 0;1 2 0;0 0 3];b=[20 34 -26]';对方程A*X=b，分别用最速下降法和共轭阶梯法，求x的值，并判断是否收敛和比较两种算法的优势二、实验原理最速下降法基于这样的观察：如果实值函数在点处可微且有定义，那么函数在点沿着梯度相反的方向下降最快。
因而，如果对于γ > 0 为一个够小数值时成立，那么。
考虑到这一点，我们可以从函数F的局部极小值的初始估计出发，并考虑如下序列使得因此可得到如果顺利的话序列收敛到期望的极值。
注意每次迭代步长γ可以改变。
三、上侧的图片示例了这一过程，这里假设F定义在平面上，并且函数图像是一个碗形。
蓝色的曲线是等高线(水平集)，即函数F为常数的集合构成的曲线。
红色的箭头指向该点梯度的反方向。
（一点处的梯度方向与通过该点的等高线垂直）。
沿着梯度下降方向，将最终到达碗底，即函数F值最小的点。
共轭梯度法是求解特定线性系统的数值解的方法，其中那些矩阵为对称和正定。
共轭梯度法是一个迭代方法，所以它适用于稀疏矩阵系统，因为这些系统对于象乔莱斯基分解这样的直接方法太大了。
这种系统在数值求解偏微分方程时相当常见。
共轭梯度法也可以用于求解无约束的最优化问题。
双共轭梯度法提供了一种处理非对称矩阵情况的推广。
设我们要求解下列线性系统Ax = b,，其中n-×-n矩阵A是对称的（也即，A T = A），正定的（也即，x T Ax > 0对于所有非0向量x属于R n），并且是实系数的。
将系统的唯一解记作x*。
四、实验内容实验步骤：1、根据最速下降法和共轭梯度法的原理编程2、设置迭代次数n3、取得er的不同，研究对收敛的影响4、设置tx把每步的x的用plot（）函数画出图像，研究对收敛的影响输出结果最速下降法（1）取er=10^(-6)时，算出x的结果及迭代次数n x =2.0000 16.0000 -8.6667n=12收敛性可以看出收敛比较好（2）取er=10^(-13)时，算出x的结果及迭代次数n x = 2.0000 16.0000 -8.6667n=21收敛性（3）取er=10^(-20)时，算出x的结果及迭代次数n x = 2.0000 16.0000 -8.6667n=27收敛性共轭梯度法（1）取er=10^(-6)时，算出x的结果及迭代次数n x = 2.0000 16.0000 -8.6667n=3收敛性（2）取er=10^(-20)时，算出x的结果及迭代次数n x = 2.0000 16.0000 -8.6667n=4收敛性（3）取er=10^(-30)时，算出x的结果及迭代次数n x = 2.0000 16.0000 -8.6667n=4收敛性数据分析（比较两种算法的优势）根据实验结果可以知道，在相同的精度下，最速下降法的迭代次数要比共轭梯度法的迭代次数要多，而但是它们算出来的结果是一样的，另外从收敛性行可以看出两者都是线性收敛，但总体来说共轭梯度法要比最深下降法要好。
线性方程组的最速下降法与共轭梯度法解法___实验代码参考
一、实验原理：1．最速下降法：从某个初始点)0(X 出发，沿)(X f 在点)0(X 处的负梯度方向)0()0()0()(AX b X f r -=-∇=求得)(X f 的极小值点)1(X , 即)(min )0()0(0r X f λλ+> 然后从)1(X 出发，重复上面的过程得到)2(X 。
如此下去，得到序列{)(k X })(...)()()()1()0(k X f X f X f >>>可以证明，从任一初始点)0(X 出发， 用最速下降法所得到的序列{)(k X }均收敛于问题使X 最小化)(X f 的解，也就是方程组b AX =的解。
其收敛速度取决于11λλλλ+-n n ，其中1λ ,n λ分别为A 的最小,最大特征值。
最速下降法迭代格式：给定初值)0(X ，)(k X 按如下方法决定：())()(1)(k )()()()(k )()(X ,,)(k k k k T k k T k k k k r X Ar r r r AX b X f r λλ+=><><=-=-∇=+ 2．共轭梯度法其基本步骤是在点)(k X 处选取搜索方向)(k d , 使其与前一次的搜索方向)1(-k d 关于A 共轭，即(1)()(1),0k k k d d Ad --<>=然后从点)(k X 出发，沿方向)(k d 求得)(X f 的极小值点)1(+k X , 即)(min )()()(0)1(k d X f X f k k λλ+=>+ 如此下去, 得到序列{)(k X }。
不难求得0,)1()(>=<-k k Ad d 的解为 )()1()1()()()()1(,,k k k k k k k d Ad d d AX b X X ><>-<+=--+ 注意到)(k d 的选取不唯一，我们可取)1(1)()()(--+-∇=k k k k d X f d β由共轭的定义0,)1()(>=<-k k Ad d 可得：><><-=----)1()1()1()(1,,k k k k k Ad d Ad r β共轭梯度法的计算过程如下：第一步：取初始向量)0(X , 计算⎪⎪⎩⎪⎪⎨⎧+=><><-=-=-∇==(0)0(0)(1))0()0()0()0(0(0)(0)(0)(0)d X X ,,X )X (r d λλAd d Ad r A b f 第1+k 步：计算⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧+=><><-=+=><><-=-=-∇=+------(k)0(k)1)(k )()()()()1(1(k))()1()1()1()(1(k)(k)(k)d X X ,,r ,,X )X (r λλββk k k k k k k k k k k k k Ad d Ad r d d Ad d Adr A b f二、实验内容：%最速下降法function [x,k]=fastest(A,b,eps); x0=zeros(size(b),1);x=x0;k=0;m=1000;tol=1;while tol>=epsr=b-A*x0;q=dot(r,r)/dot(A*r,r);x=x0+q*r;k=k+1;tol=norm(x-x0);x0=x;if k>=mdisp('迭代次数太多，可能不收敛！');return;endendxk%共轭梯度法function [k,x]=gong_e(A,b)esp=input('请输入允许误差esp=');x0=input('请输入初始值x0=');k = 0 ;r0 = b-A*x0; %求出dangqian梯度while norm(r0)>espr0 = b -A*x0;k = k + 1 ;if k==1p0 = r0 ;elselamda=(r0'*r0)/(p0'*A*p0); r1 = r0 - lamda*A*p0 ;p0=r0+(r0'*r0)/(r1'*r1)*p0; x1 = x0 + lamda*p0;x0=x1;r0=r1;endendx=r0;k;end三、实验结果：A=[5 2 0;6 4 1;1 2 5];b=[10 18 -14]';eps=1.0e-6;x =-0.87507.1875-5.5000k =60。
梯度法和共轭梯度法
极小点，其中Bk{ x Nhomakorabeax k
id (i) ,i
R}
i 1
是由 d (1) , d (2) ,, d (k) 生成的子空间。特别地，当k n时，x(n1)是
f ( x)在Rn上的唯一极小点。
推论 在上述定理条件下，必有 f ( x(k1) )T d (i) 0, i 1, 2,, k。
3、共轭方向法
i
g
i
T 1
(
gi1
gi
)
d (i)T ( gi1 gi )
||
gi1 ||2 d (i)T gi
|| gi1 ||2 || gi ||2
(4)
FR算法步骤：
1. 任取初始点x(1) ，精度要求 ，令k 1。 2. 令g1 f ( x(1) ),若 || g1 || ,停止，x(1)为所求极小点；
上式两边同时左乘d jT A，则有
k
id
jT Ad i
0,
i 1
因为d 1 , d 2 ,, d k 是 k 个 A共轭的向量，所以上式可化简为
j d jT Ad j 0 .
因为d j 0,而 A是正定矩阵，所以d jT Ad j 0，
所以
j 0, j 1, 2,, k。
因此 d1 ,d 2 ,,d k 线性无关。
2. 共轭方向
定义 设 A 是 n n的对称正定矩阵，对于Rn中的两个非零向量 d1 和 d 2， 若有 d1T Ad 2 0，则称d 1和d 2关于A共轭。 设 d1 , d 2 ,,d k 是 Rn 中一组非零向量，如果它们两两关于A 共轭，即 d iT Ad j 0, i j , i , j 1, 2,, k。 则称这组方向是关于A共轭的，也称它们是一组A共轭方向。
共轭梯度法和最速下降法的混合算法
收稿日期:1998-10-15 第一作者:男,1970年生,硕士,讲师共轭梯度法和最速下降法的混合算法欧志英 严克明 王柏岩(甘肃工业大学基础科学系,兰州 730050)摘 要 将共轭梯度法与最速下降法有机地结合起来,构造了一种共轭梯度法和最速下降法的混合算法,并证明了该算法的全局收敛.混合算法既提高了共轭梯度算法的收敛速度,又解决了目标函数“性态不优”时,最速下降法难以求解的问题.同时也可以看到共轭梯度法与最速下降法仅仅是混合算法的特例.关键词 混合算法 共轭梯度法 最速下降法 全局收敛分类号 O224FR 算法是Fletcher 和Reev s 共同提出的一种共轭梯度算法[1].由于FR 算法具有二次终止性,内存需要量小,主要存储四个n 维向量(x k ,g k ,g k -1,p k -1),程序简单,具有较快的收敛速度,因此它是解决无约束优化问题,特别是大规模问题的一个重要方法.但是FR 算法的计算效果总不如人意,当‖ f (x k )‖较小时,FR 算法收敛较慢.最速下降法具有较高的局部收敛速度,但在有的情况下全局收敛速度不高,如目标函数“性态不优”时,其收敛速度就较慢.本文给出一种将共轭梯度法与最速下降法有机地结合起来的新算法混合算法,在保证相同的计算量前提下,它将进一步提高FR 算法的计算效果.1 算法混合算法的步骤如下:1)取初始点x 0∈R n ,T >0,U >0,置k =1;2)计算g k =g (x k )= f (x k );3)若g k =0,则停止计算,否则置:p k =-g k +U k -1p k -1其中,U k -1=0 (k =1)‖g k ‖2‖g k -1‖2 (k >1);4)若‖ f (x k )‖较大,取U ∈12,1;否则,取U ∈0,12;5)若d 1(d 为目标函数f (x )的Hessian 阵的条件数),取T ∈0,12;否则,取T ∈12,1;6)p k ′=-T g k +U p k ;第25卷第1期1999年3月甘 肃 工 业 大 学 学 报Jour nal of Gansu U niv er sity o f Technolog y V o l.25No.1M ar.19997)一维搜索,求λk ,使得f (x k +λk p k ′)=min 〈f (x k +λp k ′)|λ>0〉;8)置x k +1=x k +λk p k ′,k =k +1,返回步骤2).对于上述算法,给出以下几点注释:1)p 1′,p 2′,…,p k ′,…是共轭方向.从文[1]中已知FR 算法得到的p 1,p 2,…,p k ,…是共轭方向,满足:(p i ,p j )G =0 (j ≠i )因而可得知:(p i ′,p j ′)G =0 (j ≠i )即p 1′,p 2′,…,p k ′,…是共轭方向.2)p 1′,p 2′,…,p k ′,…是目标函数的下降方向.(-g k ,p k ′)=(-g k ,-T g k +U p k )=(-g k ,T g k )+(-g k ,U p k )=T ‖g k ‖2+U (-g k ,-g k +∑k -1j =1U jp j )=‖g k ‖2>02 混合算法的全局收敛引理1 设函数f (x )二阶连续可微, f (x )满足Lipschitz 条件,且∑∞k =1cos 2〈d k ,-g k 〉<+∞,其中d k 为搜索方向.则利用共轭梯度法和精确线性搜索所得到的点列{x k }必有限终止,或者下面两极限必有一为真[2]:lim k →∞f (x k )=-∞lim k →∞inf ‖ f (x k )‖=0 由引理1可知,只要夹角〈d k ,-g k 〉不为π2,则共轭梯度算法在满足引理1的条件下是全局收敛的.定理1 混合算法在满足引理1的条件下是全局收敛的.证明 由于(-g k ,p k ′)≥‖g k ‖2,因此,co s 〈-g k ,p k ′〉≠0,则〈-g k ,p k ′〉不为π2.从引理1知定理1成立.3 混合算法的计算效果分析从算法中可知,混合算法和FR 算法的计算量大致相同,但计算效果却有不小的差异.一般来说前者优于后者,为说明这个问题,举例如下[3]:考虑min x ∈R 2f (x ),其中在圆域D 内f (x )=12(x 21+x 22),在D 的边界是连续可微函数,在D 的边界上是光滑联接的.设想从D 的边界外某一点出发,经过若干次迭代后达到D 内某点x k ,并且已经确定了共轭方向p k ,现从x k 出发,分别用FR 算法和混合算法进行计算,并对它们的进展情况作出比较.‖p k ‖‖g k ‖co s θk =(p k ,-g k )=‖g k ‖2sec θk =‖p k ‖‖g k ‖·90·甘肃工业大学学报 第25卷‖g k ‖‖p k ′‖co s θk ′=(p k ′,-g k )=(T +U )‖g k ‖2sec θk ′=‖-T g k +U p k ‖(T +U )‖g k ‖1‖g k ‖<‖p k ‖‖g k ‖=sec θk 因此,θk >θk ′,并且可以证明θk =θk +1.这说明了对于FR 算法来说,即使对于性态很好的目标函数,倘若某次搜索方向很不理想,仍会使得以后的每次迭代都只能取得很小的进展.但是对于混合算法则不一样,它具有朝最速下降方向靠拢的性质.由此可见,混合算法要比FR 算法的计算效果优.同时还可看到:T =0时,混合算法就是FR 算法;U =0时,混合算法就是最速下降法.4 结论1)最速下降法和FR 算法仅是混合算法的特例.2)混合算法的搜索方向是关于目标函数f (x )的Hessian 阵共轭.3)混合算法的计算效果要比FR 算法的计算效果优,且全局收敛;但是否具有比FR 算法更高的Q 收敛速度,还需要进一步探讨.参 考 文 献1 Fletche r R,Rv ees C M.Functio n minimization by co njug ate g radients.Co mput J ,1964(7):149～1542 袁亚湘.非线性规划数值方法.上海:上海科学技术出版社,1993.57～753 邓乃扬.无约束最优化计算方法.北京:科学技术出版社,1982.87～122A mixed method of conjugate gradient methodand steepest descent methodOu Zhiying ,Yan Keming ,Wang Baiyan(Dept.o f Basic Sciences,Ga nsu U niv.of T ech.,Lanzh ou 730050)Abstract The co njugate g radient method and the steepest descent method a re combined ,and a mix ed method o f conjug ate g radient m ethod and steepest descent is created ,and its g lobal conv ergence is prov ed.The mixed m ethod raise the co nv erg ent rate o f the co njugate g radient method,and solve the problem which the steepest descent method can no t solv e in the conditio n with badly cha racteristics for the o bjectiv e function .In the m eantime ,it can be seen that the co njuga te g radient m ethod o r steepest descent metho d is a special case o f the mix ed method.Key words mixed method,co njugate g radient m ethod,steepest descent method,g loba l conv ergence ·91·第1期 欧志英等:共轭梯度法和最速下降法的混合算法。
CH5-1,2变分原理最速下降法、共轭梯度法
k 1
使得 ( xk 1 )
xx0 span ( p1 ,, pk 1 )
min
( x) .
以下求向量 pk ， 令
xk xk 1 k pk
使得
（2）
( xk )
xx0 span ( p1 ,, pk )
min
( x)
由上节（4）式可知： ( xk ) ( xk 1 k pk )
因为
APk 1 yk 1 i Api
i 1
k 1
所以如果取 pk 满足
( pk , Api ) 0, i 1,2,, k 1
（3） （A共轭）
则
( xk ) ( xk 1 ) k ( APk 1 yk 1 , pk )
k2
2
( Apk , pk ) k (b, pk )
i, j 1,, m n，则称这组向量为 A共轭，或 A正交.
三、共轭梯度法算法（理论上）
取 x0 0; 对于 k 1,2,, n ，计算
rk 1 b Axk 1.
若 rk 1 0，则 x xk 1，停止； 否则 若 k 1，则 p1 r0 否则 选择 pk 满足 ( pi , Apk ) 0(i 1,, k 1)
的第二项为零.
这时，要使 xk 极小化 ( x) ，应有 ( xk ) 0, i 1,2,, k 1 i （4） ( xk ) 0, k
上式中前 k 1个方程自然满足， 因为 xk 1已经在
span{ p1 ,, pk 1}中极小化了 ( x) .所以只要解出上
l
又因为span{r0 , r1 ,, rl }与span{ p1 , p2 ,, pl 1}的维数都
共轭梯度法课程设计
最优化方法课程设计报告题目：共轭梯度软件设计院（系）：专业：学生姓名：指导教师：题目类型：实验研究工程设计软件开发2010 年1月15 日摘要共轭梯度法最早是由Hestenes 和Stiefle （1952）提出来的，用于解正定系数矩阵的线性方程组，在这个基础上，Fletcher 和Reeves （1964）首先提出了解非线性最优化问题的共轭梯度法。
共轭梯度法是解决无约束非线性最优化问题的重要的方法之一（Conjugate gradient method to solve unconstrained nonlinear optimization problem, one of the important ways.），因为共轭梯度法是介于最速下降法与牛顿法之间的一个方法，它仅需利用一阶导数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hesse 矩阵并求逆的缺点，共轭梯度法不仅是解决大型线性方程组最有用的方法之一，也是解大型非线性最优化最有效的算法一。
共轭梯度对于无约束优化问题Mlnf(x)，x ∈R n 给出一初始值x1，算法选代产生x2，x3,x4,x 5，…．希望某一k X是目标函数解或点收敛于解，在我们这次的运筹学课程设计当中我们正是用这种方法要求解最小问题的最优解的。
关键字：共轭梯度法；无约束优化AbstractConjugate gradient method was first used by Hestenes and Stiefle (1952) put forward for the solution of positive definite coefficient matrix of linear equations, on this basis, Fletcher and Reeves (1964) first put forward about the problem of nonlinear optimization conjugate gradient method.The conjugate gradient method to solve unconstrained nonlinear optimization problems, one of the important ways, as the conjugate gradient method is between steepest descent method and Newton's method between a method, it requires the use of a first-order derivative information, But the steepest descent method to overcome the shortcomings of slow convergence, but also avoid the need to store and calculate Newton's method and the inverse Hesse matrix of the shortcomings of the conjugate gradient method is not only a large-scale linear equations to solve one of the ways the most useful, but also large-scale solution nonlinear optimization algorithm is the most effective one.Conjugate gradient for the unconstrained optimization problem Mlnf (x), x ∈ given an initial value of x1, the election algorithm is generated on behalf of the x2, x3, x4, x 5, .... Hope that is the objective function of a solution or point of convergence in the solution, in our curriculum design, operations research this is exactly what we were using this method requires the smallest solution of the problem the optimal solution.Keywords: conjugate gradient method；unconstrained nonlinear optimization目录一、共轭梯度法的概念....................................................... 错误！未定义书签。
第二章 (2)最速下降法-Newton法-共轭梯度法
取初始向量 x0 1,1T .
解：函数的梯度为
f
(
x)
2
x1 2x2 2 x1 +4x2
4
,
第1次迭代：
f
x0
4
2
,
d 0 = f
x0
4 2
,
x0
+
d
0
=
1+4 1 2
的最大和最小特征值
对于正定二次函数 f (x) 1 xTGx bT x
步长
2
最速下降法的下一个迭代点，
xk 1
xk
gkT gk gkT Ggk
gk
二、最速下降法
其中：
二、最速下降法
得最速下降法的步长
k
gkT gk gkT Ggk
从而下一个迭代点为
xk 1
xk
gkT gk gkT Ggk
+0d
0
=
1 1
+1/4
4 2
=
2 1/2
，
f
x1
1
2
,
第2次迭代：
d1= f
x1
1 2
,
x1
+
d
1
=
2+ 1/2+2
，
()=f x1+d1 =f 2+,1/2+2
1
基于Python共轭梯度法与最速下降法之间的对比
基于Python共轭梯度法与最速下降法之间的对⽐在⼀般问题的优化中，最速下降法和共轭梯度法都是⾮常有⽤的经典⽅法，但最速下降法往往以”之”字形下降，速度较慢，不能很快的达到最优值，共轭梯度法则优于最速下降法，在前⾯的某个⽂章中，我们给出了⽜顿法和最速下降法的⽐较，⽜顿法需要初值点在最优点附近，条件较为苛刻。
算法来源：《数值最优化⽅法》⾼⽴，P111我们选⽤了64维的⼆次函数来作为验证函数，具体参见上书111页。
采⽤的三种⽅法为：共轭梯度⽅法（FR格式）、共轭梯度法（PRP格式）、最速下降法# -*- coding: utf-8 -*-"""Created on Sat Oct 01 15:01:54 2016@author: zhangweiguo"""import sympy,numpyimport mathimport matplotlib.pyplot as plfrom mpl_toolkits.mplot3d import Axes3D as ax3import SD#这个⽂件⾥有最速下降法SD的⽅法，参见前⾯的博客#共轭梯度法FR、PRP两种格式def CG_FR(x0,N,E,f,f_d):X=x0;Y=[];Y_d=[];n = 1ee = f_d(x0)e=(ee[0]**2+ee[1]**2)**0.5d=-f_d(x0)Y.append(f(x0)[0,0]);Y_d.append(e)a=sympy.Symbol('a',real=True)print '第%2s次迭代：e=%f' % (n, e)while n<N and e>E:n=n+1g1=f_d(x0)f1=f(x0+a*f_d(x0))a0=sympy.solve(sympy.diff(f1[0,0],a,1))x0=x0-d*a0X=numpy.c_[X,x0];Y.append(f(x0)[0,0])ee = f_d(x0)e = math.pow(math.pow(ee[0,0],2)+math.pow(ee[1,0],2),0.5)Y_d.append(e)g2=f_d(x0)beta=(numpy.dot(g2.T,g2))/numpy.dot(g1.T,g1)d=-f_d(x0)+beta*dprint '第%2s次迭代：e=%f'%(n,e)return X,Y,Y_ddef CG_PRP(x0,N,E,f,f_d):X=x0;Y=[];Y_d=[];n = 1ee = f_d(x0)e=(ee[0]**2+ee[1]**2)**0.5d=-f_d(x0)Y.append(f(x0)[0,0]);Y_d.append(e)a=sympy.Symbol('a',real=True)print '第%2s次迭代：e=%f' % (n, e)while n<N and e>E:n=n+1g1=f_d(x0)f1=f(x0+a*f_d(x0))a0=sympy.solve(sympy.diff(f1[0,0],a,1))x0=x0-d*a0X=numpy.c_[X,x0];Y.append(f(x0)[0,0])ee = f_d(x0)e = math.pow(math.pow(ee[0,0],2)+math.pow(ee[1,0],2),0.5)Y_d.append(e)g2=f_d(x0)beta=(numpy.dot(g2.T,g2-g1))/numpy.dot(g1.T,g1)d=-f_d(x0)+beta*dprint '第%2s次迭代：e=%f'%(n,e)return X,Y,Y_dif __name__=='__main__':'''G=numpy.array([[21.0,4.0],[4.0,15.0]])#G=numpy.array([[21.0,4.0],[4.0,1.0]])b=numpy.array([[2.0],[3.0]])c=10.0x0=numpy.array([[-10.0],[100.0]])'''m=4T=6*numpy.eye(m)T[0,1]=-1;T[m-1,m-2]=-1for i in xrange(1,m-1):T[i,i+1]=-1T[i,i-1]=-1W=numpy.zeros((m**2,m**2))W[0:m,0:m]=TW[m**2-m:m**2,m**2-m:m**2]=TW[0:m,m:2*m]=-numpy.eye(m)W[m**2-m:m**2,m**2-2*m:m**2-m]=-numpy.eye(m)for i in xrange(1,m-1):W[i*m:(i+1)*m,i*m:(i+1)*m]=TW[i*m:(i+1)*m,i*m+m:(i+1)*m+m]=-numpy.eye(m)W[i*m:(i+1)*m,i*m-m:(i+1)*m-m]=-numpy.eye(m)mm=m**2mmm=m**3G=numpy.zeros((mmm,mmm))G[0:mm,0:mm]=W;G[mmm-mm:mmm,mmm-mm:mmm]=W;G[0:mm,mm:2*mm]=-numpy.eye(mm)G[mmm-mm:mmm,mmm-2*mm:mmm-mm]=-numpy.eye(mm)for i in xrange(1,m-1):G[i*mm:(i+1)*mm,i*mm:(i+1)*mm]=WG[i*mm:(i+1)*mm,i*mm-mm:(i+1)*mm-mm]=-numpy.eye(mm)G[i*mm:(i+1)*mm,i*mm+mm:(i+1)*mm+mm]=-numpy.eye(mm)x_goal=numpy.ones((mmm,1))b=-numpy.dot(G,x_goal)c=0f = lambda x: 0.5 * (numpy.dot(numpy.dot(x.T, G), x)) + numpy.dot(b.T, x) + cf_d = lambda x: numpy.dot(G, x) + bx0=x_goal+numpy.random.rand(mmm,1)*100N=100E=10**(-6)print '共轭梯度PR'X1, Y1, Y_d1=CG_FR(x0,N,E,f,f_d)print '共轭梯度PBR'X2, Y2, Y_d2=CG_PRP(x0,N,E,f,f_d)figure1=pl.figure('trend')n1=len(Y1)n2=len(Y2)x1=numpy.arange(1,n1+1)x2=numpy.arange(1,n2+1)X3, Y3, Y_d3=SD.SD(x0,N,E,f,f_d)n3=len(Y3)x3=range(1,n3+1)pl.semilogy(x3,Y3,'g*',markersize=10,label='SD:'+str(n3))pl.semilogy(x1,Y1,'r*',markersize=10,label='CG-FR:'+str(n1))pl.semilogy(x2,Y2,'b*',markersize=10,label='CG-PRP:'+str(n2))pl.legend()#图像显⽰了三种不同的⽅法各⾃迭代的次数与最优值变化情况，共轭梯度⽅法是明显优于最速下降法的 pl.xlabel('n')pl.ylabel('f(x)')pl.show()最优值变化趋势：从图中可以看出，最速下降法SD的迭代次数是最多的，在与共轭梯度（FR与PRP两种⽅法）的⽐较中，明显较差。
共轭梯度实验报告
竭诚为您提供优质文档/双击可除共轭梯度实验报告篇一：共轭梯度法实验报告数值代数实验报告一、实验名称：用共轭梯度法解线性方程组。
二、实验目的：进一步熟悉理解掌握共轭梯度法解法思路，提高matlab编程能力。
三、实验要求：已知线性方程矩阵，应用共轭梯度法在相关软件编程求解线性方程组的解。
四、实验原理：1．共轭梯度法：考虑线性方程组Ax?b的求解问题，其中A是给定的n阶对称正定矩阵，b是给定的n维向量，x是待求解的n维向量.为此,定义二次泛函?(x)?xTAx?2bTx.定理1设A对称正定,求方程组Ax?b的解，等价于求二次泛函?(x)的极小值点.定理1表明，求解线性方程组问题就转化为求二次泛函?(x)的极小值点问题.求解二次函数极小值问题，通常好像盲人下山那样，先给定一个初始向量x0，确定一个下山方向p0，沿着经过点x0而方向为p0的直线x?x0??p0找一个点x1?x0??0p0，使得对所有实数?有??x0??0p0x0??p0?，即在这条直线上x1使?(x)达到极小.然后从x1出发，再确定一个下山的方向p1，沿着直线x?x1??p1再跨出一步，即找到?1使得??x?在x2?x1??1p1达到极小：??x1??1p1x1??p1?.重复此步骤，得到一串?0,?1,?2,x?xk??pk上确定步长?k使和p0,p1,p2,，称pk为搜索方向，?k为步长.一般情况下，先在xk点找下山方向pk，再在直线??xk??kpkxk??pk?,最后求出xk?1?xk??kpk.然而对不同的搜索方向和步长，得到各种不同的算法.由此，先考虑如何确定?k.设从xk出发，已经选定下山方向pk.令fxk??pk???xk??pk?A?xk??pk??2bT?xk??pk?T??2pkApk?2?rkTpkxk?，T其中rk?b?Apk.由一元函数极值存在的必要条件有Tf2?pkApk?2rkTpk?0所确定的?即为所求步长?k，即步长确定后，即可算出此时，只要rkTpk?0，就有rkTpk.?k?TpkApkxk?1?xk??kpk.??xk?1xkxk??kpkxk?TApk?2?krkTpkk2pk?rkTpk?2即??xk?1xk?.TpkApk?0再考虑如何确定下山方向pk.易知负梯度方向是?(x)减小最快的方向，但简单分析就会发现负梯度方向只是局部最佳的下山方向，而从整体来看并非最佳.故采用新的方法寻求更好的下山方向——共轭梯度法.下面给出共轭梯度法的具体计算过程:给定初始向量x0，第一步仍选用负梯度方向为下山方向，即p0?r0，于是有r0Tr0?0?T,x1?x0??0p0,r1?b?Ax0.p0Ap0对以后各步，例如第k+1步（k?1）,下山方向不再取rk，而是在过点由向量rk和pk?1所张成的二维平面?2?{x|x?xk??rk??pk?1,?,??R}内找出使函数?下降最快的方向作为新的下山方向pk.考虑?在?2上的限制：,?(xk??rk??pk?1)?(xk??rk??pk?1)TA(xk??rk??pk?1)?2bT(xk??rk??pk?1).??计算?关于?,??2??rTAr??rTAp?rTr?,kkkk?1kk????TTT?2?rAp??p?kk?1k?1Apk?1?,r其中最后一式用到了rkpk?1，这可由的定义直接验证.令0k????0，即知?在?2内有唯一的极小值点????x?xk??0rk??0pk?1，其中?0和?0满足??0rkTArk??0rkTApk?1?rkTrk,?TT??0rkApk?1??0pk?1Apk?1?0.1由于rk?0必有?0?0，所以可取pk?作为新的下山方向.显然，这是在平面?2内可得的最佳下山方向.令?k?1?得?0?x?xk??rk??0p?0k?1?0，则可?0rkTApk?1?k?1??T.pk?1Apk?1T注：这样确定的pk满足pkApk?1?0，即pk与pk?1是相互共轭的.总结上面的讨论，可得如下的计算公式：rkTpk，xk?1?xk??kpk，?k?TpkApkrk?1?b?Axk?1，rkT?1Apk，pk?1?rk?1??kpk.?k??TpkApk在实际计算中，常将上述公式进一步简化，从而得到一个形式上更为简单而且对称的计算公式.首先来简化rk?1的计算公式：rk?1?b?Axk?1?b?A(xk??kpk)?rk??kApk.因为Apk在计算?k是已经求出，所以计算rk?1时可以不必将xk?1代入方程计算，而是从递推关系rk?1?b??kApk 得到.再来简化?k和?k的计算公式.此处需要用到关系式rkTrk?1?rkTpk?1?rkT?1pk?0,k?1,2,从而可导出1rkT?1??rkT?1rk?1,，1T?k1TTpkApk?pk?rk?rk?1??pkrk??k1Tk1?rk?rk??k?1pk?1??rkTrk..由此可得?k?krkTrkrkT?1rk?1?k?T,，?k?T..pkApkrkrk从而有求解对称正定方程组的共轭梯度法算法如下：x0?初值r0?b?Ax0；k?0whilerk?0k?k?1ifk?1p0?r0else?k?2?rkT?1rk?1rkT?2rk?2pk?1?rk?1??k?2pk?2endT?k?1?rkT?1rk?1pk?1Apk?1xk?xk?1??k?1pk?1rk?rk?1??k?1Apk?1endx?xk注：该算法每迭代一次仅需要使用系数矩阵A做一次矩阵向量积运算.定理2由共轭梯度法得到的向量组?ri?和?pi?具有如下基本性质：（1）piTrj?0，0?i?j?k;（2）riTrj?0，i?j，0?i,j?k;（3）piTApj?0，i?j，0?i,j?k;（4）span{r0,其中,rk}?span{p0,,pk}??(A,r0,k?1)，?(A,r0,k?1)?span{r0,Ar0,,Akr0}，通常称之为Krylov子空间.下面给出共轭梯度法全局最优性定理：定理3用共轭梯度法计算得到的近似解xk满足??xk??minx?:x?x0??(A,r0,k)?或xk?x*A?min?x?x*A:x?x0??(A,r0,k)?，其中xA?x*是方程组Ax?b的解，?(A,r0,k)是由所定义的Krylov 子空间.定理2表明，向量组r0,,rk和p0,,pk分别是Krylov 子空间?(A,r0,k?1)的正交基和共轭正交基.由此可知，共轭梯度法最多n步便可得到方程组的解x*.因此，理论上来讲，共轭梯度法是直接法.然而实际使用时，由于误差的出现，使rk之间的正交性很快损失，以致于其有限步终止性已不再成立.此外，在实际应用共轭梯度法时，由于一般n很大，以至于迭代o?n?次所耗费的计算时间就已经使用户无法接受了.因此，实际上将共轭梯度法作为一种迭代法使用，而且通常是rk是否已经很小及迭代次数是否已经达到最大允许的迭代次数kmax来终止迭代.从而得到解对称正定线性方程组的实用共轭梯度法，其算法如下：x?初值k?0;r?b?Ax;??rTrwhile??b2?and?k?kmax?k?k?1ifk?1p?relse;p?r??pend??Ap;pT?;x?x??pr?r;;??rTrend算法中，系数矩阵A的作用仅仅是用来由已知向量p产生向量??Ap，这不仅可以充分利用A的稀疏性，而且对某些提供矩阵A较为困难而由已知向量p产生向量??Ap又十分方便的应用问题是十分有益的。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验2 最速下降法和共轭梯度法的程序设计
一、实验目的
1、熟悉无约束优化问题的最速下降算法和共轭梯度法。
2、培养matlab 编程与上机调试能力。
二、实验课时：2个课时
三、实验准备
1、预习无约束优化问题的最速下降算法和共轭梯度法。
2、熟悉matlab 软件的基本操作及程序编写。
四、实验内容
课堂实验演示
根据最速下降法编写程序，求函数
21222121342)(min x x x x x x x f -++-=
的极小值，其中初始点为()01,1T x =
算法步骤如下：
Step1:：给出初始点0x ，和精度1;0k ε<<=；
Step2：计算()k f x ∇，如果()k f x ε∇≤，则停止迭代，输出结果；否则转step3； Step3：令下降方向()k k d f x =-∇，计算步长因子k λ使得0()min ()k k k k k f x d f x d λλλ≥+=+，令1,1k k k k x x d k k λ+=+=+，转step2。
其程序如下：
function [x,iter,val,dval] = Steepest_Descent_Method(x,eps) k = 1;
dy = grad_obj(x);
x_mat(:,1) = x;%存储每一次迭代得到的点x
while norm(dy)>eps
d = -dy; % 搜索方向
lambda = line_search(x,d);%步长
x = x + d*lambda;
k = k + 1;
x_mat(:,k) = x;
dy = grad_obj(x);
end
iter = k - 1;
val = obj(x);%目标函数在极值点处的函数值
dval = grad_obj(x);%目标函数在极值点处的梯度
%-------------------------------------------------------------- x1 = linspace(-1.2,1.2,40);
x2 = linspace(-0.2,1.2,40);
[xx,yy] = meshgrid(x1,x2);
for i = 1:length(x1)
for j = 1:length(x2)
z(i,j) = obj([xx(i,j);yy(i,j)]);
end
end
contour(xx,yy,z,10);%画出目标函数的等高线
hold on
plot(x_mat(1,:),x_mat(2,:),'-o')%画出最速下降法的迭代路径
hold off
function y = obj(x)
%目标函数
y = x(1).^2 - 2*x(1).*x(2) + 4*x(2).^2 + x(1) - 3*x(2);
function dy = grad_obj(x)
%目标函数的梯度
dy = [2*x(1) - 2*x(2) + 1; -2*x(1) + 8*x(2) - 3];
function lambda = line_search(xk,dk)
%作线搜索，求步长
%phi(lambda) = obj( xk + lambda*dk )
%d_phi(lambda) = dk'*grad_obj( xk + lambda*dk )
syms eqn lambda
eqn = dk'*grad_obj(xk+lambda*dk);
lambda = solve(eqn); %用符号计算命令solve 求方程d_phi(\lambda)=0的根 lambda = eval(lambda);%将符号计算的结果转化为数值类型
>> x = [1;1]; eps = 1.0e-6;
>> [x,iter,val,dval] = Steepest_Descent_Method(x,eps) x =[ -0.1667 0.3333]’
val = -0.5833
dval = [0.5280*1.0e-006 -0.1760*1.0e-006]’
iter = 43
共轭梯度法的计算步骤：
Step1: 给出初始点0x ，令)(00x f d -∇=，精度1;0k ε<<=；
Step2：计算()k f x ∇，如果()k f x ε∇≤，则停止迭代，输出结果；否则转step3； Step3：计算k k k k d x x λ+=+1，其中步长因子k λ使得0
()min ()k k k k k f x d f x d λλλ≥+=+，计算下降方向k k k k k d x f x f x f d 22111||
)(||||)(||)(∇∇+-∇=+++； 令1+=k k ，转step2。
其程序如下：
function [x,iter,val,dval] = Conjugate_Gradient_Method(x,eps)
k = 1;
x_mat(:,1) = x;%存储每一次迭代得到的点x
x_old = x;
dy_old = grad_obj(x_old);
d_old = -dy_old;
while norm(dy_old)>eps
lambda = line_search(x_old,d_old);%步长
x_new = x_old + lambda*d_old;
dy_new = grad_obj(x_new);
coef = norm(dy_new)/norm(dy_old);
d_new = -dy_new + coef^2*d_old; % 搜索方向
k = k + 1;
x_old = x_new;
dy_old = dy_new;
d_old = d_new;
x_mat(:,k) = x_new;
%防止死循环
if k > 100
break;
end
end
x = x_new;
iter = k - 1;
val = obj(x_new);%目标函数在极值点处的函数值
dval = grad_obj(x_new);%目标函数在极值点处的梯度
%-------------------------------------------------------------- x1 = linspace(-1.2,1.2,40);
x2 = linspace(-0.2,1.2,40);
[xx,yy] = meshgrid(x1,x2);
for i = 1:length(x1)
for j = 1:length(x2)
z(i,j) = obj([xx(i,j);yy(i,j)]);
end
end
contour(xx,yy,z,10);%画出目标函数的等高线
hold on
plot(x_mat(1,:),x_mat(2,:),'-o')%画出最速下降法的迭代路径
hold off
function y = obj(x)
%目标函数
y = x(1)^2 - 2*x(1)*x(2) + 4*x(2)^2 + x(1) - 3*x(2);
function dy = grad_obj(x)
%目标函数的梯度
dy = [2*x(1) - 2*x(2) + 1; -2*x(1) + 8*x(2) - 3];
function lambda = line_search(xk,dk)
%作线搜索，求步长
%phi(lambda) = obj( xk + lambda*dk )
%d_phi(lambda) = dk'*grad_obj( xk + lambda*dk )
syms eqn lambda
eqn = dk'*grad_obj(xk+lambda*dk);
lambda = solve(eqn); %用符号计算命令solve 求方程d_phi(\lambda)=0的根 lambda = max(eval(lambda));%将符号计算的结果转化为数值类型
课堂实验任务
编写函数文件，实现最速下降法和共轭梯度法，并分别求解下列问题
22
214)(min x x x f +=， 初始点取()01,1T x =，精度取610-=ε；
五、实验主要步骤
1、安装matlab7.0及以上版本软件；
2、编写m 文件以创建和保存各函数；
3、运行程序，保存结果；
4、撰写实验报告。
六、实验报告的撰写要求
1. 写出实验课程名称、日期；
2. 写出姓名、学号；
3. 写出实验目的、实验内容；
4. 写出实验过程及结果（程序代码及数值解），尽量给出其图象;
5. 写出心得体会；。