（计算数学专业论文）一种基于稳定的不完全分解技术预处理非对称鞍点问题.pdf

摘要 摘要 鞍点问题广泛存在于科学与工程计算中，解决此类问题要用 到预处理，目前主要有不精确块对角与块上三角阵预处理子， 构造不精确三角预处理子要用到不完全分解技术，i l u t 是其中 一种，我们给出稳定的不完全q z 分解方法事实上，这种方法 是r i f 算法的一种变换形式，不同之处是r i f 算法主要用来求解正 规方程，而我们的算法侧重于对a 进行分解，并利用这种分解构 造a 的近似我们将其与i l u t 分解法进行比较，预处理后的鞍点 矩阵的谱性将会详细介绍，随后的数值实验将显示出此分解的优 越性 关键词：鞍点问题；不完全q z 分解；，预处理子；扰动 界；o s e e n 方程 a b s t r a c t a b s t r a c t s a d d l ep o i n tp r o b l e ma r i s e si nm a n ys c i e n t i f i ca n de n g i n e e r i n g a p p l i c a t i o n s w ep r e s e n tak i n do f i n e x a c tb l o c k d i a g o n a la n db l o c kt r i a n g u l a rp r e c o n d i t i o n e r sb a s e d o nr o b u s ti n c o m p l e t eq zf a c t o r i z a t i o nf o r s a d d l ep o i n tp r o b l e m s i nf a c t ，t h i sa l g o r i t h mi sav a r i a t i o no fr i fa l g o r i t h m t h e nw ec o m p a r es u c hm e t h o d st ot h ea p p r o x i m a t i v ei n v e r s eo f ab a s e do ni l u tf a c t o r i z a t i o n t h es p e c t r a lp r o p e r t yo ft h ep r e c o n d i t i o n e dm a t r i xi ss t u d i e di nd e t a i l n u m e r i c a le x a m p l e sa r eu s e dt o i l l u s t r a t et h ee f f i c i e n c yo fs u c hp r e c o n d i t i o n e r s k e y w o r d s ： s a d d l ep o i n tp r o b l e m ；i n c o m p l e t eq zf a c t o r i z a t i o n ；p r e c o n d i t i o n e r ；p e r t u r b a t i o nb o u n d ；o s e e ne q u a t i o n 一l v 前言 - 土上_ - j 一 刖吾 在流体力学、带有限制条件的二次优化、电磁学、线性弹力学等应用领 域中经常需要求解一类问题例如，在求解s t o k e s 方程 0 ，1 j n 若我们令q = a r ，显然，q 是列正交的( 其2 一n o r m 未必为1 ) ，通过( 2 1 3 ) 我们 有这样一个事实，a = q r ，令r 一1 兰z ，得到a = q z 2 2 稳定的不完全q z 分解 现在，我们观察发现共轭格拉姆施密特过程( 2 1 2 ) 不仅得到上三角矩 阵r ，并且同时得n q 因子，考虑到预处理的需要，我们不必要对a 进行完 全的分解，只需对其不完全分解就可以其分解过程如下：假设一个舍弃 限度0 7 r i2 r i e t 甲j 舍弃n 比7 - 小的元 e n d i f e n d f o r e n d f o r 取q = q l ，q 2 ，】，r = r l ，r 2 ，r n 】以及d = d i a g ( d l ，d 2 ，如) ，或者更 进一步可得，让国= q d 一 且a = r d 一，则可以得到 a 国a 一1( 2 2 1 ) 其中国t 国i ，詹是上三角矩阵且具有正的对角元 我们算法( 2 2 ) 类似于文献中的r i f 预处理子，它们都是基于a t a 一正交化 过程不同之处在于r i f 方法隐含着当a 满秩时对4 t a 进行c h o l e s k y 分解，而 我们的方法是当a 是非奇异时，构造出现它的“正交”分解 一9 一 筇2 章稳定的不完全q z 分解 注意虽然出现t a t a 的形式，但是并不需要将a t a 直接计算出来，不完 全格拉姆施密特过程只对a 进行操作更重要的是，在不完全分解过程中避 免了停滞这是由于主元嘞= ( a r j ) t a 乃，是正的数字，并且由r e y l e i g h 商的 极值性可以得到主元而的下界： 奶= ( a 吩) t a 巧入m 讥( a t a ) l l r jl l l = 2 饥( a ) l l r ol l ； 其中入r a i n ( a t a ) 指a t a 的最小特征值，o m i n ( a ) 指a 的最小奇异值由于 选择的a 是非奇异的，所以仃磊饥( a ) 0 ，再由正交分解过程得知勺的第j 个元 值为1 ，撇u i i r j l l ；1 由于这种分解不会出现停滞现象，我们称之为稳定的 不完全分解( i q z ) 下面有几点需要说明： 1 以上分解称为i q z 分解，最i ：i a = q r 一，称为残量矩阵，其稀疏 性由7 - 控制但算法过程中不需要计算r 为了得到a - 1 的近似，而a _ 1 r q 一1 = r d 一1 矽a t 2 d 一1 易得 第3 章预处理后的谱分析 第3 章预处理后的谱分析 3 1 预处理子构造 对于广义鞍点问题，预处理子的选取有很多类，m b e n z i ，g h g o l u b ，j l i e s e n 6 】比较全面的介绍广义鞍点问题的解的方法以及预处理技术，下面简 要预处理的一些方法 当c = 0 时，m f m u r p h y ，g h g o l u b ，和a j w a t h e n 【2 3 给出了两种预处 理子 ，、，、 r = ( a 。b a 0 1 b t ) 伤= ( a 。b 篇t ) ， 假设4 和b a - 1 b r 都是可逆的，可知预处理矩阵阿1 至多有三个非零特 征值1 ，琏乎【2 3 ，r e m a r k1 】，由k r y l o v 子空间法的有限终止性，至多3 步收 敛1 2 3 ，r e m a r k3 】。而巧1 有两个特征值1 和1 2 3 ，r e m a r k4 】，它也指出可 对伤的( 2 ，2 ) 块乘上1 后的预处理阵特征值为1 ，但是这种预处理后的系 统不再是可对角化的了 i p s e n 1 9 将m u r p h y 等的工作推广到e 0 情况，类比( 3 1 1 ) ，i p s e n 给出 了预处理子 乃= a c + 肿0 。b t ) 砍= ac + 肿b t ，矿) 地， e d es t u f f e r 和j l i e s e n 推广t m u r p h y 2 3 ( 但不是i p s e n 的) 等的预处理子， 他们利用( 1 ，1 ) 块a 的近似逆，考虑这样一个分裂形式 a = m 一 其中m 是可逆的，那么可得块对角预处理子 ( 3 1 3 ) 磊= m b 胪0b t ) 川 第3 章预处理后的谱分析 假设b m 一1 b t 是可逆的注意b m 一1 b t 是矩阵 ( 警髻) 的s c h u r 4 b ，此矩阵称为约束预处理子预处理阵 ( 3 1 5 ) 咖= ( 。赢踯bm b m ， 设它的特征值为a 巧此处特征值扰动界为 k 刈 0 代表粘度系数对此方程的许多离散都会产 生形如( 0 0 4 ) 的鞍点问题其中( 1 ，1 ) 块a 是由离散c o n v e c t i o n d i f f u s i o n 项得到的一 个非对称正实矩阵；我们利用h o w a r de l m a n ，a l i s o nr a m a g e 和d a v i ds i l v e s t e r 开发 的软件包i f i s s 【2 9 来离散( 4 0 1 ) o 的c h a n n e ld o m a i n 问题，选取q 1 一局有限元进 行离散，离散网格为正方形网格1 6 1 6 稳态参数为0 2 5 ，此实验中我们分别选 取z ：0 0 1 和= 1 那么我们可以得到形如( 0 0 4 ) n = 5 7 8 且m = 2 5 6 的系数矩阵矩 阵五和雪分别通过对a 和进行i u j t 分解( 采用m a t l a b 系统自带的l u i n c m ) ，在我们 的实验中采取i q z 分解得到，由于代码未经优化，c p u 运行时间仅供参考随着7 - 的不 同选取，p 一1 4 的特征值分布通过图4 i - 4 6 显示出来，此实验也显示出当丁一0 + 时， 预处理后的矩阵的特征值集中在两簇( 1 ，0 ) 与( 1 ，0 ) 此后我们给出预处理后采取g m r e s 迭代的次数及其收敛性，这些实验也说明了 特征值分布情况影响着收敛速度 在我们的计算中，右端项选取为r a n d ( 8 3 4 ，1 ) ，初始景x 0 = 0 ，迭代终止条件 为| i b 一a , 1 1 2 1 0 一i l b l l 2 ，i t e r 是迭代次数 一1 8 一 第4 章数值实验 首先，我们选取粘度系数y = 1 时，不精确块对角阵预处理后的矩阵特征值分布概 况( 横轴为实轴，纵轴为虚轴) ： 图4 1m - t - n u t ( o ) i ，分解的不精确块对角角馕处理基于q 可0 0 1 ) 分解的不精确块对角预处曩 图4 2 基于ni r r ( m 1 ) 分解的不精确块对角预处理基于q z 。o 1 ) 分解的不精确块对角援处理 图4 3 基于i l u t ( i c _ 4 ) 分解的不精确块对角援处理基于q z 【l e _ 4 扮解的不精确块对角璜处理 一1 9 一 第4 章数值实验 当采取不精确块三角阵预处理时，两种不完全分解后预处理阵的特征值分布如 图4 4j g t t l u t ( o 0 1 ) 分解的不精确块三角预处理基于q z ( 0 脚) 分解的不精确块三角预处理 图4 5 基于ni m o 1 ) 分解的不精确块三角预处理t - t - q z ( o 0 0 1 ) 分解的不精确块三角预处理 图4 6 基于i r r ( 1 e - 4 ，分解的不精确块二角磺处理 基于q z ( 1 e 4 ，分解的不精确块二角预处理 一2 0 一 至于l = o 0 1 时的情况，按上述方法可以同样获得下列表格给出两个分解法的 数据比较，其中声指块对角预处理阵，p 指块上三角预处理阵 t a b l e 比较发现采用块对角预处理时，选择i l u t 分解优于不完全q z 法。至于块上三 角预处理，当舍弃门槛丁变得越小时，i q z 方法要优于1 l u t 一2 l 一 参考文献 参考文献 【1 】z - z b a i ，m k n g ，o ni n e x a c tp r e c o n d i t o n e r sf o rn o n s y m m e t r i cm a t r i c e s ， s i a mj s c i c o m p u t 2 6 ( 2 0 0 5 ) 1 7 1 0 1 7 2 4 【2 z 一z b a i ，s t r u c t u r ep r e c o n d i t i o n e r sf o rn o n s i n g u l a rm a t r i c e so fb l o c kt w o - b y - t w os t r u c t u r e s ，m a t h c o m p 7 5 ( 2 0 0 6 ) 7 9 1 8 1 5 【3 】z z b a i ，i s d u f f , a n da j w a t h e n ，ac l a s so fi n c o m p l e t eo r t h o g o n a lf a c - t o r i z a t i o nm e t h o d s i ：m e t h o d sa n dt h e o r i e s ，b i t 4 1 ( 2 0 0 1 ) 5 3 7 0 【4 】z z b a i ，g h g o l u b ，l - z l u ，j - ey i n ，b l o c kt r i a n g u l a ra n ds k e w - h e r m i t i a ns p l i t t i n gm e t h o d sf o rp o s i t i v ed e f i n i t el i n e a rs y s t e m s ，s i a mj s c i c o m p u t ( 2 0 0 4 ) 【5 】m b e n z i ，g h g o l u b ，ap r e c o n d i t i o n e rf o rg e n e r a l i z e ds a d d l ep o i n tp r o b - l e m s s i a mj m a t r i xa n a l a p p l 2 6 ( 2 0 0 4 ) 2 0 - 41 【6 】m b e n z i ，g h g o l u b ，j l i e s e n ，n u m e r i c a ls o l u t i o no fs a d d l ep o i n tp r o b l e m s ，a c t an u m e r 1 4 ( 2 0 0 5 ) 1 1 3 7 【7 】m 。b e n z i ，m t u m a ，ar o b u s ti n c o m p l e t ef a c t o r i z a t i o np r e c o n d i t i o n e rf o rp o s - i t i v ed e f i n i t em a t r i c e s ，n u m e r l i n e a ra l g e b r aa p p l ，1 0 ( 2 0 0 3 ) 3 8 5 - 4 0 0 【8 】m b e n z i ，m t u m a ar o b u s tp r e c o n d i t i o n e rw i t hl o wm e m o r yr e q u i r e m e n t s f o rl a r g es p a r s el e a s ts q u a r e sp r o b l e m s s i a mj s c i c o m p u t 2 5 ( 2 0 0 3 ) 4 9 9 - 5 1 2 【9 】j h b r a m b l e ，j e p a s c i a k ，a t v a s s i l e v , u z a w at y p ea l g o r i t h m sf o rn o n s y m m e t r i cs a d d l ep o i n tp r o b l e m s ，m a t h c o m p 6 9 ( 2 0 0 0 ) 6 6 7 - 6 8 9 【1o 】eb r e z z i ，m f o r t i n ，m i x e da n dh y b r i df i n i t ee l e m e n tm e t h o d s ，s p r i n g e r - v e r l a g ，n e wy o r k ，1 9 9 1 【11 】z - h c a o ，c o n s t r a i n ts c h u rc o m p l e m e n tp r e c o n d i t i o n e r sf o rn o n s y m m e t r i c s a d d l ep o i n tp r o b l e m s ，a p p l n u m e r m a t h 5 9 ( 2 0 0 9 ) 151 16 9 【12 】k e c h e n ，m a t r i xp r e c o n d i t i o n i n gt e c h n i q u e sa n da p p l i c a t i o n s ，c a m b r i d g e u n i v e r s i t yp r e s s ，c a m b r i d g e ，2 0 0 5 【13 】h c e l m a n ，d j s i l v e s t e r , a j w a t h e n ，p e r f o r m a n c ea n da n a l y s i so fs a d d l e p o i n tp r e c o n d i t i o n e r sf o rt h ed i s c r e t es t e a d y s t a t en a v i e r - s t o k e se q u a t i o n s ， n u m e r m a t h 9 0 ( 2 0 0 2 ) 6 6 5 6 8 8 一2 2 参考文献 1 4 】g h g o l u b ，x w ua n dj - yy u a n ，s o r l i k em e t h o d sf o ra u g m e n t e ds y s - t e m s ，b i t4 1 ，( 2 0 0 1 ) 7 1 - 8 5 【1 5 】g h g o l u ba n dc g r e i f ，o ns o l v i n gb l o c k - s t r u c t u r e di n d e f i n i t el i n e a rs y s t e m s ，s i a mj s c i c o m p u t 2 4 ( 2 0 0 3 ) 2 0 7 6 2 0 9 2 【16 】e h a b e ra n du m a s c h e r ，p r e c o n d i t i o n e da l l a t o n c em e t h o d sf o rl a r g e ， s p a r s ep a r a m e t e r e s t i m a t i o np r o b l e m s ，i n v e r s ep r o b l e m s17 ( 2 0 01 ) 18 4 7 18 6 4 17 】t - z h u a n g ，s - lw u ，c - xl i ，t h es p e c t r a lp r o p e r t i e so ft h eh e r m i t i a na n d s k e w h e r m i t i a ns p l i t t i n gp r e c o n d i t i o n e rf o rg e n e r a l i z e ds a d d l ep o i n tp r o b l e m s ，j c o m p u t a p p l m a t h ( 2 0 0 8 ) 【l8 】k h a y a m i ，j - e y i n ，t i t o ，g m r e sm e t h o d sf o rl e a s ts q u a r e sp r o b l e m s ，n i i t e c h n i c a lr e p o r t ( n i l - 2 0 0 7 0 0 9 e ) ，1 - 2 8 ，2 0 0 7 【19 】i i p s e n ，an o t eo np r e c o n d i t i o n i n gn o n s y m m e t r i em a t r i c e s ，s i a mj s c i c o m p u t ，2 3 ( 2 0 0 1 ) ，p p 1 0 5 0 1 0 5 1 2 0 】y 一ql i n ，y - mw e i ，f a s tc o r r e c t e du z a w am e t h o d sf o rs o l v i n gs y m m e t r i c s a d d l ep o i n tp r o b l e m s 。c a l c o l o ，s p r i n g e r - v e r l a g 4 3 ( 2 0 0 6 ) 6 5 - 8 2 【2l 】y ql i n ，y mw e i ，c o r r e c t e du z a w am e t h o d sf o rs o l v i n gl a r g en o n s y m m e t - r i cs a d d l ep o i n tp r o b l e m s ，a p p l i e dm a t h c o m p u t 1 8 3 ( 2 0 0 6 ) 1 0 1 8 11 2 0 【2 2 】l l i t t l e ，ys a a d ，l s m o c h ，b l o c kl up r e c o n d i t o n e r sf o rs y m m e t r i ca n d n o n s y m m e t r i cs a d d l ep o i n tp r o b l e m s ，s i a mj s c i c o m p u t 2 5 ( 2 0 0 3 ) 7 2 9 7 4 8 【2 3 】m em u r p h y , g h g o l u b ，a n da j w a t h e n ，an o t eo np r e c o n d i t i o n i n gf o r i n d e f i n i t el i n e a rs y s t e m s ，s i a mj s c i c o m p u t 2 1 ( 2 0 0 0 ) 1 9 6 9 1 9 7 2 【2 4 】j 一y p a n ，m k n g ，z z b a i ，n e wp r e c o n d i t i o n e r sf o rs a d d l ep o i n tp r o b l e m s ， a p p l m a t h c o m p u t 1 7 2 ( 2 0 0 6 ) 7 6 2 - 7 7 1 【2 5 】h k p a n g ，l iw e n ，b l o c kt r i a n g u l a rp r e c o n d i t o n e di t e r a t i v em e t h o d s f o rn o n s y m m e t r i cs a d d l ep o i n tp r o b l e m s ，a c t am a t h a p p l s i n i c a ，3 1 ( 2 0 0 8 ) 4 1 9 4 3 1 【2 6 】c s i e f e r t ，e d e s t u r l e r , p r e c o n d i t i o n e r sf o rg e n e r a l i z e ds a d d l ep o i n t p r o b - l e m s ，s i a mj n u m e r a n a l 4 4 ( 2 0 0 6 ) 12 7 5 12 9 6 【2 7 】g w s t e w a r ta n dj g s u n ，m a t r i xp e r t u r b a t i o nt h e o r y , a c a d e m i cp r e s s ， b o s t o n ，19 9 0 【2 8 】vs i m o n c i n i ，b l o c kt r i a n g u l a rp r e c o n d i t i o n e r sf o rs y m m e t r i cs a d d l e p o i n t p r o b l e m s ，a p p l i e dn u m e r m a t h ，4 9 ( 2 0 0 4 ) 6 3 - 8 0 一2 3 参考文献 【2 9 】d j s i l v e s t e r , h c e l m a n ，a n da r a m a g e ，i f i s s ：l n c o m p r e s s i b l e f l o wi t e r a t i v es o l u t i o ns o f t w a r e ， a v a i l a b l eo n l i n e a t ：h t t p ：w 叫w m a n c h e s t e r a c 唑! ! 塑 【3 0 】e d es t u r l e ra n dj l i e s e n ，b l o c k - d i a g o n a la n dc o n s t r a i n tp r e c o n d i t i o n e r sf o r n o n s y m m e t r i ci n d e f i n i t el i n e a rs y s t e m s p a r ti ：t h e o r y , s i a mj s c i c o m p u t 2 6 ( 2 0 0 5 ) 1 5 9 8 - 1 6 1 9 【31 】ys a a d ，m h s c h u l t z ，g m r e s ：ag e n e r a l i z e dm i n i m a l r e s i d u a la l g o r i