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

7.8K Views

May 22, 23

スライド概要

第8回 6月8日 行列計算における高速アルゴリズム1
本講義では、科学技術計算の基盤となる行列計算について、基本的なアルゴリズムの原理と高性能計算技術を紹介する。
第1回ではクリロフ部分空間に基づく連立1次方程式反復解法と固有値計算法の基礎について解説する。

シェア

またはPlayer版

埋め込む »CMSなどでJSが使えない場合

関連スライド

各ページのテキスト
1.

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

2.

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

3.

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

4.

はじめに (3/61)

5.

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

6.

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

7.

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

8.

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

9.

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

10.

行列の特異値と条件数 • 行列の特異値 ・ 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)

11.

クリロフ部分空間と射影 (10/61)

12.

クリロフ部分空間 • 定義 ・ 行列 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)

13.

クリロフ部分空間内での解の存在 • クリロフ部分空間の拡大 – 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)

14.

クリロフ部分空間内での解の存在(続き) • 定理 – 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)

15.

残差多項式 • 残差の表式 ・ 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)

16.

クリロフ部分空間の正規直交基底の生成 • 正規直交基底の必要性 ・ 近似解 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)

17.

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

18.

クリロフ部分空間の正規直交基底の生成(続き) • 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)

19.

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

20.

固有値計算のための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)

21.

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

22.

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

23.

固有値計算のための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)

24.

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

25.

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

26.

近似解の決定法 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)

27.

近似解の決定法 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)

28.

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

29.

近似解の決定法 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)

30.

残差多項式による特徴付け • 残差最小化 ・ 定数項が –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)

31.

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

32.

代表的なクリロフ部分空間法 • 残差最小化に基づく解法 ・ 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)

33.

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)

34.

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)

35.

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)

36.

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)

37.

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)

38.

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)

39.

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

40.

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

41.

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)

42.

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

43.

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

44.

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

45.

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)

46.

Bi-CG法から派生した解法 • CGS法 ・ Bi-CG法を変形し,rn = (Rn(A))2b, sn = b* となるように工夫した手法 # snTrn は不変なので,様々な係数は Bi-CG 法と同じに計算できる ・ 残差多項式が2乗されるため,収束が速くなると予想される # ただし,収束の不規則性も2乗で効くため,Bi-CG法より不安定 ・ A による乗算が2回となり,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)

47.

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

48.

前処理 (47/61)

49.

前処理の概念 • 前処理 ・ 行列 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)

50.

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

51.

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

52.

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

53.

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

54.

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

55.

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

56.

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

57.

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)

58.

数値実験結果 • 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)

59.

終わりに (58/61)

60.

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

61.

参考文献 (60/61)

62.

参考文献 • 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)