GMRes

一般化残差最少法、GMRes法

Last Update:2009年05月15日

概要

GMRes法Y.Saadらによって提案された、非対称行列を解くことができる反復法である。

GMRes法とは回目の反復において、部分空間 $x_0+\cal{K}_k$ に属する全てのベクトル $\b{x}_k$ から残差の２乗ノルム $J(\b{x})=||\b{b}-\b{A}\b{x}||_2$ を最小化するような解 $\b{x}_k$ を探す見つける方法である。つまり、

$find \quad \b{x}\in x_0+\cal{K}_k \qquad \qquad that \quad minimize \quad ||\b{b}-\b{A}\b{x}||_2$

但し、 $\cal{K}_k$ はk階のクリロフ部分空間 $span\{\b{r}_0,\b{A}\b{r}_0,\b{A}^2\b{r}_0,\cdots,\b{A}^{k-1}\b{r}_0\}$ である。

基本的な手続き

簡単のため、まず回反復を行った後に解を求める方法を説明する。この方法では、まず個の正規直交基底を作成し、直交基底の中から残差を最小にする解を求める。残差の大きさが求まるのが解の反復の後であるので、残差が十分小さくなれば途中で反復を打ち切るようなことはできない。各反復毎に残差の多きさあを求める方法は後に説明する。

Arnoldi法を用いて正規直交化されたクリロフ部分空間 $\cal{K}_m$ の正規直交基底 $\b{v}_1,\ldots,\b{v}_m$ を並べた行列を $\b{V}_m$ ,Gram-Shumitの直交化の際の係数を行列にしたものを $\b{H}_m$ とする。 $\b{H}_m$ は $\b{H}_m=0 \quad for \quad i>j+1$ を満たす、 $m+1 \times m$ のヘッセンベルグ行列である。

ヘッセンベルグ行列を用いた残差の表現

$\b{v}_1,\ldots,\b{v}_m$ はクリロフ部分空間 $\cal{K}_m$ の基底であったから、 $\b{x}_0+\cal{K}_m$ に属するベクトル $\b{x}$ は次のように表すことができる。

$\b{x}=\b{x}_0+y_1\b{v}_1+y_2\b{v}_2+\cdots+y_m\b{v}_m=\b{x}_0+\b{V}_m\b{y}$

ここで $\b{y}$ は次のようなベクトルである。

$\b{y}^T=(y_1,y_2,\ldots,y_m)$

これらを用いてベクトル $\b{x}$ に対する残差 $\b{r}$ を計算する。

$\beta=||\b{r}_0||_2$ 、 $\b{e}^T_1=(\overbrace{1,0,\ldots,0}^{m+1})$ とおくと残差 $\b{r}$ は

$\b{r}=\b{b}-\b{A}\b{x}\\ \qquad = \b{b}-A(\b{x}_0+\b{V}_m\b{y})\\ \qquad = \b{r}_0-\b{A}\b{V}_m\b{y}\\ \qquad = \beta\b{v}_1-\b{V}_m\bar{\b{H}}_m\b{y}\\ \qquad = \b{V}_m(\beta\b{e}_1-\bar{\b{H}}_m\b{y})$

となる。但し、 $\b{v}_1=\b{r}_0/||\b{r}_0||_2$ であることから $\b{r}_0=\beta\b{v}_1$ となること、また、Arnoldi法の性質から $\b{A}\b{V}_m=\b{V}_m\bar{\b{H}}_m$ となることも用いた

これらを用いると最小にすべき残差の２乗ノルムは

$J(\b{y})=||\b{b}-\b{A}\b{x}||_2=||\b{V}_m(\beta\b{e}_1-\bar{\b{H}}_m\b{y})||_2=||(\beta\b{e}_1-\bar{\b{H}}_m\b{y})||_2$

となる。ここで、 $\b{v}_1,\ldots,\b{v}_m$ が正規直交基底であることを用いた。

結局GMRes法は $||(\beta\b{e}_1-\bar{\b{H}}_m\b{y})||_2$ を最小化するベクトル $\b{y}$ を探すことと同じといえる。

GMRes法(ヘッセンベルグ行列)

$make \qquad \b{V}_m \qquad \bar{\b{H}}_m \quad with \quad Arnoldi \quad method$
$find \quad \b{y}_m\in\cal{R}^m \qquad \qquad that \quad minimize \quad ||(\beta\b{e}_1-\bar{\b{H}}_m\b{y})||_2$
$\b{x}_m=\b{x}_0+\b{V}_m\b{y}_m$

上三角行列を用いた残差の表現

ヘッセンベルグ行列にギブンス回転(Givens Rotation)の操作を繰り返すことで上三角行列に変換することができる。回転操作はノルムを変化させないので回転した後の残差のノルムを最小化すれば、回転する前の残差のノルムを最小化したことと同じになる。

m=4の場合の例

わかりやすいように m=4 の場合の $5\times4$ のヘッセンベルグ行列 $\bar{\b{H}}_4$ を例にとって上三角行列への変換を説明する。

$\bar{\b{H}}_4 = \left(\begin{array} h_{11}&&h_{12}&&h_{13}&&h_{14}\\ h_{21}&&h_{22}&&h_{23}&&h_{24}\\ &&h_{32}&&h_{33}&&h_{34}\\ && &&h_{43}&&h_{44}\\ && && &&h_{54}\end{array}\right)$

次のような $5\times5$ の行列 $\Omega_1$ を考える

$\Omega_1 = \left(\begin{array} c_1&&s_1&& && && \\ -s_1&&c_1&& && && \\ && && 1 && && \\ && && &&1 && \\ && && && &&1 \end{array}\right)$ 但し、 $s_1=\frac{h_{21}}{\sqrt{h_{11}^2+h_{21}^2}}$ , $c_1=\frac{h_{11}}{\sqrt{h_{11}^2+h_{21}^2}}$

この行列 $\Omega_1$ を $\bar{\b{H}}_4$ の左から掛けたものを $\bar{\b{H}}_4^{(1)}$ とすると

$\b{H}^{(1)}_4=\Omega_1\bar{\b{H}}_4=\left(\begin{array} h_{11}^{(1)}&&h_{12}^{(1)}&&h_{13}^{(1)}&&h_{14}^{(1)}\\ &&h_{22}^{(1)}&&h_{23}^{(1)}&&h_{24}^{(1)}\\ &&h_{32}&&h_{33}&&h_{34}\\ && &&h_{43}&&h_{44}\\ && && &&h_{54}\end{array}\right)$

となり (2,1) 成分を消去できる。

同様にして次のような $5\times5$ の行列 $\Omega_2$ を考える

$\Omega_1 = \left(\begin{array} 1&& && && && \\ &&c_2&&s_2&& && \\ &&-s_2&&c_2&& && \\ && && &&1 && \\ && && && &&1 \end{array}\right)$ 但し、 $s_2=\frac{h_{32}}{\sqrt{{$h_{22}^{(1)}$}^2+h_{32}^2}}$ , $c_2=\frac{h_{22}^{(1)}}{\sqrt{{$h_{22}^{(1)}$}^2+h_{32}^2}}$

この行列 $\Omega_2$ を $\bar{\b{H}}_4^{(1)}$ の左から掛けたものを $\bar{\b{H}}_4^{(2)}$ とすると

$\b{H}^{(2)}_4=\Omega_2\bar{\b{H}}_4^{(1)}=\Omega_2\Omega_1\bar{\b{H}}_4=\left(\begin{array} h_{11}^{(1)}&&h_{12}^{(1)}&&h_{13}^{(1)}&&h_{14}^{(1)}\\ &&h_{22}^{(2)}&&h_{23}^{(2)}&&h_{24}^{(2)}\\ && &&h_{33}^{(2)}&&h_{34}^{(2)}\\ && &&h_{43}&&h_{44}\\ && && &&h_{54}\end{array}\right)$

となり (3,2) 成分を消去できる。

この操作を繰り替えすことでヘッセンベルグ行列 $\bar{\b{H}}_4$ を上三角行列 $\bar{\b{R}}_4$ に変換することができる。

$\bar{\b{R}}_4=\bar{\b{H}}^{(4)}=\overbrace{\Omega_4\Omega_3\Omega_2\Omega_1}^{Q_4}\bar{\b{H}}_4=Q_4\bar{\b{H}}_4=\left(\begin{array} h_{11}^{(1)}&&h_{12}^{(1)}&&h_{13}^{(1)}&&h_{14}^{(1)}\\ &&h_{22}^{(2)}&&h_{23}^{(2)}&&h_{24}^{(2)}\\ && &&h_{33}^{(3)}&&h_{34}^{(3)}\\ && && &&h_{44}^{(4)}\\ && && && 0 \end{array}\right)$

ここで $Q_4=\Omega_4\Omega_3\Omega_2\Omega_1$ とした。

一般的な場合

前では m=4 の場合について具体的にヘッセンベルグ行列から上三角行列への変換について述べた。ここでは一般的な場合について述べる。

次のよな行列 $\Omega_i$ を考える

$\Omega_i=\left(\begin{array} 1 & & & & & & & \\ &\ddots& & & & & & \\ & &1& & & & & \\ & & &c_i &s_i& & & \\ & & &-s_i&c_i& & & \\ & & & & &1& & \\ & & & & & &\ddots& \\ & & & & & & &1\\ \end{array}\right) \quad\qquad\begin{array}{l} \\ \\ \\ {\leftarrow row \quad i}\\ {\leftarrow row \quad i+1 }\\ \\ \\ \\\end{array}$ 但し、 $s_i=\frac{h_{i+1,i}}{\sqrt{{$h_{ii}^{(i-1)}$}^2+h_{i+1,i}^2}}$ , $c_i=\frac{h_{ii}^{(i-1)}}{\sqrt{{$h_{ii}^{(1)}$}^2+h_{i+1,i}^2}}$

この $\Omega$ を掛け合わせた次のような行列 Q_m を考えると

$Q_m=\Omega_m\Omega_{m-1}\ldots\Omega_1$

上三角行列 $\bar{R}_m$ への変換は次のようになる。

$\bar{R}_m=\bar{\b{H}}_m^{(m)}=Q_m\bar{\b{H}}_m$

$\bar{\b{g}}_m=Q_m(\beta \b{e}_1)=(\gamma_1,\ldots,\gamma_{m+1})^T$

これらを用いると残差 $\b{r}$ は次のように変形できる。

$\b{r}=\b{b}-\b{A}\b{x}\\ \qquad = \b{V}_m(\beta\b{e}_1-\bar{\b{H}}_m\b{y})\\ \qquad = \b{V}_mQ^T_mQ_m(\beta\b{e}_1-\bar{\b{H}}_m\b{y})\\ \qquad \ \b{V}_mQ^T_m(\bar{\b{g}}_m-\bar{R}_m\b{y})$

これらを用いると最小にすべき残差の２乗ノルムは

$J(\b{y})=||\b{b}-\b{A}\b{x}||_2=||\b{V}_mQ^T_m(\bar{\b{g}}_m-\bar{R}_m\b{y})||_2=||\bar{\b{g}}_m-\bar{R}_m\b{y}||_2$

のようになる。但し、 $\b{V}_mQ^T_m$ がユニタリーであることを用いた。

これらを用いるとGMRes法のアルゴリズムは次のようになる。

GMRes法(上三角行列)

$make \qquad \b{V}_m \qquad \bar{\b{H}}_m \quad with \quad Arnoldi \quad method$
$make \quad \bar{\b{g}}_m \quad \bar{R}_m \quad from \quad {\b{H}}_m$
$find \quad \b{y}_m\in\cal{R}^m \qquad that \quad minimize \quad ||\bar{\b{g}}_m-\bar{R}_m\b{y}||_2$
$\b{x}_m=\b{x}_0+\b{V}_m\b{y}_m$

上三角行列を用いた残差ノルムの最小値を求める方法

この章は残差を最小にする $\b{y}_m$ の値を上三角行列 $\bar{\b{R}}_m$ から求める方法について述べる。

$\b{g}_m$ を $\bar{\b{g}}_m$ の m+1 番目の成分を除いた、大きさがのベクトルであるとする。また、 $\b{R}_m$ を $\bar{\b{R}}_m$ から m+1 番目の行を除いた正方行列であるとする。 $\bar{\b{R}}_m$ の m+1 番目の行は成分が全て０であるので

$||\b{b}-A\b{x}||_2 = ||\bar{\b{g}}_m-\bar{R}_m\b{y}||_2 = ||\b{g}_m- R_m\b{y} + \gamma_{m+1}\b{e}_{m+1}||_2$

が成立つ。 $\b{g}_m- R_m\b{y}$ と $\b{e}_{m+1}$ が互いに直交していることを用いると

$||\b{g}_m- R_m\b{y} + \gamma_{m+1}\b{e}_{m+1}||_2 \ge ||\gamma_{m+1}\b{e}_{m+1}||_2 = |\gamma_{m+1}|$

が成立つ。等号成立は $\b{y}=R_m^{-1}\b{g}_m$ の時である。

以上から残差ノルムを最小化するベクトル $\b{y}_m$ は次のように求めることができる。

$\b{y}_m=R_m^{-1}\b{g}_m$

このとき、残差ノルムの最小値は $|\gamma_{m+1}|$ である。またこの時の残差ベクトル $\b{r}_m$ は、

$\b{r}_m=\b{b}-\b{A}\b{x}\\ \qquad = \b{V}_mQ^T_m(\bar{\b{g}}_m-\bar{R}_m\b{y}_m) \qquad = \b{V}_mQ^T_m(\gamma_{m+1}\b{e}_{m+1})$

のようになる。

なお、

$\b{y}_m=R_m^{-1}\b{g}_m$

の計算は R_m が上三角行列であることから後退代入を用いて簡単に求めることができる。

これらを用いるとGMResのアルゴリズムは次のようになる。

GMRes法のアルゴリズム（上三角形、後退代入）

$make \qquad \b{V}_m \qquad \bar{\b{H}}_m \quad with \quad Arnoldi \quad method$
$make \quad \bar{\b{g}}_m \quad \bar{R}_m \quad from \quad {\b{H}}_m$
$Compute \quad \b{y}_m=R_m^{-1}\b{g}_m$
$\b{x}_m=\b{x}_0+\b{V}_m\b{y}_m$

逐次的に残差の大きさを求める方法

今までの方法は個の基底 $\b{V}_m$ を作ってから残差の２乗ノルムを最小化するような解 $\b{x}_m$ を求めていたが、実際には各反復毎に残差の大きさ $||\b{r}_k||_2$ を求め、必要とされる小ささよりも小さくなった段階で反復をうちきることで無駄な反復をするすることなく解を求めることができる。

ここでは各反復毎に正規直交基底や上三角形行列、残差の大きさを更新する方法について考える。

基底の更新

Arnordi法のアルゴリズムからk回目の反復において生成されるクリロフ部分空間の基底ベクトル $\b{v}_{k+1}$ はベクトル $\b{\omega}=\b{A}\b{v}_k$ をGram-Shumitの直交化法やHausholder変換を用いてベクトル列 $\b{v}_1,\ldots,\b{v}_k$ と直交化した後に正規化したものだった。

上三角行列の更新

$\b{v}_{k+1}$ を作る際の直交化の際の係数は $h_{k,1},\ldots,h_{k,k}=(\b{A}\b{v}_k,\b{v}_1),\ldots,(\b{A}\b{v}_k,\b{v}_k)$ であり、直交化後に正規化するときの係数が $h_{k,k+1}$ であった。

ここでベクトル $\b{h}_k$ を次のように置く。 $\b{h}_k=\{h_{k,1},\ldots,h_{k,k+1}\}$

この $\b{h}_k$ を用いて回目の反復におけるヘッセンベルグ行列 $\bar{\b{H}}_k$ は次のように表される。

$\bar{\b{H}}_k=[\bar{\b{H}}_{k-1}',\b{h}_k]$

但し $\bar{\b{H}}_{k-1}'$ は k-1 回目の反復におけるヘッセンベルグ行列 $\bar{\b{H}}_{k-1}$ の下に成分が全て０である行を追加した行列である。

k-1 反復における上三角行列 $\b{R}_{k-1}$ から反復における上三角行列 $\b{R}_k$ を求める方法を考える

$\b{R}_k$ の一番下の行に成分が全て０である行を追加した行列は $\bar{\b{R}}_k$ であった。この行列は次のように回転行列を用いてヘッセンベルグ行列 $\bar{\b{H}}_k$ から求めることができた。

$\bar{\b{R}}_k=\b{Q}_k\bar{\b{H}}_k=\Omega_k\Omega_{k-1}\ldots\Omega_1\bar{\b{H}}_k$

ここで $\bar{\b{H}}_k$ に先ほどの関係を代入して

$\bar{\b{R}}_k=\Omega_k\Omega_{k-1}\ldots\Omega_1\bar{\b{H}}_k=\Omega_k\Omega_{k-1}\ldots\Omega_1[\bar{\b{H}}_{k-1}',\b{h}_k]\\ \qquad \qquad = \Omega_k[\Omega_{k-1}\ldots\Omega_1\bar{\b{H}}_{k-1}',\quad \Omega_{k-1}\ldots\Omega_1\b{h}_k]\\ \qquad \qquad =\Omega_k[\bar{\b{R}}_{k-1}',\quad \tilde{\b{h}}_k]$

ここで $\bar{\b{R}}_{k-1}'$ は $\bar{\b{R}}_{k-1}$ の一番下の行に成分が全て０の行を追加した行列である。 $\Omega_k$ は一番下の大きさが２の対角成分以外は対角が単位行列である回転行列であった。

$\Omega_k=\left(\begin{array} 1 & & & & &\\ &\ddots& & & &\\ & &1& & &\\ & & & c_k&s_k&\\ & & &-s_k&c_k&\\ \end{array}\right)$

$\bar{\b{R}}_{k-1}'$ は下２行が０であるから、 $\Omega_k\bar{\b{R}}_{k-1}=\bar{\b{R}}_{k-1}$ である。

回転行列 $\Omega_k$ はベクトル $\Omega_k\tilde{\b{h}}_k$ の k+1 番目の成分が０になるように選ばれる。

つまり、 $-s_k(\tilde{\b{h}}_k)_k+c_k(\tilde{\b{h}}_k)_{k+1}=0$ を満たす。

また $\Omega_k$ はユニタリー行列であるから s_k^2+c_k^2=1 である。

つまり、 $s_k=\frac{(\tilde{\b{h}}_k)_{k+1}}{\sqrt{(\tilde{\b{h}}_k)_k^2+(\tilde{\b{h}}_k)_{k+1}^2}} \qquad \qquad c_k=\frac{(\tilde{\b{h}}_k)_k}{\sqrt{(\tilde{\b{h}}_k)_k^2+(\tilde{\b{h}}_k)_{k+1}^2}}\quad$ である。

これらを用いると

$\bar{\b{R}}_k=\Omega_k[\bar{\b{R}}_{k-1}',\quad \tilde{\b{h}}_k]=[\bar{\b{R}}_{k-1}',\quad \Omega_k\Omega_{k-1}\ldots\Omega_1\b{h}_k]$

このようにして k-1 回目の反復における上三角行列 $\b{R}_{k-1}$ から回目の反復における上三角行列 $\b{R}_k$ を求めることができる。

残差の大きさの更新

回反復目における残差の大きさは $|\gamma_{m+1}|$ と表されるのであった。但しベクトル $\b{\gamma}$ は $\bar{\b{g}}_m$ であり、 $\gamma_{m+1}$ はベクトル $\bar{\b{g}}_m$ の m+1 番目の成分であった。つまりベクトル $\b{g}$ の一番最後の成分の絶対値が残差のノルムになるのである。

ここでベクトル $\b{g}$ がどのように更新されるのかを調べる。 $\bar{\b{g}}_k=\b{Q}_k(\beta\b{e}_1)=\Omega_k\b{Q}_{k-1}(\beta\b{e}_1)=\Omega_k\bar{\b{g}}_{k-1}'$ と表される。ここで $\bar{\b{g}}_{k-1}'$ は大きさのベクトル $\bar{\b{g}}_{k-1}$ の一番最後の成分に０を付け加えることによって大きさを k+1 にしたベクトルである。

$\Omega_k$ は対角成分がから k-1 までが単位行列、成分と k+1 成分が回転行列になっている行列であるから $\bar{\b{g}}_k$ は１ステップ前のベクトルの値 $\bar{\b{g}}_{k-1}$ を元にして次のように求めることができる。

$\bar{\b{g}}_k=\Omega_k\bar{\b{g}}_{k-1}'=\Omega_k\left(\begin{array} g_{k-1,1}\\ g_{k-1,2}\\ \vdots\\ g_{k-1,k}\\ 0 \end{array}\right) = \left(\begin{array} g_{k-1,1}\\ g_{k-1,2}\\ \vdots\\ c_k g_{k-1,k}\\ -s_k g_{k-1,k} \end{array}\right) = \left(\begin{array} \gamma_1 \\ \gamma_2 \\ \vdots\\ \gamma_k \\ -s_k g_{k-1,k} \end{array}\right)$

よって回反復目における残差のノルム $\lambda_k$ は次のようになる

$\lambda_k=|-s_k g_{k-1,k}|=|s_k|| g_{k-1,k}|=|s_k|\lambda_{k-1}$

$|s_k|\le1$ であるから残差が増加することはないこと確認できる。

以上をまとめると、回目の反復ではまず、新しいクリロフ部分空間の正規直行基底 $\b{v}_{k+1}$ を作成し、正規直交基底を作った際の係数ベクトル $\b{h}_k$ に、今までの回転行列 $\Omega_{k-1}\Omega_{k-2}\cdots\Omega_1$ を掛けた後、係数のベクトルのサイズが１つ減るような新しい回転行列 $\Omega_k$ を作成して上三角行列を作る。新しい回転行列を $\b{g}$ ベクトルにかけることによって $\b{g}$ ベクトルを更新し、 $\b{g}$ ベクトルの最後の成分から現在の残差のノルムを求める。

GMRes法のアルゴリズム（逐次更新型）

1.   Compute $\b{r}_0=\b{b}-\b{A}\b{x}_0 \qquad \b{v}_1=\b{r}_0/||\b{r}_0||$
2.   For $k=1,2,\ldots,m\quad$ Do:
3.      $\b{w}=\b{A}\b{v}_k$
4.      For $i=1,2,\ldots,k\quad$ Do:
5.         $h_{k,i}=(\b{w},\b{v}_i)$
6.         $\b{w}=\b{w}-h_{k,i}\b{v}_i$
7.      End Do
8.      $h_{k,k+1}=||\b{w}||_2$
9.      $\b{v}_{k+1}=\b{w}/||\b{w}||_2$
10.     $Compute \quad \tilde{\b{h}}_k=\Omega_{k-1}\cdots\Omega_1\b{h}_k \quad and \quad make \quad \Omega_k \quad from \quad \tilde{\b{h}}_k$
11.     $Compute \quad \bar{R}_k=[\bar{R}_{k-1}',\Omega_k\tilde{\b{h}}_k] \quad and \quad \bar{\b{g}}_k = \Omega_k\bar{\b{g}}_{k-1}'$
12.     $if \quad |\bar{\b{g}}_{k,k+1}|<tolerance \quad$ Break Loop
13. End Do:
14. $Compute \quad \b{y}_k=R_k^{-1}\b{g}_k$
15. $\b{x}_k=\b{x}_0+\b{V}_k\b{y}_k$

Full-GMRes法、GMRes(m)法

Full-GMRes法は適当な反復数を決めずに残差のノルムが小さくなるまで反復を繰り返す方法である。Full-GMRes法では残差が小さくなるまで無限に反復を繰り返し、正規直交基底も収束まで無限に増え続ける。計算機上ではリソースが限られているので現実的ではない。

GMRes(m)法はリスタートつきGMRes法とも呼ばれ、反復数がに達した時点でそれまでのクリロフ部分空間の基底や上三角行列などを全て捨ててその時点での反復解を初期値として計算しなおす（リスタートする）というアルゴリズムである。このリスタトまでの反復数の最適な値を求めることは困難である。が問題依存のある値よりも小さければいくらリスタートを繰り返しても収束が困難であることが知られている。この方法は良く使われる。

サンプルプログラム

次のプログラムは次のような移流拡散問題の例

移流拡散問題の例

から生じる行列をGMRes(m)法で解くプログラムである。

DOWNLOAD	GMRes.zip

OS:WindowsXP, Windows2000
開発環境:VC++ 2005
言語:C++

結果を表示する場合は次のようにすればよい。

結果の表示

前処理つきGMRes法

右前処理

解くべき方程式 $\b{A}\b{x}=\b{b}$ を前処理行列を使って次のように変形する。

$\b{A}\b{M}^{-1}\b{u}=\b{b}, \qquad \qquad \b{x}=\b{M}^{-1}\b{u}$

これは $\b{M}^{-1}$ が $\b{A}$ の右から掛けられているので右前処理と呼ばれる。

係数行列が $\b{A}$ から $\b{A}\b{M}^{-1}$ に変わっている。この $\b{A}\b{M}^{-1}$ の性質が単位行列に近ければ元の方程式よりも早い収束が期待される。

この前処理を使ってGMResで解を求める方法は基本的にGMRes法で $\b{A}$ の部分を $\b{A}\b{M}^{-1}$ に変えて $\b{u}$ を求め、最後に１回解く事で $\b{u}$ を求めるというだけで特別な変更はない。また残差については

$\b{r}=\b{b}-\b{A}\b{x}=\b{b}-\b{A}\b{M}^{-1}\b{u}$

であるから $\b{u}$ に対する残差はそのまま $\b{x}$ に対する残差でもある。よって逐次的に残差の大きさを更新しながら $\b{u}$ を求める際に得られる残差の大きさは $\b{x}$ に対する残差の大きさでもある。よって残差の大きさが反復途中にも得られるために右前処理法は逐次的に更新する方法によく使われる。

前処理つきGMRes法のアルゴリズム（逐次更新型）

1.   Compute $\b{r}_0=\b{b}-\b{A}\b{x}_0 \qquad \b{v}_1=\b{r}_0/||\b{r}_0||$
2.   For $k=1,2,\ldots,m\quad$ Do:
3.      $\b{w}=\b{A}\b{M}^{-1}\b{v}_k$
4.      For $i=1,2,\ldots,k\quad$ Do:
5.         $h_{k,i}=(\b{w},\b{v}_i)$
6.         $\b{w}=\b{w}-h_{k,i}\b{v}_i$
7.      End Do
8.      $h_{k,k+1}=||\b{w}||_2$
9.      $\b{v}_{k+1}=\b{w}/||\b{w}||_2$
10.     $Compute \quad \tilde{\b{h}}_k=\Omega_{k-1}\cdots\Omega_1\b{h}_k \quad and \quad make \quad \Omega_k \quad from \quad \tilde{\b{h}}_k$
11.     $Compute \quad \bar{R}_k=[\bar{R}_{k-1}',\Omega_k\tilde{\b{h}}_k] \quad and \quad \bar{\b{g}}_k = \Omega_k\bar{\b{g}}_{k-1}'$
12.     $if \quad |\bar{\b{g}}_{k,k+1}|<tolerance \quad$ Break Loop
13. End Do:
14. $Compute \quad \b{y}_k=R_k^{-1}\b{g}_k$
15. $\b{u}_k=\b{u}_0+\b{V}_k\b{y}_k$
16. $\b{x}_k=\b{M}^{-1}\b{u}_k$

サンプルプログラム

次のプログラムは次のような移流拡散問題の例

移流拡散問題の例

から生じる行列を「ILU分解による右前処理つきGMRes(m)法」で解くプログラムである。

DOWNLOAD	ILU-GMRes.zip

OS:WindowsXP
開発環境:VC++ 2005
言語:C++

結果を表示する場合は次のようにすればよい。

結果の表示

複素数行列に対するGMRes法

基本的に実数行列に対するGMRes法と手続きは同じであるが、複素行列の場合ギブンス回転のやり方のみ注意が必要である。

回転行列を $\Omega$ とすると、回転行列は内積の値を変化させないので任意の複素ベクトル $\b{a},\b{b}$ について以下が成り立つ。

$(\Omega\b{a},\Omega\b{b})=(\b{a},\Omega^*\Omega\b{b})=(\b{a},\b{b})$

つまり、エルミート行列(Hermetian matrix) $\Omega^*\Omega=\b{I}$ であることがわかる。ここで回転 $\Omega_i$ が i,i+1 平面の中での回転であるとすると、 i,i+1 成分以外は全て変化しない。つまり $\Omega$ は i,i+1 に関する項以外は単位行列であることがわかる。そこで i,i+1 に関する小行列を $\Omega_i'$ とする。つまり $\Omega$ は次のような行列

$\Omega_i=\left(\begin{array} \b{I}&& && \\ &&\Omega_i'&& \\ && &&\b{I} \end{array}\right) \begin{array} \\ \leftarrow row \quad i \quad i+1 \\ \end{array}$

このような場合、 $\Omega_i$ がエルミート行列であるためには、小行列 $\Omega_i'$ もエルミート行列でなくてはならない。

また、 $\Omega_i$ はベクトル $\Omega_k'\left(\begin{array} (\tilde{\b{h}}_k)_k \\ (\tilde{\b{h}}_k)_{k+1} \end{array}\right)$ の第２成分が０となるように選ばれた。

この２つの条件つまり、エルミート行列であるという条件、ベクトル $\tilde{\b{h}}_k$ との積の最後の成分が０であるという条件からは一意に $\Omega_i'$ は定まらない。つまり、この２つの条件を満たすならどのように $\Omega_i'$ を選んでも構わない。

例えば、ここでは $\Omega_i'$ が次のような行列であると限定して考える。

$\Omega_i'=\left(\begin{array}c_k&&\bar{s}_k\\-s_k&&c_k\end{array}\right)$

但し c_k は実数である。

すると、 $\Omega_i'$ がエルミート行列であるという条件は次のように書ける。

$c_k^2+\bar{s}_k s_k=1$

また、ベクトル $\tilde{\b{h}}_k$ との積の最後の成分が０であるという条件は次のとおり、

$-s_k(\tilde{\b{h}}_k)_k+c_k(\tilde{\b{h}}_k)_{k+1}=0$

これらを解くと c_k, s_k は次のようになる。

$c_k=\frac{||(\tilde{\b{h}}_k)_k||}{\sqrt{(||\tilde{\b{h}}_k)_k^2||+||(\tilde{\b{h}}_k)_{k+1}^2||}}\quad\quad s_k=c_k\frac{(\tilde{\b{h}}_k)_{k+1}}{(\tilde{\b{h}}_k)_k}$

例えばこのような s_k,c_k をもちいた回転行列 $\Omega_i'$ を使うことで、GMRes法を使って複素行列を解くことができる。

また、複素数の場合は $h_{k,i}=(\b{w},\b{v}_i)\ne(\b{v}_i,\b{w})$ である。間違えやすくバグに陥り易いので注意が必要である。

n.umetani@gmail.com

GMRes

目次

概要

基本的な手続き

ヘッセンベルグ行列を用いた残差の表現

GMRes法(ヘッセンベルグ行列)

上三角行列を用いた残差の表現

m=4の場合の例

一般的な場合

GMRes法(上三角行列)

上三角行列を用いた残差ノルムの最小値を求める方法

GMRes法のアルゴリズム（上三角形、後退代入）

逐次的に残差の大きさを求める方法

基底の更新

上三角行列の更新

残差の大きさの更新

GMRes法のアルゴリズム（逐次更新型）

Full-GMRes法、GMRes(m)法

サンプルプログラム

前処理つきGMRes法

右前処理

前処理つきGMRes法のアルゴリズム（逐次更新型）

サンプルプログラム

複素数行列に対するGMRes法

参考にしたもの

Books

Paper