第三章 线性代数方程组的共轭梯度法

华北电力大学数理学院
School of mathematics & physics
§3 共轭梯度法
可以证明：CG产生的序列有 1、p(0) , p(1) , , p(k) , 是两两关于A共轭的； 2、(p(i) )T r(k) = 0, i = 0,1, , k −1; 3、r(0) , r(1) , , r(k) , 是两两正交的； 3、(r(k) , Ap(i) ) = 0, i ≠ k, k +1 4、r(n) = 0.
tk
=
(r (k ) )T r (k ) (r (k ) )T Ar (k )
x(k +1) = x(k ) + t k r (k )
华北电力大学数理学院
School of mathematics & physics
§2 最速下降法
function x=zuisu(a,b,x,e) n=length(b); r=b-a*x; q=r’*r; while q>e
q=0.7906. 最速下降法迭代29次，得到满足精度0.0005的解：
X=[2.9998 4.0003 -4.9999]
华北电力大学数理学院
School of mathematics & physics
§2 最速下降法
例2：用最速下降法解方程组
⎡ 4 3.5 0 ⎤⎡x1 ⎤ ⎡ 24 ⎤
⎢⎢3.5
2：在射线x = x(k ) + tp (k )上求f ( x)的极小点x(k +1) , 即 f ( x(k +1) ) = min f ( x(k ) + tp (k ) )
t>0
3： 判 断 ∇ f ( x (k +1) ) ≤ ε , 或 x(k +1 ) − x(k) ≤ ε ！
是 ！ x* ≈ x (k +1)， 停 止 ; 否 ！ k = k + 1, 转 1.
p(k+1) = r (k+1) + sk p(k)
华北电力大学数理学院
School of mathematics & physics
§3 共轭梯度法
实用共轭梯度算法：
实际计算中由于舍入误差的影r响(k) ， 之间的正交性很快就会消失，以至于有限步 内不能得到精确解，所以CG实际是一种迭代 算法。
∑ = (r(0) , p( j) ) − (r(0) − tk Ap(k ) , p( j) ) k =1
=0 ∴ r(n)与p(1)，p(2)， ，p(n)正交，必有r(n) = 0
华北电力大学数理学院
School of mathematics & physics
§3 共轭梯度法
定理1说明：取A共轭的方向的下降算法至多n步 得到n阶线性代数方程组的精确解。把这种方法 称之为共轭梯度法（CG）。
4
−
1.5⎥⎥
⎢⎢x
2
⎥ ⎥
=
⎢ ⎢
30
⎥ ⎥
⎢⎣ 0 −1.5 4 ⎥⎦⎢⎣x3 ⎥⎦ ⎢⎣− 24⎥⎦
A的特征值：0.1921 4.0000 7.8079 A的真解：x=[6 0 -6] q=0.9520
最速下降法迭代141次，得到满足精度0.0005的解： X=[5.9987 0.0015 -5.9994]
华北电力大学数理学院
School of mathematics & physics
第三章 线性代数方程组的共轭梯度法(CG)
问题：解Ax=b，其中，A为对称正定矩阵。
原理：Ax = b ⇔ min f (x)。A是对称正定矩阵。 x∈R n
其中，f(x) = 1 xT Ax − bT x。 2
注： ∵ ∇f ( x) = Ax − b, ∇ 2 f ( x) = A
华北电力大学数理学院
School of mathematics & physics
§2 最速下降法
例3：用最速下降法解方程组
⎡ 4 3.8 0 ⎤⎡x1 ⎤ ⎡ 24 ⎤
⎢⎢3.8
4
−
1.1⎥⎥
⎢ ⎢
x
2
⎥ ⎥
=
⎢ ⎢
30
⎥ ⎥
⎢⎣ 0 −1.1 4 ⎥⎦⎢⎣x3 ⎥⎦ ⎢⎣− 24⎥⎦
A的特征值：0.0440 4.0000 7.9560
如何寻找A共轭的方向呢？
华北电力大学数理学院
School of mathematics & physics
§3 共轭梯度法
开始： ∀x(0) ∈ R n , r(0) = b − Ax(0) , p(0) = r(0) ;
t0
=
(p(0) )T p(0) (p(k) )T A p(k)
, x(1)
华北电力大学数理学院
School of mathematics & physics
定理1：设p 1 , p 2 , p n 是关于n阶对称正定 矩阵A共轭的向量组，则以 p(k)为下降方向的算
法：
x(0) ∈ Rn ; r(0) = b − Ax(0); k = 1, 2,
tk
=
(r (k −1) )T p(k ) ； ( p(k) )T A p(k)
∴(r(n) , p( j) ) = (r(0) , p( j) ) − (t j Ap( j) , p( j) ) = (r (0) , p( j) ) − (r ( j-1) − r ( j) , p( j) )
= (r (0) , p( j) ) − (r ( j-1) , p( j) )
j -1
A的真解：x=[-0.51428571428573 6.85714285714288 -4.11428571428571]
q=0.9890
最速下降法迭代568次，得到满足精度0.0005的解： X=[ -0.5088 6.8513 -4.1159]
华北电力大学数理学院
School of mathematics & physics
x(k ) = x(k −1) + t k p(k )； r (k ) = b − Ax(k ) = r (k −1) − t k Ap(k )
有Ax(n) = b,即r(n) = 0.
华北电力大学数理学院
School of mathematics & physics
§3 共轭梯度法
n
∑ 证明：∵ Ax(n) = Ax(n−1) + tn Ap(n) = Ax(0) + tk Ap(k ) k =1 n ∑ r(n) = b − Ax(n) = r(0) − tk Ap(k ) k =1
华北电力大学数理学院
School of mathematics & physics
§2 最速下降法
例1：用最速下降法解方程组
⎡4 3 0 ⎤⎡x1 ⎤ ⎡ 24 ⎤
⎢⎢3
4
−
1⎥⎥
⎢⎢x
2
⎥ ⎥
Fra Baidu bibliotek
=
⎢ ⎢
30
⎥ ⎥
⎢⎣0 −1 4 ⎥⎦⎢⎣x3 ⎥⎦ ⎢⎣− 24⎥⎦
A的特征值：0.8377,4.0000,7.1623 A的真解：x=[3 4 -5]
华北电力大学数理学院
School of mathematics & physics
§2 最速下降法
下降方向 p(k) 的确定！
函数的负梯度方向是函数值下降最快的方向，因
此，取 p(k) = − f ′(x(k) ) = b − Ax(k) = r(k)为下降方向。
最速下降法： r (k ) = b − Ax(k )
∴ 当A对称正定时， ∇f ( x) = Ax − b = 0 ⇔ min f ( x)。 x∈R n
如何寻找f(x)的极小点呢？
华北电力大学数理学院
School of mathematics & physics
§1 下降法的一般理论
基本思想就是下降法:
从某一点x(0)出发, 逐步产生一串点： x(0) , x(1) , , x(k ) ,
导出 x(k +1) ：
华北电力大学数理学院
School of mathematics & physics
§1 下降法的一般理论
r (k ) = b − Ax(k )
tk
=
(r (k) )T p(k) (p(k) )T Ap(k)
x(k +1) = x(k ) +t k p (k )
t k称之为x(k)到x(k+1)的步长。
使 f (x(0) > f (x(1) ) > f (x(k) ) >
并以“最快的速度”下降到f (x)的极小值。
完成任务的关键就是确定每步的下降方向！
华北电力大学数理学院
School of mathematics & physics
§1 下降法的一般理论
具体过程： 已求得x(k)
1：在x(k )处确定一个使f ( x)下降的方向p(k )；
w=a*r; t=q/(r'*w); x=x+t*r; r=r-t*w; q=r’*r; end
华北电力大学数理学院
School of mathematics & physics
§2 最速下降法
近似解的误差估计
x(k) − x*
A
=
⎜⎜⎝⎛
λn λn
− +
λ1 λ1
⎟⎟⎠⎞k
x(0) − x*
A
注：λ1、λ n分别是A的最小和最大特征值, 收敛速度 由q = λn − λ1 决定。q越小收敛越快。 λn + λ1
华北电力大学数理学院
School of mathematics & physics
§3 共轭梯度法
由此，可使CG进一步简化为：
∀x(0) ∈ R n , r(0) = b − Ax(0) , p(0) = r(0) ;
对于 k = 0 ,1, 2 ,
tk
=
(r(k) )T p(k) (p(k) )T Ap(k)
sk
=
−
(p(k−1) )T Ar(k) (p(k−1) )T Ap(k−1)
华北电力大学数理学院
School of mathematics & physics
§3 共轭梯度法
对于 k = 1, 2 ,
现有x(k)及共轭方向p(k) ,则
tk
=
(r(k) )T p(k) (p(k) )T Ap(k)
=
(r(k) )T r(k) (p(k) )T Ap(k)
x (k+1) = x (k) + t k p(k)
r(k+1) = b − Ax(k+1) = r(k) − t k Ap(k)
sk
= − (p(k) )T Ar(k+1) (p(k) )T Ap(k)
=
(r (k+1) )T r (k+1) (r(k) )T r(k)
§2 最速下降法
小结
最速下降法有简单易行，保稀疏性等特点，但 当A的最大特征值远远大于最小特征值时收敛速度变 得非常慢。最速下降法并非最速！
共轭梯度法可使这一问题得到一定改善！
华北电力大学数理学院
School of mathematics & physics
§3 共轭梯度法
定义：向量p、q称之为关于对称正定矩阵A共轭，如 果满足 pΤAq = 0 。
=
x(0)
+t0p(0).
华北电力大学数理学院
School of mathematics & physics
§3 共轭梯度法
对于 k = 1,2,
现有x ( k )及共轭方向p ( k −1)，则
r(k) = b − Ax(k)
在r ( k )和p ( k −1) 确定的超平面上 找共轭方向
其中
p(k) = r(k) + sk p(k−1)
x (k+1) = x (k) + t k p(k) r(k+1) = b − Ax(k+1)
华北电力大学数理学院
School of mathematics & physics
§3 共轭梯度法
综合上述，有下面的共轭梯度法：
∀x(0) ∈ R n , r(0) = b − Ax(0) , p(0) = r(0) ;
对于 k = 0 ,1, 2 ,
tk
=
(r(k) )T p(k) (p(k) )T Ap(k)
x (k+1) = x (k) + t k p(k)
r(k+1) = b − Ax(k+1) = r(k) − t k Ap(k)
sk
=
− (p(k) )T Ar(k+1) (p(k) )T Ap(k)
p(k+1) = r (k+1) + sk p(k)
= 1 ( p(k ) )T Ap(k ) ⋅ t 2 + ( p(k ) )T ( Ax(k ) − b) ⋅ t 2 + 1 (x(k ) )T Ax(k ) − bT x(k ) 2
df = ( p(k ) )T Ap(k) ⋅ t + ( p(k ) )T ( Ax(k ) − b) = 0 dt
华北电力大学数理学院
School of mathematics & physics
§1 下降法的一般理论
问题 f (x(k +1) ) = min f (x(k) + tp(k) ) 的解决！ t>0
f (x(k ) + tp(k ) ) = 1 (x(k ) + tp(k ) )T A(x(k ) + tp(k ) ) − bT (x(k ) + tp(k ) ) 2