共役勾配法の解説で出てきた前処理行列の求め方について詳しく説明します。
目次
前処理行列
ある線型方程式系、
$${\bf A}{\bf x} = {\bf b} \tag{1}$$
において、行列解法の収束性を上げるため、
$${\bf P}^{-1}{\bf A}{\bf x} ={\bf P}^{-1} {\bf b} \tag{2}$$
のように変形したときの行列${\bf P}$が前処理行列です（詳細は前述の記事参照）。
(2)式を見ると、もし${\bf P}={\bf A}$であれば、${\bf x} ={\bf A}^{-1} {\bf b}$となり方程式は直接解けていることになります。つまり、究極的には${\bf P}={\bf A}$と置けばよいのです。
しかし、${\bf A}$の逆行列を求めることはメモリや計算量が多くかかり、一般的に容易ではないので現実的ではありません（そもそもこれができれば、わざわざ共役勾配法を使うまでもありません）。
そこで、何か近似的に${\bf A}$に近い行列を${\bf P}$とすることを考えます。
対角スケーリング
最も簡単な方法は行列${\bf A}$の対角成分のみを抜き出した対角行列を前処理行列とすることです。この方法は対角スケーリングと呼ばれています。
$${\bf P}=\left(\begin{array}{cccc} a_{11} & 0 & \cdots & 0 \\ 0 & a_{22} &\cdots & 0 \\ \vdots & \vdots & \ddots & 0 \\ 0 & 0 &\cdots & a_{NN} \end{array}\right) \tag{3}$$
対角行列の逆行列は容易に求めることができるので、行列計算もすぐに出来ます。
$${\bf P}^{-1}=\left(\begin{array}{cccc} 1/a_{11} & 0 & \cdots & 0 \\ 0 & 1/a_{22} &\cdots & 0 \\ \vdots & \vdots & \ddots & 0 \\ 0 & 0 &\cdots & 1/a_{NN} \end{array}\right) \tag{4}$$
ただし、この前処理行列${\bf P}$はそれほど${\bf A}$に近くはないため、行列解法の収束性が大きく改善されない場合があります。しかし、計算が簡単で計算量も少なくて済むため最初の選択肢として利用できます。
不完全コレスキー分解
次に${\bf A}$が正定値対称行列の場合に適用できる不完全コレスキー分解（Incomplete Cholesky factorization）という方法があります。
コレスキー分解（Cholesky factorization）とは対称行列${\bf A}$を下三角行列${\bf L}$とその転置行列${\bf L}^{\rm T}$の積に分解することをいいます。
$${\bf A}={\bf L}{\bf L}^{\rm T} \tag{5}$$
この場合、行列${\bf L}$の要素のうち、もとの行列${\bf A}$でゼロだった要素がゼロでなくなるfill-inと呼ばれる現象が起こります。fill-inが起きるとせっかく疎行列で少なかった非ゼロ要素が増えてしまい、計算効率やメモリ容量などが悪化します。
そこで、コレスキー分解を完全に行わず、
$${\bf A}=\bar{\bf L} \bar{\bf L}^{\rm T} +{\bf R} \tag{6}$$
のように残差${\bf R}$が残る不完全な形に分解したものを不完全コレスキー分解といいます。ここで行列${\bf A}$のゼロの要素は下三角行列$\bar{\bf L}$でもゼロになっています。fill-inがない不完全コレスキー分解はIC(0)と書かれたりします。
そして、このときの${\bf P}=\bar{\bf L} \bar{\bf L}^{\rm T}$を前処理行列としてやります。こうすると、もとの${\bf A}$行列に近い前処理行列を得ることができます。
不完全コレスキー分解の方法は行列${\bf A}$のゼロの要素が下三角行列${\bf L}$でもゼロになるようにコレスキー分解を行うことです（別の方法としてfill-inを少しだけ許すというものもあります）。
Pythonでサンプルコードを書くと次のようになります。
for k in range(n):
a[k][k] = math.sqrt(a[k][k])
for i in range(k+1, n):
if a[i][k] != 0:
a[i][k] /= a[k][k]
for j in range(k+1, n):
for i in range(j, n):
if a[i][j] != 0:
a[i][j] -= a[i][k] * a[j][k] ここで計算された行列aの下三角部分が求める$\bar{\bf L}$になっています。
※ここではサンプルのためN×Nの２次元配列で行列を定義していますが、大規模な疎行列の場合は、次の記事のようにデータの持たせ方に工夫が必要です。
不完全LU分解
不完全コレスキー分解は対称行列の場合でしたが、非対称行列の場合は不完全LU分解（Incomplete LU factorization）という方法があります。
これは(6)式と同様に不完全な形でLU分解を行い${\bf P}=\bar{\bf L} \bar{\bf U}$を前処理行列とするものです。この場合も、下三角行列$\bar{\bf L}$と上三角行列$\bar{\bf U}$は行列${\bf A}$のゼロ要素と同じ要素がゼロになっています。fill-inがないのでILU(0)と書かれます。
これもPythonでサンプルコードを示しておきます。
for k in range(n-1):
for i in range(k+1, n):
if a[i][k] != 0:
a[i][k] /= a[k][k]
for j in range(k+1, n):
if a[i][j] != 0:
a[i][j] -= a[i][k] * a[k][j] 行列aの下三角が$\bar{\bf L}$、上三角が$\bar{\bf U}$になります。
前処理行列の扱い方
(2)式を見ると前処理行列${\bf P}$の逆行列${\bf P}^{-1}$が必要になると思われがちです。ところが、「共役勾配法の種類」のアルゴリズムをよく見ると、前処理行列は何らかのベクトル${\bf y}$に対して${\bf P}^{-1}{\bf y}$の形で作用していることがわかります。これは言い換えると、
$${\bf P}{\bf x}={\bf y} \tag{7}$$
の解${\bf x}={\bf P}^{-1}{\bf y}$であることを意味します。
不完全コレスキー分解や不完全LU分解を行った利点は、${\bf P}$を上三角、下三角行列で表した点にあります。なぜならば、この形に分解すると(7)式を解くことが容易にできるからです。
$${\bf L}{\bf U}{\bf x}={\bf y} \tag{8}$$
は、
$${\bf L}{\bf z}={\bf y} \tag{9}$$
$${\bf U}{\bf x}={\bf z} \tag{10}$$
という２式に分けることができます。(9)式は前進代入、(10)式は後退代入で容易に解を求めることができます。この解を${\bf P}^{-1}{\bf y}$として計算を進めていけばよいのです。
まとめ
前処理行列の求め方について解説しました。ここであげた以外にも様々な前処理行列があり、解く行列や解法により使い分けられています。当サイトの流体解析ツールCATCFDzeroではILU(0)の前処理行列を使っています。