ハウスホルダー変換に基づく直交化法の最近の進展 : 並列計算・高性能計算の観点から (科学技術計算における理論と応用の新展開)

ハウスホルダー変換に基づく直交化法の最近の進展

∼並列計算高性能計算の観点から∼

山本有作

神戸大学大学院システム情報学研究科計算科学専攻/JST-CREST

1

はじめに

任意の$m\cross n$行列$A(m\geq n)$

は，

$m\cross m$直交行列$Q$ と $m\cross n$上三角行列$R$_を用いて$A=QR$ と分解で

きる．また，

$Q$の最初の$n$列からなる行列を $\tilde{Q},$ _$R$の最初の

$n$

行からなる行列を左とするとき，

$A=Q^{-}\tilde{R}$

という分解も可能である．本稿では前者を完全

QR分解 (full QRdecomposition), 後者を部分QR分解

(thin _QRdecomposition)

_{と呼ぶことにする．}

$A$

の列ベクトルが独立な場合，部分

QR 分解において髭の

対角成分が正であるという条件を付けると，

$\tilde{Q}$, $\tilde{R}$

は一意に定まる．一方，完全

QR

分解においては，

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

_{分解は，行列}

$\tilde{Q}$

に着目すると，

$A$

の列ベクトルを正規直交化する手続きと見なせる．一方，行列

$R$

に着目すると，直交行列により

$A$を上三角化する手続きと見なせる．どちらの目的のためにも，QR_分解 は数値計算において頻繁に利用される．QR分解を行うには多くの手法があるが，本発表では，数値的安定 性に優れたハウスホルダー変換による手法を取り上げる． QR 分解は，用途に応じて様々な条件での利用が考えられる．具体的には， (i) 直交行列 $Q$ (または $\tilde{Q}$) を陽的に求める必要がある力1, (ii) 直交行列$Q$ (あるいは$Q^{T}$) を，後で他の行列に作用させる必要があるか， (iii) $A$_{の列ベクトルが最初から全部与えられる力 1, あるいは直交化の過程で 1 本ずつ与えられる力 1,} などの条件が考えられる．本発表では，実際の数値計算アルゴリズムでよく出てくる組台せとして，次のよ うな条件下でのQR分解を考える．

(A) $A$

は最初から全部与えられ，

$\tilde{Q}$

を陽的に求める．

_$Q$

を後で他の行列に作用させる必要はない． (B) $A$

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

$Q,\tilde{Q}$

を陽的に求める必要はない．

$Q^{T}$を後で他の行列に作用させる． (C) $A$の列ベクトルは直交化の過程で

1

本ずつ与えられ，$\tilde{Q}$を陽的に求める．_$Q$を後で他の行列に作用さ せる必要はない． (A)は標準的なベクトルの正規直交化である．たとえば第一原理分子動力学法では，各時間ステップで波動関 数を直交化するため，このタイプの直交化を行う．また，ブロック版のクリロフ部分空間法による連立1次方

程式解法でも，このタイプの直交化が用いられる．(B)

は，たとえばブロック行列

$C=[A, B](B\in R^{m\cross l})$

において，$A$_を QR_{分解した後，}$Q^{T}B$ _{を計算する問題に相当する．この計算は，ブロック化された行列の} QR 分解や，固有値計算のための帯行列化などで現れる．(C) は，GMRES

_{法など，アーノルディ過程に基}

づく行列計算で典型的に現れる．アーノルディ過程では，反復ごとにクリロフ部分空間が拡張されて新しい

基底ベクトルが生成され，それを今までに求めた正規直交基底に対して直交化する．これがちょうど

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

2

タイプ

(A):

TSQR

_{法とその後退誤差解析}

タイプ (A) _{の計算は，ハウスホルダーに基づく標準的な}QR_{分解アルゴリズムにより実行できる．ただし，}

このアルゴリズムを並列化すると，ハウスホルダー変換ごとに2回のプロセッサ間同期 (あるいはプロセッ

サ間通信) が必要となり，全体で $2n$回の同期が必要となって，オーバーヘッドが大きい．これに対して

Langou

は，行列

$A$_を上下に2_{つの部分行列}$A_{1},$$A_{2}\in R^{(m/2)\cross n}$に分割してそれぞれの QR分解を計算し，

それらの結果を合成して$A$_の QR分解を求める方法を提案した [1]. _これを TSQR (Tall and SkinnyQR)

アルゴリズムまたはAllReduceアルゴリズムと呼ぶ．アルゴリズムは次のように書ける．

$A=\{\begin{array}{l}A_{1}A_{2}\end{array}\}$ , $A_{1}=Q_{1}R_{1}$, $A_{2}=Q_{2}R_{2}$, (1)

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

.

(2)

ここで，式

(1)

_{の第 2 式，第 3 式，式}

(2) の第2式はそれぞれ$A_{1},$ $A_{2}$, 互のハウスホルダー QR分解であ

り，

$Q_{i}\in R^{(m/2)\cross(m/2)},$ $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}$

である．式

(1), (2) より，

$A$は次のように書ける． $A=\{\begin{array}{ll}Q_{1} OO Q_{2}\end{array}\}\overline{Q}R^{-}$

.

(3) 右辺の最初の2つの行列は直交行列であるから，これは$A$_の QR分解を与える． $\overline{A}$ は，第

1

行から第$n$行，第$m/2+1$行から第$m/2+n$行がそれぞれ上三角行列，それ以外の要素が すべてゼロの行列である．この非ゼロ構造を考慮すると，$\overline{A}$ のQR分解は実質的に $(2n)\cross n$行列のQR分 解となる．したがって，$m\gg n$_{の場合，この部分の演算量は小さく，演算量の大部分は} $A_{1},$ $A_{2}$ の QR分 解で占められる．これらのQR分解は並列に計算できるため，2 個のプロセッサによる大粒度並列化が可 能となる．また，$A_{1},$ $A_{2}$の QR分解に対しても TSQR法を再帰的に適用することで，さらに多くのプロ セッサを計算に利用できる． 以上のように，TSQR法は並列計算に極めて適した手法であり，多くの応用が期待されている [2]. し かし，TSQR法では，プロセッサ数が$2^{L}$ 個のとき，ハウスホルダー変換を $L$回再帰的に行う必要があり， 丸め誤差が増大しないかが心配になる．これに対して我々は，TSQR法の後退誤差解析を行い，TSQR法 を有限精度で実行して得られる QR分解の因子$\hat{Q},\hat{R}$が次の式を満たすことを示した [4].

$A+\Delta A=Q\hat{R}$, _{$\Vert\Delta A||_{F}\simeq n\gamma_{c}.\mathscr{X}+Ln\gamma_{c\cdot 2n}$}, (4) $\hat{Q}=Q+\Delta Q$, $||\Delta Q||_{F}\simeq(n\gamma_{c}1^{m_{L}}+Ln\gamma_{c\cdot 2n})\sqrt{n}$

.

(5)

ここで，$Q$ はある厳密に直交な行列 (定義は [4]参照) であり，$\Delta A,$ $\Delta Q$ は誤差行列である．また，

$\gamma_{cm}=\frac{cmu}{1-cmu}$ (6)

であり，$u$はマシンイプシロン，$c$はある正の定数である．式(4), (5)はそれぞれ，TSQR法が後退誤差の

意味で安定であること，計算された$\hat{Q}$

が直交行列に極めて近いことを示す．なお，これらの式で$L=0$ と

すると，よく知られたハウスホルダー QR分解の後退誤差解析の結果 [3]が得られる．これらの式で注目す

べきは，$L$があまり大きくない範囲では，誤差限界 $||\Delta A||_{F},$ $||\Delta Q||_{F}$が$L$に関して減少関数となっている

ことである．すなわち，TSQR法は，$L$が比較的小さい範囲では，従来のハウスホルダー QR分解よりむ しろ高精度となりうる．この結果は数値実験でも確認されており，TSQR法の応用を拡げるための基盤に なると考えられる．

3

タイプ

(B): TSQR

法で生成されるコンパクト

WY

表現の合成

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

この場合，

の QR

分解は比較的複雑な計算であるから，マルチコア上で行い，

の計算は，単純な計

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

場合，適当な $m\cross n$行列$Y,$ $n\cross n$下三角行列$T$を用いて$Q^{T}=I_{m}-YTY^{T}$ ($I_{m}$は$m$次単位行列) と表

現できることが知られており [6], これをコンパクト WY

表現と呼ぶ．コンパクト

WY表現を作成してお

けば，$Q^{T}B$の計算は，ハウスホルダー変換を1個ずつ作用させるのと同じ $4mnl$_{の計算量で，しかも行列}

乗算のみを使って計算できるので，高性能計算の観点から有利である．そこで，TSQR 法の場合でも，適

当な $m\cross n$行列$Y,$ $n\cross n$行列 $S$を用いて，$I_{m}-YSY^{T}$の形の行列を構成することを考える．

ただし，TSQR法で得られる $Q$ を直接この形に表現しようとすると，$Y$_{の列数は一般に} $n$では収まら

ない．これは，この行列が実質的に (すなわち恒等変換としてでなく) 作用する部分空間の次元が，通常の

ハウスホルダー変換の場合と異なり，$n$ より大きくなるからである．そこで，制限を緩め，$A$を $\overline{R}$に変換

する直交行列で，適当な $Y\in R^{m\cross n},$ $S\in R^{n\cross n}$_を用いて $I_{m}-YSY^{T}$の形に書ける行列を求めることに

する．これは，完全 QR分解における $Q$ の非一意性を利用していることに相当する．

いま，

$Q$ の最初の$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}\}$ , $S=(I_{n}-\tilde{Q}_{1})^{-1}$ (7) とおく．すると，この $Y$ _と $S$_{は条件を満たす} [5]. _実際， $(I_{m}-YSY^{T})A$ $=$ $\overline{R}$ , (8) $(I_{m}-YSY^{T})^{T}(I_{m}-YSY^{T})$ $=$ $I_{m}$ (9) が成り立つことが容易にわかる．こうして，求める表現が構成できた．この書き換えは，$\tilde{Q}$ の計算のため に約$2mn^{2}$_{の余分な演算量を必要とするものの，}$Q^{T}B$_{の計算を大型の行列乗算に置き換えられることから，} ハイブリッド環境での計算には有効となる可能性がある．

4

タイプ

(C):

コンパクト

WY 表現に基づく逐次添加型の直交化

ベクトル$a_{1},$ $a_{2},$_$\ldots,$$a_{n}\in R^{m}(m\geq n)$ に対し，それらを正規直交化したベクトル$q_{1},$ $q_{2},$ $\ldots,$ $q_{n}$ を求め

る問題を考える．ただし，ベクトル $a_{i}(2\leq i\leq n)$ _{は最初から与えられるのではなく，$q_{1},$ $\ldots,$}$q_{i-1}$がわ

かって初めて計算できるとする．これを，ベクトル逐次添加型の直交化と呼ぶ．

ベクトル逐次添加型の直交化では，修正グラムシュミット法 (MGS法) _{がよく使われる．しかし，MGS}

法は計算に逐次性があることに加え，行列$A=[a_{1}, a_{2}, \ldots, a_{n}]$の条件数を$\kappa(A)$ とするとき，ql,_{$q_{2},$ $\ldots,$ $q_{n}$}

の直交性が$\kappa(A)$ に比例して悪化するという問題点がある．一方，古典的グラムシュミット法 (CGS法) を

2

回繰り返す直交化法も提案されている．この手法は，並列性および $\{q_{i}\}$の直交性に関しては優れてい るが，条件$O(u\kappa(A))<1$が成り立たないと適用できないという問題点がある [2].

これに対して，

Walker

[7]

_は，

GMRES

法においてグラムシュミット法の代わりにハウスホルダー変換 に基づく直交化を用いるアルゴリズムを提案している．この手法は，$\kappa(A)$ によらずに高い直交性を確保で

きる点に特徴がある．さらに

[7]

_{では，複数のハウスホルダー変換を合成することで，直交化演算を行列ベ}

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

41

ハウスホルダー変換に基づくベクトル逐次添加型の直交化

準備として，ハウスホルダー変換に基づくベクトル逐次添加型の直交化

[9]を Algorithml

_{に示す．ここで，}

$House_{i}(x)$ _{は，ベクトル}$x$ に対し，その最初の$i-1$個の要素を不変にし，第$i+1$要素以下を $0$にするハ

ウスホルダー変換$P_{i}=I$-tiyi$y_{*}^{T}$ ($I$は単位行列) を求める操作である．また，

ej

は単位行列の第$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}$ の第$i$列であること，および$P_{i+1},$

$\ldots,$$P_{n}$

による乗算が$e_{i}$ を不変にすることを用いた．このアルゴリズムでは，$a_{i}’$ および $q_{i}$の生成において，複数

のハウスホルダー変換を順次適用していくため，MGS法と同じ逐次性がある．

4.2

コンパクト

WY

表現

一方，複数のハウスホルダー変換

$P_{k}=I-t_{k}y_{k}y_{k}^{T}(1\leq k\leq i)$

_{の作用をまとめて行う方法として，コン}

パクト WY表現[6]がある．コンパクト WY 表現では，$Y_{1}=[y_{1}],$ $T_{1}=[t_{1}]$ として，$n\cross k$行列$Y_{k},$ $k\cross k$

下三角行列$T_{k}$ を次のように帰納的に定義する．

$\ovalbox{\tt\small REJECT}$ _$=$ $[Y_{k-1}y_{k}]$, (10) $T_{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}$ の積は次のように表現できる． $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

Algorithm2 は，[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}$, $|\Delta v_{i}|\leq\gamma_{cm}|v_{i}|$ (13)

を満たす $v_{i}$ の近似ベクトルを魁とする．また，$P_{i}=I-vv^{T}$ とし，$U_{i}=P_{i}\cdots P_{2}P_{1}$ とする．いま，

$\hat{v}_{1},\hat{v}_{2},$$\ldots,\hat{v}_{r}$を用いた (近似的な) ハウスホルダー変換を浮動小数点演算により行列$A\in R^{m\cross n}$に作用さ

せて得られる結果を $\hat{A}_{r+1}$

とする．このとき，

$\hat{A}_{r+1}=U_{r}(A+\Delta A)$, $||\Delta A||_{F}\leq r\gamma_{cm}||A||_{F}$ (14)

が成り立っ．

補題2([3], Section 18.4) $v_{i},\hat{v}_{i},$ $P_{i},$ $U_{i}$

は補題

1

と同様とする．このとき，

$\{\hat{v}_{k}\}_{k=1}^{i}$ を用いて浮動小

数点演算で求めた WY表現を $\hat{U}_{i}=I+\hat{W}_{i}\hat{Y}_{i}^{T}$

とすると，

$i=1,2,$

$\ldots,$$r$について次の式力城り立つ．

$\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}$, $||\triangle U_{i}||_{F}\leq d_{1}(m,i)u$

.

(15)

ただし，$d_{1}(m, i)$ は$m,$ $i$ に関して低次の多項式オーダーで増大する関数である．

証明 まず$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}$

$\equiv$ $U_{1}+\Delta U_{1}$

.

(16)

よって，

$||\Delta U_{1}||_{F}$ $\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}$ (17)

次に，$\Delta U_{i-1}$が与えられたときに $\Delta$砿を求める漸化式を導く．まず，WY表現の定義より，次の式が

成り立っ． $\hat{U}_{i}$ $=$ $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) この最後の式は，行列$\hat{U}_{i-1}$ _に$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)

に似ているが，むしろそれより丸め誤差が入る機会が少ない．したがって，補題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$ $U_{i}+\Delta U_{i}$

.

(21)

これより，$\Delta U_{i}$ に関する漸化式が次のように得られる．

$\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}$

.

(23)

ただし，

$P_{i}P_{i-1}\cdots P_{j+1}$

のような行列の積は，最初の添字より最後の添字のほうが大きい場合は，単位行

列を表すとする．これより，

$\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}$

.

(24)

ここで，漸化式

$u_{i}= \{(i-1)\sqrt{m}+2\}\gamma_{cm}+\gamma_{cm}\sum_{j=1}^{1-1}u_{j}$ _{の解は次のように与えられる．} $u_{i}$ $=$ $\sum_{j=1}^{i}(\gamma_{cm})^{i-j+1}\{(i- 1)$〉杭$+2\}$

$=$ $\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)

これは$m,$$i$ に関して低次の多項式オーダーの関数である．一方，$||\Delta U_{1}||_{F}$ は，この漸化式において等号を

不等号に変えた式を満たすから，その値は式 (25)の最右辺以下に抑えられる．(証明終) 注意 1

_{本補題では，}

$\hat{U}_{i}=I+\hat{W}_{2}\hat{Y}_{\dot{*}}^{T}$

と定義している．すなわち，

$\hat{W}_{i},\hat{Y}_{i}$

それぞれの誤差は考えるが，そ

れらを掛け合わせたり，

$I$

に加えたりするときに生じる丸め誤差は考えない．

WY

表現は，

$I+\hat{W}_{i}\hat{Y}_{i}^{T}$ とい う行列を陽的に計算してから，それを作用させるわけではないからである．WY表現を作用させる際の誤 差については，補題 4 で考察する．

2 本補題の証明で鍵となる式は (18)

_{である．この式は}

[3]

_{で与えられている．しかし，厳密に考える}

とこの式には疑問がある．なぜなら，この式では行列

$\hat{U}_{i-1}$

を実在の行列として扱い，それに

$\hat{v}_{i}$ を用いた ハウスホルダー変換を作用させると考えて補題1を適用しているが，実際の計算はこのようには行わない

からである．

WY

表現の生成における実際の計算では，

$\hat{v}_{i}^{T}\hat{U}_{i-1}$

を求めるのに，

$\hat{U}_{i-1}$ を陽的には計算せず， $fl(\hat{v}_{i}^{T}+fl(fl(\hat{v}_{1}^{T}\hat{W}_{i-1})\hat{Y}_{i-1}^{T}))$ (26) のように計算を行う．このような計算に対して補題 1 をそのまま適用できるかどうかは，検討を要する．た だし本稿では，式(18) が正しいと考えて誤差解析を行う． 注意3 ここでは WY 表現に対する補題を示したが，同様の補題がコンパクト WY表現についても成り立 つ．以下の補題についても同様である．

補題

2

より，直ちに

Algorithm2で計算したベクトル$q_{1},q_{2},$ $\ldots,q_{n}$の直交性に関する次の命題が導ける．

系3 $q_{1},$ $q_{2},$_$\ldots,$$q_{n}$ を Algorithm

2

で計算したベクトルとする．このとき，

$m\cross n$ 列直交行列$\tilde{Q}$ および

$m\cross n$行列$\triangle\tilde{Q}$

が存在して，

$[q_{1}, q_{2}, \ldots, q_{n}]=\tilde{Q}+\Delta\tilde{Q}$, $||\Delta\tilde{Q}\Vert_{F}\leq d_{1}(m,n)u$ (27)

が成り立つ．すなわち，

$q_{1}$,q2,

..

.,$q_{n}$

の正規直交系からのずれは，

$A=[a_{1}, \ldots, a_{n}]$ の条件数によらず，

$O(u)$ である．

証明 Algorithm 2 での $q_{i}$

の計算方法より，

$q_{i}$ は $I-Y_{n}T_{n}^{T}Y_{n}^{T}$の第$i$列となる $(Y_{i},$ $T_{i}$ の非ゼロ構造よ

り，

$(I -Y_{i}T_{i}^{T}Y_{i}^{T})e_{i}=(I-Y_{n}T_{n}^{T}Y_{n}^{T})e_{i}$ _{に注意). したがって} $[q_{1}, q_{2}, \ldots, q_{n}]$_は$I-Y_{n}T_{n}^{T}Y_{n}^{T}$の部分行列 となるから，補題 2 と注意 3 より，明らかに主張が成り立っ．(証明終)

補題 2 は WY 表現の直交性に関する補題であったが，WY表現の後退安定性については次の補題が成

り立っ．

補題4 $B\in R^{m\cross l}$

とし，浮動小数点演算により

$I-\hat{W}_{i}\hat{Y}_{i}^{T}$ を $B$_{に作用させて得られる行列を} $C$ _とする．

このとき，行列$\Delta C\in R^{m\cross l}$_{が存在して，次の式が成り立つ．}

$C$ $=$ $U_{i}B+\Delta C=U_{i}(B+U_{i}^{T}\Delta C)$, (28)

$\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})$

.

(29)

ここで，鵜は補題

1

で定義した行列であり，$d_{1},$$d_{2},$ $d_{3}$は $i$ と _$m$のみに依存する定数，

$c_{1}$ は $i,$ $n,$ $l$のみに

依存する定数である．

補題 4 は，補題 2,

行列乗算に対する誤差解析結果[3], および $||\hat{W}_{i}||_{F},$ $||\hat{Y}_{i}||_{F}$に対する上界の評価から

導かれる．

$\hat{W}_{i}$

は長さ而のベクトル$\hat{v}_{1},\hat{v}_{2},$$\ldots,\hat{v}_{i}$

を並べた行列であるから，

$||\hat{W}_{i}||_{F}=O(\sqrt{n})$ は明らかで

ある．一方，

$\hat{Y}_{i}$は式(18)

により定義されるため，この式を用いてノルムを評価できる．なお，

[3]

では2 ノ

ルムに対して同様の式を導いているが，ここでは前の補題との関係から，フロベニウスノルムで統一した． 補題 4 を用いると，Algorithm 2 の後退誤差に関する次の定理が導ける．

定理5 $R\in R^{m\cross n}$

_を，第

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

要素で，それ以外の要素がゼロの

行列とする．すなわち，

$R$は浮動小数点演算を用いてAlgorithm 2_{により計算した} $A=[a_{1}, \ldots,a_{n}]$ _の $R$

因子とする．このとき，行列$\Delta A\in R^{m\cross n}$が存在して次の式が成り立つ．

$A+\Delta A=U_{n}^{T}R$, $||\Delta A||_{F}\leq d_{4}(m,n)u||A||_{F}$

.

(30)

証明

_{まず，アルゴリズムより}

$a_{j}’=(I-Y_{j-1}T_{j-1}Y_{j-1}^{T})$_aj

である．そこで，補題

4(

をコンパクト

WY表

現用に直した補題)

_{より，あるベクトル}

$\Delta aj$ と定数$d_{4}’(m,j)$が存在して，

$a_{j}’=U_{j-1}(a_{j}+\Delta a_{j})$, $\Vert\Delta a_{j}||_{2}\leq d_{4}’(m,j)u||a_{j}||_{2}$ (31)

が成り立つ．さらに，

$R$の第$j$列を $r_{j}\in R^{m}$

とすると，補題

1

より，あるベクトル

$\Delta a_{j}’$が存在して，

$r_{j}=P_{j}(a_{j}’+\Delta a_{j}’)$, $||\Delta a_{j}’||_{2}\leq\gamma_{cm}||a_{j}’||_{2}\simeq\gamma_{cm}||a_{j}||_{2}$ (32)

が成り立つ．そこで，これらの式より，

$r_{j}$ $=$ $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}’’)$, (33)

$||\Delta a_{j}’’||_{2}$ $\leq$ $(d_{4}’(m,j)u+\gamma_{cm})||a_{j}||_{2}$

$\equiv$ _{$d_{4}’’(m,j)u||a_{j}||_{2}$}

.

(34)

さらに，

$rj$ は第$j+1$

要素以降がゼロであるから，

$P_{n}P_{n-1}\cdots P_{j+1}$ をかけても変化しないことを用いると，

$r_{j}=U_{n}(a_{j}+\Delta a_{j}’’)$

.

(35)

これらを $j=1,2,$$\ldots,$$n$についてまとめて書くと，

$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)

したがって，

$d_{4}(m,n)=1 \leq j\leq n\max d_{4}’’(m,j)$ (38) とおくと，定理が得られる．(証明終)

定理 5 は，Algorithm 2が後退安定であることを意味する．

4.5

他の手法との比較

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

ここで，CGS2はCGS法を 2 回繰り返す方法，Houseはハウスホルダー変換に基づく Algorithmlの方法，

$cWY$はコンパクト WY表現を用いた Algorithm

2

の方法を表す．また，直交性は，$Q=[q_{4}, \ldots, q_{n}]$ とした

ときの行列$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}$ に対

5

おわりに

ハウスホルダー変換に基づく QR分解は，精度と安定性に優れた手法であるが，高性能計算の観点からは， これまでボトルネックと見なされることが多かった．しかし，最近ではアルゴリズムの進展により，様々な タイプの QR分解について効率的な計算ができるようになっている．本発表では

3

つのタイプの QR分解 を取り上げ，このことを示した．今後の課題としては，様々な行列計算アルゴリズムへの組み込みと性能． 精度評価，各アルゴリズムに対するより完全な理論誤差解析，計算機環境・問題サイズに応じたアルゴリズ ムの自動チューニングなどが挙げられる．

謝辞

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

References

[1] J. Langou: AllReduce algorithms: application to Householder QRfactorization, presented at Pre-conditioning2007, Toulouse, France, July 9-12,2007.

[2] J Demmel, L. Grigori,M. Hoemmen and J. Langou: Communication-optimalparallel 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: Backwarderror analysis of the AllReduce algorithm for

Householder QRdecomposition,toappearinJapan Journal

_of

IndustrialandAppliedMathematics.

[5] 山本有作: TSQRアルゴリズムで生成される Compact WY表現の合成について，日本応用数理学会

2011 年度年会講演予稿集，pp. 95-96 (2011).

[6] R. Schreiber and C. van Loan: A storage-efficientWY representationfor productsof Householder transformations, SIAM Journal

on

_Scientific

and Statistical Computing, Vol. 10, No. 1 pp. 53-57 (1989).

[7] H.Walker: Implementationof theGMRES methodusingHouseholdertransformations,SIAM Jour-nal 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

the

compactWYrepresentation, JSIAMLetters,Vol. 3, pp. S9-92 (2011).