第1回(6月22日)

第2回(6月29日)

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

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

係数行列の分類
。
密行列／疎行列／帯行列／構造を持つ行列. .対称行列／非対称行列正定値行列／その他の行列. ●
行列の分類に応じた解法(例)
•幅の狭い帯行列ならば直接解法(LU分解法)が有利.対称正定値行列ならばコレスキー分解が利用可能.対称正定値行列ならば共役勾配法(CG法)が利用可能

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

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

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

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

。
定義•行列A∈R™×m,ベクトルb∈R™に対し,次のように定義されるR™の部分空間K„(A; b)をA, bに関するn次のクリロフ部分空間と呼ぶKn(A; b) = span{b, Ab, A2b, …, An-1b}
●
クリロフ部分空間法
•  x(0)=0から始めて, x(n)∈K,(A;b) (n ≥1)となるように近似解の列{x(η)}を構成していく反復解法を,クリロフ部分空間法と呼ぶ＃  0以外の初期値から始める方法もある• K 1 (A; b)ニK 2 (A; b)ニ•••ニK n (A; b)なので,クリロフ部分空間法では,探索する空間を広げながら近似解を更新していくことになる• K„(A; b)の中で,どのような基準で近似解x (n)を定めるかにより,様々なクリロフ部分空間法がある
(11/60)

•クリロフ部分空間の拡大
• Kn(A; b) ⊆ Kn+1(A; b)は明らか• K n (A; b) = K n+1 (A; b)となると,空間はそれ以上拡大しないこれにより，解に到達できなくなることはないか？
定理
。
・ Aを正則行列とする。このとき，K n (A; b) = K n+1 (A; b)を満たす最小のnに対し, Ax = bの解x*はK„(A; b)に含まれる
証明
.
• 仮定より,  Anb = c b + c 1 Ab + • • • + c n − 1 An− 1 bを満たすc 1 , c 2 , …, c n − 1が存在する.ここで, c 0 = 0とすると, Aの正則性より, nが最小であることと矛盾が生じるので, c 0 ⁄ 0そこで,両辺をc 0で割って整理すると, A((1/c o )An-1b - (c n-1/c o )An-2b - • • • - (c 1/c o )b ) = b .これはx * ∈ K n (A; b )を意味する
(13/60)
•残差の表式• K n (A; b)における近似解をx (n) = d o b + d 1 Ab + • • • + d n −1 A n−1 bと定めたとするこのとき，残差r(n) = Ax(n) − bは次のように書けるr(n) = − b + doAb + d 1 A 2 b + • • • + d n − 1 A n b = φ n (A)bただし, φn(z) = − 1 + doz + d1z² + ▪▪▪ + dn−1z” |定数項が−1のn次多項式これを残差多項式と呼ぶ. K„(A; b)でx (η)を1つ定めるとは,残差多項式を1つ定めることと等価残差多項式は,クリロフ部分空間法の解析で重要な役割を果たす
(14/60)

.
正規直交基底の必要性
.近似解x (n)をK,(A; b)の基底ベクトルの線形結合として表す.基底として正規直交基底を用いると,解法の理論的導出のし易さでも,数値的安定性の面からも有利
正規直交基底の生成
●
考え方. •  bを正規化してq 1 = b / || b || 2を作成•   νς = Aq 1 ・ v 2をq 1に対して直交化し,正規化してq 2を作成•  ・ v3をq1, q2に対して直交化し,正規化してq3を作成これは,逐次的に生成されるベクトル{v,}に対し,グラム・シュミット法による直交化を行っていることに相当.これをArnordi過程と呼ぶ
(15/60)

• Arnoldi過程
[Arnoldi process] q1 = b/||b||2 for n = 1, 2, 3, . . . v = Aqn for j = 1 to n hjn = qj v v = v — hjnQj end for hn+1,n = ‖ v‖ 2 9n+1 = v/hn+1,n end for新たなベクトルの生成(空間の拡張)直交化正規化
Arnoldi過程の変種
.直交化にハウスホルダー変換を用いる方法もある.演算量は多くなるが,基底の直交性と並列化可能性に関して優れるY. Yamamoto and Y. Hirota, “A parallel algorithm for incremental orthogonalization based on the compact WY representation", JSIAM Letters, 3, 89-92, 2011. (16/60)
. Arnoldi分解. Arnoldi過程より, Aq iはq 1 , q 2 , …, q i+1を使って次のように書けるAq i = h 1 i q 1 + h 2 i q 2  +  • • •  + h i+1,i q i+1そこで, . Qn = [q1|q2| · · · |qn] ,
h1n h2n
hn+1,nと定義すると,次の式が成り立つAQn=Qn+1În (m×n列直交行列) ((n+1)×n行列) .これを, Aのn次のArnoldi分解と呼ぶ(17/60)
• Arnoldi分解(続き) Arnoldi分解の式AQn = Qn+1Ã„に左からQ„をかけ, Q„の各列が直交することに注意すると, QTAQn = QTQn+1Hn 1 0 : 0 : () ll Hn () h1n h2n ll ≡ H n.  (n×nヘッセンベルグ行列)これは, Aをspan{q 1 , q 2 , …, q n } = K n (A; b )に射影した行列と見ることができる# 固有値計算のためのArnoldi法の原理
(18/60)

.クリロフ部分空間の中での固有ベクトルの近似– K„(A;b)の中でAの近似固有ベクトルxを求めることとし, x = Qy (γ∈R″)とおく.また,対応する近似固有値をλとする– 残差Ax–λxがK„(A; x (⁰)に直交するという条件(Ritz-Galerkin条件)を課すと, QnT(Ax – Ax) = QnTAQny – AQnTQn = Hny – ly = 0.したがって, k×k行列に対する固有値問題Hny = lyを解けばよいこれは,ヘッセンベルグ行列に対する固有値問題となる– 近似固有値はλ, 近似固有ベクトルはx = Q,yで与えられる
(19/60)

• Lanczos過程• Q n TAQ n = H nは対称行列となるから,非ゼロ構造を考えると,
となり, H nは対称3重対角行列T nとなるこれは, Anoldi過程において計算した内積h,i ( i ≥ j+1)が自動的に0になることを意味する.すなわち, Aq,はq,とq, 1に対してのみ直交化すればよいこれを利用してArnoldi過程の計算量を減らした方法をLanczos過程と呼ぶ
(20/60)

• Lanczos過程
[Lanczos process] β0 = 0,  qo = 0,  q1 = b/||b||2 for n = 1, 2, 3, . . . v = Aqn Xn = qn v v = v — βn–19n–1 — αn9n βn = ‖ v‖ 2 Qn+1 = v/ βn end for新たなベクトルの生成(空間の拡張)直交化正規化•ただし, α n = h nn, βn = h n+1,n = h n,n+1とおいた
• Lanczos過程の特徴
q nとq n −1のみからq n+1が求まる(3項間漸化式) • 1反復あたりの演算量はnによらず一定
.クリロフ部分空間の中での固有ベクトルの近似– K„(A;b)の中でAの近似固有ベクトルxを求めることとし, x = Qy (y∈R″)とおく.また,対応する近似固有値をλとする– 残差Ax–λxがK„(A; x (⁰)に直交するという条件(Ritz-Galerkin条件)を課すと, QnT(Ax – λx) = QnTAQny – λQnTQn = Tny – λy = 0.したがって, k×k行列に対する固有値問題Tny = Ayを解けばよいこれは, 3重対角行列に対する固有値問題となる– 近似固有値はλ, 近似固有ベクトルはx = Q,yで与えられる
(22/60)

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

●
方法
.残差r(n) = Ax(n) − bが最小になるx(n)をK n (A; b)の中で選ぶ最も自然なアプローチ. # 本当は誤差e(η) = x(η) − x*を最小にしたいが, 一般には困難• GMRES法, MINRES法などがこのタイプに属する

●
方法•残差r(n) = Ax(n) − bがK,(A; b)に直交するようにx(n)を選ぶAx(n) - b  ⊥  Kn(A; b)これは,残差のうちK n (A; bに正射影した成分が0であることを示す
残差の直交性
。
. K 1 (A; b) ⊂ K 2 (A; b) ⊂ ••• ⊂ K n (A; b) (真部分空間)とするとき,残差r(0) = b, r(1), r(2), ..., r(n-1)はK,(A; b)の直交系をなす#  r(n-1)∈ Kn(A; b)を用いて, nに関する帰納法により容易に示せる
(25/60)

•  Aが対称正定値行列の場合の意味付け▪ φ(x) = (1/2) x TAx − x Tbとおくと，Ax = bの解はφ(x)の最小値ここで,  x(n)∈ K,(A; b)と制限すると, x(n) = Q,y(n) (y(n)∈ R n )と基底の線形結合で書けるので, φ(x(η) = (1/2) (Q,y(η) T AQ,y(η) = (Q,y(η) Tbこれを最小にするy(w)は,微分してQ, T(AQ,y(w) – b) = 0より求まるこれは,ちょうどAx(n) – b  ⊥ Kn(A; b)と等価•すなわち, Ritz-Galerkin法はφ(x(η)を最小化するx(η)を求めている
(26/60)

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

方法
●
・ R r r rの適当な部分空間の列L 1 CL 2 CL 3 C ・・・ (L nはn次元空間)を用意し,残差r(n) = Ax(n) − bがL,に直交するようにx(n)を選ぶAx(n) - b  L  Lnこれは,残差のうち空間L,に正射影した成分が0であることを示す•  ∑ ,としては, ♭*を適当なベクトルとし,  A Tと♭*に対するクリロフ部分空間K„(A™; b*)を使うことが多いPetrov-Galerkin法に基づく解法。 •  Bi-CG法, QMR法#  Petrov-Galerkin法は関数空間での近似にも使われる•  Petrov-Galerkin法から派生した解法:  CGS法,  Bi-CGSTAB法, GPBi-CG法
(28/60)

残差最小化.定数項が–1のn次多項式の集合をP nとするこのとき,残差最小化は,各nにおける残差多項式φ„(z)を,次の最小化問題の解になるように選ぶアプローチだと解釈できるmin φη∈Ρη ‖ φη(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) = bTφn(A)A-1φn(A)b = bTA-1φn(A)Aφn(A)A-1b = e(0)Tφn(A)Aφn(A)e(0)=  || φn(A)e(0) || Aすなわち, Ritz-Galerkin法は次の最小化問題を解いているmin φη ∈ ρη ‖ φη (Α) e (0) ‖ Α
(29/60)

•
残差最小化に基づく解法
• 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/60)
原理
K„(A; b)の中で残差r(n) = Ax(n) − bを最小にするx(n)を選ぶ. x(n)∈ Kn(A; b)よりx(n) = Qny(n) (y(n)∈ Rn)と書くと, || Ax(n) - b||2    =    || AQny(n) - b||2 =   ||Qn+1 Ãny(n) - b||2 =   ||Qn+1(Îny(n) - ||b||2e1)||2 =   || Îny(n) - || b || 2e1 || 2. • ただし,第2の等号ではArnoldi分解の式AQ n = Q n+1H nを用いた第3の等号では, q 1 = b / || b || 2であることを用いた第4の等号では, Ω n +1が列直交行列であることから, Ω n +1を除いても2ノルムが変わらないことを用いた• 最後の式の|| • || ,の中はn+1次元ベクトル, y(n)はn次元ベクトルであるから,残差の最小化はy(n)に関する最小2乗問題となる
(32/60)

最小2乗問題の解法.最小2乗問題min ‖ ñ„ν ω – ‖ b ‖ ₂ e₁ ‖ ₂を解くには, 係数行列ñ„をÃ„= U„R„ (U„∈R(n+1)×n, R„∈Rn×n)とQR分解し,連立1次方程式R y(n) = U n T || b || 2 e 1を解けばよい•  Å„は(n+1)×nヘッセンベルグ行列なので,そのQR分解はギブンス回転を用いてO(n²)の計算量で行える• さらに, ˆ β n − 1のQR分解が与えられたとき, ˆ β nのQR分解はO(n)で計算できる#  追加された第n列のみについて,ギブンス回転を作用させればよい
(33/60)

アルゴリズム
[GMRES algorithm] q1 = b/||b||2,  δ = ||b||2e1 for n = 1, 2, 3, . . . v = Aqn for j = 1 to n hjn = q, Tv,  v = v − hjnQj end for hn+1,n = || U || 2,   9n+1 = v /hn+1,n T1n = h1n for j = 2 to n Y = Cj-1Yj-1,n + Sj-1hjn Tjn = − Sj−1Tj−1,n + Cj−1hjn Tj-1,n = Y end for δ = √ r²n + h²+1,n, cn = rnn/δ, sn = hn+1,n/δ rnn = Cnrnn + Snhn+1,n bn+1 = − snbn, bn = cnbn ρ = |δn+1| (= ||Αα(n) - δ||2) if ρ is small enough break end for for j = n,n − 1,...,1 yj = (bj - ∑ n end for x =  Li=1 YiQi Arnoldi過程によるK n (A; b)の正規直交基底の生成ギブンス回転による批nのQR分解( H nのQR分解の結果を用いる) UnT || b || 2 e 1の計算後退代入によりRny(n) = UnT || b || 2 e 1を解くx(n) = Qny(n)の計算(34/60)
GMRES法の特徴
●
.計算の途中で破綻が起こることはない.残差がnととに単調に減少#非対称行列向けのクリロフ部分空間法としては,最も頑健な方法•所要メモリ量と1反復あたりの演算量がnに比例して増加# q nをq 1, q 1, …, q n –1に対して直交化する必要があるため#所要メモリ･演算量を削減する方法としてリスタートがある
• MINRES法
対称行列(必ずしも正定値でない)のための残差最小化法. .  演算量をnに依存しない定数にできる.丸め誤差に弱く, x ( )中に( κ ( A ))²のオーダーで丸め誤差が混入* * H. A. van der Vorst:  “Iterative Krylov Methods for Large Linear Systems”, Cambridge Univ. Press, 2003.
原理
.
•  Aは対称正定値行列とする•残差r(n) = Ax(n)-b ∈ Kn+1(A; b)がKn(A; b)に直交するようx(n)を選ぶ
• Lanczos過程との関係
• Lanczos過程では,基底ベクトルq n +1 ∈ K n +1 (A; b)を, K n (A; b)に直交するよう選ぶ•したがって, Lanczos過程を実行し, q n+1を適当にスケーリングすれば, CG法の残差r(n)が得られるはずである•残差r (n)の計算法が決まれば,解x (n)の計算法も決まる(p. 14) • Lanczos過程における基底ベクトルのスケーリングq"+1は3項漸化式を満たすから,それをスケーリングしたr (n)も3項漸化式を満たすはずAr(n-1) = √n-2r(n-2) + δn-1r(n-1) + εn-1r(n) r(n)がr(n-1), r(n-2)と直交するという条件より,
•さらに,両辺の各項を{b, Ab, …, Anb}で展開したときのもの係数を考えると,左辺は0,右辺のr(n-1), r(n-2), r(n)はそれぞれ–1(p. 14)であるから, cn-1 = -(Yn-2 + δn-1) •これらによりr ⁄ ⁄ ⁄の漸化式が定まる
(37/60)
• CG法のアルゴリズム(3項漸化式) –  r(n)の漸化式にr(n) = Ax(n) – bを代入してx(n)の漸化式を導出[CG method (3-term recurrence)] 7-1 = 0, r(-1) = 0, r(0) = −b, x(-1) = x(0) = 0 for n = 1, 2, 3, . . . v = Ar(n-1) δn-1 =  r(n-1)T Ar(n-1) v = v — Yn-2r(n-2) — δn-1r(n-1) €n-1 = -(Yn-2 + δn-1) x(n) = (r(n-1) — 7n-2æ(n-2) — δn-1æ(n-1))/ɛn-1 r(n) = v/en-1 Yn-1 =  r(n-1)T Ar(n) end for
アルゴリズム（連立2項漸化式）
.
[CG algorithm] xo = 0,  ro = b,  po = ro for n = 1, 2, 3, . . . αn = (rπ−1rn−1)/(pπ−1Apn−1) Xn = Xn-1 + QnPn-1 rn = rn-1 — an Apn-1 βn = (rTrn)/(rT-1rn-1) Pn = rn + βnPn−1 end forステップ幅の決定探索方向に向かってx nを動かす残差の計算次の探索方向の決定最小化問題の解法としての解釈• 関数φ(x) = (1/2) x TAx – x Tbを最小化P n-1が第nステップでの探索方向・ p n − 1方向に進んで関数値が極小になるように,ステップ幅をα,を決定
• CG法の性質
• K,(A; b) = span{ b, Ab, A 2 b, …, A n-1b} = span{ x 1 , x 2 ,…, x n } = span{r 0, r 1,..., r n-1} = span{p 0, p 1,..., p n-1}残差の直交性: r n Tr j = 0   (j < n)探索方向のA-共役性: p n TAp j = 0   (j < n) # 最小化問題の解法と見たとき,各方向p,,への探索の繰り返しにより,大域的な最小値が求まることを保証
• CG法の特徴
計算の途中で破綻が起こることはない誤差のAノルム‖ e n ‖ a = ‖ x n − x * ‖ aがnとともに単調に減少所要メモリ量と1反復あたりの演算量はnに依存しない# 3項間漸化式により計算が行えるため
CG法の収束性
●
.誤差のAノルム‖ e, ‖ λ =の減少について,次の評価式が成り立つ
ただし, κはAの条件数
●
証明の概略
.残差多項式で表現した‖ e n ‖ aの最小性より, Li=1 a²λ j (φ(λ j ))² ||¢(A)eo|| A | En | A = intdePr - ≤ inf φ EPn   max  | φ ( λ ) | | E0 | A | E0 | A i=1 ažλjただし, P,は定数項が1のn次多項式, Λ(A)はAの固有値の集合ここで，特にφ ( λ )をチェビシェフ多項式に取って最右辺を評価すると上の定理が得られる
(41/60)
.
• CG法の収束性(続き)
・ Aの相異なる固有値の個数がm*の場合, CG法は高々m*回で収束する#  前ページのφ ( λ )として,これらm*個の固有値で0となるm*次多項式が取れるから, n = m*のとき,最右辺は0となる.固有値が少数のクラスタ一に固まっている場合も,収束は速い# これらのクラスタ一上で小さな値を取る低次の多項式が作れるから
前処理
●
•固有値を1個あるいは少数個のクラスターに集中させることができれば, CG法の収束は速くなる•この目的のため, Ax = bを別の(対称正定値な係数行列を持つ)連立1次方程式に変形して解く前処理(後述)が行われる
(42/60)
原理。 ・方程式[  O    A ] [ x ] _ [ b ]に対して形式的にCG法を適用する
_ _ _
。
アルゴリズム[Bi-CG algorithm] x0 = 0,  p0 = r0 = b,  q0 = s0 = arbitrary for n = 1, 2, 3, . . . αn = (s n −1 rn−1)/(g n−1 Apn−1) Xn = Xn-1 + QnPn-1 rn = rn-1 — anApn-1 Sn = Sn-1 — anA' qn-1 βn = (s Trn)/(s T-1 Tn-1) Pn = rn + BnPn-1 Qn = Sn + BnQn-1 end for
●
Breakdown
・ r nとs nとが直交すると,アルゴリズム中で分母が0になり,計算を続行できなくなる。これは，双直交基底を生成できないことに相当する。これを第1種のbreakdownと呼ぶー方, T,,のLU分解において,ピボットが0になることに相当するbreakdownもある。これを第2種のbreakdownと呼ぶこれらは.いずれもアルゴリズム上の工夫によって回避可能
Bi-CG法の特徴
。
.残差は必ずしもnととに単調に減少しない.所要メモリ量と1反復あたりの演算量はnに依存しない# GMRES法と比べた場合の大きな長所# 3項間漸化式により計算が行えるため・ AとA Tの両方による乗算が必要#反復当たりの行列ベクトル積の回数はGMRES法の2倍
(44/60)

• CGS法
•  Bi-CG法を変形し, r,, = (R,(A))²b, s,, = b*となるように工夫した手法# s,Tr,は不変なので,様々な係数はBi-CG法と同じに計算できる・残差多項式が2乗されるため,収束が速くなると予想される#ただし,収束の不規則性も2乗で効くため, Bi-CG法より不安定・ Aによる乗算が2回となり, A Tによる乗算は不要
Bi-CGSTAB法
●
. S,(z) = (1 − ω,z)(1 − ω,z) • • (1 − ω,−,z)という形のn次多項式S,(z) (安定化多項式)を用いて, r, = S,(A)R(A)b, s, = b*となるように工夫した手法#　この場合も,様々な係数はBi-CG法と同じに計算できる•各ステップで局所的に残差を最小化するようにω,を決めることで, CGS法よりも収束が安定化される
(45/60)

GPBi-CG法
●
•  Bi-CGSTAB法における安定化多項式S„(z)として, 3項漸化式により生成される多項式を用いた方法一般に, Bi-CGSTAB法よりもさらに優れた収束性を持つ.
QMR法
●
. Bi-CG法において,残差を擬似的に最小化するようにした方式.一般に,残差の収束はBi-CG法よりもスムーズ(振動が少ない) #収束は, Bi-CG法に比べて必ずしも速くない#しかし,残差の振動が少ないほうが丸め誤差の点で有利
前処理
●
.行列Kが，行列Aを何らかの意味で近似しているとするこのとき,方程式K-¹Ax = K-¹bを考えると,その係数行列K-¹AはAと比べて単位行列に近くなっていると考えられる。これにより，クリロフ部分空間法の収束はより速くなると期待できるこのような方程式の変形を前処理(左前処理)と呼ぶ.各反復では, Apの代わりにK-¹Apを計算する
右前処理,左前処理,
。
.右前処理:   AK-¹ y = b,  x = K-¹ y .両側前処理:   K 1 −1AK 2 −1y = K 1 −1b,  x = K 2 −1y 。
前処理行列Kに求められる性質
. K-¹Aがなるべく単位行列に近いこと. Kの生成が容易であること. K-¹xの計算が容易であること
(48/60)

フィルインのレベル
.
– 元々の非ゼロ要素a i lに対して,そのレベルをl(a i i ) = 0と定義する－フィルインによりa ij≠0となったとき, a ijのレベルを次式で定義するl(aij) = max(l(aik), l(akj)) + 1 	レベル0 □レベル1 ■レベル2 □レベル3 ■レベル4 □レベル5
(49/60)
•不完全LU分解
•行列AのLU分解の際に,フィルインを無視して近似的にA ≠ LUと分解し,  K = LUとおく。これをILU(0)前処理と呼ぶ#  1段階のみのフィルインを許すILU(1)前処理, 2段階のフィルイン(フィルインによって生じたフィルイン)までを許すILU(1)前処理も定義される•対称正定値行列の場合は.不完全コレスキー分解型前処理(IC(0)前処理)が使われる•幅広いクラスの行列に対して効果がある一般に並列化は困難. #　K-¹xの計算で,前進消去L-¹x，後退代入U-¹yが必要となるため並列化向きの不完全LU分解についても,多くの研究がある# •フィルインの絶対値の大きさによって無視するかどうかを決めるILUT前処理, ひ-1の列ノルムの大きさによって無視するかどうかを決めるILUC前処理もあり,悪条件の問題で有効
(50/60)
• 近似逆行列前処理(SPAI; SParse Approximate Inverse) ・行列Aそのものでなく, A-¹を直接近似する行列Mを求める•具体的には, Mの非零要素位置を適当に決めた上で,最小化問題min M ‖ I – AM ‖ Fを最小2乗法により解く•不完全LU分解に比べ,並列化が容易•  ILU分解と同程度の効果を出すには,一般により多くの非零が必要
対角スケーリング
.
. .対角行列D 1 , D 2によりA → D 1 -1 AD 2 -1と変換し,対角要素を1とする簡単で並列化も容易だが,問題によっては有効
行と列の並べ替え
●
.適当な並べ替えを行った後に不完全LU分解を行うと,収束性が改善することがある60
(51/60)
波束シミュレーションへの応用
•波束シミュレーションの基礎方程式–時間依存シュレディンガー方程式きゆ∂t = iH ψ – H:エルミート演算子(ハミルトニアン) –以下では簡単のため, Hがψ, tに依存しない場合を考える•空間離散化－差分法dx , = iHx,   H ∈ Rm×m,   x ∈ Rm以下ではこちらを扱うdt – Ritz-Galerkin法(有限要素法/基底関数による展開) □ B qx = iHx,  H,B ∈ R m×m,  x ∈ R m （B：重なり行列） dt
(53/60)

導出● −  (1,1)次パデ近似:  exp(iHΔt) ÷ (I − (1/2)iHΔt) − 1(I + (1/2)iHΔt )標準的な導出法:時刻t + (1/2)Δtでの中心差分
●特徴–時間発展のために連立1次方程式を解くことが必要－時間に関して2次精度－波動関数のノルムを厳密に保存1 + iek At人気が出来
•時間発展演算子の固有値の絶対値は| 1 — ick ∆ t
シュレディンガー方程式の時間発展の計算に最適
(54/60)

方程式.
係数行列の性質。 –大規模･疎行列･非エルミート–  ¿倍するとA + σ (A:正定値エルミート行列, σ:複素数)の形数値解法● – 直接法:ガウスの消去法(帯行列向け解法,スパースソルバなど) – 反復法: クリロフ部分空間法(GMRES法, BiCG法など)係数行列の構造を利用した効率的なアルゴリズムの設計を考える
(55/60)

•  Aのエルミート性の利用
–  GMRES法は, クリロフ部分空間K„(A; b) = span{b, Ab, …, A”−¹b} (n = 1, 2, …)の中で残差を最小にする解を探索– クリロフ部分空間のシフト不変性Kn(A+σΙ; b) = Kn(A; b)より, エルミ一ト行列Aに対するクリロフ部分空間が使える–ランチョス法と同様,直近の3本の基底ベクトルのみを用いてK„(A+σI; b)の正規直交基底を生成可能
得られる解法の特徴
.
3項間漸化式に基づく解法→演算量・メモリ量の負担軽減残差最小性を持つ→スムーズな収束性|
●
問題点
－前処理を行うと,シフト不変性が崩れてしまう－前処理なしで,どの程度の収束性が得られるかが鍵
(56/60)
。本発表では,まずクリロフ部分空間とその基底の生成法を説明し,方程式Ax = bをクリロフ部分空間に射影する3種類の方法を述べたこれに基づき, GMRES法, CG法, Bi-CG法の3つのアルゴ.リズムと,その変種を説明した●また,収束を加速させるための前処理法について,代表的な手法を説明した●実際の問題にクリロフ部分空間法を適用するには，係数行列の分類を検討した後, LisやPETScなどのパッケージを使って様々な解法･前処理を試してみるのがよい
(58/60)