Drmač と Veselić の⼿法

Drmač らによって実装された LAPACK における⽚側ヤコビ法の前処理 [2, Section 5] は QR 前処理 を基にしているが，計算したい値（左特異ベクトル・右特異ベクトルを計算するかどうか）や，QR 前 処理がどれほどうまくいったか（⾏スケーリングした条件数が⼗分⼩さくなったかどうか）に応じて 処理を切り替える．具体的には，implicit Cholesky 法の反復回数を変更し，ごくまれに 1 回だけのも の，通常は 2 回，そして最⼤ 3 回まで⾏う．また，その結果⾏列の条件数がよくなっていることを⽤

いて，本来計算が必要な値を⾏列⽅程式から算出する．そこでここでは，左と右の両⽅の特異ベクト ルを計算する場合に選択される 2 段階の前処理と 3 段階の前処理についてのみ説明する．

2 段階の前処理

2 段階の前処理でははじめ，QR 前処理の⼿順通り

𝐴𝑃 = 𝑄 𝑅 (3.52)

を計算する．このとき𝑅の⾏スケーリングした条件数が⼗分良ければ，次の反復ではピボッティング を⾏わず，𝑅 = 𝑄 𝑅 を計算し，𝑅 に対して⽚側ヤコビ法を適⽤する．ただし，QR 分解と転置をま とめて LQ 分解

𝑅 = 𝐿𝑄 (3.53)

を計算する．

この後，𝐿 = 𝑈Σ𝑉 を⽚側ヤコビ法によって計算すれば，𝐴の SVD は

𝐴 = 𝑄 𝑅𝑃 = (𝑄 𝑈) Σ 𝑉 𝑄 𝑃 (3.54)

を得る．これには 2 回の⾏列積（正確には Householder 変換の列として𝑄 や𝑄 が保存されるため，

Householder 変換の適⽤）が含まれており，演算量が多い．そこで式 (3.53) より

𝑅 = 𝑈Σ𝑉 𝑄 (3.55)

を式変形し，

𝑄 𝑉 = 𝑅 𝑈Σ (3.56)

から𝑄 𝑉を計算できる．これは Householder 変換の適⽤の演算量の主要項2𝑛 に対して主要項とし て𝑛 の演算量で計算できるため，より⾼速である．また，Householder 変換の適⽤は複雑な計算とな るため，Level-3 BLAS で定義された演算である，三⾓⾏列による⾏列⽅程式の計算 xTRSM はより⾼

性能であると期待される．さらに⽚側ヤコビ法の中での𝑉の更新を省くことができる．これは⽚側ヤ コビ法の演算量を にする効果があり，⾼速化につながる．

この前処理は反復回数の削減のためには⼗分強⼒であり，また，次に⽰す 3 段階の前処理と⽐べて 演算量が少ないため，5 章で⽰す⾼性能計算向けの⽚側ブロックヤコビ法実装においてはこちらのも ののみを実装し，利⽤した．ただし，4 章の実験においては，2 段階の前処理と 3 段階の前処理を⾃動

で切り替える LAPACK のルーチンを利⽤しており，2 段階だけのものを実装した⾼性能計算機向け実 装との⽐較をしている．

3 段階の前処理 1 度⽬の前処理

𝐴𝑃 = 𝑄 𝑅 (3.57)

によって⼗分𝑅 の⾏スケーリングした条件数が⼩さくならなければ，さらに 2 回の implicit Cholesky 法の反復を繰り返す．つまり，

𝑅 𝑃 = 𝑄 𝑅 , (3.58)

𝑅 = 𝐿𝑄 . (3.59)

この後，𝐿 = 𝑈Σ𝑉 を計算する．そこで

𝐴 = 𝑄 𝑅 𝑃 = 𝑄 𝑃 𝑅 𝑄 𝑃 = 𝑄 𝑃 𝑄 𝐿 𝑄 𝑃 = (𝑄 𝑃 𝑄 𝑉)Σ(𝑈 𝑄 𝑃 ) (3.60) によって𝐴の SVD が計算できる．

ここでも 2 段階の前処理と同様に𝑉の計算を省くことができる．すなわち式 (3.58) より

𝑅 = 𝑈Σ𝑉 𝑄 , (3.61)

すなわち

𝑄 𝑉 = 𝑅 𝑈Σ. (3.62)

この後，𝑄 𝑃 (𝑄 𝑉)と𝑃 𝑄 𝑈を計算すればよいため，合計で 2 回の Householder 変換の適⽤と 1 回の

⾏列⽅程式の計算によって計算できる．

第 4 章

列ブロックペアの直交化⼿法の誤差解析

⽚側ブロックヤコビ法の主な計算は，あるステップにおける⾏列の列ブロックペア𝑋 = 𝐴 𝐴 の SVD を計算し，列直交化することである．いま𝑋の SVD を𝑋 = 𝑈 Σ 𝑉 と書けば，𝑉 を右側からか けることで𝑌 ≡ 𝑋𝑉 = 𝑈 Σ となるため，列直交化できる．この計算の特徴は 2 つあり，第⼀に𝑋が 縦⻑⾏列であることで，簡単のためすべての列ブロックが同じ幅を持っていると仮定すると，縦の⻑

さ対横の⻑さの⽐が𝑤 ∶ 2となる．第⼆に，𝑋の列スケーリングした条件数が多くの場合，⼩さくなる ことである．これは，もともとの⾏列𝐴が QR 前処理された⾏列であることが主な理由である．場合 によっては𝐴の列スケーリングした条件数が反復の間に上昇することがあるが，Demmel らの実験 [1]

によれば，この上昇幅は⼩さく，また反復に従って急速に1に収束していくことが⽰されている．𝑋 はある反復における𝐴の列ブロックペアであるから，𝑋の列スケーリングした条件数はその反復にお ける𝐴のものよりも⼩さくなる．そこで，𝑋の列スケーリングした条件数が⼩さいものと想定する．

𝑋の右特異ベクトル𝑉 によって列直交化することは単なる SVD と⾏列積の計算であるが，有限精 度の演算によって⾏う場合には問題があることが Drmač [92, Eq. 40] によって⽰されている．Drmač が解析している誤差は列直交化で⽣じる誤差の⼀部に過ぎないが，彼の結果は誤差が直交性に及ぼす 影響について重要なことを⽰唆する．いま計算によって求めた𝑋の右特異ベクトル ̂𝑉 が真の右特異 ベクトル𝑉 から誤差𝛿𝑉 だけずれており， ̂𝑉 = 𝑉 + 𝛿𝑉 と書けるとする．またある正則な対⾓⾏列𝐷

によって𝑋 = 𝑋 𝐷と書き，𝑋 はすべての列ベクトルがノルム1になっているとする．このとき𝑋 ̂𝑉 に

ついて次の⽅程式が成り⽴つ：

𝑋 ̂𝑉 = 𝑋𝑉 + 𝑋𝛿𝑉 = (𝑈 + 𝑋 𝛿𝐹)Σ , (4.1)

𝛿𝐹 = 𝐷𝛿𝑉 Σ . (4.2)

これは，𝛿𝑉 のノルムが⼩さい場合でも，直交⾏列𝑈 からのずれ𝑋’𝛿𝐹は⼤きくなりうることを⽰し ている．これは𝑋を⼗分に列直交化できないことにつながり，⽚側ブロックヤコビ法全体の収束が滞 ることにつながる．これは，Drmač [92, Eq. 41] によって後退誤差が⼩さくなることが⽰されているこ ととは対照的である．

しかし，上記は粗い解析であり，⾏列𝑋の構造を考慮してより詳しい解析を⾏うことで，直交性の 誤差に対してより⼩さな上界を求められる可能性がある．そこで，本研究では列ブロックペアの直交 化⼿法の 1 つである Hari らの V2 ⼿法 [101, Section 4.1]（以降から単に V2 ⼿法と呼ぶ）について，理 論的な誤差解析を⾏い，いくつかの条件が成り⽴つならば，直交性の誤差が⼩さいことを⽰す．V2 ⼿ 法は，Cholesky QR 法という誤差の⼤きな⼿法を内部的に利⽤するにも関わらず，実験においては⾼

精度に計算できている．そこでこの章では，まず V2 ⼿法についてその⼿順を⽰した後に，理論的誤差 解析の結果を⽰す．また，Hari の V2 ⼿法を⽤いて⽚側ブロックヤコビ法を実装したときの，実際の 𝑌の直交性や，SVD 全体の誤差を実験的に調べる．この章の内容は関連論⽂ [a] において簡略に⽰し た証明を詳細化し，⾮並列向けプログラムの実験結果を⽰したものである．