参考文献

1. Wendland, H. (2017). Numerical Linear Algebra: An Introduction (Cambridge Texts in Applied Mathematics). Cambridge: Cambridge University Press. doi:10.1017/9781316544938
2. 森 正武, 共立数学講座12 数値解析 (第2版), 共立出版 (2002)

前処理とは

前処理 (preconditioning) とは，ある線形連立方程式，Ax = b をより簡単に解ける同じ解を持つ方程式，By = c に変換することである．

もう少し具体的に

Ax = b を

ただし，P_R y = x として解くことである．

P_L や P_R をうまく選んで，反復回数を少なく解を得ようとすることである．

良い前処理とは

極端な例を挙げると，前処理によって，P_L A = I とできれば，何もしなくて良くなるわけです．

したがって，単位行列に近くなるように，前処理を施すことが，良い前処理だということになります．

共役勾配法では，固有値が縮退していると反復回数が少なくなります．単位行列は1つの固有値しか持たないので，この点で，固有値の縮退している数が多いほど，単位行列に近いと言えます．

より数学的に厳密なことが知りたい方は，参考文献を参照してください．

前処理付き共役勾配法

今回は共役勾配法に前処理をつけることを考えます．

そのため A は正定値対称行列として，前処理後の形も，正定値対称行列であるものを探します．すなわち，

となる，行列 S を求め，共役勾配法を実行します．

共役勾配法のアルゴリズム

ここで，共役勾配法のアルゴリズムを書いておきます．

これは，Ax=b を満たす解 x を求めるアルゴリズムです．

Input: x_0 and \epsilon
Set p_0 = r_0 = b - A x_0
while || r_j || > \epsilon do;

1. \alpha_j \coloneqq \frac{|| r_j ||_2^2}{\langle Ap_j, p_j\rangle_2}
2. x_{j+1} \coloneqq x_j + \alpha_j p_j
3. r_{j+1} \coloneqq r_j - \alpha_j Ap_j
4. \beta_{j+1} \coloneqq \frac{||r_{j+1}||_2^2}{||r_j||_2^2}
5. p_{j+1} \coloneqq r_{j+1} + \beta_{j+1} p_j
6. j \coloneqq j+1

done

前処理付き共役勾配法のアルゴリズム（導出）

行列 S が与えられたとき，連立方程式

の共役勾配法は，先述のアルゴリズムを適用すると，

Input: y_0 and \epsilon
Set p_0 = r_0 = S^\mathrm{T} (b - AS y_0)
while || r_j || > \epsilon do;

1. \alpha_j \coloneqq \frac{|| r_j ||_2^2}{\langle S^\mathrm{T}AS p_j, p_j\rangle_2}
2. y_{j+1} \coloneqq y_j + \alpha_j p_j
3. r_{j+1} \coloneqq r_j - \alpha_j S^\mathrm{T}AS p_j
4. \beta_{j+1} \coloneqq \frac{||r_{j+1}||_2^2}{||r_j||_2^2}
5. p_{j+1} \coloneqq r_{j+1} + \beta_{j+1} p_j
6. j \coloneqq j+1

done

となる．この反復法によって得られた y に対し，x = Sy とすることで，求めたい解 x が求まる．

これをもう少し変形して，直接解 x を求めるアルゴリズムに書き換えておこう．

変形

以下のように記号を定義する．

今，r^\prime_j = (S^\mathrm{T})^{-1} r_j = (S^\mathrm{T})^{-1} S^{-1} z_j = P^\mathrm{T} z_j である．
これらを使って，

となることがわかる．したがって，

となる．

前処理付き共役勾配法のアルゴリズム

したがって，前処理付き共役勾配法のアルゴリズムは以下のようになる．
（注意：上までは \prime がついていたものの \prime を省略している．）

Input: x_0 and \epsilon
Set r_0 = b - A x_0 and p_0 = z_0 = Pr_0
while || r_j || > \epsilon do;

1. \alpha_j \coloneqq \frac{ \langle r_j, z_j \rangle_2}{\langle Ap_j, p_j \rangle_2}
2. x_{j+1} \coloneqq x_j + \alpha_j p_j
3. r_{j+1} \coloneqq r_j - \alpha_j Ap_j
4. \beta_j \coloneqq \frac{\langle r_{j+1}, z_{j+1} \rangle_2}{\langle r_j, z_j \rangle_2}
5. z_{j+1} \coloneqq P r_{j+1}
6. p_{j+1} \coloneqq z_{j+1} + \beta_{j+1} p_j
7. j \coloneqq j+1

done

P として何を選べば良いか

これは，様々なものがあるが，ここでは，不完全コレスキー分解による前処理を紹介する．

正定値対称行列 A はコレスキー分解ができる：

ここで，L は下三角行列で，D は対角行列である．
D の対角成分の平方根を対角成分とする対角行列を D^{1/2} と表すと，

となり，すなわち，

となる． I は単位行列である．

これを利用する．もし，コレスキー分解が速やかに実行できて，L と D を得ることができるとすると，
前処理行列 S として，S = ((LD^{1/2})^\mathrm{T})^{-1} を選べば，

と求めることができる．

しかし，コレスキー分解をすることと，連立方程式 Ax = b を解くことは同じことである．
一方で，これと似たようなことをすれば，良い前処理行列が見つかりそうである．

したがって，コレスキー分解に近いことをして，前処理行列 S を選び，できる限り単位行列に近くなるように前処理を行いたい．

不完全コレスキー分解

上記の通り，コレスキー分解のようなことをして，前処理行列 S を探したい．

ここでは，不完全コレスキー分解を紹介する．

行列 A の要素が0である添字の集合を \mathcal{M} とすると，
行列 A の不完全コレスキー分解とは，コレスキー分解を行う時に，

1. \mathrm{M} に属する添字の場合には，コレスキー分解の計算を行わず，その要素を0とする．
2. \mathrm{M} に属さない添字の場合には，通常通りコレスキー分解の計算を行い，要素の値を決定する．

と言う手順で，下三角行列 \tilde{L} と対角行列 \tilde{D} を決定する方法のことである．

この方法では，A が疎行列の場合には，分解が速やかに(Aの0でない要素の個数程度の処理で)実行できる．

一方で，コレスキー分解ではないので，得られた分解，

は A ではない．しかしとてもおしいものではあると期待できる．

不完全コレスキー分解前処理付き共役勾配法

不完全コレスキー分解によって得られる行列 \tilde{L}，\tilde{D} をつかって前処理行列を構成するのが，不完全コレスキー分解前処理であり，これを用いた前処理付き共役勾配法が，不完全コレスキー分解前処理付き共役勾配法である．

すなわち，

とする．これはすなわち，

と選ぶことに対応する．

不完全コレスキー分解前処理付き共役勾配法のアルゴリズム

不完全コレスキー分解前処理付き共役勾配法のアルゴリズムをまとめておく．

Input: x_0 and \epsilon
Set r_0 = b - A x_0 and p_0 = z_0 = (\tilde{L} \tilde{D} \tilde{L}^\mathrm{T})^{-1} r_0
while || r_j || > \epsilon do;

1. \alpha_j \coloneqq \frac{ \langle r_j, z_j \rangle_2}{\langle Ap_j, p_j \rangle_2}
2. x_{j+1} \coloneqq x_j + \alpha_j p_j
3. r_{j+1} \coloneqq r_j - \alpha_j Ap_j
4. z_{j+1} \coloneqq (\tilde{L} \tilde{D} \tilde{L}^\mathrm{T})^{-1} r_{j+1}
5. \beta_{j+1} \coloneqq \frac{\langle r_{j+1}, z_{j+1} \rangle_2}{\langle r_j, z_j \rangle_2}
6. p_{j+1} \coloneqq z_{j+1} + \beta_{j+1} p_j
7. j \coloneqq j+1

done

ここで，z_j = (\tilde{L} \tilde{D} \tilde{L}^\mathrm{T})^{-1} r_j \,\, (j = 0, 1, \dots) の計算は，\tilde{L}\tilde{D}^{1/2}(\tilde{L}\tilde{D}^{1/2})^\mathrm{T}) z_j = r_j としておいて，

と2段階で解く．第1式は前進代入，第2式は後退代入で解くことができる．

コレスキー分解

最後に付録としてコレスキー分解のアルゴリズムを示す．

以下は，正確には修正コレスキー分解と呼ばれる．

まずは導出．

コレスキー分解とは，正定値対称行列 A を

と下三角行列 L と対角行列 D を用いて分解することであった．

ここで，下三角行列 L = (l_{ij})_{ij} は，

これは，S = LD^{1/2} とすれば，

となって，S は下三角行列，S^\mathrm{T} は上三角行列であるから，LU分解の対称行列の場合の特殊な表記と見ることができる．

LU分解は，もしできるのであれば一意に定まるので，コレスキー分解も一意である．
実際，A = SS^\mathrm{T} に対して，A = (a_{ij})_{ij}, S = (s_{ij})_{ij} とすると，i \leq j に対して，

となるが，これは，

となって，これによって，S の各成分 s_{ij} を逐次的に求めることができる．

いま，S と L, D の関係は，L = (l_{ij})_{ij}, D = \mathrm{diag}(d_1, \dots, d_n) とすると，

であるから，結局，

となる．

コレスキー分解のアルゴリズム

コレスキー分解のアルゴリズムは以下のようになる．

Input: A \in \mathbb{R}^{n\times n} symmetric and positive definite

1. l_{11} = a_{11}
2. for i =2 to n do;
   1. for j = 1 to i-1 do;
      1. l_{ij} \coloneqq a_{ij}
      2. for k = 1 to j-1 do;
         • l_{ij} \coloneqq l_{ij} - l_{ik}l_{jk}
      3. l_{ij} \coloneqq l_{ij}/d_j
   2. d_i \coloneqq a_{ii}
   3. for k=1 to i-1 do;
      • d_i \coloneqq d_i - l_{ik} l_{ik}
   4. d_i \coloneqq \sqrt{d_i}

不完全コレスキー分解

不完全コレスキー分解は，a_{ij} = 0 となる時は，上記の処理を行わず，l_{ij} = 0 とすることである．