3 反復法の基礎

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

3.1 反復法とは

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

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

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

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

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

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

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

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

$$\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)}$は，

$$\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)は

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

となる．この式(15)から式(13)を引くと， となる．

$$\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)}$とすると，

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

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

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

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

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

$$=\boldsymbol{S}^{-1}\boldsymbol{T}\boldsymbol{S}^{-1}\boldsymbol{T}\boldsymbol{S}^{-1}\boldsymbol{T}\boldsymbol{e}^{(k-2)}$$

$$=\boldsymbol{S}^{-1}\boldsymbol{T}\boldsymbol{S}^{-1}\boldsymbol{... ...^{-1}\boldsymbol{T}\cdots \boldsymbol{S}^{-1}\boldsymbol{T}\boldsymbol{e}^{(0)}$$

$$=\left(\boldsymbol{S}^{-1}\boldsymbol{T}\right)^k\boldsymbol{e}^{(0)}$$ (18)

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

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

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

$$\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
著者: 山本昌志