3 反復法の基礎

ここの説明は,文献 [1]を参考にした.これは線形代数を実際に どのように使うか--を詳細に述べた教科書で工学系の学生は一度は読んでもらいたい. 初めて私が線形代数の講義を受けたとき,あまりにも抽象的で,さっぱりわからなかった. その後,この教科書を読むことにより,なるほど線形代数は便利なものであるとやっとわ かった.

3.1 反復法とは

さて,いままで学習した直接法はしつこく計算すれば,必ず解が求まる.しかし,大きな 連立方程式を計算するには不向きである.なぜならば,ガウス・ジョルダン法の計算回数 は,方程式の次元$ n$の三乗に比例するため,大きな行列ではとたんに計算時間が必要に なるからである.

実用的なプログラムでは,非常に大きな連立方程式を計算しなくてはならない.たとえば, 私の研究室での計算でも10万元くらいは計算している.これをガウス・ジョルダン法で計 算すると膨大な時間が必要となり,現実的ではない.そこで,これよりは格段に計算の速 い反復法を用いている.ここでは,その反復法を簡単に説明する.

当然ここでも,連立方程式

$\displaystyle \boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}$ (11)

を満たす $ \boldsymbol{x}$を数値計算で求める.反復法の理論を考えるために,この連立方程式の 真の解 $ \boldsymbol{x}$とする.$ n$回目の反復計算によりで求められたものを $ \boldsymbol{x}^{(n)}$とする. そして,反復の計算回数を増やして,

$\displaystyle \lim_{n\rightarrow\infty}\boldsymbol{x}^{(n)}=\boldsymbol{x}$ (12)

になったとする.反復の計算方法を上手に選ぶと,真の解に収束させることができる.こ のように反復計算を行い真の解に収束させる方法を反復法と言う.

どのようにして反復計算をするのか?.例えば,行列 $ \boldsymbol{A}$ $ \boldsymbol{S}-\boldsymbol{T}$ と分解するだけで,反復計算の式を作成することができる.

$\displaystyle \boldsymbol{S}\boldsymbol{x}^{(k+1)}=\boldsymbol{T}\boldsymbol{x}^{(k)}+\boldsymbol{b}$ (13)

ここで, $ \boldsymbol{x}^{(k)}$ $ \boldsymbol{\alpha}$に収束するとする.すると,式 (13)と式(11)を比べれば, $ \boldsymbol{\alpha}$ $ \boldsymbol{x}$は等しいことがわかる.すなわち,式(13)で元の方程式 (11)を表した場合, $ \boldsymbol{x}^{(k)}$が収束すれば,必ず真の解 $ \boldsymbol{x}$に収束するのである.別の解に収束することはなく,真の解に収束するか,発散 するかのいずれかである.振動することはないのか?.それはよい質問である.興味があ る人が調べてみてほしい.

言うまでもないと思うが,式(13)をつかって,$ k$番めの近似解 $ \boldsymbol{x}^{(k)}$から$ k+1$番めの近似解 $ \boldsymbol{x}^{(k+1)}$は,

$\displaystyle \boldsymbol{x}^{(k+1)}=\boldsymbol{S}^{-1}\left(\boldsymbol{b}+\boldsymbol{T}\boldsymbol{x}^{(k)}\right)$ (14)

の計算により求める.この式の中には係数行列 $ \boldsymbol{A}$と非同次項の情報は入っており, 情報の過不足はない--ことに注意が必要である.ある意味ではこれは連立方程式の解の 公式と考えることもできる.もちろん,この計算のためには初期値$ x^{(0)}$は必要で, それはプログラマーあるいはユーザー適当に決めなくてはならない.

3.2 解の収束の条件

先の説明で,式(13)を使った反復法の場合, $ \boldsymbol{x}^{(k)}$の収束が重要 であることがわかった.ここでは,これが収束する条件を示す.

真の解の場合,式(13)は

$\displaystyle \boldsymbol{S}\boldsymbol{x}=\boldsymbol{T}\boldsymbol{x}+\boldsymbol{b}$ (15)

となる.この式(15)から式(13)を引くと, となる.

$\displaystyle \boldsymbol{S}(\boldsymbol{x}-\boldsymbol{x}^{(k+1)})=\boldsymbol{T}(\boldsymbol{x}-\boldsymbol{x}^{(k)})$ (16)

となる.ここで, $ \boldsymbol{x}-\boldsymbol{x}^{(k+1)}$ $ \boldsymbol{x}-\boldsymbol{x}^{(k)}$は,真の解からの差,すな わち,誤差を示している.$ k$回目の計算の誤差を $ \boldsymbol{e}^{(k)}$とすると,

$\displaystyle \boldsymbol{e}^{(k+1)}=\boldsymbol{S}^{-1}\boldsymbol{T}\boldsymbol{e}^{(k)}$ (17)

と表すことができる.この誤差ベクトル $ \boldsymbol{e}^{(k)}$がゼロに収束すれば,ハッピーなのだ.

ハッピーになるための条件を探すために,計算の最初の誤差を $ \boldsymbol{e}^{(0)}$とする.すると,

$\displaystyle \boldsymbol{e}^{(k+1)}$ $\displaystyle =\boldsymbol{S}^{-1}\boldsymbol{T}\boldsymbol{e}^{(k)}$    
  $\displaystyle =\boldsymbol{S}^{-1}\boldsymbol{T}\boldsymbol{S}^{-1}\boldsymbol{T}\boldsymbol{e}^{(k-1)}$    
  $\displaystyle =\boldsymbol{S}^{-1}\boldsymbol{T}\boldsymbol{S}^{-1}\boldsymbol{T}\boldsymbol{S}^{-1}\boldsymbol{T}\boldsymbol{e}^{(k-2)}$    
  $\displaystyle =\boldsymbol{S}^{-1}\boldsymbol{T}\boldsymbol{S}^{-1}\boldsymbol{...
...^{-1}\boldsymbol{T}\cdots \boldsymbol{S}^{-1}\boldsymbol{T}\boldsymbol{e}^{(0)}$    
  $\displaystyle =\left(\boldsymbol{S}^{-1}\boldsymbol{T}\right)^k\boldsymbol{e}^{(0)}$ (18)

となる.この式の右辺には,やっかいそうな行列の$ k$乗の計算がある.しかし, 2.3 節で得た結果を利用するとその計算も簡単である.行列 $ \boldsymbol{S}^{-1}\boldsymbol{T}$の固有値と固有ベクトルで作る行列を,$ \Lambda$ $ \boldsymbol{X}$とする と,式(18)は

$\displaystyle \boldsymbol{e}^{(k+1)}=\boldsymbol{X}\Lambda^k\boldsymbol{X}^{-1}\boldsymbol{e}^{(0)}$ (19)

となる.明らかに,計算回数$ k$を増やしていくと,誤差のベクトル $ \boldsymbol{e}^{(k)}$$ \Lambda^k$に依存する.これは,

$\displaystyle \Lambda^k=\left[ \begin{array}{@{\,}ccccc@{\,}} \lambda_1^k & & &...
...ash{\Huge$0$}}\quad} & & \ddots & \\ & & & & \lambda_n^k \\ \end{array} \right]$ (20)

となるので, $ k\rightarrow\infty$の場合,誤差 $ \boldsymbol{e}^{(k)}$がゼロに収束するためには,固有 値すべてが $ \vert\lambda_i\vert<1$でなくてはならない.そして,収束の速度は,最大の固有値 $ \max\vert\lambda_i\vert$に依存する.この絶対値が最大の固有値をスペクトル半径と言う.

ここで言いたいのは,連立方程式を式(13)の反復法で計算する場合, 結果が真の値に収束するためには,行列 $ \boldsymbol{S}^{-1}\boldsymbol{T}$の最大固有値の絶対値が1以 下でなくてはならないと言うことである.

最大固有値が1以下になる行列の条件を探すことは難しい.また,予め行列 $ \boldsymbol{S}^{-1}\boldsymbol{T}$の最大固有値を計算することも考えられるが,それもかなりの計算 量が必要で,反復法を使って計算時間を短縮するメリットが無くなってしまう.このよう なことから,反復法はとりあえず試してみて,発散するようであれば他の方法に切り替え るのが良いだろう.後で述べるSOR法の加速緩和係数$ \omega$を1以下にするという方法も ある.


ホームページ: Yamamoto's laboratory
著者: 山本昌志
Yamamoto Masashi
平成18年11月12日


no counter