数理解析研究所講究録
第 1791 巻 2012 年 57-65

57

ハウスホルダー変換に基づく直交化法の最近の進展〜並列計算高性能計算の観点から〜
山本 有作
神戸大学大学院システム情報学研究科計算科学専攻/JST-CREST

1

はじめに

任意の $m\cross n$ 行列 $A(m\geq n)$ は，$m\cross m$ 直交行列 $Q$ と $m\cross n$ 上三角行列 $R$ を用いて $A=QR$ と分解で$Q$ の最初の$R$ の最初の列からなる行列をきる．また，
行からなる行列を左とするとき，という分解も可能である．本稿では前者を完全 QR 分解 (full QR decomposition), 後者を部分 QR 分解(thin QR decomposition) と呼ぶことにする．
$A$の列ベクトルが独立な場合，部分 QR 分解において髭の$\tilde{Q},$$n$対角成分が正であるという条件を付けると， ,$\tilde{Q}$$A=Q^{-}\tilde{R}$$n$$\tilde{R}$は一意に定まる．一方，完全 QR 分解においては， は$R$

一意に定まるが， は一意には定まらず．後半の $m-n$ 本の列間で直交変換を行う自由度がある．
$Q$

QR 分解は，行列に着目すると， の列ベクトルを正規直交化する手続きと見なせる．一方，行列に着目すると，直交行列により を上三角化する手続きと見なせる．どちらの目的のためにも，QR 分解
は数値計算において頻繁に利用される．QR 分解を行うには多くの手法があるが，本発表では，数値的安定$A$$\tilde{Q}$

$R$

$A$

性に優れたハウスホルダー変換による手法を取り上げる．
QR 分解は，用途に応じて様々な条件での利用が考えられる．具体的には，(i) 直交行列

$Q$

(または

(ii) 直交行列

$Q$

(あるいは

(iii)

$A$

$\tilde{Q}$) を陽的に求める必要がある力 1,$Q^{T}$)

を，後で他の行列に作用させる必要があるか，の列ベクトルが最初から全部与えられる力 1, あるいは直交化の過程で 1 本ずつ与えられる力 1,などの条件が考えられる．本発表では，実際の数値計算アルゴリズムでよく出てくる組台せとして，次のよ
うな条件下での QR 分解を考える．
(A)

$A$

は最初から全部与えられ， を陽的に求める． を後で他の行列に作用させる必要はない．

(B)

$A$

は最初から全部与えられる．

$Q$$\tilde{Q}$$Q,\tilde{Q}$を陽的に求める必要はない．

$Q^{T}$

を後で他の行列に作用させる．

$Q$ を後で他の行列に作用さ(C) $A$ の列ベクトルは直交化の過程で 1 本ずつ与えられ， を陽的に求める．
せる必要はない．
$\tilde{Q}$(A) は標準的なベクトルの正規直交化である．たとえば第一原理分子動力学法では，各時間ステップで波動関数を直交化するため，このタイプの直交化を行う．また，ブロック版のクリロフ部分空間法による連立 1 次方
程式解法でも，このタイプの直交化が用いられる．(B) は，たとえばブロック行列 $C=[A, B](B\in R^{m\cross l})$$Q^{T}B$
において， を QR 分解した後，を計算する問題に相当する．この計算は，ブロック化された行列のQR 分解や，固有値計算のための帯行列化などで現れる．(C) は，GMRES 法など，アーノルディ過程に基づく行列計算で典型的に現れる．アーノルディ過程では，反復ごとにクリロフ部分空間が拡張されて新しい$A$

基底ベクトルが生成され，それを今までに求めた正規直交基底に対して直交化する．これがちょうど (C) のタイプの直交化となる．このタイプの直交化は，固有ベクトル計算のための逆反復法や MR3 法において，固有値が密集する場合にも現れる．

近年では，問題の大規模化に伴い，並列計算機や GPU(グラフィックスプロセッサ) など，高性能計算技術の利用が重要となっている．従来，QR 分解あるいは直交化は，並列化が難しい演算と見なされ，高性能計算におけるボトルネックと見なされることが多かった．しかし最近では，新しいアルゴリズムの開発などにより，様々なタイプの QR 分解において，効率的な計算が可能となっている．そこで本発表では，上記(A), (B), (C) の 3 つのタイプの QR 分解に対し，高性能計算の観点から，最近のトピックスを紹介する．

58
2

タイプ (A): TSQR 法とその後退誤差解析

タイプ (A) の計算は，ハウスホルダーに基づく標準的な QR 分解アルゴリズムにより実行できる．ただし，
このアルゴリズムを並列化すると，ハウスホルダー変換ごとに 2 回のプロセッサ間同期 (あるいはプロセッ
サ間通信) が必要となり，全体で $2n$ 回の同期が必要となって，オーバーヘッドが大きい．これに対して$A_{2}\in R^{(m/2)\cross n}$ に分割してそれぞれのを上下に 2 つの部分行列
QR 分解を計算し，
$A$
(Tall
これを
それらの結果を合成して の QR 分解を求める方法を提案した [1].
TSQR
and Skinny QR)
アルゴリズムまたは AllReduce アルゴリズムと呼ぶ．アルゴリズムは次のように書ける．

Langou は，行列$A$

$A_{1},$$A=\{\begin{array}{l}A_{1}A_{2}\end{array}\}$$\overline{A}=\{\begin{array}{l}R_{1}R_{2}\end{array}\}$,$A_{1}=Q_{1}R_{1}$,$\overline{A}=\overline{Q}\overline{R}$,$A_{2}=Q_{2}R_{2}$

,

.

(2)

ここで，式 (1) の第 2 式，第 3 式，式 (2) の第 2 式はそれぞれ$Q_{i}\in R^{(m/2)\cross(m/2)},$
り，$A_{1},$$A_{2}$, 互のハウスホルダー QR 分解であ$R_{i}\in R^{(m/2)\cross n}(i=1,2),\overline{Q}\in R^{m\cross m},\overline{R}\in R^{m\cross n}$

$A$

(1)

である．式 (1), (2) より，は次のように書ける．
$A=\{\begin{array}{ll}Q_{1} OO Q_{2}\end{array}\}\overline{Q}R^{-}$

右辺の最初の 2 つの行列は直交行列であるから，これは

$A$

.

(3)

の QR 分解を与える．

は，第 1 行から第すべてゼロの行列である．この非ゼロ構造を考慮すると， の QR 分解は実質的に $(2n)\cross n$ 行列の QR 分$\overline{A}$$n$行，第 $m/2+1$ 行から第 $m/2+n$ 行がそれぞれ上三角行列，それ以外の要素が$\overline{A}$の QR 分$m\gg n$ の場合，この部分の演算量は小さく，演算量の大部分は解となる．したがって，
解で占められる．これらの QR 分解は並列に計算できるため，2 個のプロセッサによる大粒度並列化が可の QR 分解に対しても TSQR 法を再帰的に適用することで，さらに多くのプロ能となる．また，$A_{1},$$A_{1},$$A_{2}$$A_{2}$

セッサを計算に利用できる．

以上のように，TSQR 法は並列計算に極めて適した手法であり，多くの応用が期待されている [2]. しかし，TSQR 法では，プロセッサ数が個のとき，ハウスホルダー変換を 回再帰的に行う必要があり，丸め誤差が増大しないかが心配になる．これに対して我々は，TSQR 法の後退誤差解析を行い，TSQR 法
が次の式を満たすことを示した [4].
を有限精度で実行して得られる QR 分解の因子$L$$2^{L}$$\hat{Q},\hat{R}$$A+\Delta A=Q\hat{R}$$\hat{Q}=Q+\Delta Q$,,$\Vert\Delta A||_{F}\simeq n\gamma_{c}.\mathscr{X}+Ln\gamma_{c\cdot 2n}$,$||\Delta Q||_{F}\simeq(n\gamma_{c}1^{m_{L}}+Ln\gamma_{c\cdot 2n})\sqrt{n}$ここで， はある厳密に直交な行列 (定義は [4] 参照) であり，$\Delta A,$$Q$$\Delta Q$

(4)

.

(5)

は誤差行列である．また，(6)$\gamma_{cm}=\frac{cmu}{1-cmu}$

であり， はマシンイプシロン， はある正の定数である．式 (4), (5) はそれぞれ，TSQR 法が後退誤差の意味で安定であること，計算された が直交行列に極めて近いことを示す．なお，これらの式で $L=0$ とすると，よく知られたハウスホルダー QR 分解の後退誤差解析の結果 [3] が得られる．これらの式で注目すがに関して減少関数となっている
べきは， があまり大きくない範囲では，誤差限界ことである．すなわち，TSQR 法は， が比較的小さい範囲では，従来のハウスホルダー QR 分解よりむ
しろ高精度となりうる．この結果は数値実験でも確認されており，TSQR 法の応用を拡げるための基盤に
$u$

$c$

$\hat{Q}$$L$$||\Delta A||_{F},$$||\Delta Q||_{F}$

$L$

$L$

なると考えられる．

3

タイプ (B): TSQR 法で生成されるコンパクト WY 表現の合成

本節では，TSQR 法で計算した

$Q$

を使って，

$m\cross l$$Q^{T}B$
行列 $B$ に対し，を計算する問題を扱う．計算機環境として，最近普及してきたマルチコアプロセッサと GPU のハイブリッド環境を用いることを考える．

59
この場合， の QR 分解は比較的複雑な計算であるから，マルチコア上で行い，$Q^{T}B$$A$の計算は，単純な計算であるから，GPU 上で行うのが自然である．しかし，前者に TSQR 法を用いた場合，得られる は，$L+1$ 個の直交行列を掛け合わせた形になり ( $L=1$ の場合が式 (3)), さTSQR 法の再帰構造に対応して，らに，各直交行列は多数の細かい非ゼロブロックを持つ行列となる．このような行列に対する $Q^{T}B$ の計算
では，ベクトル長が短くなり，GPU などのアクセラレータの性能を引き出しにくい．
そこで， を 1 つの大きな行列として表現し直すことを考える [5]. 従来のハウスホルダー QR 分解の$Q$$Q^{T}$

は $m$ 次単位行列) と表
行列 $Y,$ $n\cross n$ 下三角行列 $T$ を用いて $Q^{T}=I_{m}-YTY^{T}$ (現できることが知られており [6], これをコンパクト WY 表現と呼ぶ．コンパクト WY 表現を作成してお場合，適当な$I_{m}$$m\cross n$$Q^{T}B$ の計算は，ハウスホルダー変換を 1 個ずつ作用させるのと同じ $4mnl$ の計算量で，しかも行列
けば，乗算のみを使って計算できるので，高性能計算の観点から有利である．そこで，TSQR 法の場合でも，適$I_{m}-YSY^{T}$ の形の行列を構成することを考える．
当な $m\cross n$ 行列 $Y,$ $n\cross n$ 行列 を用いて，$Y$ の列数は一般に
$Q$
では収まら
法で得られるただし，TSQR
を直接この形に表現しようとすると，(
すなわち恒等変換としてでなく)ない．これは，この行列が実質的に作用する部分空間の次元が，通常の
$A$ を
に変換
ハウスホルダー変換の場合と異なり， より大きくなるからである．そこで，制限を緩め，する直交行列で，適当な $Y\in R^{m\cross n},$ $S\in R^{n\cross n}$ を用いて $I_{m}-YSY^{T}$ の形に書ける行列を求めることにする．これは，完全 QR 分解における $Q$ の非一意性を利用していることに相当する．
$Q$ の最初のとし，列からなる行列を
いま，
$S$

$n$

$\overline{R}$

$n$

$n$

$\tilde{Q}=[\tilde{Q}_{1}^{T}\tilde{Q}_{2}^{T}]^{T}(\tilde{Q}_{1}\in R^{n\cross n},\tilde{Q}_{2}\in R^{(m-n)\cross n})$$Y=\{\begin{array}{l}\tilde{Q}_{1}-I_{n}\tilde{Q}_{2}\end{array}\}$とおく．すると，この

$Y$

,

(7)

$S=(I_{n}-\tilde{Q}_{1})^{-1}$と $S$ は条件を満たす [5]. 実際，$(I_{m}-YSY^{T})A$$=$$\overline{R}$,(8)$(I_{m}-YSY^{T})^{T}(I_{m}-YSY^{T})$$=$$I_{m}$(9)が成り立つことが容易にわかる．こうして，求める表現が構成できた．この書き換えは， の計算のため$Q^{T}B$ の計算を大型の行列乗算に置き換えられることから，に約 $2mn^{2}$ の余分な演算量を必要とするものの，$\tilde{Q}$

ハイブリッド環境での計算には有効となる可能性がある．

4タイプ (C):

コンパクト WY 表現に基づく逐次添加型の直交化$a_{n}\in R^{m}(m\geq n)$ に対し，それらを正規直交化したベクトルる問題を考える．ただし，ベクトル $a_{i}(2\leq i\leq n)$ は最初から与えられるのではなく，

ベクトル$a_{1},$ $a_{2},$$q_{n}$$q_{1},$ $q_{2},$$\ldots,$$\ldots,$$q_{1},$を求め$\ldots,$ $q_{i-1}$

がわ

かって初めて計算できるとする．これを，ベクトル逐次添加型の直交化と呼ぶ．
ベクトル逐次添加型の直交化では，修正グラムシュミット法 (MGS 法) がよく使われる．しかし，MGSとするとき，ql,
法は計算に逐次性があることに加え，行列 $A=[a_{1}, a_{2}, \ldots, a_{n}]$ の条件数を$\kappa(A)$$q_{2},$ $\ldots,$ $q_{n}$に比例して悪化するという問題点がある．一方，古典的グラムシュミット法 (CGS 法)の直交性に関しては優れていを 2 回繰り返す直交化法も提案されている．この手法は，並列性およびの直交性が$\kappa(A)$$\{q_{i}\}$

るが，条件 $O(u\kappa(A))<1$ が成り立たないと適用できないという問題点がある [2].

これに対して，Walker [7] は，GMRES 法においてグラムシュミット法の代わりにハウスホルダー変換によらずに高い直交性を確保でに基づく直交化を用いるアルゴリズムを提案している．この手法は，きる点に特徴がある．さらに [7] では，複数のハウスホルダー変換を合成することで，直交化演算を行列ベ$\kappa(A)$

クトル積の形に変形し，並列粒度を高める方法も提案されている．しかし，後者については，直交性などの数値的安定性に関する理論的検討は行われておらず，また，並列計算機上での性能評価も行われていない．
我々は，Walker の後者の方法をハウスホルダー変換のコンパクト WY 表現を用いて書き直し，これが通常のハウスホルダー変換と同等の数値的安定性を持つことを示した [8]. また，この方法が GMRES 法に限らず一般のベクトル逐次添加型の直交化に利用できることを指摘するとともに，並列計算機上で実際に高い加速率を達成できることを示した．そのため，この方法は，クリロフ部分空間法による行列指数関数$\exp(C)x$の計算，固有値が密集する場合の逆反復法，$MR^{3}$法など，様々なアルゴリズムにおいて活用できると考えられる．以下，この方法について説明し，数値的安定性の解析を行う．

60
41

ハウスホルダー変換に基づくベクトル逐次添加型の直交化準備として，ハウスホルダー変換に基づくベクトル逐次添加型の直交化 [9] を Algorithml に示す．ここで，$House_{i}(x)$は，ベクトル$x$に対し，その最初の $i-1$ 個の要素を不変にし，第 $i+1$ 要素以下を にするハ
( は単位行列) を求める操作である．また，ej は単位行列の第 列ベ$0$

ウスホルダー変換 $P_{i}=I$ -tiyi
クトルを表す．

$I$$y_{*}^{T}$$i$[Algorithm 1]
do $i=1,n$Generate$a_{\dot{\triangleleft}}$from$q_{1},q_{2},$ $\ldots,q_{i-1}$.

$a_{i}’=P_{2-1}\cdots P_{2}P_{1}$ 趣$H_{i}=House_{i}(a’.\cdot)$$q_{i}=P_{1}P_{2}\cdots P_{*}e_{i}$end do
ここで $q$: の計算では，QR 分解の定義より$q_{i}$

が行列$P_{1}P_{2}\cdots P_{n}$の第 列であること，および $P_{i+1},$$i$による乗算がを不変にすることを用いた．このアルゴリズムでは， およびのハウスホルダー変換を順次適用していくため，MGS 法と同じ逐次性がある．
$a_{i}’$$e_{i}$

4.2$q_{i}$$\ldots,$$P_{n}$の生成において，複数コンパクト WY 表現

一方，複数のハウスホルダー変換 $P_{k}=I-t_{k}y_{k}y_{k}^{T}(1\leq k\leq i)$ の作用をまとめて行う方法として，コン$Y_{1}=[y_{1}],$
パクト WY 表現 [6] がある．コンパクト WY 表現では，

下三角行列$T_{k}$$\ovalbox{\tt\small REJECT}$$T_{k}$

すると，$T_{1}=[t_{1}]$として，$n\cross k$

行列$Y_{k},$$k\cross k$を次のように帰納的に定義する．

$P_{i}\cdots P_{2}P_{1}$

,

(10)

$=$

$[Y_{k-1}y_{k}]$$=$$\{\begin{array}{ll}T_{k-1} 0-t_{k}y_{k}^{T}Y_{k-1}T_{k-1} t_{k}\end{array}\}$

.

(11)

の積は次のように表現できる．
$P_{i}\cdots P_{2}P_{1}=I-Y_{i}T_{i}Y_{:}^{T}$

.

(12)

これにより，複数のハウスホルダー変換を順次適用する計算を，行列ベクトル積として実行できる．

43

コンパクト WY 表現に基づくベクトル逐次添加型の直交化

コンパクト WY 表現を使うと，Algorithml を行列ベクトル積を用いて書き直すことができる．アルゴリ
ズムは次のようになる．

[Algorithm 2]

do $i=1,n$Generate$a_{\dot{4}}$from$q_{1},$ $q_{2},$$\ldots,$ $q_{i-1}$$a_{i}’=(I-Y_{i-1}T_{1-1}Y_{i-1}^{T})*$$(t_{i}, y_{i})=House_{i}(a_{i}’)$$Y_{i}=[Y_{i-1}y_{i}]$$T_{i}=\{\begin{array}{ll}T_{i-1} 0-t_{i}y_{i}^{T}Y_{i-1}T_{i-1} t_{i}\end{array}\}$$q_{i}=(I-Y_{i}T_{i}^{T}Y_{i}^{T})e_{i}$end do.

61
Algorithm 2 は，[7] における直交化部分の計算を GMRES 法に特化しない形に切り出し，コンパクト WYから表現を用いて書き直したものに相当する．
を計算する方法は任意であり，ベクトル$q_{1},$ $q_{2},$$\ldots,$ $q_{i-1}$$a_{i}$

逐次添加型の直交化を用いる任意のアルゴリズムに適用できる．主要な演算はすべて行列ベクトル積の形なので，MGS 法や Algorithm 1 に比べて粒度の大きい並列化が可能となる．

4.4

数値的安定性

本項では，Algorithm 2 について数値的安定性の解析を行う．以下，ベクトル値を要素とするベクトルを

$|x|$

$x$

に対し，その各要素の絶対$|x|\leq|y|$ のような不等式は，要素ごとの不等式を表すとすで表す．また，る．まず準備として，ハウスホルダー変換および WY 表現に対する既存の誤差解析の結果を述べる．なお，
WY 表現に対する結果は [3] で与えられているが，証明や誤差限界における定数の詳細が省かれているため，
ここでは証明も含めて詳しく述べる．

補題 1([3], Lemma 18.3) $v_{i}(1\leq i\leq r)$ を正確なハウスホルダー変換のベクトルとし，条件$\hat{v}_{i}=v_{i}+\Delta v_{i}$を満たす$\hat{v}_{1},\hat{v}_{2},$$v_{i}$,(13)$|\Delta v_{i}|\leq\gamma_{cm}|v_{i}|$$U_{i}=P_{i}\cdots P_{2}P_{1}$ とする．いま，$P_{i}=I-vv^{T}$ とし，の近似ベクトルを魁とする．また，$\ldots,\hat{v}_{r}$を用いた (近似的な) ハウスホルダー変換を浮動小数点演算により行列 $A\in R^{m\cross n}$ に作用させて得られる結果をとする．このとき，$\hat{A}_{r+1}$$\hat{A}_{r+1}=U_{r}(A+\Delta A)$

,

(14)

$||\Delta A||_{F}\leq r\gamma_{cm}||A||_{F}$が成り立っ．

補題 2([3], Section 18.4)$v_{i},\hat{v}_{i},$数点演算で求めた WY 表現を$P_{i},$$U_{i}$

を用いて浮動小は補題 1 と同様とする．このとき，について次の式力城り立つ．
$i=1,2,$とすると，$\{\hat{v}_{k}\}_{k=1}^{i}$$\hat{U}_{i}=I+\hat{W}_{i}\hat{Y}_{i}^{T}$$\ldots,$$\hat{U}_{i}=I+\hat{W}_{i}\hat{Y}_{i}^{T}=P_{i}\cdots P_{2}P_{1}+\Delta\hat{U}_{i}=U_{i}+\Delta U_{i}$$d_{1}(m, i)$ は $m,$
ただし，

証明

$i$

,

$r$

.

$||\triangle U_{i}||_{F}\leq d_{1}(m,i)u$(15)

に関して低次の多項式オーダーで増大する関数である．

まず $i=1$ のときを考えると，$\hat{U}_{1}$$=$$I-$ vv $T=I-$ $(v_{i}+ Avi)(v_{i}+\Delta v_{i})^{T}$$=$$U_{1}+v_{i}\Delta v_{i}^{T}+\Delta v_{i}v_{i}^{T}+\Delta v_{i}\triangle v_{i}^{T}$$U_{1}+\Delta U_{1}$$\equiv$

.

(16)

よって，
$||\Delta U_{1}||_{F}$次に，$\Delta U_{i-1}$$\leq$$||v_{1}\Delta v_{1}^{T}||_{F}+\Vert\Delta v_{1}v_{1}^{T}||_{F}+||\Delta v_{1}\Delta v_{1}^{T}\Vert_{F}$$=$$2||v_{1}||_{2}||\Delta v_{1}||_{2}+||\Delta v_{1}||_{2}^{2}$$\simeq$$2\gamma_{cm}$が与えられたときに$\Delta$(17)砿を求める漸化式を導く．まず，WY 表現の定義より，次の式が成り立っ．
$\hat{U}_{i}$

この最後の式は，行列$\hat{U}_{i-1}$$=$$I+\hat{W}_{i}\hat{Y}_{i}^{T}=I+[\hat{W}_{i-1}-\hat{v}_{i}][fl(\hat{v}_{ii-1}^{T})\hat{Y}_{i\frac{T}{\hat{U}}1}]$$=$$I+\hat{W}_{i-1}\hat{Y}_{i-1}^{T}-\hat{v}_{i}fl(\hat{v}_{i}^{T}\hat{U}_{i-1})$$=$$\hat{U}_{i-1}-\hat{v}_{i}fl(\hat{v}_{i}^{T}\hat{U}_{i-1})$

.

(18)

に $I-\hat{v}\hat{v}^{T}$ を浮動小数点演算により作用させた場合の結果の式$fl(\hat{U}_{i-1}-fl(\hat{v}_{i}fl(\hat{v}_{i}^{T}\hat{U}_{i-1})))$

(19)

62
に似ているが，むしろそれより丸め誤差が入る機会が少ない．したがって，補題 1 の誤差評価がそのまま使える．補題 1 において $r=1,$ $A=\hat{U}_{i-1}$ とすると，最後の式は次のように書ける．
$\hat{U}_{i-1}-\hat{v}_{i}fl(\hat{v}_{1}^{T}\hat{U}_{i-1})=P_{\dot{*}}(\hat{U}_{i-1}+\Delta\hat{U}_{*-1})$,$||\Delta\hat{U}_{i-1}\Vert_{F}\leq\gamma_{cm}||\hat{U}_{i-1}||_{F}$

.

(20)

これを式 (18) に代入し，式 (15) で $i-1$ とした式を代入すると，$\hat{U}_{i}$$=$$P_{1}(\hat{U}_{i-1}+\Delta\hat{U}_{i-1})$$=$$P_{i}(U_{i-1}+\Delta U_{i-1}+\Delta\hat{U}_{i-1})$$=$$U_{i}+P_{i}(\Delta U_{i-1}+\Delta\hat{U}_{i-1})$$\equiv$

これより，$\Delta U_{i}$$U_{i}+\Delta U_{i}$

.

(21)

に関する漸化式が次のように得られる．
$\Delta U_{i}=P_{i}(\Delta U_{i-1}+\Delta\hat{U}_{i-1})$

.

(22)

これを解くと，次の式が得られる．
$\Delta U_{i}=\sum_{j=1}^{i-1}P_{i}P_{i-1}\cdots P_{j+1}\Delta\hat{U}_{j}+P_{i}P_{i-1}\cdots P_{2}\Delta U_{1}$ただし，$P_{i}P_{i-1}\cdots P_{j+1}$

.

(23)

のような行列の積は，最初の添字より最後の添字のほうが大きい場合は，単位行列を表すとする．これより，$\Vert\Delta U_{i}||_{F}$ここで，漸化式$\leq$$\sum_{j=1}^{i-1}||P_{i}P_{i-1}\cdots P_{j+1}\Delta\hat{U}_{j}\Vert_{F}+||P_{1}P_{i-1}\cdots P_{2}\Delta U_{1}||_{F}$$=$$\sum_{j=1}^{i-1}||\Delta\hat{U}_{j}||_{F}+||\Delta U_{1}||_{F}$$\leq$$\gamma_{cm}\sum_{j=1}^{i-1}||\hat{U}_{j}||_{F}+||\Delta U_{1}\Vert_{F}$$\leq$$\gamma_{cm}\sum_{j=1}^{i-1}||U_{j}||_{F}+\gamma_{cm}\sum_{j=1}^{i-1}||\Delta U_{j}||_{F}+||\Delta U_{1}||_{F}$$\simeq$$\{(i-1)\sqrt{m}+2\}\gamma_{cm}+\gamma_{cm}\sum_{j=1}^{i-1}||\Delta U_{j}||_{F}$$u_{i}= \{(i-1)\sqrt{m}+2\}\gamma_{cm}+\gamma_{cm}\sum_{j=1}^{1-1}u_{j}$$u_{i}$.

(24)の解は次のように与えられる．

〉杭 $+2\}$$=$$\sum_{j=1}^{i}(\gamma_{cm})^{i-j+1}\{(i- 1)$$=$$\gamma_{cm}\{\sqrt{m}\cdot\frac{\gamma_{cm}^{i}-i\gamma_{cm}+i-1}{(1-\gamma_{cm})^{2}}+2\cdot\frac{1-\gamma_{cm}^{i}}{1-\gamma_{cm}}\}$

.

(25)

に関して低次の多項式オーダーの関数である．一方，は，この漸化式において等号を不等号に変えた式を満たすから，その値は式 (25) の最右辺以下に抑えられる．(証明終)
これは $m,$$||\Delta U_{1}||_{F}$$i$

本補題では，
それぞれの誤差は考えるが，そと定義している．すなわち，といれらを掛け合わせたり， に加えたりするときに生じる丸め誤差は考えない．WY 表現は，
注意 1$\hat{U}_{i}=I+\hat{W}_{2}\hat{Y}_{\dot{*}}^{T}$$I$$\hat{W}_{i},\hat{Y}_{i}$$I+\hat{W}_{i}\hat{Y}_{i}^{T}$

う行列を陽的に計算してから，それを作用させるわけではないからである．WY 表現を作用させる際の誤差については，補題 4 で考察する．

63
注意 2

本補題の証明で鍵となる式は (18) である．この式は [3] で与えられている．しかし，厳密に考えるとこの式には疑問がある．なぜなら，この式では行列$\hat{U}_{i-1}$

を実在の行列として扱い，それに$\hat{v}_{i}$を用いた

ハウスホルダー変換を作用させると考えて補題 1 を適用しているが，実際の計算はこのようには行わないからである．WY 表現の生成における実際の計算では，$\hat{v}_{i}^{T}\hat{U}_{i-1}$を求めるのに，$\hat{U}_{i-1}$を陽的には計算せず，(26)$fl(\hat{v}_{i}^{T}+fl(fl(\hat{v}_{1}^{T}\hat{W}_{i-1})\hat{Y}_{i-1}^{T}))$のように計算を行う．このような計算に対して補題 1 をそのまま適用できるかどうかは，検討を要する．ただし本稿では，式 (18) が正しいと考えて誤差解析を行う．
注意 3

ここでは WY 表現に対する補題を示したが，同様の補題がコンパクト WY 表現についても成り立
つ．以下の補題についても同様である．

補題 2 より，直ちに Algorithm 2 で計算したベクトル系3$q_{1},$ $q_{2},$$m\cross n$

行列$\ldots,$$\triangle\tilde{Q}$$q_{n}$の直交性に関する次の命題が導ける．

$q_{1},q_{2},$ $\ldots,q_{n}$を Algorithm 2 で計算したベクトルとする．このとき，$m\cross n$ 列直交行列$\tilde{Q}$およびが存在して，$[q_{1}, q_{2}, \ldots, q_{n}]=\tilde{Q}+\Delta\tilde{Q}$が成り立つ．すなわち， , q2, .. .,$q_{1}$$q_{n}$,(27)$||\Delta\tilde{Q}\Vert_{F}\leq d_{1}(m,n)u$の正規直交系からのずれは，$A=[a_{1}, \ldots, a_{n}]$ の条件数によらず，$O(u)$ である．

Algorithm 2 での証明

$q_{i}$

り，の計算方法より， は$q_{i}$$(I -Y_{i}T_{i}^{T}Y_{i}^{T})e_{i}=(I-Y_{n}T_{n}^{T}Y_{n}^{T})e_{i}$$I-Y_{n}T_{n}^{T}Y_{n}^{T}$に注意). したがって

の第

$i$

列となる

$[q_{1}, q_{2}, \ldots, q_{n}]$$(Y_{i},$$T_{i}$の非ゼロ構造よは $I-Y_{n}T_{n}^{T}Y_{n}^{T}$ の部分行列となるから，補題 2 と注意 3 より，明らかに主張が成り立っ．(証明終)補題 2 は WY 表現の直交性に関する補題であったが，WY 表現の後退安定性については次の補題が成
り立っ．
補題 4

$B\in R^{m\cross l}$

このとき，行列とし，浮動小数点演算により$\Delta C\in R^{m\cross l}$$I-\hat{W}_{i}\hat{Y}_{i}^{T}$

を $B$ に作用させて得られる行列を $C$ とする．

が存在して，次の式が成り立つ．
(28)

,

$C$

$=$

$U_{i}B+\Delta C=U_{i}(B+U_{i}^{T}\Delta C)$$\Vert\Delta C||_{F}$$\leq$$[1+d_{1}(i,m)+d_{2}(i, m)d_{3}(i,m)(1+c_{1}(i, m, l)+c_{1}(m, i, l))]u||B||_{F}+O(u^{2})$ここで，鵜は補題 1 で定義した行列であり，$d_{1},$ $d_{2},$ $d_{3}$は$i$と $m$ のみに依存する定数， は $i,$$c_{1}$

$n,$

.
$l$

(29)
のみに

依存する定数である．

補題 4 は，補題 2, 行列乗算に対する誤差解析結果 [3], および$||\hat{W}_{i}||_{F},$$||\hat{Y}_{i}||_{F}$に対する上界の評価からは明らかで導かれる． は長さ而のベクトルを並べた行列であるから，ある．一方， は式 (18) により定義されるため，この式を用いてノルムを評価できる．なお，[3] では 2 ノ$\hat{W}_{i}$$\hat{v}_{1},\hat{v}_{2},$$||\hat{W}_{i}||_{F}=O(\sqrt{n})$$\ldots,\hat{v}_{i}$$\hat{Y}_{i}$ルムに対して同様の式を導いているが，ここでは前の補題との関係から，フロベニウスノルムで統一した．
補題 4 を用いると，Algorithm 2 の後退誤差に関する次の定理が導ける．
定理 5の第 要素で，それ以外の要素がゼロの行列とする．すなわち， は浮動小数点演算を用いて Algorithm 2 により計算した $A=[a_{1}, \ldots,a_{n}]$ の $R$因子とする．このとき，行列 $\Delta A\in R^{m\cross n}$ が存在して次の式が成り立つ．
$R\in R^{m\cross n}$

を，第 $(i,j)$ 要素 $(i\leq j)$ が$(I-t_{j}y_{j}y_{j}^{T})a_{j}’$

$i$

$R$

$A+\Delta A=U_{n}^{T}R$,$||\Delta A||_{F}\leq d_{4}(m,n)u||A||_{F}$$d_{4}(m, n)$ はここで， は補題 1 で定義した行列であり，$U_{n}$$m$.

と $n$ のみに依存する定数である．

(30)

64
証明

まず，アルゴリズムより現用に直した補題)$a_{j}’=(I-Y_{j-1}T_{j-1}Y_{j-1}^{T})$より，あるベクトル$\Delta aj$$a_{j}’=U_{j-1}(a_{j}+\Delta a_{j})$が成り立つ．さらに， の第
$R$

$j$

aj である．そこで，補題 4(をコンパクト WY 表と定数 $d_{4}’(m,j)$ が存在して，,列を $r_{j}\in R^{m}$ とすると，補題 1 より，あるベクトル$r_{j}=P_{j}(a_{j}’+\Delta a_{j}’)$

,

(31)

$\Vert\Delta a_{j}||_{2}\leq d_{4}’(m,j)u||a_{j}||_{2}$$\Delta a_{j}’$が存在して，(32)$||\Delta a_{j}’||_{2}\leq\gamma_{cm}||a_{j}’||_{2}\simeq\gamma_{cm}||a_{j}||_{2}$が成り立つ．そこで，これらの式より，$r_{j}$$||\Delta a_{j}’’||_{2}$$=$$P_{j}(U_{j-1}(a_{j}+\Delta a_{j})+\Delta a_{j}’)$$=$$U_{j}(a_{j}+(\Delta a_{j}+U_{j-1}^{T}\Delta a_{j}’))$$\equiv$$U_{j}(a_{j}+\Delta a_{j}’’)$$\leq$$(d_{4}’(m,j)u+\gamma_{cm})||a_{j}||_{2}$$\equiv$,$d_{4}’’(m,j)u||a_{j}||_{2}$

(33)

.

(34)

さらに， は第 $j+1$ 要素以降がゼロであるから，$P_{n}P_{n-1}\cdots P_{j+1}$$rj$$r_{j}=U_{n}(a_{j}+\Delta a_{j}’’)$

これらを $j=1,2,$$\ldots,$$n$をかけても変化しないことを用いると，

.

(35)

についてまとめて書くと，$R=U_{n}(A+\Delta A)$ ,$\Delta A=[\Delta a_{1}’’, \Delta a_{2}’’, \ldots, \Delta a_{n}’’]$

.

(36)

ここで，
$|| \Delta A||_{F}\leq\{\max_{1\leq j\leq n}d_{4}’’(m,j)\}u\sqrt{\sum_{j--1}^{n}||\Delta a_{j}’’||_{2}^{2}}=\{\max_{1\leq j\leq n}d_{4}’’(m,j)\}u||\Delta A||_{F}$

.

(37)

したがって，(38)$d_{4}(m,n)=1 \leq j\leq n\max d_{4}’’(m,j)$とおくと，定理が得られる．(証明終)定理 5 は，Algorithm 2 が後退安定であることを意味する．

4.5

他の手法との比較

コンパクト WY 表現に基づくベクトル逐次添加型の直交化手法を他の手法と比較した結果を表に示す．

ここで，CGS2 は CGS 法を 2 回繰り返す方法，House はハウスホルダー変換に基づく Algorithml の方法，
$Q=[q_{4}, \ldots, q_{n}]$ とした
WY 表現を用いた Algorithm 2 の方法を表す．また，直交性は，$cWY$ はコンパクトときの行列 $Q^{T}Q-I_{n}$ のノルムを表す．条件とは，この直交性を実現するために必要な，行列 $A=[a_{1}, \ldots, a_{n}]$に関する条件である．MGS, CGS2, House に関する結果は [2] によった．
表より，コンパクト WY 表現に基づく手法は，演算量は修正グラム・シュミット法に比べて多いものの，に対CGS2, Algorithm 1 とは同等であり，並列化のための同期回数が少なく，かつ，任意の$a_{1},$ $a_{2},$して適用できるという長所を持つ手法と位置付けられる．

$\ldots,a_{n}$

65
5

おわりに

ハウスホルダー変換に基づく QR 分解は，精度と安定性に優れた手法であるが，高性能計算の観点からは，

これまでボトルネックと見なされることが多かった．しかし，最近ではアルゴリズムの進展により，様々なタイプの QR 分解について効率的な計算ができるようになっている．本発表では 3 つのタイプの QR 分解を取り上げ，このことを示した．今後の課題としては，様々な行列計算アルゴリズムへの組み込みと性能．
精度評価，各アルゴリズムに対するより完全な理論誤差解析，計算機環境・問題サイズに応じたアルゴリズ
ムの自動チューニングなどが挙げられる．

謝辞
研究集会「科学技術計算における理論と応用の新展開」において発表の機会を下さった山本野人先生，小林
健太先生に感謝いたします．また，有益なコメントを下さった同研究集会の参加者の皆様に感謝いたします．本研究は科学研究費補助金の補助を受けている．

References[1] J. Langou: AllReduce algorithms: application to Householder QR factorization, presented at Preconditioning 2007, Toulouse, France, July 9-12, 2007.
[2] J Demmel, L. Grigori, M. Hoemmen and J. Langou: Communication-optimal parallel and sequential
QR and LU factorizations, LAPACK Working Notes, No. 204, 2008.
[3] N. Higham: Accumcy and Stability of Numerical Algorithms, SIAM, Philadelphia, 2002.

[4] D. Mori, Y. Yamamoto and S. -L. Zhang: Backward error analysis of the AllReduce algorithm forHouseholder QR decomposition, to appear in Japan Journal of Industrial and Applied Mathematics.
[5] 山本有作: TSQR アルゴリズムで生成される Compact WY 表現の合成について，日本応用数理学会
2011 年度年会講演予稿集，pp. 95-96 (2011).
[6] R. Schreiber and C. van Loan: A storage-efficient WY representation for products of Householdertransformations, SIAM Journal on Scientific and Statistical Computing, Vol. 10, No. 1 pp. 53-57(1989).

[7] H. Walker: Implementation of the GMRES method using Householder transformations, SIAM Journal on Scientific and Statistical Computing, Vol. 9, No. 1, pp. 152-163 (1988).

[8] Y. Yamamoto and Y. Hirota: A parallel algorithm for incremental orthogonalization based on thecompact WY representation, JSIAM Letters, Vol. 3, pp. S9-92 (2011).
[9] G. Golub and C. van Loan: Matrix Computations, Johns Hopkins Univ. Press, 1996.