第8回 配信講義　計算科学技術特論A （2023）

行列計算における高速アルゴリズム - 大規模連立1次方程式の反復解法 2023年 6月 8日 計算科学技術特論A 電気通信大学 大学院情報理工学研究科 情報・ネットワーク工学専攻 山本 有作

本講義の構成 • 第1回（6月8日） 「大規模連立1次方程式の反復解法」 – 代表的な行列計算である連立1次方程式の解法について，現在 最も広く利用されているクリロフ部分空間法の基礎を紹介する． • 第2回（6月15日） 「ポストペタ時代に向けた線形計算アルゴリズムの課題と 研究動向」 – 様々な行列計算について，ポストペタ計算機上での実装における 課題と，その解決のための取り組みを紹介する． – また，行列計算における量子アニーリングマシンの活用の可能性 についても簡単に触れる． （1/61）

目次 1. はじめに 2. クリロフ部分空間と射影 3. 代表的なクリロフ部分空間法 4. 前処理 5. 波束シミュレーション向けの解法 6. 終わりに （2/61）

はじめに （3/61）

大規模連立1次方程式 • 今回の講義で扱う問題 ・ 連立1次方程式 Ax = b ・ ただし，A∈Rm×m， A は正則， b∈Rm （右辺ベクトル）， x∈Rm （未知ベクトル） ・ 特に A が大規模な疎行列の場合を考える • 応用 ・ 差分法・有限要素法による偏微分方程式の解法 ・ 連立常微分方程式の解法（回路シミュレーションなど） ・ 連立非線形方程式の解法 ・ 他の線形計算の部品（shift-and invert Lanczos 法，櫻井-杉浦法） （4/61）

係数行列の分類と数値解法 • 係数行列の分類 ・ 密行列／疎行列／帯行列／構造を持つ行列 ・ 対称行列／非対称行列 ・ 正定値行列／その他の行列 • 行列の分類に応じた解法（例） ・ 幅の狭い帯行列ならば直接解法（LU分解法）が有利 ・ 対称正定値行列ならばコレスキー分解が利用可能 ・ 対称正定値行列ならば共役勾配法（CG法）が利用可能 参考文献： ・ 櫻井鉄也, 松尾宇泰, 片桐孝洋編： 「数値線形代数の数理とHPC」, 共立出版, 2018. ・ 張紹良編著： 「計算科学のための基本数理アルゴリズム」, 共立出版, 2019. （5/61）

直接解法と反復解法 • 直接解法 ・ 行列を A = LU とLU分解し，Ly = b，Ux = y を順に解く ・ A が疎行列の場合，ゼロでない要素のみに対して演算を行うことで， 演算量と記憶領域を削減可能（疎行列LU分解） ・ A が対称正定値行列の場合，コレスキー分解 A = LLT を用いることで， 記憶領域と演算量を半分にできる • 特徴 ・ 有限回の演算で解が求まる（所要時間がほぼ予測可能） ・ 悪条件の問題に強い ・ 行列の非零要素数に比べて、数倍～数十倍の記憶容量が必要なこと が多い（フィルインの問題） ・ 反復解法に比べると，演算量大 ・ 1回LU分解を行えば， b のみが異なる方程式を少ないコストで解ける （6/61）

直接解法と反復解法（続き） • 反復解法 ・ 真の解 x* に収束するベクトル列 x(0), x(1), x(2), … を順次生成していく 方法 • 特徴 ・ 反復回数が事前にはわからない（ことが多い） ・ 収束しない場合がある ・ 行列 A を変形せず，行列ベクトル積 Ap の形でのみ用いるため，行列 の記憶領域は A の非ゼロ要素数の分のみ # ただしベクトルの記憶領域は別に必要 ・ うまく収束する場合，演算量は直接解法に比べて遥かに小さい （7/61）

定常反復法と非定常反復法 • 定常反復法 ・ 線形の漸化式 x(n+1) = Cx(n) + d （C：定数行列，d：定数ベクトル）により 近似解を更新していく反復解法 ・ 上記の反復が収束し，収束先が真の解 x* になるように C，d を選ぶ ・ 収束次数は1次。収束率は C の絶対値最大の固有値に依存 ・ ガウス-ザイデル法，SOR法，ADI法，マルチグリッド法など • 非定常反復法 ・ x(n) → x(n+1) の漸化式が，x(n) に関して非線形・非定常な反復解法 ・ 一般に，定常反復法よりも収束が速い ・ クリロフ部分空間法，チェビシェフ準反復法など （8/61）

行列の特異値と条件数 • 行列の特異値 ・ A∈Rm×m に対し，ATA は対称非負定値行列であり，その m 個の固 有値はすべて 0 以上の実数となる ・ これらの平方根を大きい順に s1(A) ≧ s2(A) ≧ ・・・ ≧ sm(A) ≧0と 書き，Aの特異値と呼ぶ ・ 行列 A が正則 ⇔ sm(A) > 0 • 行列の条件数 ・ 正則行列 A∈Rm×m に対し，k(A) = s1(A)/sm(A) を A の条件数と呼ぶ ・ 単位行列の条件数は 1 ・ 条件数の大きな連立1次方程式は，次の2つの意味で解きにくい ～ # 残差 r(n) = Ax(n) – b が小さくても，x(n) が真の解から遠い場合がある # クリロフ部分空間法が収束しにくい （9/61）

クリロフ部分空間と射影 （10/61）

クリロフ部分空間 • 定義 ・ 行列 A∈Rm×m，ベクトル b∈Rm に対し，次のように定義される Rm の 部分空間 Kn(A; b) を A，b に関する n 次のクリロフ部分空間と呼ぶ Kn(A; b) = span{b, Ab, A2b, …, An–1b} • クリロフ部分空間法 ・ x(0) = 0 から始めて，x(n)∈Kn(A; b) （n ≧1）となるように近似解の列 {x(n)} を構成していく反復解法を，クリロフ部分空間法と呼ぶ # 0 以外の初期値から始める方法もある ・ K1(A; b) ⊆ K2(A; b) ⊆ ・・・ ⊆ Kn(A; b) なので，クリロフ部分空間法で は，探索する空間を広げながら近似解を更新していくことになる ・ Kn(A; b) の中で，どのような基準で近似解 x(n) を定めるかにより，様々 なクリロフ部分空間法がある （11/61）

クリロフ部分空間内での解の存在 • クリロフ部分空間の拡大 – Kk(A; b) ⊆ Kk+1(A; b) は明らか – あるところで拡大が止まって Kk(A; b) = Kk+1(A; b) となるはず． • 疑問 – 初めて Kk(A; b) = Kk+1(A; b) となるのはどんなときか？ – 拡大が止まることで，真の解 x* に到達できない可能性は？ KK44 = K5 K3 x3 K2 x2 K1 x1 x4 x* x0 = 0 （12/61）

クリロフ部分空間内での解の存在（続き） • 定理 – A を正則行列とする。このとき，Kk(A; b) = Kk+1(A; b) を満たす最小の k に対し，Ax = b の解 x* は Kk(A; b) に含まれる． • 証明 – 仮定より， Akb = c0b + c1Ab + ・・・ + ck–1Ak–1b を満たす c1, c2, …, ck–1 が存在する。ここで，c0 = 0 とすると，A の正則性より，k が最小である ことと矛盾が生じるので，c0≠0。そこで，両辺を c0 で割って整理すると， A((1/c0)Ak–1b – (ck–1/c0)Ak–2b – ・・・ – (c1/c0)b) = b。これは x*∈Kk(A; b) を意味する． 部分空間の拡大が止まって真の解 x* に到達できないという心配は不要． （13/61）

残差多項式 • 残差の表式 ・ Kn(A; b) における近似解を x(n) = d0b + d1Ab + ・・・ + dn–1An–1b と定め たとする ・ このとき，残差 r(n) = Ax(n) – b は次のように書ける r(n) = – b + d0Ab + d1A2b + ・・・ + dn–1Anb = jn(A)b ただし， jn(z) = – 1 + d0z + d1z2 + ・・・ + dn–1zn 定数項が –1の n 次多項式 これを残差多項式と呼ぶ ・ Kn(A; b) で x(n) を1つ定めるとは，残差多項式を1つ定めることと等価 ・ 残差多項式は，クリロフ部分空間法の解析で重要な役割を果たす （14/61）

クリロフ部分空間の正規直交基底の生成 • 正規直交基底の必要性 ・ 近似解 x(n) を Kn(A; b) の基底ベクトルの線形結合として表す ・ 基底として正規直交基底を用いると，解法の理論的導出のし易さでも， 数値的安定性の面からも有利 • 正規直交基底の生成 ・ 考え方 ・ ・ ・ ・ ・ b を正規化して q1 = b /∥b∥2 を作成 v2 = Aq1 v2 を q1 に対して直交化し，正規化して q2 を作成 v3 = Aq2 v3 を q1, q2 に対して直交化し，正規化して q3 を作成 ・ これは，逐次的に生成されるベクトル {vn} に対し，グラム・シュミット 法による直交化を行っていることに相当 ・ これを Arnordi 過程と呼ぶ （15/61）

クリロフ部分空間の正規直交基底の生成（続き） • Arnoldi 過程 新たなベクトルの生成（空間の拡張） 直交化 正規化 （16/61）

クリロフ部分空間の正規直交基底の生成（続き） • Arnoldi 分解 ・ Arnoldi 過程 より，Aqi は q1, q2, …, qi+1 を使って次のように書ける Aqi = h1i q1 + h2i q2 + ・・・ + hi+1,i qi+1 ・ そこで， （m×n 列直交行列） （(n+1)×n 行列） と定義すると，次の式が成り立つ ～ AQn = Qn+1Hn ・ これを，A の n 次の Arnoldi 分解と呼ぶ （17/61）

クリロフ部分空間の正規直交基底の生成（続き） • Arnoldi 分解（続き） ～ ・ Arnoldi 分解の式 AQn = Qn+1Hn に左から QnT をかけ，Qn の各列が直 交することに注意すると， （n×n ヘッセンベルグ行列） ・ これは，A を span{q1, q2, …, qn} = Kn(A; b) に射影した行列と見ること ができる # 固有値計算のための Arnoldi 法の原理 （18/61）

固有値計算のためのArnoldi法 • クリロフ部分空間の中での固有ベクトルの近似 – Kn(A; b) の中で A の近似固有ベクトル x を求めることとし，x = Qy （y∈Rn）とおく．また，対応する近似固有値を l とする – 残差 Ax–lx が Kn(A; x(0)) に直交するという条件（Ritz-Galerkin 条件） を課すと， QnT(Ax – lx) = QnTAQny – lQnTQn = Hny – ly = 0. したがって，k×k 行列に対する小型の固有値問題 Hn y = l y を解けばよい これは，ヘッセンベルグ行列に対する固有値問題となる – 近似固有値は l，近似固有ベクトルは x = Qny で与えられる （19/61）

A が対称行列の場合 • Lanczos 過程 ・ QnTAQn = Hn は対称行列となるから，非ゼロ構造を考えると， となり，Hn は対称3重対角行列 Tn となる ・ これは，Anoldi 過程において計算した内積 hji （i ≧ j+1）が自動的に 0 になることを意味する ・ すなわち，Aqi は qi と qi–1 に対してのみ直交化すればよい ・ これを利用して Arnoldi 過程の計算量を減らした方法を Lanczos 過 程と呼ぶ （20/61）

A が対称行列の場合（続き） • Lanczos 過程 新たなベクトルの生成（空間の拡張） 直交化 正規化 ・ ただし，an = hnn, bn = hn+1,n = hn,n+1 とおいた • Lanczos 過程の特徴 ・ qn と qn–1 のみから qn+1 が求まる（3項間漸化式） ・ 1反復あたりの演算量は n によらず一定 （21/61）

固有値計算のためのLanczos法 • クリロフ部分空間の中での固有ベクトルの近似 – Kn(A; b) の中で A の近似固有ベクトル x を求めることとし，x = Qy （y∈Rn）とおく．また，対応する近似固有値を l とする – 残差 Ax–lx が Kn(A; x(0)) に直交するという条件（Ritz-Galerkin 条件） を課すと， QnT(Ax – lx) = QnTAQny – lQnTQn = Tny – ly = 0. したがって，k×k 行列に対する小型の固有値問題 Tn y = l y を解けばよい これは，3重対角行列に対する固有値問題となる – 近似固有値は l，近似固有ベクトルは x = Qny で与えられる （22/61）

クリロフ部分空間内での近似解の決定 • 近似解の選び方 ・ 近似解 x(n) としてクリロフ部分空間 Kn(A; b) 内のどんなベクトルを選 ぶかは，自由度がある ・ 主に，残差最小化， Ritz-Galerkin 法，Petrov-Galerkin 法の3種類 の原理が使われる ・ どの原理を使うかによって，クリロフ部分空間法の様々な変種が生ま れる （23/61）

近似解の決定法 I： 残差最小化 • 方法 ・ 残差 r(n) = Ax(n) – b が最小になる x(n) を Kn(A; b) の中で選ぶ ・ 最も自然なアプローチ # 本当は誤差 e(n) = x(n) – x* を最小にしたいが，一般には困難 ・ GMRES法，MINRES法などがこのタイプに属する （24/61）

近似解の決定法 II： Ritz-Galerkin法 • 方法 ・ 残差 r(n) = Ax(n) – b が Kn(A; b) に直交するように x(n) を選ぶ Ax(n) – b ⊥ Kn(A; b) ・ これは，残差のうち Kn(A; b) に正射影した成分が 0 であることを示す • 残差の直交性 ・ K1(A; b) ⊂ K2(A; b) ⊂ ・・・ ⊂ Kn(A; b) （真部分空間）とするとき，残 差 r(0) = b, r(1), r(2), …, r(n–1) は Kn(A; b) の直交系をなす # r(n–1)∈Kn(A; b) を用いて，n に関する帰納法により容易に示せる （25/61）

近似解の決定法 II： Ritz-Galerkin法（続き） • A が対称正定値行列の場合の意味付け ・ f (x) = (1/2) xTAx – xTb とおくと， Ax = b の解は f (x) の最小値 ・ ここで， x(n)∈Kn(A; b) と制限すると，x(n) = Qny(n) （y(n)∈Rn）と基底の 線形結合で書けるので，f (x(n)) = (1/2) (Qny(n))T AQny(n) – (Qny(n))Tb ・ これを最小にする y(n) は，微分して QnT(AQny(n) – b) = 0 より求まる ・ これは，ちょうど Ax(n) – b ⊥ Kn(A; b) と等価 ・ すなわち，Ritz-Galerkin法はクリロフ部分空間 Kn(A; b)の中で f (x(n)) を最小化する x(n) を求めている （26/61）

近似解の決定法 II： Ritz-Galerkin法（続き） • A が対称正定値行列の場合の意味付け（続き） ・ いま，誤差 e(n) = x(n) – x* の A ノルムを考えると， ・ したがって，Ritz-Galerkin法は，誤差の A ノルムを最小化する x(n) を 求めていると考えることもできる • Ritz-Galerkin 法に基づく解法 ・ 対称正定値行列向け： CG法（共役勾配法） ・ 非対称行列向け： FOM法 # Ritz-Galerkin 法は関数空間での近似（有限要素法など）にも使われる （27/61）

近似解の決定法 III： Petrov-Galerkin法 • 方法 ・ Rm の適当な部分空間の列 L1⊂L2⊂L3⊂ ・・・ （Ln は n 次元空間）を 用意し，残差 r(n) = Ax(n) – b が Ln に直交するように x(n) を選ぶ Ax(n) – b ⊥ Ln ・ これは，残差のうち空間 Ln に正射影した成分が 0 であることを示す ・ Ln としては，b* を適当なベクトルとし，AT と b* に対するクリロフ部分 空間 Kn(AT; b*) を使うことが多い • Petrov-Galerkin 法に基づく解法 ・ Bi-CG法，QMR法 # Petrov-Galerkin 法は関数空間での近似にも使われる ・ Petrov-Galerkin 法から派生した解法： CGS法，Bi-CGSTAB法， GPBi-CG法 （28/61）

残差多項式による特徴付け • 残差最小化 ・ 定数項が –1 の n 次多項式の集合を Pn とする ・ このとき，残差最小化は，各 n における残差多項式 jn(z) を，次の最 小化問題の解になるように選ぶアプローチだと解釈できる min jn ∈Pn∥jn(A)b∥2 • Ritz-Galerkin法（A が対称正定値行列の場合） ・ x(n) の誤差を e(n) = x(n) – x* とすると，Ae(n) = Ax(n) – b = r(n) ・ よって， e(n)TAe(n) = r(n)TA–1r(n) = bTjn(A)A–1jn(A)b = bTA–1jn(A)Ajn(A)A–1b = e(0)Tjn(A)Ajn(A)e(0) = ∥jn(A)e(0)∥A ・ すなわち，Ritz-Galerkin 法は次の最小化問題を解いている min jn ∈Pn∥jn(A)e(0)∥A （29/61）

代表的なクリロフ部分空間法 （30/61）

代表的なクリロフ部分空間法 • 残差最小化に基づく解法 ・ GMRES法（Generalized Minimum Residual，一般化最小残差法） ・ MINRES法（Minimum Residual，最小残差法） • Ritz-Galerkin 法に基づく解法 ・ CG法（Conjugate Gradient，共役勾配法） ・ FOM法（Full Orthogonalizaion Method，完全直交化法） • Petrov-Galerkin 法に基づく解法 ・ ・ ・ ・ Bi-CG法（Bi-Conjugate Gradient，双共役勾配法） QMR法（Quasi Minimal Residual，準最小残差法） CGS法（Conjugate Gradient Squared，自乗共役勾配法） Bi-CGSTAB法（Stabilized Bi-CG，安定化双共役勾配法） （31/61）

GMRES法 • 原理 ・ Kn(A; b) の中で残差 r(n) = Ax(n) – b を最小にする x(n) を 選ぶ ・ x(n)∈Kn(A; b) より x(n) = Qny(n) （y(n)∈Rn）と書くと， ～ ・ ただし，第2の等号では Arnoldi 分解の式 AQn = Qn+1Hn を用いた ・ 第3の等号では，q1 = b /∥b∥2 であることを用いた ・ 第4の等号では，Qn+1 が列直交行列であることから，Qn+1 を除いても 2ノルムが変わらないことを用いた ・ 最後の式の∥・∥2 の中は n+1 次元ベクトル，y(n) は n 次元ベクトル であるから，残差の最小化は y(n) に関する最小2乗問題となる （32/61）

GMRES法（続き） • 最小2乗問題の解法 ～ ～ ・ 最小2乗問題 min∥Hny(n) –∥b∥2 e1∥2 を解くには，係数行列 Hn を ～ Hn= Un Rn （Un∈R(n+1)×n, Rn∈Rn×n）とQR分解し，連立1次方程式 Rny(n) = UnT∥b∥2 e1 を解けばよい ～ ・ Hn は (n+1)×n ヘッセンベルグ行列なので，そのQR分解はギブンス 回転を用いて O(n2) の計算量で行える ～ ～ ・ さらに，Hn–1 のQR分解が与えられたとき，Hn のQR分解は O(n) で計 算できる # 追加された第 n 列のみについて，ギブンス回転を作用させればよい （33/61）

GMRES法（続き） • アルゴリズム Arnoldi 過程による Kn(A; b) の 正規直交基底の生成 ～ ギブンス回転による Hn のQR分解 ～ （Hn のQR分解の結果を用いる） UnT∥b∥2 e1 の計算 後退代入により Rny(n) = UnT∥b∥2 e1 を解く x(n) = Qny(n) の計算 （34/61）

GMRES法（続き） • GMRES法の特徴 ・ 計算の途中で破綻が起こることはない ・ 残差が n とともに単調に減少 # 非対称行列向けのクリロフ部分空間法としては，最も頑健な方法 ・ 所要メモリ量と1反復あたりの演算量が n に比例して増加 # qn を q1, q1, …, qn–1 に対して直交化する必要があるため # 所要メモリ・演算量を削減する方法としてリスタートがある • MINRES法 ・ 対称行列（必ずしも正定値でない）のための残差最小化法 ～ ・ Hn が3重対角行列になることを利用し，所要メモリ量と1反復あたりの 演算量を n に依存しない定数にできる ・ 丸め誤差に弱く，x(n) 中に (k(A))2 のオーダーで丸め誤差が混入* * H. A. van der Vorst： “Iterative Krylov Methods for Large Linear Systems”, Cambridge Univ. Press, 2003. （35/61）

CG法 • 原理 ・ A は対称正定値行列とする ・ 残差 r(n) = Ax(n)–b ∈ Kn+1(A; b)が Kn(A; b) に直交するよう x(n) を選ぶ # Ritz-Galerkin 法 • Lanczos過程との関係 ・ Lanczos過程では，基底ベクトル qn+1∈ Kn+1(A; b) を，Kn(A; b) に直交 するよう選ぶ ・ したがって，Lanczos過程を実行し，qn+1 を適当にスケーリングすれば， CG法の残差 r(n) が得られるはずである ・ 残差 r(n) の計算法が決まれば，解 x(n) の計算法も決まる（p. 14） （36/61）

CG法（続き） • Lanczos過程における基底ベクトルのスケーリング ・ qn+1 は3項漸化式を満たすから，それをスケーリングした r(n) も3項漸 化式を満たすはず ・ r(n) が r(n–1)，r(n–2) と直交するという条件より， ・ さらに，両辺の各項を {b, Ab, …, Anb} で展開したときの b の係数を 考えると，左辺は 0，右辺の r(n–1)，r(n–2)，r(n) はそれぞれ –1（p. 14）で あるから， ・ これらにより r(n) の漸化式が定まる （37/61）

CG法（続き） • CG法のアルゴリズム（3項漸化式） – r(n) の漸化式に r(n) = Ax(n) – b を代入して x(n) の漸化式を導出 （38/61）

CG法（続き） • アルゴリズム（連立2項漸化式） ステップ幅の決定 探索方向に向かって xn を動かす 残差の計算 次の探索方向の決定 • 最小化問題の解法としての解釈 ・ 関数 f (x) = (1/2) xTAx – xTb を最小化 ・ pn–1 が第 n ステップでの探索方向 ・ pn–1 方向に進んで関数値が極小になるように，ステップ幅をan を決定 （39/61）

CG法（続き） • CG法の性質 ・ Kn(A; b) = span{b, Ab, A2b, …, An–1b} = span{x1, x2,…, xn} = span{r0, r1,…, rn–1} = span{p0, p1,…, pn–1} ・ 残差の直交性： rnTrj = 0 (j < n) ・ 探索方向の A-共役性： pnTApj = 0 (j < n) # 最小化問題の解法と見たとき，各方向 pn への探索の繰り返しにより，大 域的な最小値が求まることを保証 • CG法の特徴 ・ 計算の途中で破綻が起こることはない ～ ・ 誤差の A ノルム∥en∥A =∥xn – x*∥A が n とともに単調に減少 ・ 所要メモリ量と1反復あたりの演算量は n に依存しない # 3項間漸化式により計算が行えるため （40/61）

CG法（続き） • CG法の収束性 ・ 誤差の A ノルム∥en∥A = の減少について，次の評価式が成り立つ ただし，k は A の条件数 • 証明の概略 ・ 残差多項式で表現した∥en∥A の最小性より， ～ ただし，Pn は定数項が1の n 次多項式，L(A) は A の固有値の集合 ・ ここで，特に f (l) をチェビシェフ多項式に取って最右辺を評価すると 上の定理が得られる （41/61）

CG法（続き） • CG法の収束性（続き） ・ A の相異なる固有値の個数が m* の場合，CG法は高々m* 回で収束 する # 前ページの f (l) として，これら m* 個の固有値で 0 となる m* 次多項式 が取れるから，n = m* のとき，最右辺は 0 となる ・ 固有値が少数のクラスターに固まっている場合も，収束は速い # これらのクラスター上で小さな値を取る低次の多項式が作れるから • 前処理 ・ 固有値を1個あるいは少数個のクラスターに集中させることができれ ば，CG法の収束は速くなる ・ この目的のため，Ax = b を別の（対称正定値な係数行列を持つ）連 立1次方程式に変形して解く前処理（後述）が行われる （42/61）

Bi-CG法 • 原理 ・ 方程式 に対して形式的にCG法を適用する • アルゴリズム （43/61）

Bi-CG法（続き） • Breakdown ・ rn と sn とが直交すると，アルゴリズム中で分母が0になり，計算を続 行できなくなる。これは，双直交基底を生成できないことに相当する。 これを第1種の breakdown と呼ぶ ・ 一方，Tn のLU分解において，ピボットが0になることに相当する breakdown もある。これを第2種の breakdown と呼ぶ ・ これらは，いずれもアルゴリズム上の工夫によって回避可能 • Bi-CG法の特徴 ・ 残差は必ずしも n とともに単調に減少しない ・ 所要メモリ量と1反復あたりの演算量は n に依存しない # GMRES法と比べた場合の大きな長所 # 3項間漸化式により計算が行えるため ・ A と AT の両方による乗算が必要 # 反復当たりの行列ベクトル積の回数はGMRES法の2倍 （44/61）

Bi-CG法から派生した解法 • CGS法 ・ Bi-CG法を変形し，rn = (Rn(A))2b, sn = b* となるように工夫した手法 # snTrn は不変なので，様々な係数は Bi-CG 法と同じに計算できる ・ 残差多項式が2乗されるため，収束が速くなると予想される # ただし，収束の不規則性も2乗で効くため，Bi-CG法より不安定 ・ A による乗算が２回となり，AT による乗算は不要 • Bi-CGSTAB法 ・ Sn(z) = (1 – w0z)(1 – w1z)・・・(1 – wn–1z) という形の n 次多項式 Sn(z) （安定化多項式）を用いて，rn = Sn(A)R(A)b, sn = b* となるように工夫し た手法 # この場合も，様々な係数は Bi-CG 法と同じに計算できる ・ 各ステップで局所的に残差を最小化するように wn を決めることで， CGS法よりも収束が安定化される （45/61）

Bi-CG法から派生した解法 • GPBi-CG法 ・ Bi-CGSTAB法における安定化多項式 Sn(z) として，３項漸化式によ り生成される多項式を用いた方法 ・ 一般に，Bi-CGSTAB 法よりもさらに優れた収束性を持つ • QMR法 ・ Bi-CG法において，残差を擬似的に最小化するようにした方式 ・ 一般に，残差の収束は Bi-CG 法よりもスムーズ（振動が少ない） # 収束は，Bi-CG法に比べて必ずしも速くない # しかし，残差の振動が少ないほうが丸め誤差の点で有利 （46/61）

前処理 （47/61）

前処理の概念 • 前処理 ・ 行列 K が，行列 A を何らかの意味で近似しているとする ・ このとき，方程式 K–1Ax = K–1b を考えると，その係数行列 K–1A は A と比べて単位行列に近くなっていると考えられる。これにより，クリロ フ部分空間法の収束はより速くなると期待できる ・ このような方程式の変形を前処理（左前処理）と呼ぶ ・ 各反復では，Ap の代わりに K–1Ap を計算する • 右前処理，左前処理， ・ 右前処理： AK–1 y = b， x = K–1 y ・ 両側前処理： K1–1AK2–1 y = K1–1b， x = K2–1 y • 前処理行列 K に求められる性質 ・ K–1A がなるべく単位行列に近いこと ・ K の生成が容易であること ・ K–1x の計算が容易であること （48/61）

疎行列に対するLU分解とフィルイン • フィルインのレベル – 元々の非ゼロ要素 aij に対して，そのレベルを l(aij) = 0 と定義する – フィルインにより aij≠0 となったとき，aij のレベルを次式で定義する l(aij) = max(l(aik), l(akj)) + 1 レベル 0 レベル 1 レベル 2 レベル 3 レベル 4 レベル 5 （49/61）

様々な前処理 • 不完全LU分解 ・ 行列 A のLU分解の際に，フィルインを無視して近似的に A≒LU と 分解し，K = LU とおく。これを ILU(0) 前処理と呼ぶ # １段階のみのフィルインを許す ILU(1) 前処理，２段階のフィルイン（フィ ルインによって生じたフィルイン）までを許すILU(1) 前処理も定義される ・ 対称正定値行列の場合は，不完全コレスキー分解型前処理（IC(0) 前処理）が使われる ・ 幅広いクラスの行列に対して効果がある ・ 一般に並列化は困難 # K–1x の計算で，前進消去 L–1x ，後退代入 U–1y が必要となるため # 並列化向きの不完全LU分解についても，多くの研究がある ・ フィルインの絶対値の大きさによって無視するかどうかを決める ILUT前処理，U–1 の列ノルムの大きさによって無視するかどうかを決 める ILUC 前処理もあり，悪条件の問題で有効 （50/61）

様々な前処理（続き） • 近似逆行列前処理（SPAI; SParse Approximate Inverse） ・ 行列 A そのものでなく，A–1 を直接近似する行列 M を求める ・ 具体的には，M の非零要素位置を適当に決めた上で，最小化問題 minM∥I – AM∥F を最小２乗法により解く ・ 不完全LU分解に比べ，並列化が容易 ・ ILU 分解と同程度の効果を出すには，一般により多くの非零が必要 • 対角スケーリング ・ 対角行列 D1, D2により A → D1–1AD2–1 と変換し，対角要素を１とする ・ 簡単で並列化も容易だが，問題によっては有効 • 行と列の並べ替え ・ 適当な並べ替えを行った後に不完全LU分解を行うと，収束性が改善 することがある （51/61）

波束シミュレーション への応用 （52/61）

時間依存シュレーディンガー方程式とその離散化 • 波束シミュレーションの基礎方程式 – 時間依存シュレディンガー方程式 – H： エルミート演算子（ハミルトニアン） – 以下では簡単のため，H が , t に依存しない場合を考える • 空間離散化 – 差分法 以下ではこちらを扱う – Ritz-Galerkin法（有限要素法 / 基底関数による展開） （B ： 重なり行列） （53/61）

クランク・ニコルソン法 • 導出 – (1,1)次パデ近似： exp(iHDt) ≒ (I – (1/2)iHDt) –1(I + (1/2)iHDt ) – 標準的な導出法： 時刻 t + (1/2)Dt での中心差分 • 特徴 – 時間発展のために連立1次方程式を解くことが必要 – 時間に関して2次精度 – 波動関数のノルムを厳密に保存 • 時間発展演算子の固有値の絶対値は シュレディンガー方程式の時間発展の計算に最適 （54/61）

クランク・ニコルソン法に現れる連立1次方程式 • 方程式 • 係数行列の性質 – 大規模・疎行列・非エルミート – i 倍すると A + sI （A： 正定値エルミート行列，s ： 複素数）の形 • 数値解法 – 直接法： ガウスの消去法（帯行列向け解法，スパースソルバなど） – 反復法： クリロフ部分空間法（GMRES法，BiCG法など） 係数行列の構造を利用した効率的なアルゴリズムの設計を考える （55/61）

A + sI に対するGMRES法 • A のエルミート性の利用 – GMRES法は，クリロフ部分空間 Kn(A; b) = span{b, Ab, …, An–1b} （n = 1, 2, …）の中で残差を最小にする解を探索 – クリロフ部分空間のシフト不変性 Kn(A+sI; b) = Kn(A; b) より，エルミ ート行列 A に対するクリロフ部分空間が使える – ランチョス法と同様，直近の3本の基底ベクトルのみを用いて Kn(A+sI; b) の正規直交基底を生成可能 • 得られる解法の特徴 – 3項間漸化式に基づく解法 → 演算量・メモリ量の負担軽減 – 残差最小性を持つ → スムーズな収束性 • 問題点 – 前処理を行うと，シフト不変性が崩れてしまう – 前処理なしで，どの程度の収束性が得られるかが鍵 （56/61）

数値実験結果 • GMRES法は COCG法（Bi-CG法系統の解法）に比べてスムーズな収束 を実現．計算量は同程度． H. Seito, T. Hoshi and Y. Yamamoto: On Using the Shifted Minimal Residual Method for Quantum-mechanical Wave Packet Simulation, JSIAM Letters, Vol. 11, pp. 13-16 (2019). （57/61）

終わりに （58/61）

終わりに • 本発表では，まずクリロフ部分空間とその基底の生成法を説 明し，方程式 Ax = b をクリロフ部分空間に射影する３種類 の方法を述べた • これに基づき，GMRES法，CG法，Bi-CG法の３つのアルゴ リズムと，その変種を説明した • また，収束を加速させるための前処理法について，代表的な 手法を説明した • 実際の問題にクリロフ部分空間法を適用するには，係数行 列の分類を検討した後，Lis や PETScなどのパッケージを 使って様々な解法・前処理を試してみるのがよい （59/61）

参考文献 （60/61）

参考文献 • L. N. Trefethen and D. Bau III: “Numerical Linear Algebra”, SIAM, Philadelphia, 1997. • H. A. van der Vorst: “Iterative Krylov Methods for Large Linear Systems”, Cambridge University Press, Cambridge, 2003. • Y. Saad: “Iterative Methods for Sparse Linear Systems”, 2nd Ed., SIAM, Philadelphia, 2003. • R. Barrett et al.: “Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods”, SIAM, Philadelphia, 1994. • 藤野清次, 張紹良: “反復法の数理”, 朝倉書店, 1996. • 杉原正顯, 室田一雄: “線形計算の数理”, 岩波書店, 2009. • 櫻井鉄也, 松尾宇泰, 片桐孝洋： “数値線形代数の数理とHPC”, 共立出版, 2018. • 張紹良編著： “計算科学のための基本数理アルゴリズム”, 共立出版, 2019. （61/61）