预条件I+S+R下的AOR迭代方法

新的L-矩阵线性方程组的预条件AOR迭代法李园;韩海山【摘要】在预条件矩阵Pa=(I+Sa)和Pαβ=(I+Sαβ)的基础上提出一个新的预条件矩阵为Pαβ=(I+Sαβ)的预条件AOR迭代法,建立了新的预条件AOR迭代法与经典的AOR迭代法的比较定理,数值试验表明预条件AOR迭代法更为有效.【期刊名称】《湖北民族学院学报（自然科学版）》【年(卷),期】2012(030)001【总页数】6页(P39-44)【关键词】计算数学;预条件;线性方程组;AOR迭代法;谱半径;L-矩阵【作者】李园;韩海山【作者单位】内蒙古民族大学数学学院,内蒙古通辽028043;内蒙古民族大学数学学院,内蒙古通辽028043【正文语种】中文【中图分类】O241.6本文考虑如下线性方程组：Ax=b，(1)其中:A∈Rn×n为非奇异矩阵，b∈Rn为已知向量，x∈Rn为未知量.预条件方法通常是找到合适的预条件算子P，将Ax=b等价地转化为如下预条件线性方程组:PAx=Pb，其中P∈Rn为非奇异预条件矩阵.文献[1]中给出了预条件矩阵为Pα=I+Sα的预条件AOR迭代法，其中:α是参数.文献[2]中改进了上述迭代法，给出了预条件矩阵为Pαβ=I+Sαβ的预条件AOR迭代法，该预条件也是文献[3]中预条件的推广，其中:(2)α，β是参数，当β=0时，Sαβ=Sα.本文给出了预条件为的预条件AOR迭代法，其中:α，β，σ=(σ1,σ2,…,σn-1)T是参数，当σ=(σ1,σ2,…,σn-1)T=(0,0,…,0)T时，本文建立了新的预条件AOR迭代法与文献[1]和文献[2]以及经典的AOR迭代法的比较定理.通过比较定理，得出本文提出的预条件方法比文献[1]和文献[2]以及经典的AOR迭代法更有效.为方便起见，令A=I-L-U，其中I是单位矩阵，-L和-U分别是矩阵A的严格下三角和严格上三角矩阵，则线性方程组(1)相应的预条件AOR方法的迭代矩阵为：Lγω=(I-γL)-1[(1-ω)I+(ω-γ)L+ωU]，(3)其中：ω和γ是参数，ω≠0.记本文讨论的预条件线性系统为:(4)其中:令其中Sαβ为式(2)，若记SαβU=DS+ES,RL=DR+FR，其中DS和ES分别为SαβU的对角和严格下三角矩阵，DR和FR分别为RL的对角和严格下三角矩阵，则有:I-L+Sαβ-(U-R+RU)-DS-ES-DR-FR=(I-DS-DR)-(L+ES+FR-Sαβ)-(U-R+RU)=其中:则预条件AOR方法的迭代矩阵为：1 预备知识为方便起见，本文给出如下记号.设A=(aij)∈Rn×n为n×n实矩阵.diag(A)为矩阵A的对角元素aii(i=1,2,…,n)构成的n×n对角矩阵.对于任意矩阵A=(aij)，B=(bij)∈Rn×n，称A≥B，如果对所有i,j=1,2,…,n，成立aij≥bij.若矩阵A的每一个元素aij≥0，i,j=1,2,…,n，则称矩阵矩阵A是非负矩阵，记作A≥0.称A-B≥0当且仅当A≥B.对于n维向量也有类似的定义，ρ(·)表示矩阵的谱半径.定义1[3] 矩阵A=(aij)∈Rn×n，1)若aij≤0，i≠j，i,j=1,2,…,n，则称矩阵A为Z-矩阵；2)若A∈Z且aii≥0，i=1,2,…,n，则称A为L-矩阵；3)若A∈Z非奇异且A-1≥0，则称A为M-矩阵.定义2[4] 设矩阵M∈Rn×n非奇异，矩阵分裂A=M-N称为：1)收敛的，如果ρ(M-1N)<1；2)正则分裂，如果M-1≥0且N≥0；3)弱正则分裂，如果M-1≥0且M-1N≥0；4)M-分裂，如果M是M-矩阵且N≥0.定义3[5] 矩阵A=(aij)∈Rn×n称为可约的，如果存在n×n阶置换矩阵P，使得:其中:A11是r×r阶矩阵，A22是(n-r)×(n-r)阶矩阵，1rn.如果矩阵A不是可约的，则称A不可约.引理1[5](Perron-Frobenius) 若矩阵A=(aij)∈Rn×n为非负不可约矩阵，则:1)ρ(A)为矩阵A的一个正特征值；2)对于ρ(A)，相应地存在一个正的特征向量x>0；3)ρ(A)是矩阵A的一个单特征值；4)ρ(A)随矩阵A的任一元素增加而增加.引理2[6] 如果矩阵A是一个L-矩阵，则A是M-矩阵当且仅当存在一个正的向量x>0，使得Ax>0.引理3[7] 设A=(aij)∈Rn×n为非负矩阵，则：1)若存在x≥0，x≠0满足Ax≥αx，则ρ(A)≥α，进一步地，若Ax>αx，则ρ(A)>α；2)若存在x≥0，x≠0满足Axβx，则ρ(A)≤β，进一步地，若Ax<βx，则ρ(A)<β；3)若A不可约并且有0≠αx≤Ax≤βx，αx≠Ax和Ax≠βx对某一非负向量x成立，则α<ρ(A)<β且x是一个正向量.引理4[8] 设A=M-N是M-分裂，则ρ(M-1N)<1当且仅当A是非奇异M-矩阵.2 主要结论本文的证明需要用到下面的定理：定理1 设矩阵A∈Rn×n是一个M-矩阵，且A的元素满足0<a1nan1<1，1-σiai,i+1ai+1,i>0(i=1,2,…,n-1)，则也是M-矩阵.证明设则有：因为A是M-矩阵，aij≤0，i≠j，且aii=1，所以当i=1,2,…,n-1；j=1,2,…,n时，aij-σiai,i+1ai+1,j≤0.又因为可以得到所以：；j=2,…,n-1，并且由已知1-σiai,i+1ai+1,i>0，i=j=1,2,…,n-1，因此也是L-矩阵.因为A是M-矩阵，由引理2知，存在向量x>0，满足Ax>0，所以再根据引理2得是M-矩阵.根据定理1，可以建立下面的比较定理：定理2 设A∈Rn×n为非奇异的L-矩阵，Lγω和分别是线性方程组(1)和(4)相应的AOR方法的迭代矩阵.若0≤γ≤ω≤1(ω≠0,γ≠1)且：α>1，β∈∩，σi∈[0,1](i=1,2,…,n-1)，矩阵A的元素满足0<a1nan1<1，1-σiai,i+1ai+1,i>0(i=1,2,…,n-1)，则有：(i)若ρ(Lγω)<1，则(ii)若A为不可约矩阵，那么或证明 (i)令：因为A为非奇异L-矩阵，且0≤γ≤ω≤1(ω≠0,γ≠1)，即是非奇异M-矩阵，F≥0，所以A=E-F是M-分裂，由引理4知，ρ(Lγω)<1，又因为A为非奇异M-矩阵，由定理1知也是非奇异的M-矩阵.又因为所以因此(ii)设A=I-L-U是不可约矩阵，因为：Lγω=(I-γL)-1[(1-ω)I+(ω-γ)L+ωU]=(1-ω)I+ω(1-γ)L+ωU+T，其中：T=(I-γL)-1γL[ω(1-γ)L+ωU]≥0，可知Lγω是不可约矩阵，由引理1，存在一个向量x>0，使得Lγωx=λx，其中λ=ρ(Lγω)，则有：[(1-ω)I+(ω-γ)L+ωU]x=λ(I-γL)x，上式等价于： [(1-ω-λ)I+(ω-γ-λγ)L+ωU]x=0和：(λ-1)(I-γL)x=ω(L+U-I)x，于是：(ω-γ+λγ)(L+ES+FR-Sαβ)+ω(U-R+RU)]x=ωU-(1-ω-λ)(DS+DR)+(ω-γ+λγ)(ES+FR-Sαβ)+ω(RU-R)]x=(-γ)-1[(λ-1)(DS+DR)+γ(λ-1)(ES+FR-Sαβ)+ω(DS+DR+ES+FR-Sαβ+RU-R)]x=(-γ)-1[(λ-1)(DS+DR)+γ(λ-1)(ES+FR-Sαβ)+ω(R+Sαβ)(L+U-I)]x=(-γ)-1[(λ-1)(DS+DR)+γ(λ-1)(ES+FR-Sαβ)+(λ-1)(R+Sαβ)(I-γL)]x=(-γ)-1[(λ-1)(DS+DR)+γ(λ-1)(ES+FR-Sαβ-RL)+(λ-1)(R+Sαβ)]x=(-γ)-1[(λ-1)DS+(λ-1)DR+γ(λ-1)(ES-DR-Sαβ)+(λ-1)(R+Sαβ)]x=(-γ)-1[(λ-1)DS+(λ-1)DR+γ(λ-1)(ES-Sαβ)-γ(λ-1)DR+(λ-1)(R+Sαβ)]x=γ(λ-1)(ES-Sαβ)+(λ-1)(R+Sαβ)]x=(-γ)-1(λ-1)[DS+(1-γ)DR+γES+(1-γ)Sαβ+R]x若λ<1，则即由引理3，得若λ>1，则即由引理3，得由上述定理可得如下推论：当ω=γ时，AOR迭代法即为超松弛(SOR)迭代法，相应地有如下结论：推论1 设A∈Rn×n为非奇异的L-矩阵，Lω和分别是线性方程组(1)和(4)相应的SOR方法的迭代矩阵.若0<ω<1且α>矩阵A的元素满足0<a1nan1<1，1-σiai,i+1ai+1,i>0(i=1,2,…,n-1)，则有：(i)若ρ(Lω)<1，则(ii)若A为不可约矩阵，那么或当ω=1，γ=0时，AOR迭代法即为雅可比(Jacobi)迭代法，相应地有如下结论：推论2 设A∈Rn×n为非奇异的L-矩阵，L和分别是线性方程组(1)和(4)相应的Jacobi方法的迭代矩阵.若α>矩阵A的元素满足0<a1nan1<1，1-σiai,i+1ai+1,i>0(i=1,2,…,n-1)，则有：(i)若ρ(Lγω)<1，则(ii)若A为不可约矩阵，那么ρ(Lγω)<1或3 数值举例运用文献[2]中的例子，将本文的新算法与已知算法进行比较，来说明新算法更为有效.文献[2]中线性方程组(1)的系数矩阵为：表1 几种预条件AOR迭代法迭代矩阵谱半径的比较Tab.1 The comparison of the iterative matrices′ spectral rad ius for preconditioned AOR iterative methodsωγαβρ(Lγω)ρ(L～γω)ρ(L′γω)ρ(L＾γω)0．90．82-0．081．28881．30171．30911．37060．70．63-0．121．18811．19311．19951．23330．60．54-0．161．14871．15151．15801．18260．50．45-0．21．11491．11651．12171．1402表2 几种预条件SOR迭代法迭代矩阵谱半径的比较Tab.2 The comparison of the iterative matrices′ spectral radius for preconditioned SOR iterative methodsωγαβρ(Lγω)ρ(L～γω)ρ(L′γω)ρ(L＾γω)0．90．92-0．071．31851．33351．34101．41520．60．63-0．121．16121．16551．17101．2000表3 几种预条件Jacobi迭代法迭代矩阵谱半径的比较Tab.3 The comparison of the iterat ive matrices′ spectral radius for preconditioned Jacobi iterative methodsωγαβρ(Lγω)ρ(L～γω)ρ(L′γω)ρ(L＾γω)102-0．081．17671．18241．18561．2067103-0．151．17671．18041．18641．2075用和分别表示经典的预条件AOR迭代法、文献[1-2]以及本文所提出的预条件AOR 迭代法中迭代矩阵的谱半径.取σ=(σ1,σ2,σ3,σ4,σ5)T=(0.9,0.8,0.7,0.6,0.5)T，表1给出了当ω和γ取不同值时新的预条件AOR迭代法与文献[1]和文献[2]以及经典的AOR迭代法的迭代矩阵谱半径之间的大小关系.表2和表3分别给出了新的预条件AOR迭代法中当ω和γ取特定值时相应的预条件SOR迭代法、预条件Jacobi迭代法与文献[1]和文献[2]以及经典的预条件SOR迭代法、预条件Jacobi 迭代法的迭代矩阵谱半径之间的大小关系.由表1～3，当ω和γ取不同值时，可以看出本文所提出的预条件AOR迭代法、预条件SOR迭代法和Jacobi迭代法的收敛速度更快，这也正好验证了本文的结论.参考文献:[1] Li Y T，Li C X，Wu S L.Improvements of preconditioned AOR iterative methods forL-matrices[J].J Comput Appl Math，2007，206：656-665. [2] Wang H J，Li Y T.A New preconditioned AOR iterative methods for L-matrices[J]. J Comput Appl Math，2009，229：47-53.[3] Li Y T，Wang H J，Zhang C Y. New preconditioned AOR iterative methods for Linear System[J]. SEAMS Bull Math，2007，31：295-306. [4] Young D M. Iterative Solution of Large Linear Systems[M].London：Academic Press，NY，1971.[5] Varga R S. Matrix Iterative Analysis[M].NY：Prentice-Hall，Englewood Cliffs，1962.[6] 陈景良，陈向晖.特殊矩阵[M].北京：清华大学出版社，2001.[7] Gunawardena A D，Jain S K，Snyder L.Modified Iterative Method for Consintent Linear System[J].Linear Algebra Appl，1991，154/156：123-143.[8] Andreas F，Daniel B S.H-splitting and two-stage iterativemethods[J].Numer Math，1992，63：345-356.。

CG算法的预处理技术：、为什么要对A进行预处理：其收敛速度依赖于对称正定阵A的特征值分布特征值如何影响收敛性：特征值分布在较小的范围内，从而加速CG的收敛性特征值和特征向量的定义是什么？（见笔记本以及收藏的网页）求解特征值和特征向量的方法：Davidson方法：Davidson 方法是用矩阵( D - θI)- 1( A - θI) 产生子空间，这里D 是A 的对角元所组成的对角矩阵。

θ是由Rayleigh-Ritz 过程所得到的A的近似特征值。

什么是子空间法：Krylov子空间叠代法是用来求解形如Ax=b 的方程，A是一个n*n 的矩阵，当n充分大时，直接计算变得非常困难，而Krylov方法则巧妙地将其变为Kxi+1=Kxi+b-Axi 的迭代形式来求解。

这里的K(来源于作者俄国人Nikolai Krylov姓氏的首字母)是一个构造出来的接近于A的矩阵，而迭代形式的算法的妙处在于，它将复杂问题化简为阶段性的易于计算的子步骤。

如何取正定矩阵Mk为：Span是什么？:设x_(1,)...,x_m∈V ,称它们的线性组合∑_(i=1)^m?〖k_i x_i \|k_i∈K,i=1,2...m〗为向量x_(1,)...,x_m的生成子空间，也称为由x_(1,)...,x_m张成的子空间。

记为L（x_(1,)...,x_m），也可以记为Span（x_(1,)...,x_m）什么是Jacobi迭代法：什么是G_S迭代法：请见PPT《迭代法求解线性方程组》什么是SOR迭代法：什么是收敛速度：什么是可约矩阵与不可约矩阵？：不可约矩阵（irreducible matrix）和可约矩阵（reducible matrix）两个相对的概念。

定义1：对于n 阶方阵A 而言，如果存在一个排列阵P 使得P'AP 为一个分块上三角阵，我们就称矩阵A 是可约的；否则称矩阵A 是不可约的。

定义2：对于n 阶方阵A=(aij) 而言，如果指标集{1，2，...，n} 能够被划分成两个不相交的非空指标集J 和K，使得对任意的j∈J 和任意的k∈K 都有ajk=0, 则称矩阵 A 是可约的；否则称矩阵A 是不可约的。

牛顿迭代法解超越方程
牛顿迭代法是一种非常常用的数值计算方法，可以用来求解各种方程。

在解超越方程时，牛顿迭代法也是一种非常有效的方法。

超越方程是指含有超越函数的方程，比如三角函数、指数函数、对数函数等等。

这种方程一般很难直接求解，需要借助数值计算方法求解。

牛顿迭代法通过不断逼近函数零点来求解方程。

具体来说，就是从一个初始值开始，根据函数的导数和函数值，不断更新当前的值，直到满足一定的精度要求为止。

在解超越方程时，我们需要先将其转化为等价的形式，比如通过移项将方程变为 f(x) = 0 的形式。

然后选择一个初始值，根据牛顿迭代公式不断更新当前值，直到满足精度要求。

需要注意的是，在使用牛顿迭代法时，初始值的选择非常重要。

如果初始值选择得不好，可能会导致迭代过程不收敛，甚至发散。

因此，我们需要根据具体的方程性质和数值计算经验来选择初始值。

总之，牛顿迭代法是解超越方程的一种非常有效的数值计算方法。

通过不断逼近函数零点，可以在一定程度上克服超越函数的难题，为我们提供了一个可靠的求解方案。

- 1 -。

基于Matlab的数值分析实验的设计r——以矩阵的奇异值分解为例杨海涛【摘要】基于Matlab的数学实验,目的是为了提高普通高校学生的数学兴趣和动手能力,为培养解决实际数学问题的能力而设计的实验方法和步骤.利用矩阵的奇异值分解svd函数的数值算法设计了一种在实验课上用svd算法仿真的实验方案.【期刊名称】《内蒙古民族大学学报（自然科学版）》【年(卷),期】2016(031)006【总页数】4页(P465-468)【关键词】数学实验;数值计算;Matlab;奇异值分解【作者】杨海涛【作者单位】内蒙古民族大学数学学院,内蒙古通辽 028043【正文语种】中文【中图分类】O241.4利用计算机计算解决实际问题时需要建立数学模型，同时也需要相应数学软件的支持，为了更好的掌握数学软件这一有利的工具，本文通过Matlab中的矩阵的奇异值分解〔1〕svd函数设计了如何学习和掌握该函数的实验方法，从而可以进一步学习其它函数.目前，Matlab在计算仿真和图形图像的压缩处理〔2～4〕有着广泛的应用，希望这个实验方法能够给正在学习数学软件的学生起到抛砖引玉的作用.1.1 理论基础任意复长方矩阵的奇异值分解定理如下（Autonne-Eckart-Young定理）：〔5〕设A为一个m×n的实矩阵（m≥n），则存在正交矩阵V和U，使得：其中∑=diag(σ1,σ2,…,σr)，而且σ1≥σ2≥…≥σr≥0.奇异值分解在统计分析、信号与图像处理、系统理论和控制中有着广泛的应用. 1.2 SVD的数值算法〔6〕求解一般矩阵奇异值问题最常用的算法可分为两大类：一类是QR分解的方法；另一类是Jacobi旋转的方法.奇异值分解的QR分解法分为两个阶段：一是矩阵的二重对角化，即利用Householder变换把矩阵A变换成二重对角矩阵.这个过程就是追赶（chasing）过程〔7〕.二是QR分解，即保持二重对角矩阵形式不变，利用正交变换使对角线元素逐渐减小，使得矩阵接近对角矩阵.1.2.1 svd函数.对于matlab中的svd函数不是开源的，但通过调研得知该函数的算法就是采用propack中的lansvd函数的算法.其算法过程如下：对于任意的矩阵Am×n可以利用lanczos过程将其双对角化，即存在m阶正交矩阵Q和n阶正交矩阵W满足：其中为一双对角阵，且αi, βi∈R (i=1, 2, …, n)，这里假设了m＞n.上面的分解是逐步的利用Householder左右变换得到的.算法中利用了已知函数完成的双对角化过程.这里如果我们已经知道Bn的svd分解为其中Vn+1为n+1阶正交矩阵，的广义对角矩阵Un为n阶正交矩阵.则有令V=Qn+1Vn+1，其前r列为对应的左奇异向量，令U=WUn，其前r列为对应的右奇异向量.因此算法的关键在于计算矩阵Bn的svd分解.因为Bn是一个双对角矩阵，它的奇异值分解计算比较简单，运算量较小，利用Givens变换，算法复杂度仅为O(n).在我们这个程序中，直接调用了库函数svd 计算Bn的奇异值.1.2.2 svds函数.对于matlab中的svds函数，算法流程如下：（1）对于m×n阶矩阵A，构造(m+n)×(m+n)阶矩阵，其中O是稀疏矩阵. （2）利用eigs函数计算出B矩阵的特征值矩阵D，以及D对应的正交阵W. （3）根据输入参数，按要求输出最大特征值、最接近某参数的特征值或对特征值进行排序等.（4）根据设定的阀值，对矩阵进行筛选得到U矩阵，对矩阵筛选得到V矩阵，对D矩阵进行筛选得到S矩阵.（5）根据S中元素的大小对U、V、S进行排序.（6）输出结果.有如下J-W定理（Jordan-Wielandt定理）：若σ1≥σ2≥…≥σp-1≥σp是矩阵Am×n的奇异值（其中，p=min{m,n}），则上述B矩阵具有2p个非零特征值-σ1,…,-σp,σp,…,σ1和|m-n|个零特征值，与非零特征值（±σj）相对应的特征向量为若m≠n，另有特征向量分别对应于m＞n或m＜n.svds函数的基本思想就是构造出这样的矩阵B，并且求出其特征值和特征向量，然后根据函数重构的条件（所需奇异值的数量、最大最小等），求出所需的奇异值和对应的向量.2.1 矩阵生成生成10个不同的m×n的矩阵，这里取m=120 , n=50，对这10个矩阵进行svd分解.矩阵生成的方法是先做出矩阵的奇异值，然后用QR分解做出10个120阶正交阵和10个50阶正交阵，再利用奇异值分解公式生成10个矩阵A.10个矩阵的奇异值（均为25个）如表1所示.2.2 仿真结果2.2.1 计算时间svd和svds两种方法对十个矩阵进行奇异值分解消耗的时间如表2，单位是秒.可见，svd函数运行时间一般要小于svds函数运行时间.2.2.2 最小奇异值计算准确率.通过对比最小特征值的计算准确度，将每种算法特征值的最小值与矩阵最开始生成时的最小特征值进行比较.比较结果两种方法对10个矩阵的svd分解运算最小特征值的误差均几乎为0，都在10-15次方的量级左右，见表3.其中delta_min_svd和delta_min_svds是matlab程序中两个变量值.这个误差已经和截断误差非常接近了，因此可以认为两种方法都能较好的进行奇异值计算.2.2.3 左右奇异值向量计算准确率.由于svd分解后U、V矩阵的不唯一性，这里我们只看生成的U、V的正交性是不是依然比较完好.通过计算每一列和它的转置相乘减去单位阵也都在10-15这样的量级，因此svd和svds对于正交性的保持都是非常好的.综上所述，通过上述设计的方法和步骤，完全可以掌握matlab中的M-文件和.mat文件的使用方法，并且能熟练掌握奇异值分解的理论知识和具体计算方法，才能为以后svd应用打下基础，这样的实验方法才有实际的意义，才是真正的综合设计实验.本文的宗旨在于引导学生如何用数学软件完成一个实验设计，同时掌握这个实验的目的、相关内容实验的方法和手段.〔1〕尹芳黎，杨雁莹，王传栋，等.矩阵奇异值分解及其在高维数据处理中的应用〔J〕.数学的实践与认识，2011，45（15）: 171-176.〔2〕胡乡峰，卫金茂.基于奇异值分解（SVD）的图像压缩〔J〕.东北师大学报:自然科学版，2006，38（3）:36-39.〔3〕曾超，张陈，徐永利.基于奇异值分解的图像压缩及其matlab实现〔J〕.科技信息，2010，14:484.〔4〕向培素.一种自适应AP算法的matlab实现〔J〕.西南民族大学学报:自然科学版，2014，40（16）:877-882.〔5〕张贤达.矩阵分析与应用〔M〕.北京：清华大学出版社，2004.358-360. 〔6〕姜健飞，胡良剑，唐俭.数值分析及其matlab实验〔M〕.北京：科学出版社，2004.6，219-224.〔7〕山晓东，李园.一类不可约L-矩阵的预条件AOR迭代法〔J〕.湖北民族学院学报:自然科学版，2014，32（2）:128-132.。