第三章 线性代数方程组的共轭梯度法
共轭梯度法公式
共轭梯度法公式
共轭梯度法是一种用于求解线性方程组的迭代算法。
其主要思想是通过利用前一次迭代的信息来加速当前迭代的速度，从而减少迭代次数和计算量。
共轭梯度法公式包括以下几个步骤：
1. 初始化：设初始解为x0，残量b0为Ax0-b，共轭方向d0=b0。
2. 迭代求解：对于第k次迭代，计算步长αk，使得xk+1=xk+αkd，其中d是共轭方向，满足dTkAd=0，即d是A的共轭向量。
3. 更新残量：计算新的残量bk+1=Axk+1-b，如果bk+1小于预设精度，则停止迭代。
4. 更新共轭方向：计算新的共轭方向dk+1=bk+1+βkdk，其中βk=(bk+1)Tbk+1/(bk)Tbk，保证dk+1与之前的共轭方向都是A的共轭向量。
5. 重复迭代，直到满足收敛条件，返回最终解xk+1。
共轭梯度法是一种高效的求解大型线性方程组的方法，尤其适用于稀疏矩阵和对称正定矩阵。
公式简单易懂，容易实现，且具有较快的收敛速度。
- 1 -。
共轭梯度法步骤
共轭梯度法步骤共轭梯度法是一种求解线性方程组的迭代算法，它以高效稳定的特点而广受欢迎。
以下是共轭梯度法的步骤：步骤1：初始化首先，我们需要有一个初始向量x0和一个初始残量r0=b-Ax0。
其中，A为系数矩阵，b为常数向量。
步骤2：计算方向向量令d0=r0，表示第一次迭代的方向向量。
步骤3：计算步进长度令α0=(r0·r0)/(d0·Ad0)，其中·表示向量的点积。
α0表示迭代过程中每个方向向量的步进长度。
步骤4：更新解向量令x1=x0+α0d0，表示迭代后的解向量。
步骤5：计算新残量令r1=r0-α0Ad0。
步骤6：判断终止条件如果r1的范数小于预设阈值，或者迭代次数达到预设次数，终止迭代。
否则，进入下一次迭代。
步骤7：更新方向向量令β1=(r1·r1)/(r0·r0)，表示更新方向向量的轴线。
步骤8：计算新方向向量令d1=r1+β1d0，表示新的迭代方向向量。
步骤9：计算新的步进长度令α1=(r1·r1)/(d1·Ad1)。
步骤10：更新解向量令x2=x1+α1d1。
步骤11：更新残量令r2=r1-α1Ad1。
步骤12：重复步骤6至11，直至满足终止条件。
总结起来，共轭梯度法的步骤主要包括初始化、计算方向向量、计算步进长度、更新解向量、计算新残量、判断终止条件、更新方向向量、计算新的步进长度、更新解向量和更新残量等。
该算法迭代次数较少，收敛速度快，适用于大规模线性方程组的求解。
共轭梯度法例题详解
共轭梯度法例题详解共轭梯度法是一种用于求解线性方程组的迭代方法。
它的基本思想是利用一组共轭的搜索方向来逐步逼近线性方程组的解。
下面我将从多个角度详细解释共轭梯度法的例题。
首先，让我们考虑一个简单的例子。
假设我们要求解一个二维线性方程组 Ax = b，其中 A 是一个对称正定矩阵，x 和 b 是列向量。
共轭梯度法的步骤如下：1. 初始化 x0 和 r0，其中 x0 是初始解向量，r0 是初始残差向量，r0 = b Ax0。
2. 初始化搜索方向 p0 = r0。
3. 迭代计算：a. 计算步长 alpha = (r_k^T r_k) / (p_k^T A p_k)，其中 k 表示第 k 次迭代。
b. 更新解向量，x_(k+1) = x_k + alpha p_k。
c. 更新残差向量，r_(k+1) = r_k alpha A p_k。
d. 计算 beta = (r_(k+1)^T r_(k+1)) / (r_k^T r_k)。
e. 更新搜索方向，p_(k+1) = r_(k+1) + beta p_k。
4. 重复步骤 3 直到满足收敛条件。
这是共轭梯度法的基本算法。
下面我们通过一个具体的例子来说明。
假设我们要求解以下线性方程组：3x + 2y = 9。
2x + 5y = -4。
将其转化为矩阵形式，Ax = b，其中。
A = [[3, 2],。
[2, 5]],。
x = [[x],。
[y]],。
b = [[9],。
[-4]].首先，我们需要判断矩阵 A 是否是对称正定矩阵。
对于本例中的 A，它是对称矩阵且特征值为正，因此满足条件。
接下来，我们进行共轭梯度法的迭代计算。
假设初始解向量 x0 = [[0], [0]]，初始残差向量 r0 = b Ax0 = [[9], [-4]]。
初始化搜索方向 p0 = r0。
第一次迭代：计算步长 alpha = (r0^T r0) / (p0^T A p0) = (81 + 16) / (9 + 20) = 97 / 29 ≈ 3.34。
共轭梯度法(讲稿)3.
• 一、共轭梯度法的适用范围 • 二、等价极小值问题 • 三、极小化迭代法基本步骤 • 四、共轭梯度法
一、共轭梯度法的适用范围
• 1、CG法适用于求解大散射体的问题也可以解谐振问题 • 2、与SIT法比较，都可以避免矩阵求逆，但SIT法收敛较慢，有时不 一定收敛，而CG法则能保证收敛，误差小，贮存量较SIT大一些，且 其初始值可任意选定。 • 3、最速下降法反映的目标函数的一种局部性质,从局部看, 最速下降 方向是目标函数值下降最快的方向,选择这样的方向进行搜索是有利 的. • 4、但从全局来看,由于锯齿现象的影响, 即使向着极小点移近不太大 的距离,也要经历不小的”弯路”,因此收敛速度大为减慢.
解 设初始点为U ( 0) (1,1)T ，U (u1 , u 2 , u3 ...un )T 2u1 F(u 1 , u 2 ) 8u 2 (1,1)T 得， 2 F(U ( 0 ) ) , F(U ( 0 ) ) 8.24621 8 p ( 0 ) F (U ( 0 ) ) (2,8)T U(1) U ( 0 ) t0 p ( 0 ) , 其中t0由 min F (U ( 0 ) tp ( 0 ) ) min[( 1 2t ) 2 4(1 8t ) 2 ] dF (U ( 0) tp ( 0 ) ) 利用必要条件 4(1 2t ) 64(1 8t ) 520t 68 0 得t 0 0.13077 dt 1 2 0.73846 U (1) 0 . 13077 1 8 0.04616 F (U (1) ) (1.47692 ,0.36923 )T , p (1) F (U (1) ), F(U (1) ) 1.52237 U( 2 ) U (1) t1 p (1)
共轭梯度法
10
11
12
13
14
由定理5.18，CG法最多迭代n次得到精确解
一方面，因舍入误差，n次迭代不一定得到精确解 另一方面，对大型方程组，n次迭代的工作量仍太大 实际上，通过比较当前残量rk的范数和初始残量r0的 范数，来判断是否终止迭代
15
16
int main(){ double A[2][2]={{6,3}, {3,2}}; double b[2]={0,-1}, x[2]={0,0}, r[2], r_p[2], p[2], Ap[2]; double alpha, beta; int i,j,k; for(i=0; i<2; i++) p[i] = r[i] = b[i] - (A[i][0]*x[0] + A[i][1]*x[1]); for(k=1; k<=2; k++){ for(i=0; i<2; i++) Ap[i] = A[i][0]*p[0] + A[i][1]*p[1]; alpha = (r[0]*r[0] + r[1]*r[1]) / (Ap[0]*p[0] + Ap[1]*p[1]); for(i=0; i<2; i++) x[i] = x[i] + alpha*p[i]; for(i=0; i<2; i++){ r_p[i] = r[i]; r[i] = b[i] - (A[i][0]*x[0] + A[i][1]*x[1]); } beta = (r[0]*r[0] + r[1]*r[1]) / (r_p[0]*r_p[0] + r_p[1]*r_p[1]); for(i=0; i<2; i++) p[i] = r[i] + beta*p[i]; printf("k=%2d, ", k); for(i=0; i<2; i++) printf("x%d=%6.2f, ", i, x[i]); printf("\n"); } }
计算方法——共轭梯度法求解线性方程组
(2)
(k)
共轭梯度法中关键的两点是迭代格式(2)中最佳步长k 和搜索方向 d
(k) (k)
的确定。其
中k 可以通过一元函数f(x +d )的极小化来求得，其表达式为公式(3)；取 d (0) = r(0) = b-Ax(0)，则 d(k+1) = r(k+1) +kd(k)，要求 d(k+1)满足 (d(k+1) , Ad(k)) = 0，可得k 的表达 式(4)。
计算方法上机报告
计算方法上机报告
1 共轭梯度法求解线性方程组
1.1 算法原理及程序框图 当线性方程组 Ax = b 的系数矩阵 A 是对称正定矩阵是，可以采用共轭梯度法对该 方程组进行求解，可以证明，式(1)所示的 n 元二次函数 1 f ( x ) x T Ax bT x (1) 2 取得极小值点 x*是方程 Ax = b 的解。共轭梯度法是把求解线性方程组的问题转化为求 解一个与之等价的二次函数极小化的问题。从任意给定的初始点出发，沿一组关于矩 阵 A 的共轭方向进行线性搜索，在无舍入误差的假定下，最多迭代 n 次（其中 n 为矩 阵 A 的阶数） ，就可求得二次函数的极小点，也就求得线性方程组 Ax = b 的解。其迭 代格式为公式(2)。
4
1 1
1 0 ，b 2 1 0 1 1 2
计算方法上机报告
矩阵 A 的阶数取 200 进行求解。 由于该线性方程组的系数矩阵阶数比较大，而且具有一定的规律，因此首先用 matlab 编程将系数矩阵、右端项以及阶数保存在 D 盘根目录的三个文件中（生成系数 矩阵， 右端项以及阶数的程序见附录 2） ， 然后运行共轭梯度法程序进行方程组的求解。 最终的运行结果为图 4 和图 5。程序运行之后 D 盘根目录下生成的文件如图 6 所示。
共轭梯度法课件
4.3共轭梯度法4.3.1共轭方向法定义4.3.1设A 是n ×n 对称正定矩阵，d 1，d 2，是n 维非零矢量，如果d 1T Ad 2=0则称d 1和d 2是A-共轭的，简称共轭的设d 1，d 2，...,d m 是R n 中一组非零向量，如果d i T Ad j =0,i ≠j ，j,i=1,2,...,k则d 1,d 2,...,d m 是A-共轭的，简称共轭的，也称它们是一组A 共个方向定理4.3.3设x 0∈Rn 是任意初始点，对于极小化二次函数min f(x)=1/2 x T Ax-b T x 共轭方向法至多经n 步精确线性搜索终止；且每一x i+1都是f(x)在x 0和方向d 1，d 2,....,di, 所张成的线性流形｛|x x=x 0+,0j i j j da ∑=j a ∀｝中的极小点。
4.3.4共轭梯度法共轭梯度法是一个典型的共轭方向法，他的每一个搜索方向是相互共轭的，而这些搜索方向d k 仅仅是负梯度方向-g k 与上一次迭代的搜索方向d k-1组合。
因此，存储量小，计算方便。
定理4.3.6对于正定二次函数，采用精确线性搜索的共轭梯度法在m ≦n 步后终止，且对1≦i≦n成立下列关系式：d i T Ad j=0,j=0,1,...,i-1,g i T Ag j=0,j=0,1-1,d i T Ag i= - g i T g I[g0,g1,...,g i]=[g0,Ag0,,...,A i g0][d0,d1,...,d i]=[g0,Ag0,,...,A i g0]其中[g0,g1,...,g i]和[d0,d1,...,d i]分别表示g0,g1,...,g i及d0,d1,...,d i张成的子空间，[g0,Ag0,,...,A i g0]表示g0的i阶Krylov子空间。
定理4.3.9（FR共轭梯度法的总体收敛性定理）假定f R n R在有界水平集L={x R n|f（x）≦f（x0）｝上连续可微，且有下界，那么采用精确线性搜索的F-R共轭梯度法产生的序列｛x k｝至少有一个聚点是驻点，即1当｛x k｝是有穷数列时，其最后一个点是f（x）的驻点；2当｛x k｝是无穷数列时，它必有聚点，且任一聚点都是f（x）的驻点。
利用共轭梯度法求解线性方程组
利用共轭梯度法求解线性方程组翟莹1012205052在自然科学和工程技术中很多问题的解决常常归结为解线性方程组，而这些方程组的系数矩阵大致可分为两种:低阶稠密矩阵和大型稀疏矩阵。
而求解方程组的方法通常为直接法和迭代法。
直接法用于较低阶方程组的求解，效率较高；迭代法更适用于高阶方程组的求解，常用的经典迭代法有高斯-赛德尔迭代法和雅各比迭代法，但收敛效率较低；共轭梯度法（CG）以较高的收敛速度而经常被采用。
从理论上讲，一个n阶方程组最多迭代n 步就可求出精确解。
1 直接法直接法就是经过有限步算术运算，无需迭代可直接求得方程组精确解的方法。
但实际计算中由于舍入误差的存在和影响，这种方法也只能得到线性方程组的近似解，该方法是求解低阶稠密矩阵方程组的有效方法。
如Cramer法则，Gauss消元法及其变形（LU分解法、Cholesky 分解法、QR分解法）等。
Matlab中，用矩阵除法“/”或“\”直接求解线性方程组（见附录一），它是一个内部包含着许许多多的自适应算法，对超定方程用最小二乘法求解；对欠定方程因为它的解不唯一，Matlab给出所有解中范数最小的一个特解；对于三对角阵方程组，采用追赶法求解。
2 迭代法迭代法就是用某种极限过程去逐步逼近线性方程组精确解的方法。
该方法具有对计算机的存贮单元需求少，程序设计简单、原始系数矩阵在计算过程中不变等优点，是求解大型稀疏矩阵方程组的重要方法。
迭代法不是用有限步运算求精确解，而是通过迭代产生近似解逼近精确解。
如Jacobi迭代、Gauss-Seidel迭代、SOR迭代法等。
在科学研究和大型工程设计中提出的求解问题，经离散后，常常归结为求解形如Ax = b 的大型线性方程组，此时系数矩阵A和b是通过测量或其它方法得到的，但是在多数情况下方程可能是病态的，即A的条件数非常大。
此时，我们仍然用Matlab中的常规算法，得到的解则可能不是真解。
通常情况下，对系数矩阵A和方程的右端b作微小的扰动，再用上述方法求解，扰动后方程组的解与原方程组的解相差甚远。
共轭梯度法
3. 算法的 MATLAB 实现
在 MATLAB 中编程实现的共轭梯度法函数为： min GETD 功能：用共轭梯度法求解多维函数的极值。
调用格式： [ x, min f ] min GETD( f , x0, var, eps) 其中， f ：目标函数；
x0 ：初始点； ：自变量向量； var x ：目标函数取最小值时的自变量值；
（6） 若 k 1 n ，令 x(0) x( n ) ，转步骤（3） ，否则转步骤（7） ；
（7） 令 p
( k 1)
f ( x
( k 1)
) k p
(k )
， k
f ( x ( k 1) ) f ( x )
(k )
2 2
， 置 k k 1， 转步骤 （4） 。
t 0
(k ) ( k ) （4） 用 一 维 搜 索 法 求 t k ， 使 得 f ( x kt p ) m i n f
( ) 5； x( k 1 ) x k( ) kt p k，转步骤
k( ) x
) 令， tk p( ，
（5） 若 f ( x
( k 1)
) ，停止，极小值点为 x ( k 1) ，否则转步骤（6） ；
1 f ( X ) d1 + x1 -x1* 2 1 =d1 + x1 -x1* 2
2
*
2
1 d1 a x1 -x1* 2
2b x1 -x1*
x -x c x -x
2
2b x1 -x1*
共轭梯度方法
共轭梯度方法（Conjugate Gradient Method）是求解线性方程组的一种迭代算法。
该方法适用于求解大型稀疏的对称正定线性方程组，可以显著减少计算量和存储空间。
该方法的主要思想是利用共轭方向（Conjugate Directions）的性质，在有限次迭代中求解方程组的解。
共轭梯度方法的基本步骤如下：
选取一个初值$x_0$，并令$r_0=b-Ax_0$，其中$b$ 为方程组的右端向量，$A$ 为系数矩阵。
计算一个共轭方向$p_0=r_0$，即$p_0$ 与$r_0$ 正交，并满足$Ap_0 \neq 0$。
对于$k=0,1,2,\ldots$，执行以下操作：
a. 计算$\alpha_k=\frac{r_k^Tr_k}{p_k^TAp_k}$。
b. 更新解向量$x_{k+1}=x_k+\alpha_kp_k$。
c. 计算残差向量$r_{k+1}=r_k-\alpha_kAp_k$。
d. 计算$\beta_k=\frac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k}$。
e. 更新共轭方向$p_{k+1}=r_{k+1}+\beta_kp_k$，即$p_{k+1}$ 与$p_k$ 具有共轭性。
如果残差向量$r_k$ 较小，则停止迭代，输出解向量$x_k$。
共轭梯度方法具有收敛速度快、存储空间小等优点，但对于非对称和非正定的线性方程组，该方法可能不收敛。
同时，该方法也有一些变体，如预处理共轭梯度法、共轭残差法等，可以更好地解决不同类型的线性方程组求解问题。
线性方程组的共轭梯度法30页PPT文档
W(1)
W(2)
沿两个相互正交的方向，进行精确一维搜索， 即可得到最优解（二维情形）
罗林开
模式识别与智能系统研究所
2020/6/27
6
n维情形: 沿相互正交的n个方向,进行精确 一维搜索,至多迭代n次,即可得到正定二次 函数
f (w) 1 wT w rT w 2
的最优解w* r.
f (x) f (x*)f (x*)T(xx*)1(xx*)T A(xx*) 2
f (x*)1(xx*)T A(xx*) 2
因此f (x)的等高面是一簇超椭球面.
罗林开
模式识别与智能系统研究所
2020/6/27
9
x(3) =x* d(2)
x(1) d(1) x(2)
从x(1)出发,沿d(1)作一维搜索,与等高面相切于x(2),记d(2) x* x(2), 再沿d(2)作一维搜索,即可得到最优解x*.由于(d(1))Tf (x(2)) 0, 因此(d(1))T Ad(2) (d(1))T A(x* x(2)) (d(1))T(Ax* Ax(2))
f (w ) f (w *) f (w *)(w w*)
1 (w w *)T I (w w *) 2
f ( w * ) 1 || w w * ||2 . 2
因 此 f ( w )的 等 高模式识别与智能系统研究所
2020/6/27
5
正交方向及其性质
( d ) (1) T A d ( 2 ) 0 则 称 这 两 个 方 向 关 于 A 共 轭 ,或 称 它 们 关 于 A 正 交 . 若 R n中 的 k 个 方 向 d (1) , d ( 2 ) ,L , d ( k ) , 它 们 两 两 关 于 A 共 轭 ,即
共轭梯度法求解方程组
共轭梯度法是一种常用的迭代方法，用于求解线性方程组Ax = b。
它适用于对称正定矩阵的情况，可以高效地求解大规模的线性方程组。
下面是使用共轭梯度法求解方程组的一般步骤：1. 初始化：选择一个初始解x0 和初始残差r0 = b - Ax0，设置初始搜索方向d0 = r0。
2. 迭代计算：进行迭代计算，直到满足停止准则（如残差的大小或迭代次数达到一定阈值）为止。
a. 计算步长αk = (rk^T rk) / (dk^T A dk)，其中rk = b - A xk 是当前的残差。
b. 更新解xk+1 = xk + αk dk。
c. 计算新的残差rk+1 = rk - αk A dk。
d. 计算新的搜索方向dk+1 = rk+1 + (rk+1^T rk+1) / (rk^T rk) dk。
e. 更新迭代次数k = k + 1。
3. 输出解：当满足停止准则时，输出最终的解x。
需要注意的是，共轭梯度法的效率和收敛速度与矩阵的条件数有关。
对于病态矩阵或条件数较大的情况，可能需要进行预处理或使用其他更适合的求解方法。
此外，共轭梯度法还可以应用于非线性方程组的求解，采用牛顿法等方法来迭代求解。
在实际应用中，可以使用现有的数值计算库或软件来实现共轭梯度法，以提高计算的效率和精度。
线性方程组的共轭梯度法
迭代过程
计算方程组的雅可比矩阵A和右端项b，得到线性方程组Ax=b。 计算初始残差r0=b-Ax0。 进行迭代，对于k=0,1,2,...,max_iter，执行以下步骤
迭代过程
01
1. 计算搜索方向pk=-Ak^T。
02
2. 在搜索方向pk上进行线搜索，找到步长λk，使得 Axk+1=b-λk*r^k最小化。
感谢观看
THANKS
定义
线性方程组是由一组线性方程组成的 数学模型，其中包含未知数和已知数。
分类
根据方程的系数矩阵和常数项矩阵， 线性方程组可以分为多种类型，如超 定方程组、欠定方程组和恰定方程组。
线性方程组的求解方法
直接法
通过消元或迭代等方法将方程组化为标准形式，然后 求解。
迭代法
通过不断迭代更新解的近似值，逐步逼近方程的解。
在金融工程中的应用
投资组合优化
共轭梯度法可以用于求解投资组合优化问题 ，以最大化投资收益或最小化风险。
期权定价
在期权定价模型中，共轭梯度法可以用于求解 Black-Scholes方程，以得到期权的合理价格。
风险管理
在风险管理方面，共轭梯度法可以用于求解 风险评估模型中的最优化问题，以评估和管 理金融风险。
解效率。
02
常用的预处理方法包括对角占优预处理、不完全LU
分解预处理等。
03
预处理技术可以消除原始方程组中的病态条件，降低
数值误差的放大效应。
自适应步长调整策略
自适应步长调整策略可以根据上 一步的搜索结果动态调整步长， 提高算法的稳定性和收敛速度。
常见的自适应步长调整策略包括 Armijo线搜索、Goldstein线搜
科学计算
第三章 线性代数方程组的共轭梯度法
对于 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)
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.
=
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)
华北电力大学数理学院
School of mathematics & physics
定理1：设p 1 , p 2 , p n 是关于n阶对称正定 矩阵A共轭的向量组，则以 p(k)为下降方向的算
线性代数方程组中的预处理共轭梯度法ppt课件
4 预处理方法
• 预处理方法
•
•
•
•
•
取预优矩阵（预处理矩阵）为A的一个小带宽部分（如三对角或对角线
部分）
矩阵分裂，尤其是线性稳定迭代中的矩阵A的分裂构造预处理矩阵
通过A的各种近似分解得到预处理矩阵（如不完全分解）
通过矩阵A的多项式构造预处理矩阵
子结构，区域分裂，EBE预处理途径等等
4 预处理方法
R中的非零元素远远少于A。R需要显式存储
•
•
•
•
•
第二种高效实现分裂
A是对称的
= 0 − − E ，-E严格下三角，D0是对角线矩阵
M = ( − ) −1 ( − E )，D不完全等于D0
SSOR-PCG
谢 谢！
4 预处理方法
• 预处理方法
•
•
•
不完全分解预处理方法：
主要有不完全LU分解、不完全Cholesky分解以及相应的块不完全分解等
设 是A的一个近似分解
= + ，L是下三角矩阵，R称为剩余矩阵，M= 为预处理矩阵。
此方法与CG迭代方法相结合，就形成了不完全Cholesky分解预处理共轭
1 引入
• 线性代数方程组的解法
• 直接法：高斯消去法，分解法
• 迭代法：
•
•
古典迭代法：Jacobi，Gauss-Seidel，SOR，SSOR
现代迭代法：（投影方法，子空间法）
•
•
•
•
正交化的误差投影型Krylov：FOM，IOM，DIOM
对称情形误差投影型Krylov：Lanczos，CG，PCG
• 预处理方法
•
对角预处理法：
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
定理1：设p 1 , p 2 ,
pn
有Ax
(n)
= b,即r
(n)
= 0.
华北电力大学数理学院
School of mathematics & physics
§3 共轭梯度法
证明： Ax ∵ r
(n)
= Ax
( n −1)
+ tn Ap =r
(0)
(n)
= Ax
n
(0)
+ ∑ tk Ap ( k )
华北电力大学数理学院
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 ) + s k p ( k −1) 其中 (p ( k −1) ) T Ar ( k ) s k = − ( k −1) T ( k −1) (p ) Ap
华北电力大学数理学院
School of mathematics & physics
§2 最速下降法
例1：用最速下降法解方程组
⎡4 3 0 ⎤ ⎡ x1 ⎤ ⎡ 24 ⎤ ⎢3 4 − 1⎥ ⎢ x ⎥ = ⎢ 30 ⎥ ⎢ ⎥⎢ 2 ⎥ ⎢ ⎥ ⎢0 − 1 4 ⎥ ⎢ x 3 ⎥ ⎢− 24⎥ ⎣ ⎦⎣ ⎦ ⎣ ⎦
(p ( k ) ) T Ar ( k +1) sk = − (k) T (k) (p ) Ap p ( k +1) = r ( k +1) + s k p ( k )
华北电力大学数理学院
School of mathematics & physics
§3 共轭梯度法
可以证明：CG产生的序列有 1、 p , p ,
华北电力大学数理学院
School of mathematics & physics
§1 下降法的一般理论 具体过程： 已求得x
(k )
1：在 x ( k )处确定一个使 f ( x )下降的方向 p ( k )；
2：在射线 x = x ( k ) + tp ( k ) 上求 f ( x )的极小点 x ( k +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 = 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
(k ) (k )
(r ) p t k = (k) T (k) (p ) Ap x
( k +1)
=x
(k )
+t k p
(k )
t k 称之为x (k)到x (k +1)的步长。
华北电力大学数理学院
School of mathematics & physics
§2 最速下降法
p ( k ) 的确定！ 下降方向
A的特征值：.8377,4.0000,7.1623 0
A的真解：x=[3 q=0.7906.
4
-5]
最速下降法迭代29次，得到满足精度0.0005的解： X=[2.9998 4.0003 -4.9999]
华北电力大学数理学院
School of mathematics & physics
§2 最速下降法
( k +1)
=x
(k )
+tkr
(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 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 − λ1 ⎞ ( 0 ) ⎟ x − x* =⎜ ⎜λ +λ ⎟ ⎝ n 1⎠
k
A
注：λ1、λ n分别是A的最小和最大特征值, 收敛速度 λ n − λ1 由q = 决定。q越小收敛越快。 λ n + λ1
§3 共轭梯度法
定义：向量p、q称之为关于对称正定矩阵A共轭，如 p Τ Aq = 0 。 果满足
华北电力大学数理学院
School of mathematics & physics
是关于n阶对称正定 p ( k )为下降方向的算 矩阵A共轭的向量组，则以 法：
x (0) ∈ R n ; r (0) = b − Ax (0) ; k = 1, 2, (r ( k −1) )T p ( k ) ； t k = (k) T (k) (p ) A p x ( k ) = x ( k −1) + t k p ( k )； r ( k ) = b − Ax ( k ) = r ( k −1) − t k Ap ( k )
华北电力大学数理学院
School of mathematics & physics
§2 最速下降法 小 结
最速下降法有简单易行，保稀疏性等特点，但 当A的最大特征值远远大于最小特征值时收敛速度变 得非常慢。最速下降法并非最速！ 共轭梯度法可使这一问题得到一定改善！
华北电力大学数理学院
School of mathematics & physics
函数的负梯度方向是函数值下降最快的方向，因 此，取 p ( k ) = − f ′( x ( k ) ) = b − Ax ( k ) = r ( k )为下降方向。 最速下降法：
r
(k )
= b − Ax
(k ) T
(k ) (k )
(r ) r tk = ( k ) T ( k ) (r ) Ar x
华北电力大学数理学院
School of mathematics & physics
§3 共轭梯度法
对于 k = 1 , 2 ,
现有x 及共轭方向p , 则
(k) (k )
(r ( k ) ) T p ( k ) t k = (k) T (k) (p ) Ap x ( k +1) = x ( k ) + t k p ( k ) r ( k +1) = b − Ax ( k +1)
∵ ∇ f ( x ) = Ax − b , ∇ 2 f ( x ) = A 注： ∴ 当 A对称正定时， ∇ f ( x ) = Ax − b = 0 ⇔ min f ( x )。 n
x∈ R
如何寻找f(x)的极小点呢？
华北电力大学数理学院
School of mathematics & physics
k =1
n
(n)
= b − Ax
(n)
− ∑ tk Ap ( k )
k =1
∴ (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 ) ) = (r (0) , p ( j ) ) − (r (0) − ∑ tk Ap ( k ) , p ( j ) )
§1 下降法的一般理论 基本思想就是下降法:
从某一点x (0)出发, 逐步产生一串点： x (0) , x (1) ,
使 f ( x (0) > f ( x (1) ) > f ( x(k ) ) > 并以“最快的速度”下降到f ( x)的极小值。
, x(k ) ,
完成任务的关键就是确定每步的下降方向！
如何寻找A共轭的方向呢？
华北电力大学数理学院
School of mathematics & physics
§3 共轭梯度法
∀x ( 0 ) ∈ R n , r ( 0) = b − Ax ( 0 ) , p ( 0) = r ( 0) ; 开始：
(p ( 0) ) T p ( 0) t 0 = (k) T , x (1) = x ( 0 ) + t 0 p ( 0 ) . (p ) A p (k)
A的特征值：0.0440
4.0000
7.9560
A的真解：x=[-0.51428571428573 6.85714285714288 -4.11428571428571] q=0.9890 最速下降法迭代568次，得到满足精度0.0005的解： X=[ -0.5088 6.8513 -4.1159]
例2：用最速下降法解方程组
3.5 0 ⎤ ⎡ x1 ⎤ ⎡ 24 ⎤ ⎡4 ⎢3.5 4 − 1.5⎥ ⎢ x 2 ⎥ = ⎢ 30 ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ 0 − 1.5 4 ⎥ ⎢ x 3 ⎥ ⎢− 24⎥ ⎣ ⎦⎣ ⎦ ⎣ ⎦