5.3 連立一次方程式 (クリロフ部分空間法)

n × n 正則行列 A と非ゼロベクトル u の積により生成されるベクトルで張られる空間 \(K_{k+1}(A; u)\) をクリロフ部分空間とよぶ.
\[
K_{k+1}(A; u) = span(u, Au, … , A^k u)
\] クリロフ部分空間を生成して連立一次方程式の近似解を求める解法をクリロフ部分空間法という. 現在用いられている非定常反復解法の多くはクリロフ部分空間法に属する反復法である.

5.3.1 アーノルディ過程

\(K_{k+1}(A; u)\) の正規直交基底 \(V_{k+1} = { v_1, v_2, \dots, v_{k+1} }\) を \(u, Au, \dots, A^k u\) からグラム・シュミットの直交化により求めることを考える. この手順をアーノルディ過程と呼ぶ.

• \(v_1 = u / \|u\|\) とし, \(k = 1, 2, \dots\) について以下を繰り返す
  • \(h_{i,k} = (v_i, Av_k)\) とする \((i = 1, \dots, k)\) (ベクトル \(v_i\) と \(\tilde{v}_{k+1}\) の内積が 0 になる (直交する) ようにするため)
  • \(\tilde{v}_{k+1} = Av_k – \sum_{j=1}^k h_{j,k}v_j\)
  • \(v_{k+1} = \tilde{v}_{k+1} / h_{k+1,k}\) と正規化する. ただし, \(h_{k+1,k} = \|\tilde{v}_{k+1}\|\) とする

以上より, \(Av_i = \sum_{j=1}^{i+1} h_{j,i}v_j\) となるので, \(i = 1, \dots, k\) について行列にまとめて表すと次のようになる.
\[
\begin{align}
AV_k & = V_k H_k + h_{k+1,k}v_{k+1}e_k^T \\
& = V_{k+1}\tilde{H}_k \\
\end{align}
\] ここで, \(e_k = (0, 0, \dots, 0, 1)^T\) である.

k x k 行列 \(H_k\) は上ヘッセンベルグ形 (i ≥ j + 1 のとき (i ,j)要素が 0 になる行列) である.
\[
H_k =
\begin{pmatrix}
h_{1,1} & h_{1,2} & \cdots & h_{1,k-1} & h_{1,k} \\
h_{2,1} & h_{2,2} & \cdots & h_{2,k-1} & h_{2,k} \\
0 & h_{3,3} & \cdots & \cdot & \cdot \\
0 & 0 & \cdots & \cdot & \cdot \\
& & \ddots \\
0 & 0 & \cdots & h_{k,k-1} & h_{k,k} \\
\end{pmatrix}
\] 行列 \(\tilde{H}_k\) は \(H_k\) に 1 行加えた (k+1) x k 行列である.
\[
\tilde{H}_k =
\begin{pmatrix}
h_{1,1} & h_{1,2} & \cdots & h_{1,k-1} & h_{1,k} \\
h_{2,1} & h_{2,2} & \cdots & h_{2,k-1} & h_{2,k} \\
0 & h_{3,3} & \cdots & \cdot & \cdot \\
0 & 0 & \cdots & \cdot & \cdot \\
& & \ddots \\
0 & 0 & \cdots & h_{k,k-1} & h_{k,k} \\
0 & 0 & \cdots & 0 & h_{k+1,k} \\
\end{pmatrix}
\] なお, 行列 \(H_k\) は次のように表すことができる (上の \(AV_k = \cdots\) の式の左から \(V_k^T\) を掛ける).
\[
H_k = V_k^TAV_k
\]

5.3.2 ランチョス過程

行列 A が対称行列の場合, アーノルディ過程において \(h_{i,j} = (v_i, Av_j) = (Av_i, v_j) = 0 \space (i \ge j – 2)\) となり, \(H_k\) は三重対角行列になる.
この場合, \(V_k\) を 3 項漸化式から計算することができるようになる.
\[
v_{k+1} = (Av_k – h_{k-1,k}v_{k-1} – h_{k,k}v_k) / h_{k+1,k} \space (k = 1, 2, \dots )
\] ただし, \(h_{k+1,k} = \|Av_k – h_{k-1,k}v_{k-1} – h_{k,k}v_k\|\).

この場合を特に \(T_k = H_k, \alpha_k = h_{k,k}, \beta_k = h_{k-1,k} = h_{k,k-1}\) と書くことにすると次のようになる.
\[
T_k =
\begin{pmatrix}
\alpha_1 & \beta_2 \\
\beta_2 & \alpha_2 & \beta_3 & & & 0 \\
& \beta_3 & \alpha_3 & \beta_4 \\
& & & \ddots \\
0 & & & \beta_{k-1} & \alpha_{k-1} & \beta_k \\
& & & & \beta_k & \alpha_k \\
\end{pmatrix}
\] \(T_k\) に 1 行加えた (k+1) x k 行列 \(\tilde{T}_k\) は次のようになる.
\[
\tilde{T}_k =
\begin{pmatrix}
\alpha_1 & \beta_2 \\
\beta_2 & \alpha_2 & \beta_3 & & & 0 \\
& \beta_3 & \alpha_3 & \beta_4 \\
& & & \ddots \\
& & & \beta_{k-1} & \alpha_{k-1} & \beta_k \\
0 & & & & \beta_k & \alpha_k \\
& & & & & \beta_{k+1} \\
\end{pmatrix}
\] 対称行列用の 3 項漸化式を用いて \(V_k\) を求める以下の手順をランチョス過程とよぶ.

• \(\beta_1 = \|u\|, v_0 = 0, v_1 = u / \beta_1\) とし, 以下を \(k = 1, 2, \dots\) について繰り返す
  • \(\alpha_k = (v_k, Av_k)\) とする
  • \(v_{k+1} = Av_k – \alpha_kv_k – \beta_kv_{k-1}\)
  • \(\beta_{k+1} = \|v_{k+1}\|\) として, \(v_k+1 \gets v_{k+1}/\beta_{k+1}\) と正規化する

5.3.3 クリロフ部分空間法による連立一次方程式の解法

連立一次方程式 \(Ax = b\) において残差を \(r = b – Ax\) と表す. 初期近似解を \(x_0\) とするとき, 初期残差は \(r_0 = b – Ax_0\) である.

第 k ステップでの近似解 \(x_k\) を, それと初期近似解 \(x_0\) の差が係数行列 \(A\) と初期残差 \(r_0\) から形成されたクリロフ部分空間 \(K_k(A; r_0)\) に属するように選ぶ.
\[
x_k – x_0 \in K_k(A; r_0)
\] このとき, \(x_k\) の残差 \(r_k\) は \(K_{k+1}(A; r_0)\) に属する.
\[
\begin{align}
r_k & = b – Ax_k \\
& = b – Ax_0 – A(x_k – x_0) \\
& = r_0 – A(x_k – x_0) \in K_k+1(A; r_0) \\
\end{align}
\] クリロフ部分空間法では, まずクリロフ部分空間 \(K_k(A; r_0)\) の直交基底ベクトルをアーノルディ過程またはランチョス過程により求め, 次に \(x_k\) または \(r_k\) に一定の条件をつけることにより近似解を定める.

近似解を求めるために \(x_k\) または \(r_k\) につける条件には次のようにいくつかの方法がある.

5.3.4 対称行列のクリロフ部分空間法

クリロフ部分空間法において \(x_k – x_0\) は \(K_k(A; r_0)\) に属するから正規直交基底 \(v_1, v_2, \dots, v_k\) の線形結合で表すことができ, 近似解 \(x_k\) は次のように表される. ただし, \(y_k\) は k-ベクトルである.
\[
x_k = x_0 + V_ky_k
\] A が対称行列のときはランチョス過程により \(V_{k+1}\) と \(\tilde{T}_k\) を求めることができ \(x_k\) の残差 \(r_k\) は次のように表される.
\[
\begin{align}
r_k & = b – Ax_k \\
& = r_0 – A(x_k – x_0) \\
& = V_k+1(\beta_1 e_1 – \tilde{T}_ky_k) \\
\end{align}
\] ただし, \(e_1 = (1, 0, \dots, 0)^T\) である.

これより, \(\beta_1 e_1 – \tilde{T}_ky_k\) ができるだけ小さくなるように \(y_k\) を選べばよいが, \(y_k\) を一意に定めるためにはもうひとつ条件が必要でその定め方にはいくつかの方法があった. 対称行列の場合の主な解法と適用される条件は次のとおりである.

偏微分方程式を差分法で解く場合など多くの問題に現れる係数行列は対称でしかも正定値の場合が多いので, 共役勾配法 (CG 法) がよく使われる. 最小残差法 (MINRES 法) と共役残差法 (CR 法) は, 収束の速さは CG 法とほぼ同じで, 正定値でなくても使えるが, 計算量はやや多くなる. 対称 LQ 法 (SYMMLQ 法) は, 正定値でなくても使えるが, 他の方法に比べて収束がやや遅い傾向がある.

5.3.4.1 共役勾配 (CG) 法

(1) クリロフ部分空間法としての CG 法

リッツ・ガレルキン条件: \(r_k \perp K_k(A; r_0)\) より,
\[
\begin{align}
0 & = V_k^Tr_k \\
& = V_k^Tr_0 – V_k^TAV_ky_k \\
& = \beta_1e_1 – T_ky_k \\
\end{align}
\] これより次の連立一次方程式を解いて \(y_k\) を定めることができる.
\[
T_ky_k = \beta_1e_1
\] 求めた \(y_k\) を用いて \(x_k\) および \(r_k\) を更新する.

(2) CG 法の別の解釈

CG 法はクリロフ部分空間法の一種として上のように近似解を計算する反復式を得ることができるが, これに対して, 以下に説明する別の解釈による導出がよく知られており, 得られる反復式は等価である.

係数が正定値対称行列の連立一次方程式 Ax = b において, 次の関数 \(\phi(x)\) は x が連立一次方程式の解のときに最小値をとる.
\[
\phi(x) = \frac{1}{2}(x, Ax) – (x, b)
\] 従って, \(\phi(x)\) の最小点を求めれば方程式の解が求められることになる.

CG 法では, \(p_k\) を方向ベクトルとして, その方向上の \(\phi(x)\) の最小点を次のように探索して求める.
\[
x_{k+1} = x_k + \alpha_k p_k
\] \(\alpha_k\) は探索方向 \(p_k\) 上での修正量であり, \(\phi(x)\) が最小になるように定めると次のようになる.
\[
\alpha_k = \frac{(r_k, r_k)}{(p_k, Ap_k)}
\] このとき残差は次のようになり, 漸化式で計算できる.
\[
r_{k+1} = r_k – \alpha_kAp_k
\] 次に, 探索方向ベクトルを次のように定める. ただし, \(p_0 = r_0\) とする.
\[
p_{k+1} = r_{k+1} + \beta_k p_k
\] CG 法では \(\beta_k\) を \(p_{k+1}\) と \(Ap_k\) が直交する (\(p_{k+1}\) と \(p_k\) が A に関して共役になる) ように定める.
\[
\beta_k = -\frac{(r_{k+1}, Ap_k)}{(p_k, Ap_k)}
\] 以上を整理すると次の計算手順が得られる.

• \(x_0\) = 初期推定値, \(r_0 = b – Ax_0\), \(p_0 = r_0\) とする
• 以下を \(k = 0, 1, 2, \cdots\) について繰り返す
  • \(\alpha_k = \frac{(r_k, r_k)}{(p_k, Ap_k)}\)
  • \(x_{k+1} = x_k + \alpha_kp_k\)
  • \(r_{k+1} = r_k – \alpha_kAp_k\)
  • \(\|r_{k+1}\|/\|b\| \le\) 収束判定値 であれば終了する
  • \(\beta_k = \frac{(r_{k+1}, r_{k+1})}{(r_k, r_k)}\)
  • \(p_{k+1} = r_{k+1} + \beta_k p_k\)

ここで, \((p_k, r_k) = (r_k, r_k),\space (p_k, Ap_k) = (r_k, Ap_k)\) などの関係があるので

• \(\alpha_k = \frac{(p_k, r_k)}{(p_k, Ap_k)}\)
• \(\beta_k = -\frac{(r_{k+1}, Ap_k)}{(p_k, Ap_k)}\)

とも書ける.

なお, \(p_k\) および \(r_k\) について次式が成り立ち, クリロフ部分空間法になっている.
\[
span(p_0, p_1, \dots, p_k) = span(r_0, r_1, \dots, r_k) = K_{k+1}(A; r_0)
\]

プログラムの実装上は (1) よりもこちらの手順を使っているものが多いようである.

(3) CG 法の収束性

CG 法は反復法に分類されるが, 残差ベクトル \(r_k\) の直交性から n 回の反復で, 計算誤差がなければ, \(r_k = 0\) になる (解が求められる) という直接法的な性質も持つ.

CG 法の収束に関しては, 目的関数 \(\phi(x)\) について次の式が成り立つ (文献[3]).
\[
\phi(x_k) \le \phi(x_0) \cdot 4(\frac{\sqrt{\kappa} – 1}{\sqrt{\kappa} + 1})^{2k}
\] ここで, κ は係数行列の条件数で最大固有値と最小固有値の比である\((\kappa = \lambda_{max} / \lambda_{min})\).

すなわち, 係数行列の条件数が小さい (固有値が密集している) ほど CG 法は速く収束する. n 回よりも少ない回数で収束することも少なくない.

5.3.4.2 最小残差法 (MINRES 法)

(1) クリロフ部分空間法としての MINRES 法

CG 法は係数行列が正定値でない場合には適用できないが, その場合でも使える方法として最小残差法 (MINRES 法) が開発された.

最小残差条件: \(min \|r_k\|_2\) より, 最小二乗問題
\[
y_k = arg min \|\beta_1e_1 – \tilde{T}_ky_k\|_2
\] を解いて \(y_k\) を定める方法である.

三重対角行列 \(\tilde{T}_k\) は, 下副対角要素が 0 になるように (2 次元の) ハウスホルダー変換 (鏡映) を左から k 回施すことにより, 上三角行列 (上副対角要素が 2 つの帯行列) \(R_k\) にすることができる. k 回の変換をまとめて \(Q_k\) と表すと,
\[
Q_k\tilde{T}_k =
\begin{pmatrix}
R_k \\
0 \\
\end{pmatrix}
\] となる. また,
\[
Q_k\beta_1e_1 =
\begin{pmatrix}
f_k \\
\phi_k \\
\end{pmatrix}
\] と表すと,
\[
R_ky_k = f_k
\] より \(y_k\) を定めることができる.

求めた \(y_k\) を用いて \(x_k\) および \(r_k\) を更新する.

(2) 共役残差 (CR) 法

5.3.4.1 (2) の CG 法の手順を, 残差二乗ノルム \(\|r\|_2 = \|b – Ax\|_2\) を最小にする x を求めるように変更する.

CG 法では探索方向ベクトルを \(p_{k+1}\) と \(Ap_k\) が A に関して共役になるように定めた. これに対して, \(A^{T}A\) 共役になる, すなわち, \(Ap_{k+1}\) と \(Ap_k\) が直交するように定める. これにより, 残差二乗ノルムが効率よく減少する反復が得られる.

手順の変更点は次のとおりである.

• \(\alpha_k = \frac{(r_k, Ar_k)}{(Ap_k, Ap_k)}\)
• \(\beta_k = \frac{(r_{k+1}, Ar_{k+1})}{(r_k, Ar_k)}\)

このようにして得られる反復式は MINRES 法と等価であるが, こちらは 共役残差 (CR) 法とよばれる. 正定値対称行列に対しては MINRES 法と CR 法の結果は一致する. 不定値行列については MINRES 法のほうが有利とされる.

5.3.4.3 対称 LQ (SYMMLQ) 法

SYMMLQ 法も係数行列が正定値でない場合にも使える方法である.

最小誤差条件: \(min \|x_k – x\|_2\) より, 最小ノルム問題
\[
y_{k+1} = arg min \|y_{k+1}\| \space ただし \space \tilde{T}_k^Ty_{k+1} = \beta_1e_1
\] を解いて \(y_{k+1}\) を定める方法である.

今度は右から k 回の変換を施すことにより, \(\tilde{T}_k^T\) を下三角行列 (下副対角要素が 2 つの帯行列) \(L_k\) に変換する.
\[
\tilde{T}_k^TQ_k =
\begin{pmatrix}
L_k 0 \\
\end{pmatrix}
\] \(y_{k+1}\) を次のように置き換える.
\[
y_{k+1} = Q_k\tilde{z}_{k+1}
\] ただし,
\[
\tilde{z}_{k+1} =
\begin{pmatrix}
z_k \\
0 \\
\end{pmatrix}
\] そうすると, 次式が成り立つ.
\[
\begin{align}
\tilde{T}_k^Ty_{k+1} & = \tilde{T}_k^T Q_k \tilde{z}_{k+1} \\
& =
\begin{pmatrix}
L_k 0 \\
\end{pmatrix}
\tilde{z}_{k+1} \\
\end{align}
\] これより, 次式を解いて \(z_k\) を求めることができる.
\[
L_k z_k = \beta_1 e_1
\] \(x_k\) は \(V_{k+1}y_{k+1}\) より求められる.

5.3.4.4 対称行列のクリロフ部分空間法の比較

5.2.4 のポアソン方程式の 5 点差分近似を例として, CG 法, MINRES 法, CR 法, および, SYMMLQ 法を比較する.

数値実験 (3)

ここでは分点を多くとり, 31×31 分割の 900 元連立一次方程式 (非ゼロ要素数 = 4380) とした. 横軸は反復回数, 縦軸は残差ノルム (初期残差ノルムで正規化したもの) である. 収束条件は, 残差の相対誤差 < \(10^{-10}\) とした. 比較のために SOR法 (ω = 1.82 (最適値)) の結果も示した.

CG 法と MINRES 法はそれぞれ 64 回, 63 回の反復で収束した. CR 法は MINRES 法と全く同じ結果であった. 図では上書きされて MINRES 法のプロットが見えない. SYMMLQ 法はやや遅くて 70 回かかった. SOR 法は 127 回必要とした.

数値実験 (4)

CG 法は A が正定値対称行列という条件のもとに導かれた. これに対して MINRES 法, CR 法, および, SYMMLQ 法はこの条件なしに使えることを特長とする.

そこで, 係数行列 A が対称だが不定値な連立一次方程式を解いて比較する.

この例では 64 元連立一次方程式とし, 係数行列は要素がランダムな値の帯行列 (帯幅上下それぞれ 4) とした. ただし, その固有値の絶対値を 0.01, 0.02, 0.03, …, 0.64 とし, そのうちランダムに 30% の符号を負にして不定値行列になるようにした.

横軸は反復回数, 縦軸は残差ノルム (初期残差ノルムで正規化したもの) である. 収束条件は, 残差の相対誤差 < \(10^{-10}\) とした.

CG 法では正定値でないことが検出されても, 警告を出して計算を継続するようにした. これは, この例のように正定値でなくても CG 法は実際には収束することがあるためである.

図では N (A の行および列数) 回繰り返したあたりからどの解法も急速に収束しているように見える. CG 法は残差が増減し不安定ながらも収束した. MINRES 法と CR 法は残差が単調減少して安定に収束したが, 終盤で MINRES 法の方がやや早く抜け出している. SYMMLQ 法は MINRES 法および CR 法よりやや不安定に見える.

5.3.5 非対称行列のクリロフ部分空間法

非対称行列のクリロフ部分空間法はアーノルディ過程による解法とランチョス過程による解法の２種類に大別される.

非対称行列のクリロフ部分空間の正規直交基底は本来アーノルディ過程により求める必要がある. すなわち, 漸化式 (ランチョス過程) を使うことができず, 直交基底を陽に構成する必要がある. これは, 直交基底を保存しておくためメモリー量と計算量が反復毎に増えていくことを意味する. そのため, 実装上は一定回数ごとにリスタートする (初期化し直す) か, 保存する直交基底数を一定数に制限して古いものを破棄する (切り捨てる) かしてメモリー量と計算量が際限なく増大しないようにする必要がある. これに属する代表的な解法には最小残差 (GMRES) 法がある.

ところで, 正規直交基底の代用として双直交基底を使うことにすれば非対称行列であっても短い漸化式 (ランチョス過程) により直交基底を求めることができる. これを使用すればメモリー量と計算量が反復毎に増えることはなくアーノルディ過程による解法より有利である. ただし収束性についてよくわかっていない点もあり注意する必要がある. 双直交基底を使う解法では係数行列の転置行列に関する乗算や前処理計算が必要になる. そのため, 係数行列を陽に求めることができない場合には適用できない. これに属する代表的な解法は双共役勾配 (BICG) 法である.

ランチョス過程による解法であっても転置行列に関する計算が不要なように改善した解法も開発されている. 代表的な解法としては安定化双共役勾配 (BICGSTAB) 法があげられる.

いずれにせよ多くの解法が開発されているが決定的なものはなく問題の特性に応じて適切なものを選択する必要がある.

主要な解法を下表に整理する.

5.3.5.1 アーノルディ過程による解法

5.3.5.1.1 完全直交化法 (FOM (Full orthogonalization method))

クリロフ部分空間法において近似解 \(x_k\) は次のようにアーノルディ過程により求めた正規直交基底 \(V_k = (v_1, v_2, \dots, v_k)\) の線形結合で表された.
\[
x_k = x_0 + V_k y_k
\] ただし, \(y_k\) は k-ベクトルである.
\[
\begin{align}
r_k & = b – Ax_k \\
& = r_0 – AV_ky_k \\
\end{align}
\] より次が成り立つ.
\[
\begin{align}
V_k^Tr_k & = V_k^Tr_0 – V_k^TAV_ky_k \\
& = V_k^T(\|r_0\|v_1) – H_ky_k \\
& = \|r_0\|e_1 – H_ky_k \\
\end{align}
\] ただし, \(e_1 = (1, 0, \dots, 0)^T\) である.

ここで直交条件: \(r_k \perp K_k(A; r_0)\) を適用すると, \(V_k^Tr_k = 0\) となるから, 次の (k x k 密行列の) 連立一次方程式を解いて \(y_k\) を求めることができる.
\[
H_ky_k = \|r_0\| e_1
\] \(H_k\) は上ヘッセンベルグ行列であるからピボットの選択を行わないことにすれば LU 分解で容易に解くことができる. ただし, \(H_k\) の対角要素に 0 が現れるとブレークダウン (計算が継続できなくなること) が発生する.

このような解法を完全直交化法 (FOM) という.

FOM ではアーノルディ過程を用いるため, ヘッセンベルグ行列 \(H_k\) と正規直交基底 \(V_k = (v_1, v_2, \dots, v_k)\) を格納しておく必要がある. 反復ごとに k が大きくなっていくがメモリ量と計算量の観点から限界があるため, リスタートと呼ばれる手法が用いられる. すなわち, k の最大値 m を決めておき k = m になったところでそのときの近似解を初期近似解に設定し直してやり直す (\(x_0 = x_m, k = 0\) とする) ことにする. これは, リスタート付き FOM と呼ばれ FOM(m) 法と表すことがある. m をリスタートパラメータとよぶ.

もう一つの方法として, アーノルディ過程で \(v_{k+1}\) を計算するときに \(\sum_{j=1}^k\) とするところをあらかじめ値 m を定めておき最新の m 個の値だけを使うようにする, すなわち \(\sum_{j=k-m+1}^k\) とする方法が考えられる (切り捨て型). これにより得られる基底は m 個の間だけで直交する不完全なものであるが, 構わずに使って FOM と同様に計算を行う方法を不完全直交化法 (DIOM (Direct incomplete orthogonalization method)) という. この方法によるとヘッセンベルグ行列であった \(H_k\) が帯幅 m + 1 の帯行列になり, 解の更新には \(H_k\) および \(V_k\) の最新の m 列だけがあれば計算できるようになり, リスタートが不要になるメリットがある.

5.3.5.1.2 一般化最小残差 (GMRES (Generalized minimum residual)) 法

FOM の直交条件の代わりに最小残差条件: \(min \|r_k\|\) を適用する. すなわち, 最小二乗問題
\[
y_k = arg min y_k \|(\|r_0\| e_1 – \tilde{H_k}y_k)\|
\] を解いて \(y_k\) を定める方法を一般化最小残差法 (GMRES) という (ここでは \(H_k\) に 1 行加えた \(\tilde{H_k}\) を用いる).

上ヘッセンベルグ行列 \(\tilde{H_k}\) は, 下副対角要素が 0 になるようにギブンス変換を左から k 回施すことにより, 上三角行列 \(R_k\) にすることができる. k 回の変換をまとめて \(Q_k\) と表すと,
\[
Q_k\tilde{H_k} =
\begin{pmatrix}
R_k \\
0 \\
\end{pmatrix}
\] となる. また,
\[
Q_k \|r_0\| e_1 =
\begin{pmatrix}
f_k \\
\phi_k \\
\end{pmatrix}
\] と表すと,
\[
R_k y_k = f_k
\] より \(y_k\) を定めることができる. \(\phi_k\) は残差となる.

GMRES 法もアーノルディ過程を用いるため FOM と同様にリスタートが必要で, 実用的な解法はリスタート付き GMRES 法と呼ばれ GMRES(m) 法と表すことがある.

GMRES 法に前処理行列 M によって前処理を行うことを考える. 単純にこれを行うと \(M^{-1} V_k\) を何度か計算する必要がある. そこで, アーノルディ過程の中で計算した \(M^{-1} V_k\) を保存しておき計算量を減らす方法が考えられた. これを Flexible GMRES (FGMRES) 法という. FGMRES 法は GMRES 法と同じ結果を与え, 計算量が減る代わりに必要メモリー量が増える.

FOM に対する DIOM と同様に, GMRES 法に対して最新の m 個の値だけを使うようにする方法 (切り捨て型) を Direct quasi-minimum residual (DQGMRES) 法と呼び, リスタートが不要になるメリットがある.

5.3.5.1.3 一般化共役残差 (CGR (Generalized conjugate residual)) 法

GMRES 法と同様にクリロフ部分空間において最小残差条件を適用する. ただし, GMRES 法ではアーノルディ過程により直交基底を生成したが, GCR 法では \(W = A^TA\) に関する直交基底を用いる. すなわち, アーノルディ過程でユークリッドノルム \(\|x\|\) の代わりに \(\|x\|_W = \|x\|_{A^TA} = (x, x)_{A^TA} = (Ax, Ax)\) を用いる (GMRES 法は \(W = I\) の場合にあたると解釈できる).

残差は次のように表される.
\[
\begin{align}
r_k & = r_0 – AV_ky_k \\
& = r_0 – \sum_{i=1}^k y_k(i)Av_i
\end{align}
\] ただし, \(y_k(i)\) はベクトル \(y_k\) の i 番目の成分を表す.

最小残差条件: \(min \|r_k\|_2\) を適用すると, \((Av_1, Av_2, \dots, Av_k)\) が正規直交系であるため次式が成り立つ.
\[
y_k(i) = (Av_i, r_0)
\] 従って, \(y_k(i)\) は k に依存しないことがわかるので \(y(i)\) と書き直すと残差は次のように表される.
\[
r_k = r_0 – \sum_{i=1}^k y(i)Av_i
\] これより次の漸化式が得られる.
\[
r_k = r_{k-1} – y(k)Av_k
\] GMRES 法では近似解を求めるのに QR 分解を用いて最小二乗問題を解いたが, GCR 法では上の漸化式を用いて近似解を求めることができるので計算が簡単になる.

GCR 法でも GMRES 法と同様にベクトル列を保存しておく必要があるため実用上はリスタートが必要である. リスタート付きの GCR 法は GCR(m) 法とよばれる.

リスタートを行わずに最新の m 個の値だけを使う切り捨て型の方法もあり, これは Orthomin(m) 法とよばれる.

5.3.5.2 ランチョス過程による解法

5.3.5.2.1 非対称行列のランチョス過程

行列 A が非対称の場合であっても, A に関するクリロフ部分空間 \(K_{k+1}(A; u)\) と \(A^T\) に関するもうひとつのクリロフ部分空間 \(K_{k+1}(A^T; \tilde{u})\) の双直交系 \(V_{k+1} = (v_1, v_2, \dots, v_{k+1})\) と \(W_{k+1} = (w_1, w_2, \dots, w_{k+1})\) のペアを漸化式により求めることができる.
\[
\begin{align}
& AV_k = V_kT_k + \gamma_{k+1} v_k e_k^T = V_{k+1}\tilde{T_k} \\
& A^T W_k = W_k T_k^T + \beta_{k+1} w_k e_k^T = W_{k+1}\tilde{T_k}^T \\
\end{align}
\] ただし, \(e_k = (0, 0, \dots, 0, 1)^T\) である.

漸化式の係数を表す k x k 行列 \(T_k\) は次のとおりである.
\[
T_k =
\begin{pmatrix}
\alpha_1 & \beta_2 \\
\gamma_2 & \alpha_2 & \beta_3 & & 0 \\
& \gamma_3 & \alpha_3 & \beta_4 \\
& & \ddots & \ddots & \ddots \\
& 0 & & \gamma_{k-1} & \alpha_{k-1} & \beta_k \\
& & & & \gamma_k & \alpha_k \\
\end{pmatrix}
\] \(\tilde{T_k}\) はこれに 1 行加えた (k+1) x k 行列で次のようになる.
\[
\tilde{T_k} =
\begin{pmatrix}
\alpha_1 & \beta_2 \\
\gamma_2 & \alpha_2 & \beta_3 & & 0 \\
& \gamma_3 & \alpha_3 & \beta_4 \\
& & \ddots & \ddots & \ddots \\
& & & \gamma_{k-1} & \alpha_{k-1} & \beta_k \\
& 0 & & & \gamma_k & \alpha_k \\
& & & & & \gamma_{k+1} \\
\end{pmatrix}
\] また, \(\tilde{T_k}^T\) は \(T_k^T\) に \(\beta_{k+1}e_k^T\) の 1 行を加えた (k+1) x k 行列を表す.

2 つのクリロフ部分空間は双直交するので次式が成り立つ.
\[
V_k^TW_k = I
\] 非対称行列 \(A\) と \(A^T\) に関する双直交するクリロフ部分空間の双直交基底のペアを求める以下の手順を両側ランチョス過程という.

• \(u\) と \(\tilde{u}\) について \((u, \tilde{u}) \neq 0\) とする
• \(\beta_1 = \gamma_1 = 0, v_0 = w_0 = 0, v_1 = u/\|u\|, w_1 = \tilde{u}/(\tilde{u}, v_1)\) とし, 以下を \(k = 1, 2, \dots\) について繰り返す
  • \(\alpha_k = (w_k, Av_k)\)
  • \(v_{k+1} = Av_k – \alpha_kv_k – \beta_kv_{k-1}\)
  • \(w_{k+1} = A^Tw_k – \alpha_kw_k – \gamma_kw_{k-1}\)
  • \(\gamma_{k+1} = \|v_{k+1}\|\) として, \(v_{k+1} := v_{k+1}/\gamma_{k+1}\) と正規化する
  • \(\beta_{k+1} = (v_{k+1}, w_{k+1})\) として, \(w_{k+1} := w_{k+1}/\beta_{k+1}\) と正規化する

この手順中, \(\gamma_{k+1} = 0\) になったときには解が求められているので停止する. \(\beta_{k+1} = 0\) の場合, 割り算ができなくなるので計算が継続できなくなる (ブレークダウンが発生する). その場合, \(w_{k+1} = 0\) ならば初期値を変えて計算し直す必要がある.

5.3.5.2.2 双共役勾配 (BICG (Bi-conjugate gradient)) 法

FOM や対称行列の CG 法では \(r_k\) がクリロフ部分空間 \(K_k(A^T; r_0)\) と直交するという条件を付けた. 代わりに, BICG 法では非対称行列の 2 つの双直交するクリロフ部分空間 \(K_k(A; r_0)\) および \(K_k(A^T; \tilde{r_0})\) について次のような双直交条件を付ける.
\[
r_k \perp K_k(A^T; \tilde{r_0})
\] これより, 次の関係が得られる.
\[
\begin{align}
0 & = W_k^Tr_k \\
& = W_k^Tr_0 – W_k^TAV_ky_k \\
& = \|r_0\| e_1 – T_ky_k \\
\end{align}
\] これより次の連立一次方程式を解いて \(y_k\) を定め, それを用いて \(x_k\) の近似解を求めることができる.
\[
\begin{align}
& T_ky_k = \|r_0\|e_1 \\
& x_k = x_0 + V_ky_k \\
\end{align}
\] BICG 法では三重対角行列 \(T_k\) を LU 分解して \(y_k\) を求め, それを用いて \(x_k\) および \(r_k\) を更新する. また, 双対系の残差ベクトル \(\tilde{r_k}\) も同時に更新する.

BICG 法では \(A^T\) に関する計算, すなわち, \(A^Tx\) (乗算) と \(M^Tx = b\) の解 (前処理) が必要になる. 係数行列 \(A\) が陽に与えられている場合には問題ないが, そうでない特別の場合には適用が難しいことがあるかもしれない.

\(A\) が対称行列 (\(A = A^T\)）で \(\tilde{r_0} = r_0\) のとき, BICG 法は CG 法に帰着する.

BICG 法はクリロフ部分空間法の直交条件が満たされているわけではないので理論的な収束性はよくわかっていないが, 行列の種類によっては GMRES 法と同等の反復回数で収束するとされる.

5.3.5.2.3 積型反復解法

\(\tilde{r_0}, \tilde{r_1}, \dots, \tilde{r_k}\) を BICG 法の残差列とする. \(H_0(A), H_1(A), \dots, H_k(A)\) を多項式列とする. ただし, \(H_k(A)\) は k 次多項式を表す.

このとき, 残差が \(H_0(A)\tilde{r_0}, H_1(A)\tilde{r_1}, \dots, H_k(A)\tilde{r_k}\) となるようにして収束の加速を図ることを考える. このような解法を積型反復解法という.

BICG 法では残差を \(\tilde{r_k} = R_k(A)r_0\) と表すことができる. ただし, \(R_k(A)\) は k 次多項式である. そうすると, 積型反復解法の残差 \(r_k\) は次のように表される.
\[
\begin{align}
r_k & = H_k(A)R_k(A)r_0 \\
& = b – Ax_k \\
\end{align}
\] このとき次の関係が成り立つ.
\[
x_{k+1} – x_k = A^{-1}(H_{k+1}(A)R_{k+1}(A) – H_k(A)R_k(A))r_0
\] この漸化式では \(H_k(0)R_k(0) = 1\) であるとする. BICG 法では \(R_k(0) = 1\) であるから, この条件は \(H_k(0) = 1\) に帰着される.

積型反復解法では BICG 法と異なり \(A^T\) に関する計算 (乗算と前処理の計算) は不要である.

(1) 二乗共役勾配 (CGS (Conjugate gradient squared)) 法

BICG 法では \(r_0\) に \(R_k(A)\) を掛けることにより残差は減少していくのだから \(R_k(A)^2r_0\) とする, すなわち, 積型反復解法で \(H_k(A) = R_k(A)\) とすると, 2 倍の速さで収束しそうである. これを CGS 法とよぶ. 実際には x に対する反復ごとの補正が大きすぎると計算が不安定になったり発散することがあるためいつも 2 倍の速さで収束するわけではない.

(2) 安定化双共役勾配 (BICGSTAB (Bi-conjugate gradient stabilized)) 法

CGS 法の収束は不安定になることがあるため, BICGSTAB 法ではそれを改善するために \(H_k(A) = Q_k(A)\) とする. \(Q_k(A)\) は k 次の安定化多項式と呼ばれ次の漸化式で表される.
\[
Q_{k+1}(A) = (1 – \omega_kA)Q_k(A), \space Q_0(A) = 1
\] スカラーパラメータ \(\omega_k\) は残差ノルム \(\|r_{k+1}\|\) を最小にするように決められる.

(3) 積型双共役勾配 (GPBICG (Generalized product bi-conjugate gradient)) 法

GPBICG 法では安定化多項式 \(H_k(A)\) を次の漸化式で与える.
\[
\begin{align}
& H_0(A) = 1, \space G_0(A) = \omega_0 \\
& H_{k+1}(A) = H_k(A) – AG_k(A) \\
& G_{k+1}(A) = \omega_{k+1}H_{k+1}(A) + \eta_{k+1}G_k(A) \\
\end{align}
\] スカラーパラメータ \(\omega_k\), \(\eta_k\) は残差ノルム \(\|r_{k+1}\|\) を最小にするように決められる. \(\eta_k = 0\) とすると BICGSTAB 法に帰着する.

(4) BICGSTAB2 法

偶数回目の反復で BICGSTAB 法のパラメータ (すなわち, \(\eta_k = 0\)) を適用し, 奇数回目の反復では GPBICG 法のパラメータを使う方法である.

5.3.5.2.4 疑似最小残差 (QMR (Quasi minimum residual)) 法

\(r_k = r_0 – AV_ky_k\) より次式が得られる.
\[
r_k = V_k+1(\|r_0\|e_1 – \tilde{T_k}y_k)
\] 従って, \(r_k\) のノルムは次式を満たす.
\[
\|r_k\| \le \|V_{k+1}\| \|(\|r_0\|) e_1 – \tilde{T_k}y_k\|
\] QMR 法ではこれの第 2 項を最小化するように次の最小二乗問題を解いて \(y_k\) を定める.
\[
y_k = arg min \|(\|r_0\|)e_1 – \tilde{T_k}y_k\|
\] 最小二乗解は \(\tilde{T_k}\) をギブンス変換を施して QR 分解することにより求められる.

QMR 法については理論的な収束性が解析されており, GMRES 法に匹敵することがわかっている.

BICG 法では三重対角行列 \(T_k\) が特異になるとアルゴリズムが停止するが, QMR 法では最小二乗問題を解くため解を求めることができる. さらに, ランチョス過程のブレークダウンに対応するために先読み (Look-Ahead) 技法を組み込んだ実装もある. 先読み技法を組み込んでいなくても, QMR 法は BICG 法よりも安定に収束する傾向がある.

なお, もともとの QMR 法はランチョス過程の 3 項漸化式を使うが, 等価な 2 項漸化式のペアを使う方法もあり, 少し精度がよくなるとされる.

QMR 法では BICG 法と同様に \(A^T\) に関する計算, すなわち, \(A^Tx\) (乗算) と \(M^Tx = b\) の解 (前処理) が必要であるが, それを不要にしたのが Transpose free QMR (TFQMR) 法で, CGS 法における双直交条件を疑似最小残差条件に変更した解法である.

5.3.5.3 非対称行列のクリロフ部分空間法の比較

5.3.5.3.1 ポアソン方程式の 5 点差分近似 (非対称)

対称行列の数値実験で使用したポアソン方程式の 5 点差分近似の例を少し変形して非対称行列の係数にする.

(1) ポアソン方程式の例 (非対称)

5.2.1.4 の対称行列用の例題を少し変形して, 次のような 2 次元正方形領域 (0 ≤ x ≤ 1, 0 ≤ y ≤ 1) におけるポアソン方程式を考える.
\[
-\Delta u – \partial u/\partial x = f
\] ここで, \(\Delta u = \nabla \cdot(\nabla u) = \partial^2u/\partial x^2 + \partial^2u/\partial y^2\) である.

境界条件は次のとおりとする.
\[
\begin{align}
& u(x, 0) = 0, \space u(x, 1) = 0 \space (0 \le x \le 1) \\
& u(0, y) = 0, \space u(1, y) = 0 \space (0 \le y \le 1) \\
\end{align}
\] 追加された左辺第 2 項は移流項とよばれる.

(2) 5 点差分近似による離散化

移流項を中心差分を用いて次のように表す.
\[
-\partial u/\partial x = \frac{u(x-h, y) – u(x+h, y)}{2h} = \frac{u_{i-1,j} – u_{i+1,j}}{2h}
\] これより次の差分方程式が得られる.
\[
-u_{i-1,j} – (1 – \frac{h}{2}u_{i,j-1} + 4u_{i,j} – (1 + \frac{h}{2}u_{i+1,j} – u_{i,j+1} = h^2 f_{i,j}
\] こうして得られた連立一次方程式の係数行列は上下の副対角要素の値が異なる非対称行列になる.

(3) 数値実験 (5) 差分方程式を解く

非対称行列用のポアソン方程式の例題の差分方程式を解く. ここでは N = 31 (900 元連立一次方程式, 非ゼロ要素数 = 4380) とし, 解が (1, 1, …, 1) となるように右辺を設定した.

比較する解法は FOM, GMRES, GCR, BICG, CGS, BICGSTAB, GPBICG, BICGSTAB2, QMR, TFQMR である. 初期値を (0, 0, …, 0) としたときの収束の様子を次に示す. 横軸は反復回数, 縦軸は残差ノルム (右辺ベクトルのノルムで正規化したもの) である. GMRES, FOM, GCR ではリスタートや切り捨てが起こらないようにパラメータ M を十分大きくとった.

この例では GMRES と GCR は同じ結果になったので図では重なって GCR しか見えない.

アーノルディ過程を用いた GMRES は非対称行列用の基本解法といえるが, 安定して収束しているのがわかる. 同じくアーノルディ過程を用いた FOM と GCR も同様である.

BICG はこれらよりも計算が容易なランチョス過程を用いながらも同程度の収束性を狙ったものであるが, 狙いどおりおおむね GMRES に匹敵する収束性を示している. ただし, BICG はやや不安定な挙動を示している (残差が単調減少していない). それに対して, QMR は BICG の安定化を狙ったものであるが, BICG と収束の速さが同等でかつより安定して収束しているのがわかる.

以上のアーノルディ過程を用いた解法 (GMRES, FOM, GCR) とランチョス過程を用いながらも同等の収束性を狙った BICG および QMR のグループを基本の収束速度とすると, 積型反復法はそれよりも速い2つのグループにはっきりと分かれた.

積型反復法の BICG に対する特長は, 転置に関する計算が不要なことの他に速く収束することが期待できることである. CGS は BICG よりも速いがかなり不安定なことがわかる. その安定化を狙った BICGSTAB, GPBICG, BICGSTAB2 は安定性が増しているだけでなく CGS よりさらに速いグループを形成している. TFQMR と CGS の関係は QMR と BICG の関係に同じなので, TFQMR は CGS とほぼ同じ収束速度だが安定化されているのがわかる.

この例題は比較的解きやすい問題であるためそれぞれの解法の特徴をよく表していると思われるが, 実際の問題では問題によって適切な解法が異なる.

5.3.5.3.2 GMRES におけるリスタート・パラメータ

数値実験 (6) リスタートあるいは切り捨ての影響

非対称行列用のポアソン方程式の例題 (N = 31) の差分方程式を解く.

GMRES (リスタート) および DQGMRES (切り捨て) においてパラメータ m を 5, 10, 15 と変化させて, リスタートあるいは切り捨てなしの場合と比較する.

ここで, 無印の GMRES と DQGMRES (リスタートあるいは切り捨てなし (m を十分に大きくとった)) は同じ結果を与えるので図では重なって DQGMRES しか表示されていない. 反復回数で表示したこの図では DQGMRES の方が有利に見えるが, 1 回の反復当たりの計算量は GMRES の方が少なくなるので実際のところは計算速度を計測してみる必要がある. DQGMRES では誤差が蓄積されるのか m が小さいと計算精度の限界が現れるようである.

5.3.6 前処理付きクリロフ部分空間法

CG 法の収束定理からわかるように, 条件数が小さくなるように方程式を変換してから CG 法を適用できれば収束を速くできることが期待される. これを前処理といい, CG 法に限らずクリロフ部分空間法を適用する際には一般的な方法となっている. 実際に, 前処理をしない場合に比べて大幅に収束が速くなることが多い.

n × n 行列 M を行列 A の前処理行列とする. M は次のように行列 A を近似できるように選ぶ.
\[
M = M_1 M_2 \simeq A
\] 次式が成り立つ.
\[
M_1^{-1} A M_2^{-1} \simeq I
\] これを使って, 次のように係数を変換する.
\[
\begin{align}
& \tilde{A} = M_1^{-1} A M_2^{-1} \\
& \tilde{x} = M_2 x \\
& \tilde{b} = M_1^{-1} b \\
\end{align}
\] こうして得られた方程式は, 係数行列 A がより単位行列 I に近づいた解きやすい方程式になっている.
\[
\tilde{A} \tilde{x} = \tilde{b} \space (\tilde{A} \simeq I)
\] このようにして変換された方程式にクリロフ部分空間法を適用することを前処理付きクリロフ部分空間法という.

なお, \(M_1 \ne I\) かつ \(M_2 \ne I\) の場合を分離前処理 (両側前処理), \(M_1 = I\) の場合を右前処理, \(M_2 = I\) の場合を左前処理という.
\[
\begin{align}
& AM^{-1}(Mx) = b \space (右前処理) \\
& M^{-1}Ax = M^{-1}b \space (左前処理) \\
\end{align}
\] 分離前処理, 右前処理、左前処理のうちどの前処理を使用しても同じ前処理行列であれば収束性は大きくは変わらない.

A が対称行列の場合には \(\tilde{A}\) も対称になるように \(M_1 = M_2^T\) となるような前処理行列を使った分離前処理を使用することが多い. A が非対称行列の場合には残差の大きさが変わらない右前処理を使うことが多い.

5.3.6.1 前処理付き CG 法

対称な前処理行列 \(M\) について \(M_1 = M_2 = M^{1/2}\) と表すことにする. これを用いて \(\tilde{A}\), \(\tilde{b}\) および \(\tilde{x}\) を次のように定義する.
\[
\begin{align}
& \tilde{A} = M^{-1/2} A M^{-1/2} \\
& \tilde{x} = M^{1/2} x \\
& \tilde{b} = M^{-1/2} b \\
\end{align}
\] このようにすると, \(Ax = b\) と等価な方程式 \(\tilde{A} \tilde{x} = \tilde{b}\) が得られる. これに通常と同じ解法を適用すればよい.

CG 法の場合にはこれを使って計算手順を

• \(\alpha_k = \frac{(\tilde{r}_k, \tilde{r}_k)}{(\tilde{p}_k, \tilde{A}\tilde{p}_k)}\)
• \(\tilde{x}_{k+1} = \tilde{x}_k + \alpha_k\tilde{p}_k\)
• \(\cdots\)

などと書くことができるが, \(\tilde{r}_k = M^{-1/2}r_k, \tilde{p}_k = M^{1/2}p_k\) などの関係を使って書き換えると次の手順が得られる. \(M^{1/2}\) の項が消えて \(M^{-1}y\) を求めるという計算 (= 方程式 \(Mx = y\) を解いて \(x\) を求める) を加えただけになる.

前処理付きCG法の計算手順

• \(x_0\) = 初期推定値, \(r_0 = b – Ax_0\), \(p_0 = M^{-1}r_0\) とする
• 以下を \(k = 0, 1, 2, \dots\) について繰り返す
  • \(\alpha_k = (M^{-1}r_k, r_k)/(p_k, Ap_k)\)
  • \(x_{k+1} = x_k + \alpha_kp_k\)
  • \(r_{k+1} = r_k – \alpha_kAp_k\)
  • \(\|r_{k+1}\|/\|b\| \le\) 収束判定値 であれば終了する
  • \(\beta_k = (M^{-1}r_{k+1}, r_{k+1})/(M^{-1}r_k, r_k)\)
  • \(p_{k+1} = M^{-1}r_{k+1} + \beta_k p_k\)

他のクリロフ部分空間法の場合も同様にして前処理付きにすることができる.

前処理には種々の方法が考えられるが, 収束が速くなっても前処理の計算に時間がかかっては効果がないため, ある程度シンプルなものが求められる.

5.3.7 対称行列用の前処理

対称行列用の代表的な前処理としては不完全分解型の不完全コレスキー分解と行列分離型の SSOR があげられる. 他に点ヤコビ前処理が有効なこともある.

5.3.7.1 不完全コレスキー分解

CG 法の前処理としては不完全コレスキー分解 (incomplete Cholesky decomposition) がよく使われる. 不完全コレスキー分解による前処理付きの CG 法を ICCG 法とよぶ.

正定値対称行列を係数とする連立一次方程式の直接解法ではコレスキー分解 \(A = LDL^T\) が用いられる. これをそのまま疎行列に適用すると, もともとゼロ要素だったところがゼロでなくなり分解結果は非ゼロ要素数が増えたものになってしまい (フィルインが生じて), 反復法の特長が失われてしまう. これを避けるために, 特定部分だけ計算して残りを強制的にゼロにして非ゼロ要素を極力増やさないようにするのが不完全コレスキー分解である.

このような不完全な分解では, 正しい解を求めることができなくても, 前処理の効果が期待できる.

不完全コレスキー分解を式で表すと次のようになる.
\[
A = LDL^T + R \space (R \simeq 0)
\] ここで, R は強制的に 0 にしてしまった部分に起因する A との差分である. この分解 \(LDL^T\) は完全ではないものの, R が小さければかなり本物のコレスキー分解に近いことが期待できる.

前処理行列 M として不完全コレスキー分解を使用することにする.
\[
M = M_1M_2 = (LD^{1/2})(LD^{1/2})^T
\] そうすると, \(M_1^{-1}AM_2^{-1}\) が単位行列に近いものになり収束が速くなることが期待できる.

どの部分を 0 にしてしまうかはいろいろな方法が考えられ, それにより ICCG 法もいくつかのバリエーションがある. 適用する問題 (対象とする行列の構造) ごとに適切なやり方が工夫されている. 以下に代表的な例を示す.

(1) フィルインなしの方法 (IC(0))

最も素朴なやり方として, もともと 0 だったところは 0 のままにしてフィルインをなくす方法が考えられる. 行列の構造に依存せず汎用的に使える方法である.

(2) 規則格子における ICCG 法

規則的な長方形格子上で構成した差分法の場合, 行列は規則的な構造を持つ. 2次元で m + 1 等分した格子の場合, 対角要素 \(a_{i,i}\) の他には隣 \(a_{i,i-1}, a_{i,i+1}\) と対角要素から m 離れた \(a_{i,i-m}, a_{i,i+m}\) だけが非ゼロ要素となる. 定常反復法の例題として取り上げたポアソン方程式の 5 点差分近似の係数行列はその一例である.

この構造に着目して, コレスキー分解の対象を \(a_{i,j} (j = i, i\pm 1, i\pm m)\) として他は 0 としたものを ICCG(1, 1) 法という.

ICCG(1, 1) 法に対してグスタフソンの修正 \((x_{i-m+1}, x_{i+m-1}\) の項を対角要素に修正して精度を上げる方法) を加えたものを MICCG(1, 1) 法 (Modified ICCG method) といい, ICCG(1, 1) 法より収束が速い.

さらに, 1 つ内側の要素も対象, すなわち, \(a_{i,j} (j = i, i\pm 1, i\pm (m-1), i\pm m)\) を対象とする方法を ICCG(1, 2) 法という. これにグスタフソンの修正を加えたものが MICCG(1, 2) 法である.

5.3.7.2 SSOR (Symmetric successive-over-relaxation) 前処理

定常反復法で用いられる反復行列 M を前処理行列として用いる方法を行列分離型前処理という. その中でも対称行列に対する SOR 法の反復行列を用いるのが SSOR 前処理である.

対称な係数行列 A を次のように分解する.
\[
A = L + D + L^T
\] ここで, \(L\) は対角要素が 0 の下三角行列, \(D\) は対角行列, \(L^T\) は対角要素が 0 の上三角行列である. これはコレスキー分解とは違って A を 3 つの部分に分けただけ (足し算) なので簡単に求められる. これを使って前処理行列を次のように定義する.
\[
M(\omega) = \frac{1}{2 – \omega}(\frac{D}{\omega} + L)(\frac{D}{\omega})^{-1}(\frac{D}{\omega} + L)^T
\] この方法は計算が簡単であるが効果が期待できる. ただし, ω は SOR 法と同じように加速パラメータであるが最適な値を前もって知ることは難しいため経験的に決めることになる.

5.3.7.3 点ヤコビ (Point Jacobi) 前処理

前処理行列 M として係数行列の対角行列 (対角要素だけからなる正方行列) D を用いる.
\[
M = D
\] 最も簡単な前処理であるが, 問題によっては効果がある.

5.3.7.4 対称行列用の前処理の比較

数値実験 (7)

再びポアソン方程式の 5 点差分近似を例として, CG法, ICCG法 (フィルインなし(IC(0))), ICCG(1, 1) 法, ICCG(1, 2) 法, MICCG(1, 1) 法, MICCG(1, 2) 法, および SSOR 法 (ω = 1.7) により比較する.

31×31 分割の 900 元連立一次方程式 (非ゼロ要素数 = 4380) を使用した. 横軸は反復回数, 縦軸は残差ノルム (右辺ベクトルのノルムで正規化したもの) である. 収束条件は, 残差の相対誤差 < \(10^{-10}\) とした. なお, この例では係数行列の対角要素がそろっているので点ヤコビ前処理は効果がない.

前処理を行うことにより, 反復回数を尺度とした場合, 倍以上速くなっている. この場合, IC(0) と ICCG(1, 1) は同じ結果になる.

同じ前処理は MINRES 法および SYMMLQ 法にも適用でき, CG 法と同様の傾向を示した.

数値実験 (8)

次に, 同じ例題で SSOR 前処理付き CG 法において, ω を変化させたときにどうなるか確認する.

この例では ω の値にそれほど敏感ではない.

5.3.8 非対称行列用の前処理

非対称行列用の代表的な前処理としては不完全分解型の不完全 LU 分解, 行列分離型の SSOR, および点ヤコビ前処理があげられる.

5.3.8.1 不完全 LU 分解 (ILU) による前処理

非対称行列を係数とする連立一次方程式の直接解法では LU 分解 A = LU が用いられる. これをそのまま疎行列に適用すると, もともとゼロ要素だったところがゼロでなくなり分解結果は非ゼロ要素数が増えたものになってしまう (フィルインが生じる). これを避けるために, 特定部分だけ計算して残りを強制的にゼロにして非ゼロ要素を極力増やさないようにするのが不完全 LU 分解 (ILU 分解)である.

ILU 分解は対称行列用の不完全コレスキー分解に相当する非対称行列用の前処理である.

ILU 分解を式で表すと次のようになる.
\[
A = LU + R \space (R \simeq 0)
\] ここで, R は強制的に 0 にしてしまった部分に起因する A との差分である. この LU 分解は完全ではないものの, R が小さければかなり本物の LU 分解に近いことが期待できる.

前処理行列 M として ILU 分解を使用することにする.
\[
M = LU
\]

どの部分を 0 にしてしまうかはいろいろな方法が考えられ, いくつかのバリエーションがある. 適用する問題 (対象とする行列の構造) ごとに適切なやり方が工夫されている. 以下に代表的な例を示す.

(1) フィルインなしの方法 (ILU(0))

最も素朴なやり方として, もともと 0 だったところは 0 のままにしてフィルインをなくす方法が考えられる. 行列の構造に依存せず汎用的に使える方法である. 強制的に 0 にした分の値を対角要素に補正を加える方法を修正 ILU (MILU) 前処理とよび, 行列の構造によっては有効である.

(2) フィルインを許す方法 (レベル指定) (ILU(p))

係数行列 A のそれぞれの要素 \(a_{ij}\) にレベル \(L_{ij}\) を定義する. 要素の値が小さければレベルは大きく, 要素の値が大きければレベルは小さくなるものとする. 最初は \(a_{ij} \ne 0\) であれば \(L_{ij} = 0\), \(a_{ij} = 0\) であれば \(L_{ij} = \infty\) とする. ILU 分解の途中で \(a_{ij} = a_{ij} – a_{ik} a_{kj}\) という計算を行うが, このとき次のように \(a_{ij}\) のレベルを更新する.
\[
L_{ij} = min(L_{ij}, L_{ik} + L_{kj})
\] このようにすると最初に 0 でなかった要素のレベルは分解終了時にも 0 のままである.

ILU(p) 前処理では分解終了時に p より大きいレベルの要素を 0 にする. p = 0 とすると ILU(0) に一致する.

(3) フィルインを許す方法 (しきい値指定) (ILUT)

ILU(p) ではおおむね大きい要素を残し小さい要素を捨てる方針であるが, 基準にしているのは行列の構造であり要素の値は見ていない. そこで要素の値が大きいかどうかを直接見て判断する方法が考えられ、それには種々のやり方が考えられる.

例えば, 基本的には τ より大きな値の要素を残すことにするが, 最低でももともと 0 でなかった要素は残し (すなわち, 最低でも ILU(0) に相当), 各行あたりの残す要素数の上限をもともと 0 でなかった要素 + 大きい方から p 個とする (フィルインが多くなり過ぎないように制限する). これは ILUT(τ, p) とよばれる.

ILU 分解は対角要素が 0 だと失敗する. また, 条件が悪い行列では不安定になる. そこで, 行交換を行う方法が考えられる. これは ILUTP とよばれる (P は pivoting を表す).

5.3.8.2 SSOR (Symmetric successive-over-relaxation) 前処理

定常反復法で用いられる反復行列 M を前処理行列として用いる方法を行列分離型前処理という. その中でも対称行列に対する SOR 法の反復行列を用いるのが SSOR 前処理である. この方法は対称行列用の前処理系として使用されるが, 非対称行列にも使用できる.

非対称行列 A を次のように分解する.
\[
A = L + D + U
\] ここで, L は対角要素が 0 の下三角行列, D は対角要素だけからなる行列 (対角行列), U は対角要素が 0 の上三角行列である. これは LU 分解とは違って A を 3 つの部分に分けただけ (足し算) なので簡単に求められる. これを使って前処理行列を次のように定義する.
\[
M(\omega) = (I + \omega LD^{-1}) \frac{1}{\omega (2 – \omega)} D(I + \omega D^{-1}U)
\] この方法は計算が簡単であるが効果が期待できる. ただし, ω は SOR 法と同じように加速パラメータであるが最適な値を前もって知ることは難しいため経験値を使うことになる.

5.3.8.3 点ヤコビ (Point Jacobi) 前処理

前処理行列 M として係数行列の対角行列 (対角要素だけからなる正方行列) D を用いる.
\[
M = D
\] 最も簡単な前処理であるが, 問題によっては効果がある.

5.3.8.4 非対称行列用の前処理の比較

数値実験 (9) SSOR 前処理および ILU(p) 前処理の効果

非対称行列用のポアソン方程式の例題 (N = 31) の差分方程式を解く.

GMRES における SSOR 前処理および ILU(p) 前処理の効果を比較する.

前処理の効果は大きく, 収束に 100 回かかっていたのが 16 ～ 37 回で収束するようになる. ただし, SSOR 前処理は比較的計算が簡単なのに対し ILU 前処理では分解にそれなりの時間がかかるため, トータルの計算時間を分析してみる必要がある.

対称行列の場合も同様であるが, 疎行列の連立一次方程式では解法の選択も重要であるが, 収束性が前処理に大きく依存するところがあるので問題に応じて適切な前処理を選択することが重要である.

数値実験 (10) 前処理方式の比較

非対称行列用の前処理では分離前処理, 右前処理, および, 左前処理があるが, これらの収束性にどの程度の違いがあるか確認する.

GMRES 法に ILU(0) 前処理を適用して, 非対称行列用のポアソン方程式の例題 (N = 31) の差分方程式を解く.

分離前処理では \(M = M_1 M_2 = LU\), 右前処理では \(M_2 = LU\), 左前処理では \(M_1 = LU\) とする.

この例ではほとんど差がないことが確認された. 一般的には残差の大きさが変わらない右前処理を使う.