细长矩阵的块正交化方法

矩阵＾厶ｘ晰是非奇异的，再增加计算Ｒ＝＾ｆ一１，它就是ＱＲ
分解所得到的上三角矩阵忌一。因为在很多科学计算中
９１
表１ Ｑ矩阵列正交偏离性
Ｎ ２
２＊＊３
２＊＊５
只需要得到列正交矩阵Ｇ。，所以算法中没有计算上三角 矩阵Ｒ。。
细长矩阵块正交化算法中（１）和（７）能从经典格拉姆一 施密特正交化方法中找到相应的计算部分，所以下面仅简 单估计（２）～（６）的计算复杂性。基于Ｇｉｖｅｎ变换ＱＲ分解 来求解线性方程组式（１１）总共需要２（ｋ一１）３＋１ｌ（ｋ一 １）２／２—１３（ｋ一１）／２个浮点操作，算法中（４）的计算需要 ３（ｋ＋１）ｋ／２＋１个浮点操作，（５）需要点个浮点操作，所以增 加的浮点操作个数大约为：
关键词：格拉姆一施密特正交化方法；冗余计算；并行计算
Ｋｅｙ ｗｏｒｄｓ：Ｇｒａｍ－Ｓｃｈｍｉｄｔ ｏｒｔｈｏｇｏｎａｌｉｚａｔｉｏｎ ｐｒｏｃｅｓｓ！ｒｅｄｕｎｄａｎｔ ｃｏｍｐｕｔａｔｉｏｎ；ｐａｒａｌｌｅｌ ｃｏｍｐｕｔｉｎｇ ｄｏｉ：１０．３９６９／ｊ．ｉｓｓｎ．１００７—１３０Ｘ．２０１０．０４．０２５
３块正交化算法
３．１算法推导
从经典格拉姆一施密特正交化过程‘ｓ３可以看出，对任意
恳一１，２，…，ｍ，均有：
Ｉ
ａ¨ ·．Ｉ 一２厶ｙｒｆ‰．田口·¨ ·ｉ
（‘１ｌ’）
ｆ覃１
和
ｑ¨一∑叩“
（２）
其中，口ｉ是标量。 因为生成的矩阵Ｑ，渤是列正交的，对任意忌＝１，２，
…，ｍ，有： （ｑ“，ｑ．．＾）＝（ｇ．．¨ｑ．．ｉ）＝Ｏ，ｆ＝１，２，…，ｋ一１（３） 由线性代数基本知识，对任意五＝１，２’．．·，ｍ，均有；
针对目前的大规模并行计算，研究细长矩阵列向量正
交化的高效并行算法是必要的。原有的格拉姆一施密特正 交化方法的并行算法仅将串行算法进行简单并行，未充分 考虑通讯开销的影响，并行性能受限。由于它们计算时频 繁地计算向量的内积和范数，这就涉及到向量本身的全局 计算，它带来大量的全局通讯和全局同步。随着并行度的 增大，一方面全局通讯和同步的开销呈非线性增长，另一方 面计算向量内积和范数所需要的全局通讯和同步的次数并 未减少，因而全局通讯和同步会成为大规模并行的瓶颈。 并行算法的可扩展性会受到全局通讯和同步的限制。提高 算法的可扩展性就需要尽量减少全局通讯和全局同步的次 数。文献［２～６］提出的正交化方法基本上首先都是计算矩 阵的转置与矩阵的乘积，这样在并行计算时就可将向量的 内积一起计算，避免计算向量内积带来的大量的全局通讯 和同步，但它们在最后将各向垦规范化时计算量较大而且 需要一次全局通讯和同步。本文提出了一种与经典格拉姆 一施密特正交化方法在数学上等价的细长矩阵列向量正交
ｆｏｒ Ｓｋｉｎｎｙ Ｍａｔｒｉｃｅｓ
宋君强。龚西平．张理论。赵文涛，吴建平
ｊｍ－ｑｔ啦，ＧＯＮＧ ＳＯＮＧ
Ｘａ－ｐｉＩｌｇ。ＺＩｔＡＮＧ Ｌｉ－ｌｕｎ。ＺＨＡＯ Ｗｅｎ－ｔａｏ。ＷＵ Ｊｉａｎ－ｐｉｎｇ
（国防科学技术大学计算机学院。湖南长沙４１００７３）
ｆＳｃｌ）ｏｏｌ ｏｆ Ｃｏｍｐｕｔｅｒ Ｓｃｉｅｎｃｅ。Ｎａｔｉｏｎａｌ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｄｅｆｅｎｓｅ Ｔｅｃｈｎｏｌｏｇｙ，Ｃｈａｎｇｓｈａ ４１００７３。Ｃｈｉｎａ）
摘要：本文提出了一种与经典格拉姆一施密特正交化方法得到相同理论结果的细长矩阵的正交化方法。该方法在
增加部分冗余计算的情况下，将经典格拉姆一施密特正交化方法中的向量内积计算转换为可同时计算，使之更适合于并行
wk.baidu.com
计算。数值实验表明，该方法是正确的。
Ａｂｓｔｒａｃｔ：Ｔｈｉｓ ｐａｐｅｒ ｐｒｏｐｃｌｓｅｓ ａｎ ｏｒｔｈｏｎｏｒｍａｌｉｚａｔｉｏｎ ｍｅｔｈｏｄ ｗｈｉｃｈ ｃａｎ ｇｅｔ ｔｈｅ ｓａｍｅ ｔｈｅｏｒｅｔｉｃａｌ ｏｒｔｈｏｎｏｒｍａｌ ｒｅｓｕｌｔｓ
（ｇ．，ｉ，ｑ．，Ｉ）＝（口．．＾，ｑ．．ｉ）＝０
（８）
和
（口．，ｉ，ｑ．．１）＝（口．．Ｉ＇ａ．，ｉ）＝０
（９）
万方数据
将式（６）中的否¨代入式（９），有：
ｋ＇－Ｉ
（口．．Ⅲ．‘。ｔ＋∑町“ｊ）＝０
Ｊ暑１
将式（１０）表示成矩阵向量形式为：
（口．．１酊．．１）
缸．．１’口．．２）… Ｑ．。ｌ’口．州）
（ｎ．．２皿．．１） Ｑ．ｔ２ｍ．．２）…Ｑ．．２啦．Ｈ）
ｉ
ｉ
ｉ
（口．。卜ｌ，ａ．。１）（口．，卜１．ａ．．２）…（口．．卜ｌ，ａ。．卜１）
（１０）
巨
ａ。．·）
ａ。．ｔ）
求解线性方程组（１１）就可得到否¨与式（２）类似的形
式：
＾
～¨＝∑ｏ：ｉａｑ ··＾２厶 ｏ：ｉａ。·．·ｉ
（Ｌ１ｌ２ｚ）Ｊ
其中，劬，０１２，…，口卜１是线性方程组（１１）的解，ｍ＝１。 将否¨进行归一化，最关键的是求它的范数，实际上
·收稿日期：２００８－０７．０４；修订日期：２００９—０３—２０ 基金项目：国家科技支撵计划资助项目（２００６ＢＡＣ０２１３００）；国家自然科学基金资助项目（４０５０５０２３，６０８０３０３９）
作壹茸俞：宋君强（１９６２一）．男．湖南宁乡人，研究员，博士生导师，研究方向为科学工程计算和并行算法；龚西平，博士。研究方向为 糟复ｉｆ蔓：熟值算选塑大型软丝！堂理论！焦■副研蜜晏，研究方向为大规模并行计算；赵文涛，博士，副研究员，研究方向为流体力 堂和许行计算；吴建平，博士。副研究员，研究方向为科学一［＝程计算和并行算法。
通讯地址：４１００７３湖南省长沙市国防科学技术大学计算机学院软件所；Ｔｅｌ：１３６０７４８９７８３；Ｅ－ｍａｉｌ：Ｊｕｎｑｉａｎｇ＠ｎｕｄｔ．ｅｄｕ．ｃｎ Ａｄｄｒｅｓｓ：Ｉｎｓｔｉｔｕｔｅ ｏｆ Ｓｏｆｔｗａｒｅ，Ｓｃｈｏｏｌ ｏｆ Ｃｏｍｐｕｔｅｒ Ｓｃｉｅｎｃｅ，Ｎａｔｉｏｎａｌ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｄｅｆｅｎｓｅ Ｔｅｃｈｎｏｌｏｇｙ，Ｃｈａｎｇｓｈａ。ＨＯＲＢＳ ４１００７３， Ｐ．比Ｃｈｉｎａ
ｓｐａｎ（ｑ．。１，ｑ．．２，…，ｑ．．·）２ ｓｐａｎ（ｎ．．１，口。．２，…，口．．＾） （４）
根据式（１）和式（３），对任意ｋ＝ｌ，２，…，ｍ，有：
（口“，ｑ．．★）＝（口．．Ｉ，ａ．．‘）＝０，ｉ＝１，２，…，奄一１（５） 自然，对任意屉＝１．２＇．．·，ｍ，ｑ¨的计算过程可以分 为两步，首先计算：
４．１数值实验环境 本数值实验都是在ＨＰ机群上完成的，它由７个节点
（每个节点４ＣＰＵ）组成，每个ＣＰＵ为Ｉｔａｎｉｕｍ２ １．６ＧＨｚ。 每个节点物理内存为４ＧＢ，ｓｗａｐ交换空间为８ＧＢ。测试的 细长矩阵的元素均是随机生成的，所有计算均采用双精度。 ４．２数值实验结果
块正交化方法产生的矩阵么。的列正交偏离性用
ｋ－－Ｉ
孔＝“＾＋∑ｏｔｌａ－．ｉ，五＝ｌ，２＇．．·，仇 （６）
然后再将蚕¨进行归一化，就可得到ｑ¨：
ｑ．．ｔ一否．．ｔ／６孑．．１ ０
（７）
对任意忌＝１，２＇．．·，ｍ，ｑ¨是否¨的归一化向量，所
以；“与ｑ¨具有与ｑ“（ｉ一１，２，…，矗一１）和ａ“（ｉ＝
１．２，…，惫一１）相同的如式（３）和式（５）的正交性质，即：
可先求（云¨，面¨），有：
卢＝（蚕¨，蚕“）＝（∑叩“，∑ｍ口“）＝
ｉ＝ｌ
ｌ＝１
Ｉ
卜ｌ ＾
∑ａ｝（口¨，眠ｔ）＋∑∑２ａ“眠』（口．．朋¨）（１３）
ｉ＝１
ｉｆＩｊ＝州
于是计算：
口ｉ＝ａ，／妇。ｉ＝１，２，…，ｋ
（１４）
就可得到式（２）中的系数，也就是说式（２）中的系数和式
（１２）中的系数只差一个因子１／晒。至此可一次算出 ａ“（ｉ＝１，２，…，优）和ａｑ（．『＝ｉ，ｉ＋１，…，ｍ）的内积，然 后对任何惫＝１，２，…，ｍ根据式（１１）～式（１４）得到式（２）中 的系数，再根据矩阵与上三角矩阵乘得到所需要的标准正
有很多方法比如Ｈｏｕｓｅｈｏｌｄｅｒ变换法、Ｇｉｖｅｎｓ变换法 以及格拉姆一施密特正交化方法来实现矩阵Ａ的ＱＲ分解， 从而得到矩阵Ｑ，渤以及上三角矩阵Ｒ。砌。Ｒ表示成列向 量的形式为Ｒ＝Ｅｒ¨，．．．２…ｒ．，。］，其中ｒ．ｔｉ＝ ［，１．ｔ ｒ２．ｆ．”ｒ卅，ｔ］ｒ。在这里仅考虑其正交化得到的列向 量的标准正交基矩阵Ｑ帕。。下面提出的方法称为块正交化 方法。
．苎Ｌ
＞：（（２（七一１）３＋１１（ｋ一１）２／２—１３（ｈ一１）／２）＋
ｋ。‘＝一ｌ （３（愚＋１）ｋ／２＋１）＋志）＝
２ｌ＇
＞：２（ｋ一１）３＋７七２—１５ｋ＋１３一
ｔ。。＝一ｌ
打１４／２＋优３—５Ｍ一１３ｍ／２ 当ｍ较小时候，算法１中（２）～（６）的浮点操作数目较 小，计算时间基本可以忽略。
４数值实验
ＣＮ４３—１２５８／ＴＰ ＩＳＳＮ １００７—１３０Ｘ
计算机工程与科学
ＣＯＭＰＵＴＥＲ ＥＮＧＩＮＥＥＲＩＮＧ＆ＳＣＩＥＮＣＥ
文章编号：１００７－１３０Ｘ（２０１０）０４－００９０－０３
２０１０年第３２卷第４期 Ｖ０１．３２，Ｎｏ．４。２０１０
细长矩阵的块正交化方法’
Ａ Ｂｌｏｃｋ Ｏｒｔｈｏｇｏｎａｌｉｚａｔｉｏｎ Ｐｒｏｃｅｄｕｒｅ
９０
万方数据
化方法。虽然该方法在计算量上与经典格拉姆一施密特正 交化方法相当，但不同的向量内积可以同时计算‘¨，这样它 就更适合于并行计算。
２问题引入
Ａ，一∈辩加是细长矩阵。即ｍ＜＜挖，并且它的列向 量是线性无关的。计算列向量的标准正交基矩阵Ｑ一，其 中Ｑ丁Ｑ＝Ｌ。。将矩阵Ａ表示成列向量的形式为Ａ＝ Ｅａ．．１ ａ．．２…ａ．．。］，其中ａ．．ｆ＝［口¨口２．ｆ…“．ｉ］了’。同 理，Ｑ表示成列向量的形式为Ｑ＝［ｑ．。１ ｑ¨…ｑ．．。］，其 中ｑ¨一Ｅｑｌｔｆ ｑ２．ｆ…铂。ｉ］７。
交基矩阵Ｑ一。 ３．２算法及复杂性
算法块正交化方法
（１）计算矩阵ｃ：，。一（Ａ。）ｒＡ脚
（２）Ｆｏｒｋ＝１。２。…。ｍ
（３）
如果ｋ＞１，求解线性方程ｃｌ；卜ｌｔｌＩ／ｒ－ｌｆｆｌ。卜ｌ
＝一Ｃｌ。卜ｌｏ；ｍ ２ １；
（４）
，Ｆ 算 范 数 ｐ ＝ 。∑Ⅲ 口 口 Ｃ ＋ Ｈ∑斛 。∑州 ２ 啦 Ｃ
卢＝１／扫； （５）ＭｈＩ一皿ｌ。Ｉ，
ｗｉｔｈ ｔｈｅ ｔｒａｄｉｔｉｏｎａｌ Ｇｒａｒｎ－Ｓｅｈｍｉｄｔ ｏｒｔｈｏ－ｇｏｎａｌｉｚａｔｉｏｎ ｐｒｏｃｅｓｓ ｆｏｒ ｓｋｉｎｎｙ（１０ｎｇ）ｍａｔｒｉｃｅｓ．Ｔｈｅ ｍｅｔｈｏｄ ｌｅｔｓ ｔｈｅ ｉｎｎｅｒ ｐｒｏｄｕｃｔｓ
ｂｅ ｃｏｍｐｕｔｅｄ ａｔ ｔｈｅ ｓａｍｅ ｔｉｍｅ ａｓ ｔｈｅ ｔｒａｄｉｔｉｏｎａｌ Ｇｒａｍ－Ｓｃｈｍｉｄｔ ｏｒｔｈｏｇｏｎａｌｉｚａｔｉｏｎ ｐｒｏｃｅｓｓ。ａｎｄ ｉｓ ｍｏｒｅ ｓｕｉｔａｂｌｅ ｆｏｒ ｐａｒａｌｌｅｌ ｃｏｍｐｕｔｉｎｇ ｂａｓｅｄ ｏｎ ８０ｍｅ ｒｅｄｕｎｄａｎｔ ｃｏｍｐｕｔａｔｉｏｎｓ．Ｔｈｅ ｅｘｐｅｒｉｍｅｎｔｓ ｃｏｎｆｉｒｍ ｔｈｅ ｃｏｒｒｅｃｔｎｅｓｓ ｏｆ ｔｈｅ ｍｅｔｈｏｄ．
中图分类号：ＴＰ３０１
文献标识码：Ａ
１ 引言
计算线性无关向量的标准正交基是很多科学计算应用 的一个摹本问题￡ｌ】。行维数远远超过列维数的细长矩阵在 很多科学计算问题中出现，比如统计分析中观测量个数远 大于变鼍的个数时，用较少的向量来张成求解的子空间等。 计算细长矩阵列向量的标准正交基通常采用基于ＱＲ分解 的传统方法【２￣６］。ＱＲ分解主要通过Ｈｏｕｓｅｈｏｌｄｅｒ变换法、 Ｇｉｖｅｎｓ变换法和格拉姆一施密特正交化方法得到。在没有 舍入误差的情况下，格拉姆一施密特正交化方法得到的矩阵 Ｑ和Ｒ与Ｇｉｖｅｎ或Ｈｏｕｓｅｈｏｌｄｅｒ变换法得到的是一致的 （不计符号差别）。格拉姆一施密特正交化方法在数值上不 是很稳定，但它的良好计算特性使其在实际应用中被经常 使用。
Ｑ＆。Ｑ一非对角最大元素的值来表示。表１是不同的押和 ｍ时，Ｑ§。Ｑ胁非对角最大元素的值。表ｌ中的值基本上
都是计算机的双精度的舍入误差所容许的值，说明产生的 矩阵么×埘在数值上是列正交的。表ｌ结果表明，本文提到 的正交化方法是正确的。
５结束语
该算法是正确有效的。
（６） Ｅｎｄｄｏ
（７）计算么跏＝４：：Ｍ砌 算法得到细长矩阵Ａ，。的列向量的标准正交基矩阵
Ｑ。，实际上如果在没有舍入误差的情况下，它就是Ａ一
的ＯＲ分解的列正交矩阵Ｑ。但是，该块正交化方法没得 到ＱＲ分解的上三角矩阵Ｒ湘，而是得到一个上三角矩阵
‰，它们满足Ｑ洳＝＾×Ｍ。糯，并且‰是Ｒ，一的
逆。根据线性代数知识，如果Ａ．Ⅺ，．是列满秩的。则上三角