线性方程组的最速下降法与共轭梯度法

共轭梯度法

一 共轭梯度法原理

对于线性方程组A b X =，即：

1111221n 12112222n 21122nn n n n n n n

a x a x a x

b a x a x a x b a x a x a x b +++=⎧⎪+++=⎪⎨⎪⎪++

+=

⎩ （1）

其中，()=ij n n a ⨯A 为对称正定矩阵，()1i n b b ⨯=，如何熟练地运

用最速下降法与共轭梯度法的求解线性方程组。

在求解线性方程组之前，首先用内积将问题转化为函数问题。 1 最速下降法

最速下降法是一种运用梯度与极值的性质，综合数值计算方法寻找局部极值。

基本思想：任一点的负梯度方向是函数值在该点下降最快的方向。将n 维问题转化为一系列沿负梯度方向用一维搜索方法寻优的问题，利用负梯度作为搜索方向，故称最速下降法。 具体步骤：

1、搜索方向：()k k d f x =-∇，即最速下降方向。

2、搜索步长：k λ取最优步长，即满足：

()min ()k k k k k f x d f x d λ

λλ+=+

1 给定初始点0n x R ∈，允许误差0ε≥，令1k =。

2 计算搜索方向()k k d f x =-∇。

3 若k d ε≤，则k x 为所求的极值点，否则，求解最优步长k λ，使得()min ()k k k k k f x d f x d λ

λλ+=+。

4 令1k k k k x x d λ+=+，1k k =+

最速下降方向是反映了目标函数的局部性质，它只是局部目标函数值下降最快的方向。 2 共轭梯度法

对于1

min ()2T T f x x Ax b x =+

其中，0n x R ∈，A 是对称正定矩阵。

基本思想：将共轭性与最速下降法相结合利用已知迭代点的梯度方向构造一组共轭方向，并沿此方向搜索，求出函数的极小值。 具体步骤：

1 取初始点(0)x ，取第一次搜索方向为(0)(0)()d f x =-∇。

2 设已求得(1)

k x

+，若(1)

()0k f x +∇≠，令(1)

()()k g x f x +=∇

，则下一个

搜索方向

(1)()1k k k k d g d β++=-+ （1）

由于(1)

k d +与()

k d

关于A 共轭，所以给（1）两边同时乘以()T

k d

A ，

即：

()(1)()()()10T

T

T

k k k k k k k d d d g d d β++A =-A +A =

解得：()1()()

k T k k k T k d A g d Ad β+= （2）

3 搜索步长的确定，已知迭代点()k x ，和搜索方向()k d ，确定

步长k λ，即：()()min ()k k f x d λ

λ+

记

()()()()k k f x d φλλ=+，

令

()()()()()0k k T k f x d d φλλ'=∇+=

既有：()()()[()]0k k T k A x d b d λ++=

令 ()()()k k k g f x Ax b =∇=+ 既有：()()[]0k T k k g Ad d λ+=

解得：()

()()

T k k k k T k g d d Ad

λ=-

共轭梯度法是对最速下降法的一种改进，减少了迭代次数从而提高了程序运行效率。

二 程序框图

程序：

%%%%%%%%%%%%%%%%%%%最速下降法%%%%%%%%%%%%%%%%%%

function [x,k]=fast(A,b)

esp=input('ÇëÊäÈëÔÊÐíÎó²îesp=');

N=input('ÇëÊäÈë×î´óµü´ú´ÎÊýN=');

x0=input('ÇëÊäÈë³õʼֵx0=');

k=0;

tol=1;

while tol>=esp

r=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>=N

disp('µü´ú´ÎÊýÌ«¶à£¬¿ÉÄܲ»ÊÕÁ²£¡');

return;

end

end

x

k

%%%%%%%%%%%%%%%%%%%共轭梯度法%%%%%%%%%%%%%%%%%%%

function [k,x]=C_G(A,b)

esp=input('请输入最大误差=');

x0=input('请输入初值x0=');

k = 0 ;

r0 = b-A*x0; %Çó³ödangqianÌݶÈ

while norm(r0)>esp

r0 = b -A*x0;

k = k + 1 ;

if k==1

p0 = r0 ;

else

lamda=(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;

end

end

x=r0;

k;

end