第六节 最速下降法与共轭梯度法
最优化方法
随机梯度下降每次迭代只使用一个样本，迭代 一次计算量为n 2 ，当样本个数m很大的时候， 随机梯度下降迭代一次的速度要远高于批量梯 度下降方法。 两者的关系可以这样理解：随机 梯度下降方法以损失很小的一部分精确度和增 加一定数量的迭代次数为代价，换取了总体的 优化效率的提升。增加的迭代次数远远小于样 本的数量。
2. 牛顿法和拟牛顿法（Newton's method & Quasi-Newton Methods）
牛顿法（Newton's method） 牛顿法是一种在实数域和复数域上近似求解方程 的方法。方法使用函数 f ( x ) 的泰勒级数的前 面几项来寻找方程 f ( x ) = 0 的根。牛顿法最大 的特点就在于它的收敛速度很快。
具体步骤：
首先，选择一个接近函数 f ( x ) 零点的 x 0 ， 计算相应的 f ( x 0 ) 和切线斜率 f ' (x 0 ) （这 里 f ' 表示函数 f 的导数）。然后我们计算穿 过点 (x 0 , f (x 0 )) 并且斜率为 f '(x 0 ) 的直线 和 x 轴的交点的 x 坐标，也就是求如下方程的 解：
批量梯度下降法（Batch Gradient Descent,BGD）
（1）将J(theta)对theta求偏导，得到每个theta对应 的的梯度：
（2）由于是要最小化风险函数，所以按每个参数 theta的梯度负方向，来更新每个theta：
（3）从上面公式可以注意到，它得到的是一个全 局最优解，但是每迭代一步，都要用到训练集 所有的数据，如果m很大，那么可想而知这种 方法的迭代速度会相当的慢。所以，这就引入 了另外一种方法——随机梯度下降。 对于批量梯度下降法，样本个数m，x为n维向 量，一次迭代需要把m个样本全部带入计算， 迭代一次计算量为m*n 2 。
最速下降法和共轭梯度法
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)
最速下降法与共轭梯度法
最速下降法与共轭梯度法1.最速下降法当方程组Ax = b (1)中的A为对称正定矩阵时，方程组Ax=b的解正好是二次函数(2)的唯一极小值点。
求解方程组(1)的问题等价与求(3)问题。
求解问题(3)的最简单的方法是所谓最速下降法，即从某个初始点x(0)出发，沿φ(x)在点x(0)处的负梯度方向(4)(称为搜索方向)求得φ(x)的极小值点x(1) , 即(5)然后从x(1)出发，重复上面的过程得到x(2)。
如此下去,得到序列{x(k) }(6)可以证明，从任一初始点x(0)出发，用最速下降法所得到的序列{x(k)}均收敛于问题(3)的解，也就是方程组(1)的解。
其收敛速度取决于其中λ1 ,λn分别为A的最小,最大特征值。
最速下降法迭代格式：给定初值x(0) ,x(k)按如下方法决定例8 用最速下降法求解对称正定方程组解：过程如图所示。
2.共轭梯度法共轭梯度法简称CG(Conjugate Gradient),其基本步骤是在点x(k)处选取搜索方向d(k) , 使其与前一次的搜索方向d(k-1)关于A共轭，即<d(k) ,Ad(k-1)> = 0 k=1,2, (7)然后从点x(k)出发，沿方向d(k)求得φ(x)的极小值点x(k+1) , 即(8)如此下去, 得到序列{x(k)}。
不难求得(7)的解为注意到d(k)的选取不唯一，我们可取d(k) = -▽φ(x(k) )+βd(k-1) , 由共轭的k-1定义(7)可得共轭梯度法的计算过程如下：第一步：去初始向量x(0) , 计算第k+1步(k=1,2,…)：计算例8 用共轭梯度法求解对称正定方程组解迭代过程如图所示故x2 = (1,1)T就是φ(x)的最小点，也就是索求方程的解。
共轭梯度法原理
共轭梯度法原理共轭梯度法是线性代数中一种求解稀疏矩阵下的大规模线性方程组的方法。
它的优点是它具有迭代次数小的特点，同时也能得到比较准确的解。
共轭梯度法基于梯度下降法，但是避免了梯度下降法的缺点。
梯度下降法每一次迭代都需计算整个向量的梯度，这在高纬度数据中非常复杂，同时使用步长较大时又容易产生来回震荡的现象。
共轭梯度法的优点是在每一次迭代都会寻找一个与上次方向不同的方向，这点可以消除梯度下降法的缺陷。
令A为若干个线性矩阵的乘积，如果我们要解矩阵方程Ax=b，其中b是已知向量，求解的x向量是未知向量。
首先，我们可以用梯度下降法求出一个初值向量x0，称之为迭代初始值。
然后，我们可以利用高斯打乘法和高斯消元得到向量P，并设向量R0=Ax0-b，然后再设P逆矩阵为Pt。
共轭梯度法迭代的主要步骤如下：1. 根据目标函数和梯度函数确定初值x0；2. 初始化残差向量r0=b-Ax0，并设置迭代数k=0；3. 设置搜索方向向量p0=r0；4. 使用一维搜索方法（如Armijo步长准则）寻找最优步长αk；5. 将搜索方向向量移动到新的位置：xk+1=xk+αkp，计算新的残差rk+1=rk+αkAp；6. 如果rk+1=0或者迭代次数已达到设定值，则停止迭代；否则：7. 确定下一个搜索方向pk+1，并重复步骤4-6直到满足收敛条件。
共轭梯度法迭代的优点是每一步都在原解空间的一个线性子空间中求解，因此不需要保存全部解向量，从而大大减少了计算量和存储空间，特别适用于高纬度的线性方程组的求解。
在参数优化、图像处理和物理模拟等领域中广泛应用。
在参数优化中，共轭梯度法可以加速大规模函数的梯度计算，从而加快参数搜索的速度；在图像处理中，共轭梯度法常用于求解稀疏线性系统，例如图像压缩、图像去噪和图像恢复等；在物理模拟中，共轭梯度法可以用于求解偏微分方程组、流体力学问题和计算机视觉等领域。
李庆扬-数值分析第五版第6章习题答案(20130819)
试考察解此方程组的雅可比迭代法及高斯-赛德尔迭代法的收敛性。 雅可比迭代的收敛条件是
( J ) ( D 1 ( L U )) 1
高斯赛德尔迭代法收敛条件是
(G ) (( D L) 1U ) 1
因此只需要求响应的谱半径即可。 本题仅解 a)，b）的解法类似。 解：
3.设线性方程组
a11 x1 a12 x2 b1 a11 , a12 0 a21 x1 a22 x2 b2
证明解此方程的雅可比迭代法与高斯赛德尔迭代法同时收敛或发散， 并求两种方 法收敛速度之比。 解：
a A 11 a21
则
a12 a22
5. 何谓矩阵 A 严格对角占优？何谓 A 不可约？ P190， 如果 A 的元素满足
aij aij ，i=1,2,3….
j 1 j i
n
称 A 为严格对角占优。 P190 设 A (aij )nn (n 2) ，如果存在置换矩阵 P 使得
A PT AP 11 0
x ( k 1) x ( k )
10 4 时迭代终止。
2 1 5 （a）由系数矩阵 1 4 2 为严格对角占优矩阵可知，使用雅可比、高斯 2 3 10
赛德尔迭代法求解此方程组均收敛。[精确解为 x1 4, x 2 3, x3 2 ] （b）使用雅可比迭代法：
2．给出迭代法 x ( k 1) Bx (k ) f 收敛的充分条件、误差估计及其收敛速度。 迭代矩阵收敛的条件是谱半径 ( B0 ) 1 。其误差估计为
1 k
(k) Bk (0)
R ( B) ln B k 迭代法的平均收敛速度为 k
梯度下降法、牛顿迭代法、共轭梯度法
梯度下降法、牛顿迭代法、共轭梯度法（参见：神经网络->PGM-ANN-2009-C09性能优化）优化的目的是求出目标函数的最大值点或者最小值点，这里讨论的是迭代的方法梯度下降法首先，给定一个初始猜测值 ，然后按照等式k k k k ΡαΧ+=X +1 （1）或kk k k k P =X -X =∆X +α)(1 （2）逐步修改猜测。
这里向量 kP 代表一个搜索方向，一个大于零的纯量kα 为学习速度，它确定了学习步长。
当用 k k k k ΡαΧ+=X +1 进行最优点迭代时，函数应该在每次迭代时都减小，即)()(1k k F F X <X +考虑（3）的)(X F 在k X 的一阶泰勒级数展开：kTk k k k k g F F F ∆X +X ≈∆X +X =X +)()()(1（4）其中，Tk g 为在旧猜测值k X 处的梯度kF g k X =X X ∇≡)( （5） 要使)()(1k k F F X <X +只需要（4）中右端第二项小于0，即<P =∆X k T kk k T k g g α （6）选择较小的正数k α。
这就隐含0<k Tk P g 。
满足0<k Tk P g 的任意向量成为一个下降方向。
如果沿着此方向取足够小步长，函数一定递减。
并且，最速下降的情况发生在k T k P g 最小的时候，容易知道，当k k -g P =时k Tk P g 最小，此时，方向向量与梯度方向相反。
在（1）式中，令k k -g P =，则有k k k k g αΧ-=X +1 （7）对于式（7）中学习速率k α的选取通常有两种方法：一种是选择固定的学习速率k α，另一种方法是使基于学习速率k α的性能指数或目标函数)(1k +X F 在每次迭代中最小化，即沿着梯度反方向实现最小化：k k k k g X X α-=+1。
注意：1、对于较小的学习速度最速下降轨迹的路径总是与轮廓线正交，这是因为梯度与轮廓线总是正交的。
优化方法 2016
优化方法上机大作业上机大作业Ⅰ：编写程序求解下述问题 min x f(x) = (1−x1)2 + 100(x2–x12)2. 初始点取 x0 = 0, 精度取ε=1e−4,步长由 Armijo 线搜索生成, 方向分别由下列方法生成:1 最速下降法2 牛顿法3 BFGS 方法4 共轭梯度法1.最速下降法源程序如下：function x_star = steepest(x0,eps)gk = grad(x0);res = norm(gk);k = 0;while res > eps && k<=1000dk = -gk;ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk);while f1 > f0 + 0.1*ak*slopeak = ak/2;xk = x0 + ak*dk;f1 = fun(xk);endk = k+1;x0 = xk;gk = grad(xk);res = norm(gk);fprintf('--The %d-th iter, the residual is %f\n',k,res); endx_star = xk;endfunction f = fun(x)f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2;endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);end运行结果：>> x0=[0,0]';eps=1e-4eps =1.0000e-004>> xk=steepest(x0,eps)--The 1-th iter, the residual is 3.271712--The 2-th iter, the residual is 2.504194--The 3-th iter, the residual is 2.073282 ……………………………--The 998-th iter, the residual is 0.449447 --The 999-th iter, the residual is 0.449447 --The 1000-th iter, the residual is 0.449447 --The 1001-th iter, the residual is 0.449447 xk =0.63690.40382.牛顿法源程序如下：function x_star = newton(x0,eps)gk = grad(x0);bk = [grad2(x0)]^(-1);res = norm(gk);k = 0;while res > eps && k<=1000dk=-bk*gk;xk=x0+dk;k = k+1;x0 = xk;gk = grad(xk);bk = [grad2(xk)]^(-1);res = norm(gk);fprintf('--The %d-th iter, the residual is %f\n',k,res); endx_star = xk;endfunction f = fun(x)f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2;endfunction g = grad2(x)g = zeros(2,2);g(1,1)=2+400*(3*x(1)^2-x(2));g(1,2)=-400*x(1);g(2,1)=-400*x(1);g(2,2)=200;endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);end运行结果：>> x0=[0,0]';eps=1e-4;>> xk=newton(x0,eps)--The 1-th iter, the residual is 447.213595 --The 2-th iter, the residual is 0.000000 xk =1.00001.00003.BFGS方法源程序如下：function x_star = bfgs(x0,eps)g0 = grad(x0);gk=g0;res = norm(gk);Hk=eye(2);k = 0;while res > eps && k<=1000dk = -Hk*gk;ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk);while f1 > f0 + 0.1*ak*slopeak = ak/2;xk = x0 + ak*dk;f1 = fun(xk);endk = k+1;fa0=xk-x0;x0 = xk;g0=gk;gk = grad(xk);y0=gk-g0;Hk=((eye(2)-fa0*(y0)')/((fa0)'*(y0)))*((eye(2)-(y0)*(fa0)')/((fa0)'*(y0)))+(fa0*(fa0)')/((fa0)'*(y0)); res = norm(gk);fprintf('--The %d-th iter, the residual is %f\n',k,res); endx_star = xk;endfunction f=fun(x)f=(1-x(1))^2 + 100*(x(2)-x(1)^2)^2;endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);end运行结果：>> x0=[0,0]';>>eps=1e-4;>> xk=bfgs(x0,eps)--The 1-th iter, the residual is 3.271712--The 2-th iter, the residual is 2.381565--The 3-th iter, the residual is 3.448742 ……………………………--The 998-th iter, the residual is 0.004690 --The 999-th iter, the residual is 0.008407 --The 1000-th iter, the residual is 0.005314 --The 1001-th iter, the residual is 0.010880 xk =0.99550.99114.共轭梯度法源程序如下：f unction x_star =conj(x0,eps)gk = grad(x0);res = norm(gk);k = 0;dk = -gk;while res > eps && k<=1000ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk); while f1 > f0 + 0.1*ak*slopeak = ak/2;xk = x0 + ak*dk;f1 = fun(xk);endd0=dk;g0=gk;k=k+1;x0=xk;gk=grad(xk);f=(norm(gk)/norm(g0))^2;res=norm(gk);dk=-gk+f*d0;fprintf('--The %d-th iter, the residual is %f\n',k,res);endx_star = xk;endfunction f=fun(x)f=(1-x(1))^2+100*(x(2)-x(1)^2)^2;endfunction g=grad(x)g=zeros(2,1);g(1)=400*x(1)^3-400*x(1)*x(2)+2*x(1)-2;g(2)=-200*x(1)^2+200*x(2);end运行结果：>> x0=[0,0]';>> eps=1e-4;>> xk=Conj(x0,eps)--The 1-th iter, the residual is 3.271712--The 2-th iter, the residual is 1.380542--The 3-th iter, the residual is 4.527780--The 4-th iter, the residual is 0.850596………………………………--The 73-th iter, the residual is 0.001532--The 74-th iter, the residual is 0.000402--The 75-th iter, the residual is 0.000134--The 76-th iter, the residual is 0.000057xk =0.99990.9999上机大作业Ⅱ：编写程序利用增广拉格朗日方法求解下述问题min 4x1−x22−12s.t. 25−x12−x22 = 010x1–x12 + 10x2−x22−34 ≥ 0x 1,x2≥ 0初始点取 x0 = 0, 精度取ε= 1e−4.主程序：function [x,mu,lamda,output]=main(fun,hf,gf,dfun,dhf,dgf,x0)maxk=2000;theta=0.8; eta=2.0;k=0; ink=0;ep=1e-4;sigma=0.4;x=x0; he=feval(hf,x); gi=feval(gf,x);n=length(x); l=length(he); m=length(gi);mu=0.1*ones(l,1); lamda=0.1*ones(m,1);betak=10; betaold=10;while(betak>ep && k<maxk)[ik,x,val]=bfgs('lagrang','dlagrang',x0,fun,hf,gf,dfun,dhf,dgf,mu,lamda,sigma); ink=ink+ik;he=feval(hf,x); gi=feval(gf,x);betak=sqrt(norm(he,2)^2+norm(min(gi,lamda/sigma),2)^2);if betak>epmu=mu-sigma*he;lamda=max(0.0,lamda-sigma*gi);if(k>=2 && betak>theta*betaold)sigma=eta*sigma;endk=k+1;betaold=betak;x0=x;endf=feval(fun,x);output.fval=f;output.iter=k;output.inner_iter=ink;output.beta=betak;增广拉格朗日函数function lag=lagrang(x,fun,hf,gf,dfun,dhf,dgf,mu,lamda,sigma) f=feval(fun,x); he=feval(hf,x); gi=feval(gf,x);l=length(he); m=length(gi);lag=f; s1=0.0;for(i=1:l)lag=lag-he(i)*mu(i);s1=s1+he(i)^2;endlag=lag+0.5*sigma*s1;s2=0.0;for(i=1:m)s3=max(0.0,lamda(i)-sigma*gi(i));s2=s2+s3^2-lamda(i)^2;lag=lag+s2/(2.0*sigma);增广拉格朗日梯度函数function dlag=dlagrang(x,fun,hf,gf,dfun,dhf,dgf,mu,lamda,sigma) dlag=feval(dfun,x);he=feval(hf,x); gi=feval(gf,x);dhe=feval(dhf,x); dgi=feval(dgf,x);l=length(he); m=length(gi);for(i=1:l)dlag=dlag+(sigma*he(i)-mu(i))*dhe(:,i);endfor(i=1:m)dlag=dlag+(sigma*gi(i)-lamda(i))*dgi(:,i);end线搜索程序基于BFGS方法function [k,x,val]=bfgs(fun,gfun,x0,varargin)Max=1000;ep=1.e-4;beta=0.55; sigma1=0.4;n=length(x0); Bk=eye(n);k=0;while(k<Max)gk=feval(gfun,x0,varargin{:});if(norm(gk)<ep), break; enddk=-Bk\gk;m=0; mk=0;while(m<20)newf=feval(fun,x0+beta^m*dk,varargin{:});oldf=feval(fun,x0,varargin{:});if(newf<=oldf+sigma1*beta^m*gk'*dk)mk=m; break;endm=m+1;endx=x0+beta^mk*dk;sk=x-x0;yk=feval(gfun,x,varargin{:})-gk;if(yk'*sk>0)Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk); endk=k+1;x0=x;endval=feval(fun,x0,varargin{:});目标函数文件function f=f1(x)f=4*x(1)-x(2)^2-12;等式约束文件function he=h1(x)he=25-x(1)^2-x(2)^2;不等式约束function gi=g1(x)gi=zeros(3,1);gi(1)=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;gi(2)=x(1);gi(3)=x(2);目标函数梯度文件function g=df1(x)g=[4;-2*x(1)];等式函数梯度文件function dhe=dh1(x)dhe=[-2*x(1);-2*x(2)];不等式函数梯度文件function dgi=dg1(x)dgi=[10-2*x(1),1,0；10-2*x(2)，0,1];输入数据X0=[0;0][x,mu,lamda,output]=main('f1','h1','g1','df1','dh1','dg1',x0) 输出数据x =1.00134.8987mu =0.0158lamda =0.5542output =fval: -31.9924iter: 5inner_iter: 33beta: 8.4937e-005上机大作业Ⅲ：1.解：将目标函数改写为向量形式：x'*a*x-b*x程序代码：n=2;a=[0.5,0;0,1];b=[2 4];c=[1 1];cvx_beginvariable x(n)minimize( x'*a*x-b*x)subject toc * x <= 1x>=0cvx_end运算结果:Calling SDPT3 4.0: 7 variables, 3 equality constraintsFor improved efficiency, SDPT3 is solving the dual problem. ------------------------------------------------------------num. of constraints = 3dim. of socp var = 4, num. of socp blk = 1dim. of linear var = 3******************************************************************* SDPT3: Infeasible path-following algorithms******************************************************************* version predcorr gam expon scale_dataNT 1 0.000 1 0it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime-------------------------------------------------------------------0|0.000|0.000|8.0e-001|6.5e+000|3.1e+002| 1.000000e+001 0.000000e+000|0:0:00| chol 1 11|1.000|0.987|4.3e-007|1.5e-001|1.6e+001| 9.043148e+000 -2.714056e-001|0:0:01| chol 1 12|1.000|1.000|2.6e-007|7.6e-003|1.4e+000| 1.234938e+000 -5.011630e-002|0:0:01| chol 1 13|1.000|1.000|2.4e-007|7.6e-004|3.0e-001| 4.166959e-001 1.181563e-001| 0:0:01| chol 1 14|0.892|0.877|6.4e-008|1.6e-004|5.2e-002| 2.773022e-001 2.265122e-001| 0:0:01| chol 1 15|1.000|1.000|1.0e-008|7.6e-006|1.5e-002| 2.579468e-001 2.427203e-001| 0:0:01| chol 1 16|0.905|0.904|3.1e-009|1.4e-006|2.3e-003| 2.511936e-001 2.488619e-001| 0:0:01| chol 1 17|1.000|1.000|6.1e-009|7.7e-008|6.6e-004| 2.503336e-001 2.496718e-001| 0:0:01| chol 1 18|0.903|0.903|1.8e-009|1.5e-008|1.0e-004| 2.500507e-001 2.499497e-001| 0:0:01| chol 1 19|1.000|1.000|4.9e-010|3.5e-010|2.9e-005| 2.500143e-001 2.499857e-001| 0:0:01| chol 1 110|0.904|0.904|5.7e-011|1.3e-010|4.4e-006| 2.500022e-001 2.499978e-001| 0:0:01| chol 2 211|1.000|1.000|5.2e-013|1.1e-011|1.2e-006| 2.500006e-001 2.499994e-001| 0:0:01| chol 2 212|1.000|1.000|5.9e-013|1.0e-012|1.8e-007| 2.500001e-001 2.499999e-001| 0:0:01| chol 2 213|1.000|1.000|1.7e-012|1.0e-012|4.2e-008| 2.500000e-001 2.500000e-001| 0:0:01| chol 2 214|1.000|1.000|2.3e-012|1.0e-012|7.3e-009| 2.500000e-001 2.500000e-001| 0:0:01|stop: max(relative gap, infeasibilities) < 1.49e-008-------------------------------------------------------------------number of iterations = 14primal objective value = 2.50000004e-001dual objective value = 2.49999996e-001gap := trace(XZ) = 7.29e-009relative gap = 4.86e-009actual relative gap = 4.86e-009rel. primal infeas (scaled problem) = 2.33e-012rel. dual " " " = 1.00e-012rel. primal infeas (unscaled problem) = 0.00e+000rel. dual " " " = 0.00e+000norm(X), norm(y), norm(Z) = 3.2e+000, 1.5e+000, 1.9e+000norm(A), norm(b), norm(C) = 3.9e+000, 4.2e+000, 2.6e+000Total CPU time (secs) = 0.99CPU time per iteration = 0.07termination code = 0DIMACS: 3.3e-012 0.0e+000 1.3e-012 0.0e+000 4.9e-009 4.9e-009-------------------------------------------------------------------------------------------------------------------------------Status: SolvedOptimal value (cvx_optval): -32. 解：将目标函数改写为向量形式：程序代码：n=3;a=[-3 -1 -3];b=[2;5;6];C=[2 1 1;1 2 3;2 2 1];cvx_beginvariable x(n)minimize( a*x)subject toC * x <= bx>=0cvx_end运行结果：Calling SDPT3 4.0: 6 variables, 3 equality constraints------------------------------------------------------------num. of constraints = 3dim. of linear var = 6******************************************************************* SDPT3: Infeasible path-following algorithms******************************************************************* version predcorr gam expon scale_dataNT 1 0.000 1 0it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime-------------------------------------------------------------------0|0.000|0.000|1.1e+001|5.1e+000|6.0e+002|-7.000000e+001 0.000000e+000| 0:0:00| chol 1 11|0.912|1.000|9.4e-001|4.6e-002|6.5e+001|-5.606627e+000 -2.967567e+001| 0:0:00| chol 1 12|1.000|1.000|1.3e-007|4.6e-003|8.5e+000|-2.723981e+000 -1.113509e+001| 0:0:00| chol 1 13|1.000|0.961|2.3e-008|6.2e-004|1.8e+000|-4.348354e+000 -6.122853e+000| 0:0:00| chol 1 14|0.881|1.000|2.2e-008|4.6e-005|3.7e-001|-5.255152e+000 -5.622375e+000| 0:0:00| chol 1 15|0.995|0.962|1.6e-009|6.2e-006|1.5e-002|-5.394782e+000 -5.409213e+000| 0:0:00| chol 1 16|0.989|0.989|2.7e-010|5.2e-007|1.7e-004|-5.399940e+000 -5.400100e+000| 0:0:00| chol 1 17|0.989|0.989|5.3e-011|5.8e-009|1.8e-006|-5.399999e+000 -5.400001e+000| 0:0:00| chol 1 18|1.000|0.994|2.8e-013|4.3e-011|2.7e-008|-5.400000e+000 -5.400000e+000| 0:0:00|stop: max(relative gap, infeasibilities) < 1.49e-008-------------------------------------------------------------------number of iterations = 8primal objective value = -5.39999999e+000dual objective value = -5.40000002e+000gap := trace(XZ) = 2.66e-008relative gap = 2.26e-009actual relative gap = 2.21e-009rel. primal infeas (scaled problem) = 2.77e-013rel. dual " " " = 4.31e-011rel. primal infeas (unscaled problem) = 0.00e+000rel. dual " " " = 0.00e+000norm(X), norm(y), norm(Z) = 4.3e+000, 1.3e+000, 1.9e+000norm(A), norm(b), norm(C) = 6.7e+000, 9.1e+000, 5.4e+000Total CPU time (secs) = 0.11CPU time per iteration = 0.01termination code = 0DIMACS: 3.6e-013 0.0e+000 5.8e-011 0.0e+000 2.2e-009 2.3e-009 -------------------------------------------------------------------------------------------------------------------------------Status: SolvedOptimal value (cvx_optval): -5.4。
梯度法最速下降法
4. 令 x k 1 x k k d k , 令 k : k 1 , 转2。
2 2 例. 用最速下降法求解 : min f ( x ) x1 3 x2 , 设初始点为 x 1 ( 2 , 1 )T ,
求迭代一次后的迭代点 x 2 。
解： f ( x ) ( 2 x1 , 6 x 2 )T ,
n 次迭代必可得到最优解 。
如何选取一组共轭方向？
2. 共轭梯度法
Fletcher R eeves 共轭梯度法:
1 T x Ax bT x c 2 其中 x R n , A是对称正定矩阵， b R n， c 是常数。 min f ( x)
基本思想： 将共轭性和最速下降方 向相结合，利用已知迭 代点
向量，则这个向量组线性无关。
证明
设存在实数 1 , 2 ,, k ，使得
i 1
id 0,
T
k
i
上式两边同时左乘 d j A ，则有
i 1 k
id
k
jT
Ad i 0 ,
因为 d 1 , d 2 ,, d 是 k 个 A 共轭的向量，所以上式可化简为
jd
j
jT
Ad j 0 .
d 1 f ( x 1 ) ( 4 , 6 )T . x 1 d 1 ( 2 4 , 1 6 )T . 令 ( ) f ( x 1 d 1 ) ( 2 4 ) 2 3 ( 1 6 ) 2 ，
求解
min ( )
d (1)T Ad ( 2) 0,
即等值面上一点处的切 向量与由这一点指向极小点的向量关于A 共轭。
1 T x Ax bT x c ， 2 其中 A 是 n 阶对称正定矩阵。 d (1) , d ( 2 ) ,, d ( k ) 是 一组A共轭向量。
最速下降法与共轭梯度法实验
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. 方向选择：最速下降法每次迭代选择梯度的方向作为搜索方向，而共轭梯度法每次迭代选择一个共轭方向，不同的迭代步骤之间方向是相互垂直的。
2. 收敛速度：共轭梯度法通常比最速下降法更快地收敛到最优解。
共轭梯度法利用了问题的特殊结构，通过选择共轭方向进行迭代，可以跳过一些搜索方向，从而加快收敛速度。
3. 存储需求：最速下降法仅需要存储当前点的梯度信息，而共轭梯度法需要额外存储用于计算共轭方向的历史梯度信息，因此共轭梯度法的内存需求更高。
4. 线性搜索：最速下降法一般使用线性搜索确定每次迭代的步长，而共轭梯度法可以利用解析公式或其他更高效的非线性搜索方法来选择合适的步长。
总体来说，共轭梯度法在收敛速度和存储需求方面有优势，特别适用于求解大规模优化问题。
但最速下降法相对简单，并且对一些非二次凸函数仍然有效，因此在某些情况下也是一种有效的选择。
ch6 最速下降法与共轭梯度法
例:二维正定二次函数 二维正定二次函数
2 min f ( X ) = 2 x12 + x2 + 2 x1 x2 + x1 − x2
∇f ( X ) = (4 x1 + 2 x2 + 1, 2 x1 + 2 x2 − 1)T , 4 2 H(X ) = 2 2
二次终结性
3.2 基本原理
∇f ( X ) = AX + B ⇒ ∇f ( X ( k +1) ) − ∇f ( X ( k ) ) = A( X ( k +1) − X ( k ) ) 而X ( k +1) = X ( k ) + λk P ( k ) ∇f ( X
( k +1)
) − ∇f ( X
(k )
1.2 算法步骤 （1）给定初始近似点x(0)及精度ε>0， ）给定初始近似点 及精度 ，
即为近似极小点。 若 ∇f ( X ) ≤ ，则x (0)即为近似极小点。 ε
(0) 2
(2) 若 ∇f ( X ) > ，求步长 0，并计算 ε 求步长λ
(0)
2
X (1) = X (0) −λ∇f ( X (0) ) 0
(0)
2
X
(k +1)
=X
(k )
−λ∇f ( X ) k
(k )
如此继续， 如此继续，直至达到要求的精度为止。
具有二阶连续偏导数， 若f(x)具有二阶连续偏导数，在x(k)作 f ( X 具有二阶连续偏导数 的泰勒展开： 的泰勒展开：
(k )
−λ∇f ( X (k ) ))
f ( X (k) −λ∇f ( X (k ) )) ≈ f ( X (k ) ) −∇f ( X (k ) )Tλ∇f ( X (k) ) 1 + λ∇f ( X (k ) )T H( X (k ) ) ∇f ( X (k ) ) λ 2
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}的维数都
共轭梯度法的原理
共轭梯度法的原理共轭梯度法是一种用于求解线性方程组的迭代方法，其原理是利用共轭方向的特性来加速收敛速度。
在实际问题中，线性方程组的求解是一项非常重要的任务，经常出现在科学计算、工程设计、物理模拟等领域。
共轭梯度法通过迭代的方式，逐步逼近方程组的解，从而提高求解效率和精度。
共轭梯度法的主要思想是将待求解的线性方程组转化为一个优化问题，通过最小化目标函数来找到最优解。
具体来说，假设有一个n 维向量x和一个n×n的矩阵A，以及一个n维向量b，我们需要求解Ax=b。
共轭梯度法的目标是找到一个n维向量x，使得目标函数f(x)最小化，其中f(x)=(1/2)x^TAx-x^Tb。
共轭梯度法的基本步骤如下：1. 初始化向量x0和残差r0=b-Ax0，设置迭代次数k=0。
2. 计算搜索方向d，d=rk+(k-1)βk-1d(k-1)，其中βk-1=(rk^Trk)/(rk-1^Trk-1)。
3. 计算步长αk，αk=(rk^Trk)/(d^TAd)。
4. 更新xk=xk-1+αkd，更新残差rk=rk-1-αkAd。
5. 如果残差rk的范数小于给定的收敛条件，或达到最大迭代次数，则停止迭代；否则，返回步骤2。
共轭梯度法的优点是具有快速收敛和内存占用少的特点。
由于每次迭代都在共轭方向上前进，可以避免了梯度下降法中的zigzag现象，从而加快了收敛速度。
此外，共轭梯度法只需要存储当前迭代的向量和矩阵乘积，而不需要存储所有的迭代历史，因此内存占用较小。
然而，共轭梯度法也存在一些限制和注意事项。
首先，共轭梯度法只适用于对称正定的矩阵A，否则可能出现收敛失败的情况。
其次，共轭梯度法对初始解的选择较为敏感，不同的初始解可能导致不同的收敛速度和精度。
此外，共轭梯度法在求解大规模问题时可能需要较长的迭代时间，因此对于大规模问题，需要考虑使用并行计算的方法来加速求解过程。
共轭梯度法是一种有效的线性方程组求解方法，通过利用共轭方向的特性，可以加快收敛速度和提高求解精度。
最速下降法(sd);共轭梯度法
最速下降法(sd);共轭梯度法
最速下降法(SD)和共轭梯度法(CG)都是求解非线性优化问题中的常用算法。
最速下降法是基于梯度方向的一种搜索方法，在每一步所需找到函数在当前点的最陡方向，并沿着该方向走一步，直到达到要求的精度为止。
该方法速度快，收敛性好，但容易陷入“zigzag”现象，即由于步长过大或过小，导致序列在搜索方向上反复飞奔而收敛缓慢，同时，最速下降法对函数“弯曲性”敏感，函数梯度变化太快时收敛缓慢。
共轭梯度法是一种基于梯度方向的线性搜索方法，其优势在于快速收敛，准确性高。
其核心思想是，由于函数在一般性条件下不是QUADRATIC FUNCTION，因此，图像往往不是一个明显的"碗状"，而是一个复杂的非线性图形。
在这种情况下，最速下降法很容易落入“zigzag”现象，收敛速度慢。
而共轭梯度法可以从不同方向进行极小值点的搜索，进而明显提高收敛速度。
总之，最速下降法适用于方向比较简单的情况，而共轭梯度法适用于方向较为复杂的情况。
根据不同的情况进行选择，可以有效地提高求解的效率和精度。
基于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两种⽅法）的⽐较中，明显较差。
下降方向
戴宽平一．关于下降法的下降方向问题考虑如下方程=Ax b （1）其中n n⨯∈A R是可逆方阵，且其对角线元素全部非零，列向量n∈b R 。
由此我们建立一个求解方程组（1）的变分方程12,,ϕ()=()-()x x Ax x b （2）其中x 是n 维列向量。
我们可以得到：对于n*∈x R 满足*=Ax b 的充分必要条件是nϕϕ*∈()=m i n ()x R x x 。
这样我们就把问题转化为求解ϕ()x 的极小值问题。
为了找到ϕ()x 的极小值点*x ，我们从人一k x 出发，沿着ϕ()x 在k x 点沿着某一个制定的方向搜索下一个近似点1k +x ，使得1k ϕ+()x 在该方向上达到最小值。
在最速下降法中，我们会选择ϕ()x 在k x 点下降最快的方向k z ，即在这点的负梯度方向k r 。
对于共轭梯度法，我们的选取方式为()()00r =z ，()()()1k k k k r β-=+z z ，12k ,...=()，并且()()()10k k ,-=z Az ，即可求得()()()()()()1111k k kk k ,,β----=-z Ar zAz 。
但对于非正定矩阵A ，我们不能保证()()()1kk ,-z Az 是非零的，所以我们将用矩阵D 代替公式中的A 。
D 的对角线上的元素为对应的A 的对角线元素的绝对值，其他元素全部为零。
于是我们可以得到一个新的算法如下：()()()()()()()()()()()()()()()()()()()()()()()()00000000010001111112k k k k k k k k k k k ,,,,,,,,,,r k ,...ααββ-----⎧=-⎪=⎪⎪⎪=⎪⎪⎪⎪=+⎨⎪=-⎪⎪⎪=-⎪⎪⎪=+=⎪⎩r b Ax z r z r z Dz x x r r b Ax z Dr z Dz z z ()()()()()()()()()()1k k k k kk k k k ,,αα+⎧⎪=⎪⎨⎪⎪=+⎩z r z Dz xx z 当12k ke +-<x x 成立时，则停止计算，其中e 为终止准则。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第六节最速下降法与共轭梯度法
6.1 最速下降法
当方程组
Ax = b (1)
中的A为对称正定矩阵时，方程组Ax=b的解正好是二次函数
(2)
的唯一极小值点。
求解方程组(1)的问题等价与求
(3)
问题。
求解问题(3)的最简单的方法是所谓最速下降法，即从某个初始点x(0)出发，沿φ(x)在点x(0)处的负梯度方向
(4)
(称为搜索方向)求得φ(x)的极小值点x(1) , 即
(5)
然后从x(1)出发，重复上面的过程得到x(2)。
如此下去,得到序列{x(k) }
(6)
可以证明，从任一初始点x(0)出发，用最速下降法所得到的序列{x(k)}均收敛于问题(3)的解，也就是方程组(1)的解。
其收敛速度取决于
其中λ1 ,λn分别为A的最小,最大特征值。
最速下降法迭代格式：给定初值x(0) ,x(k)按如下方法决定
对称正定方程组
解：过程如图所示。
6.2 共轭梯度法
共轭梯度法简称CG(Conjugate Gradient),其基本步骤是在点x(k)处选取搜索方向d(k) , 使其与前一次的搜索方向d(k-1)关于A共轭，即
<d(k) ,Ad(k-1)> = 0 k=1,2, (7)
然后从点x(k)出发，沿方向d(k)求得φ(x)的极小值点x(k+1) , 即
(8)
如此下去, 得到序列{x(k)}。
不难求得(7)的解为
注意到d(k)的选取不唯一，我们可取d(k) = -▽φ(x(k4) )+βk-1 d(k-1) , 由共轭的定义(7)可得
共轭梯度法的计算过程如下：
第一步：去初始向量x(0) , 计算
第k+1步(k=1,2,…)：计算
例8 用共轭梯度法求解对称正定方程组
解
迭代过程如图所示
故x2 = (1,1)T就是φ(x)的最小点，也就是索求方程的解。