第六节 最速下降法与共轭梯度法
最优化方法
随机梯度下降每次迭代只使用一个样本，迭代 一次计算量为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共轭向量。
- 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)的最小点，也就是索求方程的解。