（计算数学专业论文）关于加权最小二乘问题几种迭代方法的比较.pdf

复旦大学硕士毕业论文 a b s t r a c t 2 i nt h i sp a p e r ，s o m ei t e r a t i v em e t h o d sf o rt h ew e i g h t e dl e a s ts q u a r e sp r o b l e ma r ed i s c u s s e d t h e a l g o r i t h m s ，p r o p e r t i e sa n dn u m e r i c a le x a m p l e sa x es h o w n t h i sp a p e ri sd i v i d e di n t ot w o p a r t s t h ef i r s tp a r ti sa b o u tg s o r ( g e n e r a l i z e ds u c c e s s i v e o v e r r e l a x z t i o n ) m e t h o da n dp c g ( p r e e o n d i t i o n e dc o n j u g a t eg r a d i e n t ) m e t h o df o rt h en o r m a l e q u a t i o n s t h ec o n n e c t i o nb e t w e e nt h e mi si n v e s t i g a t e d i ti ss h o w nt h a tt h ep c gm e t h o di s a tl e a s ta sf a s ta st h eg s o rm e t h o d n u m e r i c a le x a m p l e sd e m o n s t r a t et h a tt h ep c g m e t h o di s m u c hf a t e rt h a nt h eg s o rm e t h o d t h es e c o n dp a r ti n t r o d u c e s w g m r e s ( w e i g h t e dg m r e s ) m e t h o d a n d p w g m r e s ( p r e e o n d i t i o n e w g m r e s ) m e t h o d w ef i r s tp o i n to u tas c a l i n gi n v a r i a n tp r o p e r t yo ft h e s em e t h o d s ，t h e nw e d i s c u s st h ep e r f o r m a n c eo ft h ew g m r e s ，t h ep w g m r e sh a sn o ts og o o d p e r f o r m a n c ec o m p a r e d t op g m r e s k e y w o r d s ：l e a s ts q u a r e sp r o b l e m ；w e i g h t e dl e a s ts q u a r e sp r o b l e m ；i t e r a t i v em e t h - o d s ；g s o r ；p c g ；k r y l o vs u b s p a c e s ；a r n o l d i ；g m r e s ；w g m r e s 第一章引言 现在许多科学领域都要求解线性方程组 a z = b 这里a 霹“和b r “，m n 。我们研究的问题是找向量o r “，使得方程组成 立如果方程个数多于未知量个数，即m 凡，我们称这个方程组a x = 6 是超定的。因 为b 必须是r ”的真子空间r ( a ) 的一个元素，所以通常的超定方程组没有准确解。这一 事实启发了我们，可以考虑对某些适当选取的p ，极小化i l a x 一6 ，这就是我们所谓的 最小二乘问题。 我们定义加权最小二乘问题的形式【3 为 叫n 胪一a z i | 一( 1 2 ) 即 哑n ；( 6 一a x ) t d 。( 6 一a 。) ( 1 3 ) 这里d 是一个对称正定矩阵。 在这里我们假定a 叼“且m n ，我们令余向量y = d _ 1 ( b a x ) ，那么解加 权最小二乘问题等价于解等式 = ( ：) ( 1 4 ) 此线性方程组有着广泛的应用，我们称( 1 4 ) 为k k t ( k a r u s h k u h n t u c k e r ) 系统，更多的 介绍请参考 8 在第一章节里我们考虑线性系统 a t0 x = s ， b 0 l f = z 可 m g d 为一不表阵矩用或 复旦大学硕士毕业论文 如果c = 0 ，这个线性系统就是k k t 系统。 解这个线性系统最常用的就是迭代法。如果a 是满秩的，那么方程( 1 5 ) 有唯一解+ 和x + 。在这篇文章中，我们假设a 是满秩的，那么 y + = d - 1 ( 6 一a x + )( 1 6 ) 5 第二章对于加权最小二乘问题的g s o r 方法和p c g 方法 2 1g s o r 方法 最近在 5 】中介绍了一种广义的s o r ( g s o r ) 方法。g s o r 方法包括参数p 和预条件 p 。在 6 中讨论了如何选取最优参数p 。下面我们来看具体的g s o r 算法。 我们先来考虑一个简单的2 2 矩阵系统 卜北) = b 1 ) 对于这个系统，我们给出修正的s o r ( m s o r ) 方法 i 0 。p + 1 = “j 1 ( d 茁+ b 1 ) + ( 1 一u 。) z p z 字十1 = 叫2 ( 。1 2 + 1 1 + b 2 ) + ( 1 一u 2 ) 。扩 我们将2 x2 矩阵的系数矩阵分裂为 一a 口1 ( 2 2 ) ( 2 3 ) 适当地选择参数p l ，j d 2 ，使( 1 + p 1 ) 0 ，( 1 + 见) 0 ，因此对于系统( 21 ) 我们可以定义一 种迭代方法 这里 令 z 掣= 而1 ( 。z 岫) + 去z i z = 击( 触i ) + 6 2 ) + 惫z 垆 l 姚2 而江1 ，2 ( 2 , 4 ) ( 2 5 ) ( 2 6 ) o 阮 p 0 ，、 一 、， 0 p +m h 邓 ， i | 、 n恤陆2 o o ，、 血 纯 p o ， + 、 ， = 、 + + 晴，“”2 茁 z ，t、 0 p +m ， 邓 ， 复里盔堂矍圭坐些堡壅j 7 一 l ：1 一u ，i = 1 ，2 ( 2 7 ) 1 + p i 从( 2 6 ) 和( 27 ) ，我们可以看到( 2 2 ) 和( 2 5 ) 是相同的。因此迭代方法( 2 4 3 是基于 f 23 1 的分裂，这一点是同m s o r 方法相同的。于是对于系统( 1 5 ) ，我们给出了一个新的 迭代方法。首先我们将系统( 1 5 ) 的系数矩阵按以下方法分裂 ( d = l + p ) d p 0 ) _ ( p 言7 ) 江 这里p 一1 ，且p 是一个非奇异矩阵因此我们定义对于系统( 1 5 ) 的迭代方法为 ( 二p ) 。? ) ( ：) = ( ：) + ( 9 d 。一p a ) ( ”x ( k ) ) c 2 9 ， 因为( 2 8 ) 的分裂类似于( 2 3 ) 的分裂，我们定义迭代方法( 2 9 ) 为广义s o r ( g s o r ) 方法 下面我们给出此迭代方法的完整算法： 算法1 ：g s o r 方法 i 选择z 霉k 凡，器。冗和预条件矩阵p 2 对于k ：o ，1 ，直到收敛 ( 2 1 ) ( 2 2 ) g 美茹k = d 一1 6 十p ! 龆o r d 一1 a z 龆。r ( 1 + p ) 。峨k ：茁龆。冗+ p - 1 。( c 一舻端k ) 结束算妆 g s o r 方法的一个主要性质是它结合了预条件的技术，因此可以认为它是更有效的 另一方面，如果预条件矩阵p 正确选择，g s o r 算法不仅有很快的收敛速度，而且适合于 并行算法，对于等式 a t d 一1 a 。；a t d 一。b c ( 2 1 0 ) g s o r 方法能替代一些著名的方法，比如g s ，s o r ，s s o r p ，a o r 4 ，s a o r 2 详情请参考【6 6 。 复旦大学硕士毕业论文 2 2 1 共轭梯度( c g ) 法 2 2p c g 方法 共轭梯度法是近年来常用的一种有效迭代方法。 因为线性最t j 、_ _ - - 乘问题 i i d a x lj = r a i n 等价于求解相应的法方程组 a r a o ：a t df 2 1 1 ) 按照实际问题的大多数情况，假定a 是列满秩的，即a r ? “，因而( 2 1 1 ) 的系数矩阵 a 丁a 是对称正定的，所以凡是使用于系数矩阵为对称正定矩阵的线性方程组的迭代解法， 都可以用来求解( 1 1 ) 而共轭梯度法对求解它更为简便和有效。 这里我们假设( 1 1 ) 方程组的系数矩阵a 是对称正定的下面我们给出对于方程组( 1 1 ) 的完整的共轭梯度( c g ) 算法 1 1 】： 算法2 ：c g 方法 1 选择茁( 叭，计算r ( o ) = p ( o ) = b a x ( o ) 2 对于k = 0 ，1 ，直到收敛 结束算法 ( 2 1 )o = 。，( ，r ( 。) ( r ( m ，r ( 。) 伊丽【2 硒f 石巧 ( 2 2 )z ( 。+ 1 ) = z ( 七) + 。k p ( 七) ( 2 3 )r ( + 1 ) = r ( “) 一。k a p ( ) ( 2 - 4 ) 风= 一黼 ( 2 5 )p ( 。+ 1 ) = r ( + 1 ) + 尻p ( 七) 舻 8 糍，r 复旦大学硕士毕业论文 2 2 2 预条件共轭梯度( p c g ) 法 r e i d 证明了把共轭梯度法作为迭代算法时，它是十分有效的方法。然而该迭代在有限 精度计算中搜索向量会失去正交性，对于病态的问题，过程收敛速度极慢。 避免这个困难的一个途径是对a 进行预条件处理。这个意思是说寻求一个非奇异矩阵 c ，使得a = c - 1 a c 。的条件数比a 好，然后可把方程转化为 盈：b 并用共轭梯度法求解( 具有改善的收敛性质) ，其中茁= c _ 窑和5 = c l b 下面给出代入 后的共轭梯度法： 1 选择量( ，计算f ( o ) = o ) = 5 一五童( o ) 2 对于k = 0 ，1 ，直到收敛 结束算法 ( 2 1 ) ( 2 2 ) ( 23 ) ( 2 4 ) ( 2 5 ) a t = 斋 ( 七十1 ) = 童( + o 茸) 爿+ 1 ) = f ( ) 一a e 一1 a c 一1 ) 耻帮 疹+ 1 ) = 一+ 1 + 凤纠 我们定义 代入以上算法碍到 q = c 2 p ( 2 ) = c 一1 矿) o ( 女) ：c 一1 i ( ) z ( ) = c 一1 o ) r ( k ) ：6 一a z ( 女) ：e 再女) 9 复旦盔堂塑坐些堡塞 1 0 算法3 ：p c g 方法 1 选择z ( m ，计算r ( o ) = b a z ( 叭，p ( o ) = z ( o ) = q 一1 r ( o ) 2 对于a = 0 ，1 ，直到收敛 慨，d t = 蔫端 ( 2 2 )z ( + 1 ) = z ( ) + o p ( ) ( 2 ，3 )r ( + 1 ) = r ( ) 一。女却( 。) 结束算法 ( 2 4 ) ( 2 5 ) ( 2 6 ) z ( k + n ：q 一1 r ( 。+ 1 ) 凤= 莆 p ( + 1 ) = 2 ( + 1 ) + 凤p ( ) 这就是所谓预条件共轭梯度算法，对称正定矩阵q 称为预条件矩阵。一个好的预条件 矩阵能导出十分快的收敛性。 2 2 3 对于加权t , j , - - 乘问题的p c g 算法 我们将线性系统( 1 5 ) 写成等式形式 a t d l a x = a t d 一1 b c 定义q 为一对称正定矩阵且为a t d 1 a 的预条件，以下给出p c g 方法 算法4 ：对于加权最小二乘问题的p c g 方法 1 ，选择z 龆g ，计算r ( o ) = a t d l b a t d 一1 a z 龆g c ，p ( o ) = z ( o ) = q - l r ( o ( 2 1 2 ) 复旦大学硕士毕业论文 2 对于k = 0 ，l ，直到收敛 结束算法 ( 2 1 ) ( 2 2 ) ( 2 3 ) ( 2 4 ) ( 2 5 ) ( 2 6 ) f r ( 、z ( k ) 1 钒2 硒f 万矿噶莉 z 煳= z 龆g + d p ( ) r ( + 1 ) = r ( ) 一。女a t d 一1 a p ( ) z + 1 ) = q 一1 r 咔+ 1 ) 风= p ( 2 + 1 ) = 。( + 1 ) + f l k p ( ) p c g 算法中的迭代值z 龆g 满足以下的性质 茁龆g z g + s p a n z ( 叭，q ( 一1 ) a 丁d 一1 a z ( 0 1 ，( q ( 一1 ) a t d 一1 a ) k 一1 三c o ) = z 龆g + k k ( z ( m ，o 一1 a t d 一1 a ) 并且 ( k ) - - x * t d - - 1 a = 。；。龆。+ 州r a 。i n ，q 一) i x * - i i a d1 a a r 。一a 廿z 苫二d + k ( z t ，口圳 一 】 加权范数恻l m 定义为i i x l l m = ( x t m z ) ，m 为一对称正定矩阵。 2 3 主要结果 引理2 3 1 如果在算法中满足条件嚣。兄= d 一1 b ，z 般o r = 0 ，那么算法中的迭代向 量可g s 。r ，z 龆。r 满足： 这里 可龆。r d 一1 b + d 一1 a 凰一1 ( u ，p 一1 a t d 一1 a ) ，k 1 z 龆o r k k ( u ，p 一1 a 丁d 一1 a ) ，k 1 札= p 。1 ( c a r d 。1 6 ) ，凰( ，尸。a t d 1 a ) = 0 ) ( 2 1 3 ) 1 1 复旦大学硕士毕业论文 证明：我们用数学归纳法来证明此引理。当= 1 时，运用算法1 得到 攥o r ，( 1 ) 山g s o r d 一1 b + p 可器o r d 一1 a z 嚣o r 】( 1 + p ) = 【d 一1 b + p d 一1 b ( 1 + p ) = d 一1 b z 嚣。r + p 一1 ( c a t 攥。咒) = p 一1 ( c a 7 d 一1 6 ) = “ 因此，当= 1 时引理成立。 假设对于某个k ，以上引理成立，那么运用算法1 得到 鉴芜k = d 一1 b + p 嚣。凡一d 一1 a 。龆o r 0 + p ) = d 一1 b + p d 一1 b + p d 一1 a v l 一d 一1 a v 2 ( 1 + p ) 通过假设我们得知 v l k k l ( u ，p 一1 a t d 一1 a ) ”2 。( ，p 一1 a r d 一1 a ) 那么 l 嬲k = d 一1 b + p d 一1 b + p d 一1 a 1 一d 一1 a v 2 】( 1 + p ) d - 1 b + d 。1 a 瓤( 扎，p 1 a t d 。1 a )( 2 1 4 ) 同样通过算法1 得到 。瞄k = 。龆o r + p - 1 ( c a 丁嬲k ) = l + p 一1 ( c a t d 一1 b a t d 一1 2 ) 假设通过并且 1 k k ( u ，p 一1 a t d 一1 a ) ， 2 a 凰( 让，p 一1 a t d 一1 a ) 因为 p 一1 a t d 一1 口2p 一1 a r d 一1 4 甄( u ，p 一1 a 丁d 一1 a ) k k + 1 ( u ，p 一1 a t d 一1 a ) 1 2 复旦大学硕士毕业论文 因此 z 兽茹k = l + p 一1 ( c a t d 一1 b a t d l u 2 ) = ”l + u 一尸一1 a r d 一1 v 2 k k + 1 ( u ，p 一1 a 丁d 一1 a ) 所以当= 1 ，2 ，时，以上引理成立。口 下面给出主要的定理。 定理2 3 1 令z 龆。r 和z 龆g 分别为算法j 和算法4 第步的迭代向量，在算法中令 z 嚣。r = 0 并且口嚣d r = d b ，在算法4 中令茁垛g = 0 ，并且两个预条件p 和q 满 足：q = 一尸。那么 i i 。镍o r z + j l a r d t a i i x 龆g z + i l ，d 一，( 2 1 5 ) 因为 并且 因此 所以 证明：通过引理( 2 3 1 ) 和p c g 方法的性质我们得到 z 龆。兄k k ( u ，p 一1 a 丁d 一1 a ) z 龆g k k ( z ( 叭，q 一1 a t d 一1 a ) 媲g z + 协川一。虮( 胛，m q i 嘲n 啊。m 妒一吼r 川一 q 一1 ( a t d 一1 6 一c a r d 一1 a t 龆g ) = q - 1 ( a t d 一1 b c ) 一p 一1 ( a t d 一1 b c ) = “ k k ( u ，q 一1 a t d 一1 a ) = k k ( u ，p 一1 a 丁d 一1 a ) z 龆o r k k ( z ( 叭，q 一1 a t d - 1 4 ) z g - x * lr a r 。一- a 2 d 。( ：( 。，m 。i n ， ，。一。a ) | | z 一毋| | a t 。一t a l | z 咎；o r z + i i a t d 一。a 1 3 复旦大学硕士毕业论文 证毕。口 从以上的定理可以看出，如果不考虑每一步迭代的计算工作量和内存占有量，关于 4 t d - 1 a 范数，p c g 方法要优于g s o r 方法。对于两种算法，每一步迭代都要计算a x 和a r y ，并要解方程p y = ，和q ”= 厂。在算法1 的每一步迭代中，常量乘以向量共包含 2 m 步乘法；在算法4 的每一步迭代中，向量乘以常量加上内积乘法共需5 n 步乘法。因此 从计算量的角度来看，这两种方法是差不多的。然而p c g 方法对内存的要求要比g s o r 方法来得高，并且从并行计算的角度来看，p c g 方法每一步要计算内积，这一点使它没有 g s o r 方法来得受欢迎。从另一角度来看，g s o r 也包括一些额外的工作量，比如参数p 的选取。 此线性系统为加权最i j 、- 乘问题，因此最好对c s o r 方法和p c g 方法做直接的比较。 首先我们考虑加权最小二乘问题。以下令线性系统( 1 5 ) 中c = 0 。 定理2 3 2 在定理偿3 纠的条件下 1 | a 0 罂s o r 一6 i i d 一- i i a 圳。g b l i d 一，( 2 1 6 ) 证明：另z + 是加权最小二乘问题的唯一解，令b = b l + b 2 ，那么 d 一1 ，2 b = d 一b l + d 一b 2 这里d 一1 2 b 1 属于兄- 1 2 4 ) ，即d 一1 2 且的值域空间；d 一1 2 b 2 属于k e r ( ( d - 1 2 a ) t ) ， 即d _ 1 2 a 的零空间。因此我们得到； d 一1 2 a x + = d 一1 2 b 1 ，a x + = b 1 ( d 一1 2 a ) r d 一b 2 = a t d 一1 2 d 一1 2 b 2 = a r d b 2 = 0 ( d 一1 2 k ) t ( d - 1 2 b 1 ) = b d _ 1 6 1 = 0 对于任何。r “，存在 l i z z + i | j ，d 一， = ( z z + ) t a t d 一1 a ( z z 4 ) = ( a x a x + ) 丁d ( a z a r c + ) = ( a x b + k ) t d 一1 ( a x b + b 2 ) = ( a z 一6 ) t d 一1 ( n x b ) + b i d 一1 b 2 + ( a x 一6 ) t d 一1 b 2 + 6 ；d ( a z b ) 1 4 复旦大学硕士毕业论文 = i i a 。一b 临。一慨幅一。 我们用z 龆。咒和。龆g 来代替z 得到 通过定理( 2 3 1 ) z 龆。r z + i 晤，口一。 忙龆g z 4i l j ，d 一。 i l a z 龆。r 一6 l | 刍一。一1 1 6 。 刍一 | | a z 龆g 一6 幅一，一临一， a 龆。r 一6 i | d 一。| | a z 龆g 一6 | | d 一。 证毕口 若在线性系统( 1 5 ) 中b = 0 ，在引理( 2 3 1 ) 的条件下 饕品k e d 一1 a k i ( u ，p 一1 a r d 一1 a ) ，k 1 z 龆o r 甄( 让，p 一1 a t d 一1 a ) ，k 1 这里u = p c 。 因此，存在秽k k ( u ，p 一1 a r d 一1 a ) 使 g 墨芜k = 一d 一1 a 口 ( 2 1 7 ) 定理2 3 3 通过算法和算法4 相应地得到y 4 的近似嬲k 和9 磊，在定理俾3 j j 的 假设下，得到 l i g g 嚣k y + l i d l l g 苫占一f + i l d 证明：令w 凰( u ，p 一1 a r d 一1 a ) ，定义y = 一d a w 那么当y + = d 一1 ( b a x + ) 和b = 0 时，成立 1 5 ( y 一+ ) 7 d ( y y + ) = ( d 一1 a x + 一d 一1 a ) r d ( d 一1 a x + 一d 一1 a 叫) = ( 。+ 一w ) t a t d 一1 a ( x 一w )( 2 1 8 ) = l i x 一u i l j r d 一。a 根据( 21 7 ) ，蹦k = 一d 一1 a 秽，嚣占= 一d 一1 a x p c g ，并且毋，z 龆g ( u ，p - 1 a t d 一1 a ) 得到 复旦大学硕士毕业论文 证毕。口 = l l z 笋6 g z + i i j ，d 一。a = | | 毋一z | | j r d 一 l i z 龆g z l l j ，d 一。 = l l g 譬若占一+ j | 刍 2 4主要例子 在以下的例子中“一”表示g s o r 方法，“”表示p c g 方法。 例1 ：我们用m a t l a b ( 6 0 ) 命令：b = r a n d ( 2 0 ，1 0 ) 来取一个2 0 1 0 矩阵b ，并用 命令：b = r a n d ( 2 0 ，1 ) 得到一个有2 0 个元素的向量。那么，现在的问题是：找9 4 ，使得： 恸+ 一惦一一伽m i n ，。恸一惦一 ( 2 - 1 9 ) 令b 1 是矩阵b 的最初1 0 行，那么b 1 是一个1 0 1 0 矩阵令a = b b f l ，以上的问题 就等价于问题：找岱，使得： l | a z + 一6 l | 刍一= m i 甚ol a z 一6 l l 刍一- ( 2 2 0 ) 旷和z + 有关系：。+ = b 1 矿。实际上问题( 2 1 9 ) 和问题( 22 0 ) 是等价的，但是后者有更 好的条件数。我们把g s o r 方法( 算法1 ) 和p c g 方法( 算法4 ) 运用到问题( 2 2 0 ) 上。 p c g 方法的预条件q 为a r d _ 1 a 的对角元部分；g s o r 方法的预条件p = 一q 。我们用 m a t l a b 命令d = s p a r s e ( 1 ：2 0 ，1 ：2 0 ，b ，2 0 ，2 0 ，z o ) 产生一个对称正定矩阵d 。g s o r 方法的最优参数p 的选择根据参考文献【6 】。我们令z ( 2 ) 为算法1 和算法4 第k 步的迭代 值，那么相应的余向量为r ( ) = b a x ( 舢。以下的图1 显示了两种算法迭代次数和余向量 d _ 1 范数的比值 从图中我们可以看出，p c g 方法余向量的口_ 1 范数要小于g s o r 余向量的d 。范 数，并且p c g 方法余向量的d _ 1 范数很快就达到了它的最小值。这个例子证明了我们的 1 6 复旦大学硕士毕业论文 分析，p c g 方法的收敛速度要比g s o r 方法的收敛速度快 图 1 例2 ：这个例子我们用表格及具体的数字来说明我们的结论：p c g 方法的收敛速度要 比g s o r 方法的收敛速度快。同样用m a t l a b ( 6 0 ) 命令：b = r a n d ( 2 0 ，1 0 ) 来取一个2 0 1 0 矩阵b ，并用命令：b = r a n d ( 2 0 ，1 ) 得到一个有2 0 个元素的向量令b l 是矩阵b 的最初l o 行，那么b 。是一个1 0 1 0 矩阵。p c g 方法的预条件q 为a t d _ 1 4 的对角元部分；g s o r 方法的预条件p = - q 我们用m a t l a b 命令d = s p a r s e ( 1 ：2 0 ，1 ：2 0 ，b ，2 0 ，2 0 ，2 0 ) 产生一个对称正定矩阵d 。g s o r 方法的最优参数p 的选择根据参考文献【6 。 迭代次数g s o r 方法余向量的d _ 1 范数p c g 方法余向量的d _ 1 范数 12 6 3 8 41 4 8 5 1 22 8 1 9 91 2 4 7 6 32 0 7 3 80 7 8 3 6 1 7 复旦大学硕士毕业论文 迭代次数g s o r 方法余向量的d _ 1 范数p c g 方法余向量的d _ 1 范数 422 9 7 7 0 5 9 8 1 51 7 6 7 30 5 5 6 9 61 5 2 5 5 04 5 6 3 72 1 0 6 80 4 3 2 4 81 7 8 5 7 0 4 3 0 5 91 3 5 8 6 0 4 2 9 1 1 01 3 7 5 3 0 4 2 8 7 1 1 1 0 9 5 3 0 4 2 8 7 1 21 0 5 2 3 04 2 8 7 1 31 0 5 4 4 0 4 2 8 7 1 40 6 6 8 9 04 2 8 7 1 50 5 9 8 1 0 4 2 8 7 1 60 6 9 9 9 04 2 8 7 1 70 5 9 3 304 2 8 7 1 80 6 0 5 0 0 4 2 8 7 1 90 5 7 2 3 0 4 2 8 7 2 00 4 8 6 8 0 4 2 8 7 2 10 5 5 3 5 0 4 2 8 7 2 20 5 5 6 6 0 4 2 8 7 2 30 4 9 8 8 0 4 2 8 7 2 40 4 9 6 60 4 2 8 7 2 50 4 7 7 20 4 2 8 7 图 2 例3 ：用m a t l a b ( 60 ) 命令：b = r a n d ( 2 0 0 ，1 0 0 ) 来取一个2 0 0 x1 0 0 矩阵b ，并用命 令：b = r a n d ( 2 0 0 ，1 ) 得到一个有2 0 0 个元素的向量。令b l 是矩阵b 的最初i 0 0 行，那么 复旦大学硕士毕业论文 b 1 是一个1 0 0 x1 0 0 矩阵。p c g 方法的预条件q 为4 丁d _ 1 a 的对角元部分；g s o r 方法 的预条件p = 一q 。我们用m a t l a b 命令d = s p a r s e ( 1 ：2 0 0 ，1 ：2 0 0 ，b ，2 0 0 ，2 0 0 ，2 0 0 ) 产生一个对称正定矩阵d 。g s o r 方法的最优参数p 的选择根据参考文献6 1 。 从图中我们可以看出，p c g 方法余向量的d _ 1 范数要小于g s o r 余向量的d _ 1 范 数，并且p c g 方法余向量的d - 1 范数很快就达到了它的最小值。这个例子也证明了我们 的分析，p c g 方法的收敛速度要比g s o r 方法的收敛速度快。 图 3 例4 ：此例子中的所有命令和例3 中的命令相同。 通过具体的数据我们论证了我们的结论。 迭代次数g s o r 方法余向量的d _ 1 范数p c g 方法余向量的d _ 1 范数 13 1 8 3 7 25 6 6 0 9 55 1 2 3 3 84 8 6 1 8 1 9 复旦大学硕士毕业论文 迭代次数g s o r 方法余向量的d - 1 范数p c g 方法余向量的d - 1 范数 1 03 7 6 5 8 23 2 9 0 0 1 55 9 2 1 1 02 3 6 7 4 2 04 7 2 0 4 61 9 9 8 8 2 55 9 9 8 1 21 ，9 0 7 5 3 04 1 9 3 8 11 8 7 5 3 3 54 9 1 7 5 11 8 6 0 7 4 02 7 7 2 7 61 8 5 3 4 4 53 8 1 9 6 01 8 5 0 9 5 02 3 6 3 9 51 8 4 9 9 5 53 7 5 4 9 21 8 4 9 5 6 03 0 2 7 5 01 8 4 9 3 6 52 8 0 1 5 01 8 4 9 2 7 02 8 8 1 3 81 8 4 9 2 7 52 4 8 1 9 61 8 4 9 2 8 02 2 7 3 6 11 8 4 9 2 8 51 9 2 1 2 31 8 4 9 2 9 01 6 9 1 9 31 8 4 9 2 9 52 58 1 3 618 4 9 2 1 0 01 99 4 2 018 4 9 2 这个图表和我们的理论分析一致。 图4 例5 ：用m a t l a b ( 6 0 ) 命令：b = r a n d ( 4 0 0 ，2 0 0 ) 来取一个4 0 0 2 0 0 矩阵b ，并用命 令：b = r a n d ( 4 0 0 ，1 ) 得到一个有4 0 0 个元素的向量令b l 是矩阵b 的最初2 0 0 行，那么 b l 是一个2 0 0 2 0 0 矩阵。p c g 方法的预条件q 为4 t d _ 1 a 的对角元部分；g s o r 方法 的预条件p = 一q 。我f f nm a t l a b 命令d = s p a r s e ( 1 ：4 0 0 ，l ：4 0 0 ，b ，$ o o ，4 0 0 ，4 0 0 ) 复旦大学硕士毕业论文 产生一个对称正定矩阵d 。g s o r 方法的最优参数p 的选择根据参考文献 6 】。 从图中我们可以看出，p c g 方法余向量的d - 1 范数要小于g s o r 余向量的d _ 1 范 数，并且p c g 方法余向量的d - 1 范数很快就达到了它的最小值。这个例子也证明了我们 的分析，p c g 方法的收敛速度要比g s o r 方法的收敛速度快 例6 ：令 a = 001 0 3 010 100 000 000 000 图 5 b = 1 0 s 1 0 0 0 0 c * 2 1 复旦大学硕士毕业论文 并且 h = d i a 9 ( 151 05 01 0 01 0 4 ) ，p _ h 圹黼 d=hph p c g 方法的预条件q 为a r d _ 1 a 的对角元部分；g s o r 方法的预条件p = 一q 。我们 令z ( 2 ) 为算法1 和算法4 第k 步的迭代值，那么相应的余向量为r ( 2 ) = b a x ( m 。以下 的图显示了两种算法迭代次数和余向量d _ 范数的比值 这个例子再次证明了p c g 方法的收敛速度要比g s o r 方法的收敛速度快。 图 6 2 2 第三章w g m r e s 方法和p w g m r e s 方法 在此章节中，我们讨论线性方程组( 1 1 ) ，且令a 是一个非奇异的实矩阵。一般来说， 这样的系统是高阶，稀疏和非对称的。 3 1 加权a r n o l d i 方法 为了解这样的系统，s a a d 和s c h u l t z 提出了g m r e s 方法 1 2 。此方法运用a r n o l d i 方法来构造k r y l o v 子空间 ( a ，v 0 ) = s p a n ( v o ，a v o ，a m - 1 v o ) 的正交基v ，m = p 1 7 一，v 。 ，这里z o 是给定的初始向量，v o = b a x o 。 以下的算法给出了a r n o l d i 方法，它构造了k r y l o v 子空间的一组正交基。此算法运用 了修正g r a m s c h m i d t 算法。 算法5 ：a r n o l d i 方法 1 计算i 1 = v l t v l l 2 2 对于j = 1 ，2 ，直到m ( 2 1 )w = a f j ( 2 2 )i = 1 ，j ( 2 2 1 ) h i d = ( w ，矾) 2 ( 2 2 2 )叫= w h i j v i ( 2 3 ) h j + 1 ，j = 怕恢i fh j + l d = 0 ，s t o p ( 2 4 ) v j + l = 毒 结束算法 复旦大学硕士毕业论文 令豆。尺“为h e s s e n b e r g 矩阵，它的非零元为a r n o l d i 方法产生的值无”。我们 定义矩阵百。尺( m + 1 ) m 为 。) _ ( 无。：。三) 簖吒= k a v = v ，m + l h m 瓦= 蠕a ( 3 1 ) ( 3 2 ) ( 3 3 ) 在这篇文章中，新方法用( ，一) 口来代替( ，) 。，d 为一对角矩阵，我们称新方法为加权 a r n o l d i 方法。加权a r n o l d i 方法加快了余向量的收敛速度，构造了k r y l o v 子空间k 。( a ， ) 的一组d 正交基。 在给出加权a r n o l d i 方法之前，我们先给出d 内积和d 范数的定义。 令d = d i a y ( d 。，d 2 ，厶) 为一对角阵，我们定义d 为实向量，d = ( d 1 ，d 2 ，d n ) 1 ， 其中出 o , j = 1 ，n ，那么d 内积定义为： ( u ， ) d = v t d u = e 坠1 d u t v 。 v u ，u r “ 根据d 内积的形式来定义d 范数 u 怕= 而- 、伊瓦= 藤v u 卯 2 4 若f i d f l 2 = 元，且d 的每个元素都相等，那么此范数等于2 一范数，此方法也等于a r n o l d i 方法。 下面给出加权a r n o l d i 算法： 算法6 ：加权a r n o l d i 方法 1 计算v l = 口川口怕 o ，、 = 珂 贡性下以足满阵矩此 叁旦盔堂亟坐些迨塞 2 5 2 对于j = 1 ，2 ，直到m f 2 1 ) w = a v j ( 2 2 )i = 1 ，j ( 2 2 1 ) h i j = ( w ，仇) d ( 2 2 2 ) 叫= w h i d v i ( 2 3 ) h j + 1 ，j = 怯，i ，h j + l , j = 0 ，s t o p ( 2 a ) 帅2 老 结束算法 算法6 构造的基= v l ， 。】是d 正交的，因此 v g d v , 。= k( 3 4 ) 由加权a r n o l d i 算法构造的h e s s e n b e r g 矩阵士k 可以表示成形式 日。= 咯d - 4 ( 3 5 ) 我们定义矩阵耳。冗( ”十1 ) 。”为 玩= ( 忍。：。三) 同样，由加权a r n o l d i 算法构造的矩阵满足 a = + 1 霄。 3 2w g m r e s 方法和p w g m r e s 方法 ( 3 6 ) 如同所有的k r y l o v 方法，加权g m r e s 方法的第m 步迭代值z 。属于x o + k 。( a ， g o ) 。 我们在空间z o + 甄。( a ，r o ) 中选择x 。，使余向量的d 范数最小。实际上，它就是加权最 复旦大学硕士毕业论文 小二乘问题 。+ m k i 。n ( ，。) 1 f 6 一a 。1 l 。 ( 3 _ 7 ) 的解。 迭代值z 。可写成形式 z 。= x o + ( 3 8 ) 这里y m r ”。因此，相应的余向量r 。= b a x 。满足 r 。= b a ( x o + k 。) = t o a v m y 。= v m + 1 ( 卢e 1 一日m y 。) 这里卢= 0 怯，e 。是单位矩阵的第一列。 因为矩阵v 幺+ 1 是d 正交的，所以 l l r 。l l 刍= 8 1 + ，( 卢e ，一再3 。) ij 刍= l i 卢e - 一j 再。u m l l ； 那么问题( 3 7 ) 就可以转变为以下问题：找以下最小二乘问题的解向量 。r a 即i n 慨1 一玩圳2 ( 3 ，9 ) 以上问题的解可通过对矩阵豆。作q r 分解来得到。 加权g m r e s 方法要求存储v k ，因此m 的值受到内存的限制。为了解决w g m r e s 方法中遇到的问题，算法在m 步迭代后重新开始，我们用w g m r e s ( m ) 来表示此算法 算法7 ：w g m r e s ( m ) 算法 1 选择嚣o ，m 和扰动，计算r d = b a x o 。 2 选择对角矩阵d ，计算芦= l l r o 怕，u = r o 卢。 3 用算法6 来构造d 正交基，初始向量为”。 4 i 十算g m = g 则m i n 。慨1 一瓦圳2 ，z m 2 z 。+ ，r m = 6 4 z m 。 5 若1 眺 e ，停止否则令跏= z m ，t o = r m ，回到步骤2 。 当h j + 1 1 ，= 0 时，加权a r n o l d i 方法就在j 步停止了，也就是说此w g m r e s 算法也停 止了，那么近似解为巧= x o + k 巧1 e 1 ) 。 下面我们来看一下算法7 的性质。 复旦大学硕士毕业论文 定理3 2 1 若用e d 来代替口，口是正实数，那么算法7 的结果不改变。 证明：令d = a d ，从算法6 得到 和 我们得到 并且 1 1 地2 蕊地舢2 而”、 q 五巧= ( 西，妨) 西= 三( 叫，码) 。= ( w ，v j ) 。= 。 死= 击，_ m = 瓦，= 卢= l i r o 怯= 托怖怕= 施卢 在w g m r e s 方法中：= 。o + 孰，这里 = g p r a 。i n 慨1 一日m 北 因为声= 以卢和矗。= 一h 。，所以 因此 基 如= 狐伽躲一瓦g i i 2 = 瓶妨 ( 3 1 0 ) ( 3 1 1 ) = 时( 去) ( 俩m ) = x o - i - = 盘m 证毕口 令m 为a 的( 左) 预条件矩阵，那么我们可构造预条件k r y l o v 子空间的一组d 正交 k m ( m a ，”) = s p a n ，m - 。a v ，( m a ) 1 ”) d 正交基矩阵依旧用三【u 1 ，】来表示，( 3 4 ) 成立，并且 m 一1 a = + l 耳。 基于预条件加权a r n o l d i 方法，我们定义预条件加权g m r e s ( p w g m r e s ( m ) ) 方法 复旦大学硕士毕业论文 算法8 ：p w g m r e s ( m ) 算法 1 选择z o ，m 和扰动e ，计算7 0 = b a z o ，5 0 = m _ 1 r o 。 2 选择对角矩阵d ，计算卢= 1 5 0 i l y ，”= 5 0 8 。 3 用算法6 来构造d 正交基，开始向量为u 1 。 4 坩算2 盯g m i n l i