数值分析11(共轭梯度法)

初等变分原理最速下降法共轭梯度法数值试验算例《数值分析》11*(x,x)≥0,当且仅当x=0时等号成立;(x,y)=(y,x);(kx+ly,z)=k(x,z)+l(y,z);(x,y)=||x||2||y||2cos<x,y>。预备知识设

,则实数(x,y)=xTy=x1y1+…+xnyn*设A是n阶对称正定阵(Ax,x)≥0,当且仅当x=0时等号成立;(Ax,y)=(x,Ay)=(Ay,x)=(y,Ax);(kAx+lAy,z)=k(Ax,z)+l(Ay,z).预备知识矩阵A正定,如果对于任意非零向量x满足xTAx>0.*预备知识例如

f(x1,x2,x3)=x12x22x32梯度(多元函数的一阶导数信息):*思考:多元函数的二阶导数信息？预备知识泰勒展式(数值分析的基石):*单变量函数极值点(费马引理):*注释:费马引理的价值在于将极值问题转化为非线性方程的求解问题。*多变量函数极值点:定理4.10(初等变分原理)

设A=(aij)n×n为实对称正定矩阵,则

x是二次函数

的极小值点x是线性方程组

Ax=b的解。

证明:设

u是

Ax=b的解

Au=b对任意x∈Rn,只须证明

f(x)–f(u)≥0*

设u是

f(x)极小值点。取非零向量

x∈Rn,

对任意

t∈R,有当t=0时,g(0)=f(u)达到极小值,所以

g′(0)=0,即(Au–b,x)=0Au–b=0所以u是方程组

Ax=b

的解。*最速下降方向从初值点x(0)

出发,以负梯度方向

r

为搜索方向在

x处,梯度方向是

f(x)增长最快方向负梯度方向是

f(x)下降最快方向选择步长t0,使x(1)

=x(0)+t0r

为f(x)极小值点最速下降方向:r=–f=b

–

Ax*求解得

t0=(r0

,r0)/(Ar0,r0)为选取最佳步长

t0,令取初值点

x(0),取负梯度方向

r0=b–Ax(0)求点:

x(1)=x(0)+t0r0使得记*精确搜索步长解对称正定方程组Ax=b的最速下降算法:第一步:取初值

x(0)∈R(n),>0,计算

r0=b–Ax(0),k

0;第二步:计算

tk=(rk,rk)/(Ark,rk)

x(k+1)=x(k)+tkrk,rk+1=b–Ax(k+1)第三步:kk+1,如果

||rk||≥

,转第二步;

否则输出

x(k),结束。*注释:

最速下降算法思想简单且容易实现,是求解无约束优化问题的经典方法。最好+

最好

=

最好

?

方向(最速下降)(best

rk)步长(精确搜索)(besttk)x(k+1)=x(k)+tkrk

是否最好?*Rosenbrock方程

f(x1,x2)=(1-x1)2+100(x2-x12)2*Rosenbrock方程

f(x1,x2)=(1-x1)2+100(x2-x12)2*Barzilai-Borwein方法

方向(最速下降-最好方向)

(best

rk)步长(上一次精确搜索步长)

(besttk-1)最好的rk+上一步最好的tk-1

最好参考文献:Two-PointStepSizeGradientMethod*f(x1,x2)=100x12+x22最速下降法*f(x1,x2)=100x12+x22BB方法*最速下降法思想简单,但是收敛速度慢。本质上是因为负梯度方向函数下降快是局部性质。全局思想:局部思想:

在n维空间中,任意n个线性无关的向量构成n维空间的基。换句话说,n维空间中任意向量均可以由这组基线性表示。**

共轭梯度法的关键是构造一组两两共轭的方向(即张成n维空间的基)。巧妙的是共轭方向可以由上次搜索方向和当前的梯度方向线性组合产生。pk+1:=rk+1+beta*pk

TheBestofthe20thCentury:EditorsNameTop10Algorithms,SIAMNews现代迭代方法:子空间方法——共轭——A是n阶对称正定矩阵,非零向量

p1,p2∈Rn(Ap1,p2)=0n个向量

p0,p1,···,pn-1共轭:(Api,pj)=0(i≠j;i,j=0,1,···,n-1)非零向量p0,p1,···,pn-1∈Rnp0,p1,···,pn-1关于A共轭

p0,p1,···,pn-1线性无关两个向量

p1,p2关于A共轭:*——系数计算——更加整体地考虑搜索方向的选择,选择一组关于A共轭的向量:n个向量

p0,p1,···,pn-1共轭:(Api,pj)=0(i≠j;i,j=0,1,···,n-1)*——共轭向量组构造——思路:Gram-Schmidt正交化过程*第一步:取初值

x(0)=0,>0,计算

r0=b–Ax(0),若||

r0||≤

结束;

否则p0r0,k1,转第二步;原始共轭梯度算法第二步:计算

=(pk-1,b

)/(Apk-1,pk-1)

x(k)=x(k-1)+pk-1

;(张成k维子空间)

第三步:如果k=n,则结束;

否则计算

rk=b–Ax(k),转第四步;*第四步:如果

||rk||≤

,则结束;否则计算=–(Apj,rk)/(Apj,pj),(j=0,···,k-1)pk=rk+(

p0+···+pk-1

)

kk+1,转第二步。定理4.12A是n阶对称正定矩阵,p0,p1,···,pn-1

是关于A共轭的向量组,任取

x(0)∈Rn,计算=(pk-1,b)/(Apk-1,pk-1)x(k)=x(k–1)+pk-1

(k=1,2,···,n)则有

Ax(n)=b。共轭梯度方法理论上有限步终止,故仅被视为直接法,所以在其后的20年没有受到重视。1972年共轭梯度法被首次作为迭代法研究,很有新意。第k步迭代,p0,p1,···,pk-1张成k维子空间,故此类方法称为子空间方法。*共轭梯度算法(TheConjugateGradientAlgorithm)

1.Start:x0,r0:=b-Ax0,p0:=r0,k:=0.2.Iterate:Untilconvergencedo,(a)alpha:=(rk,rk)/(Apk,pk)(b)xk+1:=xk+alpha*

pk(c)rk+1:=rk–alpha*Apk(d)beta:=(rk+1,rk+1)/(rk,rk)(e)pk+1:=rk+1+beta*pk

k:=k+1*例1用共轭梯度法求解例1用共轭梯度法求解程序片段:MatlabCode:CG方法function[x,relres,iter]=conjgrad(A,b,x,nmax,tol)res0=norm(b-A*x);iter=0;relres=1;res=b;p=res;rho1=res'*res;whilerelres>tol&&iter<nmaxrho=rho1;q=A*p;alpha=rho/(q'*p);x=x+alpha*p;res=res-alpha*q;rho1=res'*res;beta=rho1/rho;p=res+beta*p;iter=iter+1;relres=norm(res)/res0;end*算例1Possion方程:令

h=1/(n+1),xj=jh,yj=jh(i,j=0,1,···,n+1

)记

ui,j=u(xi,yj),(i,j=0,1,···,n+1)线性方程组(i,j=1,···,n)u0,j=0,ui,0

=0,ui,n+1

=0(i,j=1,···,n)算例2

A=gallery('poisson',n);重要性质:*证明的细节:x(0)=