第３２卷第１期 数学理论与应用 ＭＡＴＨＥＭＡ ｎＣＡＬ ＴＨＥＯＲＹ ＡＮＤ ＡＰＰＬＩＣＡＴＩＯＮＳ Ｖ０ｌ＿３２ Ｎｏ．１ Ｍａｒ．２０１２ 、 ２０１２年３月 块Ｔｏｅｐｌｉｔｚ方程组的快速块Ｇａｕｓｓ—Ｓｅｉｄｅｌ迭代算法 冯月华刘成志刘仲云 （长沙理工大学数学与计算科学学院，长沙，４１０００４） 摘要本文研究块Ｔｏｅｐｌｉｔｚ方程组的块Ｇａｕｓｓ—Ｓｃｉｄｅｌ迭代算法．我们首先讨论了块三角Ｔｏｅｐｌｉｔｚ矩阵的一 些性质，然后给出了求解块三角Ｔｏｏｐｌｉｔｚ矩阵逆的快速算法，由此而得到了求解块Ｔｏｅｐｌｉｔｚ方程组的快速块 Ｇａｕｓｓ—Ｓｅｉｄｅｌ迭代算法，最后证明了当系数矩阵为对称正定和Ｈ一矩阵时该方法都收敛．数值例子验证了方 法的收敛性． 关键词块Ｔｏｅｐｌｉｔｚ 块Ｇａｕｓｓ—Ｓｅｉｄｅｌ迭代 快速算法对称正定Ｈ一阵 Ｆａｓｔ Ｂｌｏｃｋ Ｇａｕｓｓ——Ｓｅｉｄｅｌ Ｉｔｅｒａｔｉｏｎｓ ｆｏｒ Ｔｏｅｐｌｉｔｚ Ｓｙｓｔｅｍｓ Ｆｅｎｇ Ｙｕｅｈｕａ Ｌｉｕ Ｃｈｅｎｇｚｈｉ Ｌｉｕ Ｚｈｏｎｇｙｕｎ （Ｓｃｈｏｏｌ ｏｆ Ｍａｔｈｅｍａｔｉｃｓ ａｎｄ Ｃｏｍｐｕｔｉｎｇ Ｓｃｉｅｎｃｅ，Ｃｈａｎｇｓｈａ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｓｃｉｅｎｃｅ ａｎｄ Ｔｅｃｈｎｏｌｏｇｙ， Ｃｈａｎｇｓｈａ，４１０００４，Ｐ．Ｒ．Ｃｈｉｎａ） Ａｂｓｔｒａｃｔ Ｔｈｅ ｂｌｏｃｋ Ｇａｕｓｓ——Ｓｃｉｄｅｌ ｉｔｅｒａｔｉｏｎｓ ｆｏｒ ｓｏｌｖｉｎｇ ｂｌｏｃｋ Ｔｏｅｐｌｉｔｚ ｓｙｓｔｅｍｓ ｔｉｒｅ ｃｏｎｓｉｄｅｒｅｄ ｉｎ ｔｈｉｓ ｐａｐｅｒ．Ｗｅ ｉｆｓｔｒ ｄｉｓｃｕｓｓ ｓｏｍｅ ｐｒｏｐｅｒｔｉｅｓ ｏｆ ｂｌｏｃｋ ｔｒｉａｎｇｕｌａｒ Ｔｏｅｐｌｉｔｚ ｍａｔｒｉｃｅｓ。ｔｈｅｎ ｐｒｅｓｅｎｔ ｆａｓｔ ａｌｇｏｒｉｔｈｍｓ ｆｏｒ ｉｆｎｄｉｎｇ ｔｈｅ ｉｎｖｅｒｓｅｓ ｏｆ ｓｕｃｈ ｍａｔｒｉｃｅｓ．ａｎｄ ｆｕｎＩｌｅｒ ｏｂｔａｉｎ ｆａｓｔ ｂｌｏｃｋ Ｇａｕｓｓ—Ｓｅｉｄｅｌ ｉｔｅｒａｔｉｖｅ ａｌｇｏｒｉｔｈｍｓ ｆｏｒ ｂｌｏｃｋ Ｔｏｏｐｌｉｔｚ ｓｙｓｔｅｍｓ．Ｆｉｎａｌｌｙ， ｗｅ ｓｈｏｗｔｈａｔ ｏｕｒｍｅｔｈｏｄｓ ａｒｅ ｃｏｎｖｅｒｇｅｎｔｗｈｅｎｔｈｅ ｃｏｅｆｆｉｃｉｅｎｔｍａｔｒｉｃｅｓ ８ｒｅ ｓｙｍｍｅｔｒｉｃ ｐｏｓｉｔｉｖｅ ｄｅｆｉｎｉｔｅ ｏｒＨ—ｍａｔｒｉｃｅｓ． Ｓｏｍｅ ｎｕｍｅｒｉｃａｌ ｅｘａｍｐｌｅｓ ｄｅｍｏｎｓｔｒａｔｅ ｔｈｅ ｃｏｎｖｅｒｇｅｎｃｅ ｏｆ ｏｕｒ ｓｃｈｅｍｅｓ． Ｋｅｙ ｗｏｒｄｓ Ｂｌｏｃｋ Ｔｏｅｐｌｉｔｚ Ｂｌｃｋ Ｇａｕｓｓ—ｏ－Ｓｅｉｄｅｌ ｉｔｅｒａｔｉｏｎ Ｆａｓｔ ａｌｇｏｒｉｈｍ Ｓｙｍｍｅｔｒｉｔｃ ｐｏｓｉｔｉｖｅ ｄｅｆｉｎｉｔｅ Ｈ—‘ ｍａｔｒｉｃｅｓ １ 引言 本文考虑方程组 Ａ ． ：ｂ （１） 的迭代解法，其中Ａ…为ｍ 阶的块Ｔｏｅｐｌｉｔｚ矩阵． 块Ｔｏｅｐｌｉｔｚ矩阵具有如下结构 湖南省教育厅重点资助项目（０９Ａ００２［２００９］） 收稿日期：２０１２年２月７日 ２ 数学理论与应用 Ａｏ Ａ—ｌ Ａ川 ｌ 以０ ●●● ＋２ Ａ ＝ ●●● （２） ●●－ ０ 其中Ａ （１一 ≤ｉ≤ 一１）为ｍ ｘ ｍ的任意矩阵，即 ． 的每一条对角线上都是同一个矩阵； 若Ａ— ＝Ａ （１一ｎ≤ｉ≤ｎ一１），则称矩阵Ｃ 为块循环矩阵，记为Ｃ…＝Ａ…＝Ｂ—Ｃｉｒｃ（Ａ。， Ａ 一，Ａ ）． Ｔｏｅｐｌｉｔｚ矩阵在信号处理等工程问题中有着广泛的应用，其理论与结构算法被广为研 究¨ ．对于一个Ⅳ阶的Ｔｏｅｐｌｉｔｚ线性方程组，其解法主要分为直接法和迭代法 ，ｍ］，虽然 直接法的计算复杂度为０（Ｎ２），但对许多Ｔｏｅｐｌｉｔｚ方程组，直接解法是非常不稳定的，因此，迭 代法就成了另一种选择．对于迭代法，一个好的预条件子能使其复杂度降为０（ＮｌｏｇＮ），由此 引出了Ｔｏｅｐｌｉｔｚ及相关线性方程组预条件子的广泛研究． 本文将讨论块Ｔｏｅｐｌｉｔｚ方程组的Ｇａｕｓｓ—Ｓｅｉｄｅｌ型迭代法，并重点考虑块Ｔｏｅｐｌｉｔｚ矩阵 … 为对称正定矩阵和Ｈ一矩阵这两种情形．下节主要介绍一些将要用到的预备知识，第三节给 出块Ｇａｕｓｓ—Ｓｅｉｄｅｌ迭代的快速算法，最后一节给出数值算例验证了算法的有效性． ２预备知识 本文在用块Ｇａｕｓｓ—Ｓｅｉｄｅｌ迭代法解方程组时，需要求块下三角Ｔｏｅｐｌｉｔｚ矩阵的逆，下面 是 些显而易见的有关块下三角Ｔｏｅｐｌｉｔｚ的结论． 引理１记ｒ 为ｎ×ｎ的块下三角Ｔｏｅｐｌｉｔｚ矩阵的集合，则有： １）Ａ，Ｂ∈＂ｉ＂ｍｎ￣ＡＢ∈ ２）Ａ，Ｂ∈ｒ ； ＡＢ＝ＢＡ； ３）Ａ∈ｒ 且Ａ０可逆 Ａ～∈７ ； ４）Ａ Ｅ ｒｍ ，Ｂ∈ ｍ ，Ａ ＝Ｂ ，ｋ＜ｆ Ａ ’＝ －１’，ｋ＜１． 由引理１易得下面的结果： 推论２设Ａ ＝（ １］，Ａ 。∈ ， ，ｃ是块Ｔ。ｅｐｈｔｚ矩阵，则 （ ］， ：一Ａ ｃＡ Ａ： ＝ｃ３ 块Ｔｏｅｐｌｉｔｚ方程组的快速块Ｇａｕｓｓ—Ｓｅｉｄｅｌ迭代算法 ３ 定义１一个矩阵Ｂ：（ｂ ）称为：（ｉ）非负矩阵，如果所有ｂ ≥０；（ｉｉ）ｚ一矩阵，如果对 于所有ｉ≠ ，ｂ ≤０；（ｉｉｉ）Ｍ一矩阵，如果 是可逆的Ｚ一矩阵且Ｂ ≥Ｏ；（ｉｖ）Ｈ一矩阵，如 果它的比较矩阵＜Ｂ＞是Ｍ一矩阵，这里＜Ｂ＞＝（ ｉ ），其中 ：Ｉ ｂ ｆ， ｉｊ＝－Ｉ ｂｉ ， ≠ ｊ． 众所周知，一个块循环矩阵ｃ ， 可块对角化，即：Ｃ ＝， ） Ｆ（ ）其中Ｆ（ ）＝Ｆ ｔｍ，Ｆ 是ｎ阶ＤＦＴ矩阵； 表示Ｋｒｏｎｅｃｋｅｒ积；Ａ表示一个ｍｎ×ｍｎ块对角矩阵Ａ＝ Ｂｌｋｄｉａｇ（Ｆ㈨Ｃ（：，１：ｍ）），这里ｃ（：，１：ｍ）指的是矩阵Ｃ的前ｍ列，Ｂｌｋｄｉａｇ表示由括号里的 块向量生成的块对角阵．因此一个块循环矩阵与一个向量相乘可以使用Ｂ—ＦＦＴ技术，使用 Ｂ—ＦＦｒ技术不仅计算快，而且便于并行计算． ３块Ｇａｕｓｓ—Ｓｅｉｄｅｌ迭代 求解方程组Ａｘ＝ｂ的迭代格式一般可表示为 Ｍｘ川＝Ｎｘ ＋ｂ ｋ＝０，１，…， （４） 其中Ａ＝Ｍ一Ⅳ称为矩阵Ａ的分裂，如果 非奇异．众所周知，对于任给的初始向量 ，由 （４）生成的迭代序列｛ ｝收敛的充要条件是迭代矩阵Ｈ：Ｍ Ｎ的谱半径ｐ（日）＜１． 关于迭代法，下面的定义和结果是熟知的． 定义２分裂Ａ＝Ｍ—Ｎ称为矩阵Ａ的：弱正则分裂，如果Ｍ ≥０且Ｍ Ｊ】ｖ≥０；Ｐ一正 则分裂，如果 ＋Ⅳ正定． 引理３如果 对称正定，且Ａ＝Ｍ—ＮＰ一正则分裂，则ｐ（Ｍ Ｎ）＜ｌ；如果Ａ为Ｍ一矩 阵，且Ａ：Ｍ—Ｎ弱正则分裂，则ｐ（Ｍ Ｎ）＜１． 令Ａ＝Ｄ一￡一 ，其中Ｄ， 和 分别表示由 …块对角，块下三角和块上三角部分所形 成的矩阵．如果在方程（４）中取Ｍ＝Ｄ—Ｌ，Ｎ＝Ｕ，那么就得到了如下求解Ａｍ，＝ｘ＝ｂ的块 Ｇａｕｓｓ—Ｓｅｉｄｅｌ迭代格式： （Ｄ—Ｌ）ｘ ＝ＵｘＩ＋ｂ ｋ＝０，ｌ，…， （５） 从方程（５），我们知道，（Ｄ一￡） 和 都是块Ｔｏｅｐｌｉｔｚ矩阵，它们与向量的乘积可先将它们嵌 入到块循环矩阵，再与向量作乘积，由此而得到快速算法．因此，其核心是如何计算Ｄ一￡的 逆．下面我们利用引理１和推论２并参照文献［２］中关于三角Ｔｏｅｐｌｉｔｚ矩阵逆的快速算法，来 推导块下三角Ｔｏｅｐｌｉｔｚ矩阵逆的快速算法，这里我们主要应用了分合方法（Ｄｉｖｉｄｅｄ ａｎｄ Ｃｏｎ． ｑｕｅｒ）的思想． 从公式（３）知，当知道了Ａ 时，则只需要计算Ｂ 由引理ｌ的３），Ａ２ 是一个块下三角 ４ 数学理论与应用 Ｔｏｅｐｌｉｔｚ矩阵，因此，我们只需求出Ｂ 的前ｍ列，它的前ｍ列等于Ａ 与Ａ －１ ｃ积的前ｍ列， 而Ａ二 与ｃ都是块Ｔｏｅｐｌｉｔｚ矩阵，能利用Ｂ—ＦＦｒ来计算．不妨设ｎ＝２ ， 胁＝２ ，（０≤ｚ＜ ｑ）． 算法该算法可快速求块下三角Ｔｏｅｐｌｉｔｚ矩阵逆． １）ｆｕｎｃｔｉｏｎＡ１＝ｂｔｒｉｎｖ（Ａ，ｎ， ） ２）ｉｆ ｎ≤ｎ ３）利用公式（３）计算 三 的前 列，再由块Ｔｏｅｐｌｉｔｚ结构，得到Ａ二 ． ４ ｅｌｓｅ ５，将Ａ分块成Ａ＝Ａ １ ）； ６）Ａ ：ｂｔｒｉｎｖ（Ａｌｌ，ｎ／２，ｎ ）； ７）用两次Ｂ—ＦＦＴ技术，先计算Ａ： Ａ 的前ｍ列，再计算Ａ 乘以这ｍ列； ８）最后由块Ｔｏｅｐｌｉｔｚ结构，得出４～． ９）ｅｎｄ ４收敛性和数值实验 在这一节，我们先分析（５）的收敛性，然后给出实例，来检验块Ｇａｕｓｓ—Ｓｅｉｄｅｌ迭代法的收 敛性． 定理４如果（１）中Ａ…为对称正定或Ｈ一矩阵，则对于任意初始向量 。，由（５）生成的 迭代序列｛ ｝收敛． 证明 由引理３，我们只需要验证块Ｇａｕｓｓ—Ｓｅｉｄｅ１分裂满足收敛性条件即可． 当Ａ…对称正定时，容易知道由它的对角块组成的块对角矩阵Ｄ也对称正定，而（Ｄ一 ￡） ＋Ｕ＝Ｄ意味Ａ…＝（Ｄ—Ｌ）一 为Ｐ一正则分裂，收敛条件成立． 当Ａ…为Ｈ一矩阵时，＜Ａ ＞为Ｍ一矩阵，＜Ｄ—Ｌ＞为Ｍ一矩阵，ｌ Ｉ为非负矩阵， 于是有＜Ｄ—Ｌ＞ ≥０，且＜Ｄ—Ｌ＞ Ｉ Ｉ≥０，即＜Ａ…＞＝＜Ｄ—Ｌ＞－Ｉ Ｕ Ｉ为弱正 则分裂，于是又ｐ（＜Ｄ一 ＞ ｆ ｆ）＜１． 利用谱半径不等式ｐ（／３）≤ｐ（１日１）和ｐ（ ≤Ｐ（Ｃ），如果Ｃ≥Ｂ≥０以及矩阵不等式 Ｉ Ａ８ ｌ≤Ｉ Ａ Ｉｌ Ｉ，Ｉ（Ｄ— ）Ｉ－１≤＜Ｄ—Ｌ＞＿。，我ｆ门有ｐ（（Ｄ—Ｌ）一Ｕ）≤Ｐ（Ｉ（Ｄ一 ）．１ Ｉ）≤Ｐ（１（Ｄ一￡）ｌ Ｉ Ｕ Ｉ）≤ｐ（＜Ｄ—Ｌ＞ Ｉ ｕ Ｉ），因此有Ｐ（（Ｄ—Ｌ） Ｕ）≤１，定理得 证． 块Ｔｏｅｐｌｉｔｚ方程组的快速块Ｇａｕｓｓ—Ｓｅｉｄｅｌ迭代算法 ５ 在我们的实验中，块Ｔｏｅｐｌｉｔｚ矩阵Ａ中元素是随机生成的，当Ａ是对称正定时，它对角线上 的元素取值在［２０，３０］，其它元素在［０，１］；当Ａ是日一矩阵时，它的对角线上的元素取值在 ［５１ｏｏ，ｌｏ２ｏｏ］，下三角（去掉对角线）的元素取值在［一ｌ，一７］，上三角（去掉对角线）的元素 取值在［一１，０），ｂ是随机生成的向量，初始点 。是零向量，终止条件是ｌＩ 川一 Ｉ ｌ≤１０～， 指的是第ｋ步迭代的近似解，由Ｍａｔｌａｂ完成计算，块Ｇａｕｓｓ—Ｓｅｉｄｅｌ的结果如表１． 表１块Ｇａｕｓｓ—Ｓｅｉｄｅｌ迭代结果 参考文献 ［１］Ｒ．Ｈ．Ｃｈａｎ ａｎｄ Ｍ．Ｋ．Ｎｇ．Ｃｏｎｊｕｇａｔｅ ｇｒａｄｉｅｎｔ ｍｅｔｈｏｄｓ ｆｏｒ Ｔｏｅｐｌｉｔｚ ｓｙｓｔｅｍｓ，ＳＩＡＭ Ｒｅｖ．，３８：４２７—４８２， １９９６． ［２］Ｄ．Ｃｏｍｍｇｅｓ，Ｍ．Ｍｏｎｓｉｏｎ．Ｆａｓｔ Ｉｎｖｅｒｓｉｏｎ ｏｆ Ｔｒｉａｎｇｕｌａｒ Ｔｏｅｐｌｉｔｚ Ｍａｔｉｒｃｅｓ．１ＥＥＥ Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ ａｕｔｏｍａｔｉｃ ｃｏｎ— ｔｒｏｌ，１９８４，２９（３）２５０—２５１． ［３］Ａ．Ｆｒｏｍｍｅｒ，Ｄ．Ｂ．Ｓｚｙｌｄ．Ｈ—Ｓｐｌｉｔｔｉｎｇｓ ａｎｄ ｔｗｏ—ｓｔａｇｅ ｉｔｅｒａｔｉｖｅ ｍｅｔｈｏｄｓ．Ｎｕｍｅｒｉｓｃｈｅ Ｍａｔｈｅｍａｔｉｃ，１９９２， ６３：３４５—３５６． ［４］Ｇ．Ｈ．Ｇｏｌｕｂ，ｃ．Ｆ．Ｖａｎ Ｌｏａｎ．Ｍａｔｉｒｘ Ｃｏｍｐｕｔａｔｉｏｎ．Ｔｈｉｒｄ Ｅｄｉｔｉｏｎ，Ｔｈｅ Ｊｏｈｎｓ Ｈｏｐｋｉｎｓ Ｕｎｉｖｅｓｒｉｔｙ Ｐｒｅｓｓ，１９９６． ［５］Ｚｈｏｎｇ—Ｙｕｎ Ｌｉｕ，Ｈｅ—Ｂｉｎｇ Ｗｕ，Ｌｕ Ｌｉｎ．Ｔｈｅ ｔｗｏ—ｓｔａｇｅ ｉｔｅｒａｔｉｖｅ ｍｅｔｈｏｄｓ ｆｏｒ ｓｙｍｍｅｔｉｒｃ ｐｏｓｉｔｉｖｅ ｄｅｆｉｎｉｔｅ ｍａ。 ｔｉｒｃｅｓ．Ａｐｐｌｉｅｄ Ｍａｔｈｅｍａｔｉｃｓ ａｎｄ Ｃｏｍｐｕｔａｔｉｏｎ，１１４（２０００）１—１２． ［６］Ｍ．Ｋ．Ｎｇ．１ｔｅｒａｔｉｖｅ Ｍｅｔｈｏｄｓ ｆｏｒ Ｔｏｅｐｌｉｔｚ Ｓｙｓｔｅｍｓ，Ｏｘｆｏｒｄ Ｕｎｉｖｅｒｓｉｔｙ Ｐｒｅｓｓ，Ｏｘｆｏｒｄ，ＵＫ，２００４．