数值代数实验报告
一、实验名称：用共轭梯度法解线性方程组。
二、实验目的：进一步熟悉理解掌握共轭梯度法解法思路，提高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+α0p0)≤ϕ(x0+αp0)，
即在这条直线上x1使ϕ(x)达到极小.然后从x1出发，再确定一个下山的方向p1，沿着直
线x=x1+αp1再跨出一步，即找到α1使得ϕ(x)在x2=x1+α1p1达到极小：
ϕ(x1+α1p1)≤ϕ(x1+αp1).
重复此步骤，得到一串
α0,α1,α2,
x=xk+αpk上确定步长αk使
和 p0,p1,p2,
，
称pk为搜索方向，αk为步长.一般情况下，先在xk点找下山方向pk，再在直线
ϕ(xk+αkpk)≤ϕ(xk+αpk),
最后求出xk+1=xk+αkpk.然而对不同的搜索方向和步长，得到各种不同的算法.
由此，先考虑如何确定αk.设从xk出发，已经选定下山方向pk.令 f(α)=ϕ(xk+αpk)
=(xk+αpk)A(xk+αpk)-2bT(xk+αpk)
T
=α2pkApk-2αrkTpk+ϕ(xk)，
T
其中rk=b-Apk.由一元函数极值存在的必要条件有
T
f'(α)=2αpkApk-2rkTpk=0
所确定的α即为所求步长αk，即 步长确定后，即可算出
此时，只要rkTpk≠0，就有
rkTpk
. αk=T
pkApk
xk+1=xk+αkpk.
ϕ(xk+1)-ϕ(xk)=ϕ(xk+αkpk)-ϕ(xk)
T
Apk-2αkrkTpk=- =αk2pk
(rkTpk)
2
即ϕ(xk+1)
T
pkApk
再考虑如何确定下山方向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+
η0
p ξ0k-1
η0
，则可ξ0
rkTApk-1
βk-1=-T.
pk-1Apk-1
T
注：这样确定的pk满足pkApk-1=0，即pk与pk-1是相互共轭的.
总结上面的讨论，可得如下的计算公式：
rkTpk
， xk+1=xk+αkpk， αk=T
pkApk
rk+1=b-Axk+1，
rkT+1Apk
， pk+1=rk+1+βkpk. βk=-T
pkApk
在实际计算中，常将上述公式进一步简化，从而得到一个形式上更为简单而且对称的计算公式.首先来简化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,
从而可导出
1
rkT+1=-rkT+1rk+1,， 1Tαk1TT
pkApk=pk(rk-rk+1)=pkrk
ααk1Tk1
=rk(rk+βk-1pk-1)=rkTrk.
.
由此可得
αkαk
rkTrkrkT+1rk+1αk=T,， βk=T..
pkApkrkrk
从而有求解对称正定方程组的共轭梯度法算法如下： x0=初值
r0=b-Ax0；k=0 while rk≠0
k=k+1
if k=1 p0=r0
else
βk-2=rkT-1rk-1rkT-2rk-2 pk-1=rk-1+βk-2pk-2 end
T
αk-1=rkT-1rk-1pk-1Apk-1
xk=xk-1+αk-1pk-1 rk=rk-1-αk-1Apk-1 end
x=xk
注：该算法每迭代一次仅需要使用系数矩阵A做一次矩阵向量积运算. 定理2 由共轭梯度法得到的向量组{ri}和{pi}具有如下基本性质： （1）piTrj=0， 0≤i
,rk}=span{p0,,pk}=κ(A,r0,k+1)，
κ(A,r0,k+1)=span{r0,Ar0,,Akr0}，
通常称之为Krylov子空间.下面给出共轭梯度法全局最优性定理： 定理3 用共轭梯度法计算得到的近似解xk满足
ϕ(xk)=min{ϕ(x):x∈x0+κ(A,r0,k)}
或
xk-x*
A
=min{x-x*
A
:x∈x0+κ(A,r0,k)}，
其中x
A
=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;ρ=rTr
while
>εb
2
)and(k
max
)
k=k+1
if k=1 p=r else
β=ρρ;p=r+βp end
ω=Ap;α=ρpTω;x=x+αp r=r-αω;ρ=ρ;ρ=rTr end
算法中，系数矩阵A的作用仅仅是用来由已知向量p产生向量ω=Ap，这不仅可以充分利用A的稀疏性，而且对某些提供矩阵A较为困难而由已知向量p产生向量
ω=Ap又十分方便的应用问题是十分有益的。
2．算例
⎡10 1 2 3 4⎤⎡x1⎤⎡ 12⎤⎢1 9 -1 2 -3⎥⎢x⎥⎢-27⎥运用共轭梯度法求解下述方程，并解释你所观察到的结果 ⎢⎥⎢2⎥⎢⎥⎢2 -1 7 3 -5⎥⎢x3⎥=⎢ 14⎥ ⎢⎥⎢⎥⎢⎥3 2 3 12 -1-17x解：已知共轭梯度法的⎢⎥⎢4⎥⎢MATLAB⎥程序代码如下所示：
⎥⎣ 12⎥⎢⎣4 -3 -5 -1 15⎥⎦⎢⎦⎣x5⎦⎢function[x,n]=conjgrad(A,b,x0) %采用共轭梯度法求线性方程组Ax=b的解 %线性方程组的系数矩阵：A %线性方程组中的常数向量：b %迭代初始向量：x0 %线性方程组的解：x
%求出所需精度的解实际的迭代步数：n r1=b-A*x0; p=r1; n=0;
for i=1:rank(A)
if(dot(p,A*p)
alpha=dot(r1,r1)/dot(p,A*p);
x=x0+alpha*p; r2=r1-alpha*A*p; if(r2
belta=dot(r2,r2)/dot(r1,r1); p=r2+belta*p; n=n+1; end
用共轭梯度法求解，在MATLAB命令窗口中输入求解程序：
>> A=[10 1 2 3 4 2 4 2 1 1;1 9 -1 2 -3 3 -3 -1 9 2;2 -1 7 3 -5 4 -5 7 -1 3;3 2 3 12 -1 5 -1 3 2 4;4 -3 -5 -1 15 1 3 2 3 2;2 3 4 5 1 10 4 5 3 4;4 -3 -5 -1 3 4 9 1 -3 2;2 -1 7 3 2 5 1 7 -1 1;1 9 -1 2 -3 3 -3 -1 12 10;1 2 3 4 2 4 2 1 10 15];
>> b=[12;-27;14;-17;12;12;-27;14;-17;12]; >> x0=[0;0;0;0;0;0;0;0;0;0]; >> [x,n]=conjgrad(A,b,x0) 输出的计算结果为： x = 8.5669 -7.1981 -7.3479 -13.9388 1.1303 11.5188 -26.8966 -2.2388 0.0829 13.2786
输出的迭代次数为 n = 10
共轭梯度法理论上只要迭代n步，就能得出方程组的解，但是由于存在计算误差，即分数向小数转化时存在舍入误差，很难保证在第n步时得到准确解。
3． 总结
本文首先给出最小二乘问题的定义，随后从盲人下山法开始讨论了共轭梯度法的具体推导过程及其相关性质与算法.继而重点给出正则化方法的实用共轭梯度算法并举例进行检验.最后，需要说明虽然正则化方法是求一般线性方程组Ax=b，
A∈Rm⨯n(m≥n)且A列满秩的最小二乘解的一种方法且简单易行，但是也有许多不
足之处，如m>n时一般无解；ATA形成时运算量大，A中某些信息会丢失；当A病
2
态时其收敛性速度由于κ2(ATA)=κ2(A)很大变得非常之慢等，故为了避免正则化方
法的缺点，还可运用残量极小化方法或残量正交化方法等更好的方法来解决此类问题.
在实际的科学与工程问题中，常常将问题归结为一个线性方程组的求解问题，而求解线性方程组的数值解法大体上可分为直接法和迭代法两大类.直接法是指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法.因此，直接法又称为精确法.迭代法则是采取逐次逼近的方法，亦即从一个初始向量出发，按照一定的计算格式，构造一个向量的无穷序列,其极限才是方程组的精确解，只经过有限次运算得不到精确解.当线性方程组的系数矩阵为对称正定矩阵时，我们常用共轭梯度法（或简称CG法）求解，目前有关的方法与理论已经相当成熟，并且已成为求解大型稀疏线性方程组最受欢迎的一类方法.
数值代数实验报告
一、实验名称：用共轭梯度法解线性方程组。
二、实验目的：进一步熟悉理解掌握共轭梯度法解法思路，提高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+α0p0)≤ϕ(x0+αp0)，
即在这条直线上x1使ϕ(x)达到极小.然后从x1出发，再确定一个下山的方向p1，沿着直
线x=x1+αp1再跨出一步，即找到α1使得ϕ(x)在x2=x1+α1p1达到极小：
ϕ(x1+α1p1)≤ϕ(x1+αp1).
重复此步骤，得到一串
α0,α1,α2,
x=xk+αpk上确定步长αk使
和 p0,p1,p2,
，
称pk为搜索方向，αk为步长.一般情况下，先在xk点找下山方向pk，再在直线
ϕ(xk+αkpk)≤ϕ(xk+αpk),
最后求出xk+1=xk+αkpk.然而对不同的搜索方向和步长，得到各种不同的算法.
由此，先考虑如何确定αk.设从xk出发，已经选定下山方向pk.令 f(α)=ϕ(xk+αpk)
=(xk+αpk)A(xk+αpk)-2bT(xk+αpk)
T
=α2pkApk-2αrkTpk+ϕ(xk)，
T
其中rk=b-Apk.由一元函数极值存在的必要条件有
T
f'(α)=2αpkApk-2rkTpk=0
所确定的α即为所求步长αk，即 步长确定后，即可算出
此时，只要rkTpk≠0，就有
rkTpk
. αk=T
pkApk
xk+1=xk+αkpk.
ϕ(xk+1)-ϕ(xk)=ϕ(xk+αkpk)-ϕ(xk)
T
Apk-2αkrkTpk=- =αk2pk
(rkTpk)
2
即ϕ(xk+1)
T
pkApk
再考虑如何确定下山方向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+
η0
p ξ0k-1
η0
，则可ξ0
rkTApk-1
βk-1=-T.
pk-1Apk-1
T
注：这样确定的pk满足pkApk-1=0，即pk与pk-1是相互共轭的.
总结上面的讨论，可得如下的计算公式：
rkTpk
， xk+1=xk+αkpk， αk=T
pkApk
rk+1=b-Axk+1，
rkT+1Apk
， pk+1=rk+1+βkpk. βk=-T
pkApk
在实际计算中，常将上述公式进一步简化，从而得到一个形式上更为简单而且对称的计算公式.首先来简化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,
从而可导出
1
rkT+1=-rkT+1rk+1,， 1Tαk1TT
pkApk=pk(rk-rk+1)=pkrk
ααk1Tk1
=rk(rk+βk-1pk-1)=rkTrk.
.
由此可得
αkαk
rkTrkrkT+1rk+1αk=T,， βk=T..
pkApkrkrk
从而有求解对称正定方程组的共轭梯度法算法如下： x0=初值
r0=b-Ax0；k=0 while rk≠0
k=k+1
if k=1 p0=r0
else
βk-2=rkT-1rk-1rkT-2rk-2 pk-1=rk-1+βk-2pk-2 end
T
αk-1=rkT-1rk-1pk-1Apk-1
xk=xk-1+αk-1pk-1 rk=rk-1-αk-1Apk-1 end
x=xk
注：该算法每迭代一次仅需要使用系数矩阵A做一次矩阵向量积运算. 定理2 由共轭梯度法得到的向量组{ri}和{pi}具有如下基本性质： （1）piTrj=0， 0≤i
,rk}=span{p0,,pk}=κ(A,r0,k+1)，
κ(A,r0,k+1)=span{r0,Ar0,,Akr0}，
通常称之为Krylov子空间.下面给出共轭梯度法全局最优性定理： 定理3 用共轭梯度法计算得到的近似解xk满足
ϕ(xk)=min{ϕ(x):x∈x0+κ(A,r0,k)}
或
xk-x*
A
=min{x-x*
A
:x∈x0+κ(A,r0,k)}，
其中x
A
=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;ρ=rTr
while
>εb
2
)and(k
max
)
k=k+1
if k=1 p=r else
β=ρρ;p=r+βp end
ω=Ap;α=ρpTω;x=x+αp r=r-αω;ρ=ρ;ρ=rTr end
算法中，系数矩阵A的作用仅仅是用来由已知向量p产生向量ω=Ap，这不仅可以充分利用A的稀疏性，而且对某些提供矩阵A较为困难而由已知向量p产生向量
ω=Ap又十分方便的应用问题是十分有益的。
2．算例
⎡10 1 2 3 4⎤⎡x1⎤⎡ 12⎤⎢1 9 -1 2 -3⎥⎢x⎥⎢-27⎥运用共轭梯度法求解下述方程，并解释你所观察到的结果 ⎢⎥⎢2⎥⎢⎥⎢2 -1 7 3 -5⎥⎢x3⎥=⎢ 14⎥ ⎢⎥⎢⎥⎢⎥3 2 3 12 -1-17x解：已知共轭梯度法的⎢⎥⎢4⎥⎢MATLAB⎥程序代码如下所示：
⎥⎣ 12⎥⎢⎣4 -3 -5 -1 15⎥⎦⎢⎦⎣x5⎦⎢function[x,n]=conjgrad(A,b,x0) %采用共轭梯度法求线性方程组Ax=b的解 %线性方程组的系数矩阵：A %线性方程组中的常数向量：b %迭代初始向量：x0 %线性方程组的解：x
%求出所需精度的解实际的迭代步数：n r1=b-A*x0; p=r1; n=0;
for i=1:rank(A)
if(dot(p,A*p)
alpha=dot(r1,r1)/dot(p,A*p);
x=x0+alpha*p; r2=r1-alpha*A*p; if(r2
belta=dot(r2,r2)/dot(r1,r1); p=r2+belta*p; n=n+1; end
用共轭梯度法求解，在MATLAB命令窗口中输入求解程序：
>> A=[10 1 2 3 4 2 4 2 1 1;1 9 -1 2 -3 3 -3 -1 9 2;2 -1 7 3 -5 4 -5 7 -1 3;3 2 3 12 -1 5 -1 3 2 4;4 -3 -5 -1 15 1 3 2 3 2;2 3 4 5 1 10 4 5 3 4;4 -3 -5 -1 3 4 9 1 -3 2;2 -1 7 3 2 5 1 7 -1 1;1 9 -1 2 -3 3 -3 -1 12 10;1 2 3 4 2 4 2 1 10 15];
>> b=[12;-27;14;-17;12;12;-27;14;-17;12]; >> x0=[0;0;0;0;0;0;0;0;0;0]; >> [x,n]=conjgrad(A,b,x0) 输出的计算结果为： x = 8.5669 -7.1981 -7.3479 -13.9388 1.1303 11.5188 -26.8966 -2.2388 0.0829 13.2786
输出的迭代次数为 n = 10
共轭梯度法理论上只要迭代n步，就能得出方程组的解，但是由于存在计算误差，即分数向小数转化时存在舍入误差，很难保证在第n步时得到准确解。
3． 总结
本文首先给出最小二乘问题的定义，随后从盲人下山法开始讨论了共轭梯度法的具体推导过程及其相关性质与算法.继而重点给出正则化方法的实用共轭梯度算法并举例进行检验.最后，需要说明虽然正则化方法是求一般线性方程组Ax=b，
A∈Rm⨯n(m≥n)且A列满秩的最小二乘解的一种方法且简单易行，但是也有许多不
足之处，如m>n时一般无解；ATA形成时运算量大，A中某些信息会丢失；当A病
2
态时其收敛性速度由于κ2(ATA)=κ2(A)很大变得非常之慢等，故为了避免正则化方
法的缺点，还可运用残量极小化方法或残量正交化方法等更好的方法来解决此类问题.
在实际的科学与工程问题中，常常将问题归结为一个线性方程组的求解问题，而求解线性方程组的数值解法大体上可分为直接法和迭代法两大类.直接法是指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法.因此，直接法又称为精确法.迭代法则是采取逐次逼近的方法，亦即从一个初始向量出发，按照一定的计算格式，构造一个向量的无穷序列,其极限才是方程组的精确解，只经过有限次运算得不到精确解.当线性方程组的系数矩阵为对称正定矩阵时，我们常用共轭梯度法（或简称CG法）求解，目前有关的方法与理论已经相当成熟，并且已成为求解大型稀疏线性方程组最受欢迎的一类方法.