数値演算法

(8) 連立方程式を解く -2-

前回は、直接法と呼ばれる連立方程式の解法を紹介しました。これらは「消去法」を応用したものであり、処理の内容も非常に素朴で分かりやすいものでした。今回は、直接法とは異なるアルゴリズムを利用して連立方程式を解く方法を紹介したいと思います。

(注) 数式などの記法について (ドキュメントの中で使用している数式の表現方法に関する注意点です)

1) ヤコビ法 ( Jacobi method )

前章で、行列とベクトルを使って連立方程式を表現できることを紹介しました。

Ax = b

上式で、A は正方行列、xb はベクトルを表します。これを、各成分を使って表すと、次のような式になります。

|a11, a12,... a1n ||x1| = |b1|
|a21, a22,... a2n ||x2||b2|
|::...:||:||:|
|an1, an2,... ann ||xn||bn|

上式の左辺は、次のように分解することができます。

|0, a12,... a1n ||x1| + |a11,0,...0||x1|
|a21, 0,... a2n ||x2||0, a22,...0||x2|
|::...:||:||::...:||:|
|an1, an2,... 0||xn||0,0,...ann||xn|

この式は、行列を対角成分とそれ以外の成分に分解した形になっています。対角成分のみを持つ行列を「対角行列 ( Diagonal Matrix )」といい、これを D で表すと、Ax = b を次のように表すことができます。

( L + U )x + Dx = b

ここで LU はそれぞれ下三角行列と上三角行列を表します。L, U を持つ項を右辺へ移項すると、上式は

Dx = b - ( L + U )x

となり、この第 r 行目の式は次のようになります。

arrxr = br - Σk{1→r-1}( arkxk ) - Σk{r+1→n}( arkxk )

これはちょうど、r 番目の未知数 xr をその他の未知数で表現した形になっています。当然、その他の未知数も値は分からないので、この式から r 番目の未知数を求めることはできませんが、その他の未知数として適当な値を決めて計算すれば、値を得ることが可能になります。この操作を繰り返したとき、i 回目の操作で利用した未知数 xr(i)の値を使って計算した結果は xr(i+1)と表され、両者の関係式は次のようになります。

arrxr(i+1)=br - Σk{1→r-1}( arkxk(i) ) - Σk{r+1→n}( arkxk(i) )
=br - Σk{k≠r}( arkxk(i) )

ここで、xr(i) と xr(i+1) が任意の r について等しくなれば、xr(i) が求める解となります。このような、仮の解を使って実際の解を求める操作は「反復法 ( Iterative Method )」と呼ばれ、その中で、上記のような操作を「ヤコビ法 ( Jacobi Method )」と言います。

次の連立方程式を、ヤコビ法を使って実際に解いてみます。

4x +  y - 2z = 6   ... (1)
 x + 6y + 3z = -2  ... (2)
2x +  y + 9z = -7   ... (3)

それぞれの式を変形して以下のように表します。

x = {  6 - (  y - 2z ) } / 4 ... (1')
y = { -2 - (  x + 3z ) } / 6 ... (2')
z = { -7 - ( 2x +  y ) } / 9 ... (3')

はじめに、仮の解を ( x(0), y(0), z(0) ) = ( 0, 0, 0 ) とすると、

x(1)=3/2
y(1)=-1/3
z(1)=-7/9

この解を使ってさらに ( x(2), y(2), z(2) ) を求めると

x(2)=[ 6 - { (-1/3) - 2×(-7/9) } ] / 4 = 43/36
y(2)=[ -2 - { 3/2 + 3×(-7/9) } ] / 6 = -7/36
z(2)=[ -7 - { 2×(3/2) + (-1/3) } ] / 9 = -29/27

この操作を繰り返すと、次のようになります。

表 1-1. ヤコビ法による収束の様子
ix(i)y(i)z(i)
0000
11.5-0.333333-0.777778
21.19444-0.194444-1.07407
31.011570.004630-1.02160
40.9880400.008873-1.00309
50.9962380.003537-0.998328
60.999952-0.000209-0.999557
71.00027-0.000213-0.999966
81.00007-0.000063-1.00004
90.9999970.000007-1.00001
100.9999940.000005-1
110.9999990.000001-0.999999
1210.000000-1

解はだんだんとある値に収束し、最終的には ( x, y, z ) = ( 1, 0, -1 ) となる様子が上記の表から分かると思います。

さて、この操作によって、解は常に収束するのでしょうか。他の連立方程式で試してみたいと思います。

 x + 7y - 8z = 9 ... (1)
9x + 2y + 4z = 5 ... (2)
6x +  y +  z = 5 ... (3)

この方程式の解は ( x, y, z ) = ( 1, 0, -1 ) になりますが、ヤコビ法を使って処理すると次のように解が変化します。

表 1-2. ヤコビ法により解が発散する例
ix(i)y(i)z(i)
0000
192.55
231.5-48-51.5
3-67-36.25-136
4-825.25576443.25
5-4772829.624380.5
615245.6-661237.375
746592-68677.6-84856.8
8-198102-39948-210869
9-1.40731e+061.31320e+061.22857e+06
106361273.87577e+067.13067e+06
112.99150e+07-1.71239e+07-7.69253e+06
125.83271e+07-1.19232e+08-1.62366e+08

何回繰り返しても解は収束せず、それどころか発散してしまいます。このように、常に収束するとは限らないため、処理回数に対して制限を掛ける必要があります。

以下に、ヤコビ法のサンプル・プログラムを示します。

/**
   ConvergenceTest_BySum 差分の和と全体の和の比率による収束判定

   反復処理において、前の値 p と今回の値 c の差分の絶対値 | c - p | と c のそれぞれに対して和をを求める。
   一回の処理で得られる和 Δs = Σ| c - p | と S = Σc の比率 Δs / S を求め、その値が指定したしきい値以下か
   どうかを判定する。但し、S = 0 であれば判定結果は常に true となる。

   一回の反復処理で扱う値がただ一つであれば、ConvergenceTest_BySum と ConvergenceTest_ByMax の結果は一致する。
**/
template< class T >
class ConvergenceTest_BySum
{
  T e_;    // 判定のためのしきい値
  T diff_; // 差分の和
  T sum_;  // 全体の和

 public:

  /*
    しきい値 e を指定して構築
  */
  ConvergenceTest_BySum( T e ) : e_( e ) { init(); }

  /*
    求めた総和を初期化(ゼロクリア)する
  */
  void init()
  { diff_ = sum_ = T(); }

  /*
    前後の値 preVal, newVal から差分と総和を求める
  */
  void operator()( T preVal, T newVal )
  {
    diff_ += std::abs( preVal - newVal );
    sum_ += std::abs( newVal );
  }

  /*
    収束しているかを判定する
  */
  bool isConvergent()
  { return( sum_ == 0 || ( diff_ / sum_ <= e_ ) ); }
};

/**
   ConvergenceTest_ByMax 差分と値の比率の最大値による収束判定

   反復処理において、前の値 p と今回の値 c の差分 c - p と c との比率の絶対値 | ( c - p ) / c |
   を求める ( 但し c = 0 ならば | c - p | をそのまま用いる )。一回の処理で得られる全ての値の中から
   その最大値を抽出し、その値が指定したしきい値以下かどうかを判定する。

   一回の反復処理で扱う値がただ一つであれば、ConvergenceTest_BySum と ConvergenceTest_ByMax の結果は一致する。
**/
template< class T >
class ConvergenceTest_ByMax
{
  T e_;    // 判定のためのしきい値
  T max_;  // 差分

 public:

  /*
    しきい値 e を指定して構築
  */
  ConvergenceTest_ByMax( T e ) : e_( e ) { init(); }

  /*
    求めた最大値を初期化(ゼロクリア)する
  */
  void init()
  { max_ = T(); }

  /*
    前後の値 preVal, newVal から差分との比率を求め、前回の最大値と比較してより大きければ入れ替える
  */
  void operator()( T preVal, T newVal )
  {
    // 新たな値による比率
    T ratio = std::abs( preVal - newVal );
    if ( newVal != 0 ) ratio /= std::abs( newVal );

    if ( max_ < ratio )
      max_ = ratio;
  }

  /*
    収束しているかを判定する
  */
  bool isConvergent()
  { return( max_ <= e_ ); }
};

/**
   LES_Jacobi ヤコビ法による連立方程式の計算
**/
template< class T, template< class T > class Test = ConvergenceTest_BySum >
class LES_Jacobi
{
  Test< T > test_; // 収束判定
  unsigned int maxCnt_; // 収束しなかった場合の最大処理回数

 public:

  /*
    しきい値 e と最大処理回数 maxCnt を指定して構築
  */
  LES_Jacobi( T e, unsigned int maxCnt )
  : test_( e ), maxCnt_( maxCnt ) {}

  /*
    連立方程式 les を解く

    計算回数が最大処理回数 maxCnt を超えた場合は例外 ExceptionExcessLimit を投げる。
     但し、それまでに得られた近似解は連立方程式 les の解ベクトルに返される。
  */
  void operator()( LinearEquationSystem< T >& les );
};

/*
  LES_Jacobi< T, Test >::operator() : Jacobi法による連立方程式 les の計算

  計算回数が maxCnt_ を超えた場合は例外 ExceptionExcessLimit を投げる。
*/
template< class T, template< class T > class Test >
  void LES_Jacobi< T, Test >::operator()( LinearEquationSystem< T >& les )
{
  typedef typename LinearEquationSystem< T >::size_type size_type;

  assert( &les != 0 );

  std::vector< T > ans( les.size(), T() ); // 仮の解
  unsigned int cnt = 1;                    // 処理回数

  // 対角成分が 0 なら処理失敗
  for ( size_type i = 0 ; i < les.size() ; ++i )
    if ( les[i][i] == 0 ) throw ExceptionUnsolvedEquation();

  while ( cnt - 1 < maxCnt_ ) {
    std::vector< T > buff( les.size() ); // 次に計算する近似解

    // 近似解の計算
    for ( size_type i = 0 ; i < les.size() ; ++i ) {
      buff[i] = les.ans( i );
      for ( size_type j = 0 ; j < les.size() ; ++j ) {
        if ( i == j ) continue;
        buff[i] -= les[i][j] * ans[j];
      }
      buff[i] /= les[i][i];
    }

    // 収束判定
    test_.init();
    for ( typename std::vector< T >::size_type i = 0 ; i < les.size() ; ++i ) {
      test_( ans[i], buff[i] );
      ans[i] = buff[i];
    }
    if ( test_.isConvergent() ) break;
    ++cnt;
  }

  // 求めた解を代入
  for ( typename std::vector< T >::size_type i = 0 ; i < ans.size() ; ++i )
    les.ans( i ) = ans[i];

  if ( cnt - 1 >= maxCnt_ )
    throw ExceptionExcessLimit( maxCnt_ );
}

解が収束したのかを判定するための関数オブジェクトとして、前回求めた解との差の絶対値について和を求め、それを解の絶対値の和で割った商をしきい値 e と比較することで行う ConvergenceTest_BySum と、差の絶対値を解の絶対値で割った商をそれぞれ求め、その中の最大値で判定する ConvergenceTest_ByMax の二通りのやり方が用意されていて、どちらを利用するかはテンプレート引数で選択ができるようにしてあります。
ヤコビ法を利用した連立方程式の計算は LES_Jacobi が行います。テンプレート引数としては、値の型を表す T の他に、前述の収束判定用関数オブジェクトを表す Test があり、デフォルトは ConvergenceTest_BySum にしてあります。なお、ConvergenceTest_BySum もテンプレートクラスなので、テンプレート引数には "template< class T > class Test" と指定されていることに注意してください。コンストラクタの引数 e は収束判定を行うためのしきい値で、値を小さくするほど解の精度は向上します。また、maxCnt は計算を最大何回実行するかを指定するためのもので、解が収束しない場合に無限ループになるのを防ぐ役割を持っています。近似解を計算する部分 operator() は、先ほど説明した計算をそのまま実装しただけのものなので容易に理解できると思います。

対角成分が 0 の場合、処理の途中でゼロ除算が発生します。サンプル・プログラムでは対角成分が 0 の場合すぐに例外を投げていますが、他の行と入れ替えて 0 にならないようにすれば処理を継続させることができます。反復法では、対角成分の値が他の成分より大きいほど解が収束しやすくなります ( 詳細は後述します )。そのため、ピボット選択を使って対角成分が大きくなるようにすると収束までの回数が少なくなることが期待できます。しかし、ピボット選択の処理は重く、かえって反復法の処理時間が長くなってしまうので、ここでは利用していません。


2) ガウス=ザイデル法(Gauss-Seidel method)

ヤコビ法は、前に得た仮の解全てを使って次の解を求めていました。しかし、最上位の行から順に解を求めることで、前回よりも ( もし収束していれば ) 精度の良い解が部分的に得られることになるので、これを利用すれば、さらに収束を速めることが期待できます。この時の漸化式は、

arrxr(i+1) = br - Σk{1→r-1}( arkxk(i+1) ) - Σk{r+1→n}( arkxk(i) )

と表すことができます。この方法を「ガウス=ザイデル法 ( Gauss-Seidel Method )」と言います。ヤコビ法の中で例として使った以下の連立方程式をガウス=ザイデル法で処理したときの近似解の推移を、ヤコビ法の場合と比較した結果を以下に示します。

4x +  y - 2z = 6  ... (1)
 x + 6y + 3z = -2 ... (2)
2x +  y + 9z = -7 ... (3)
表 2-1. ヤコビ法とガウス=ザイデル法の収束比較
iヤコビ法ガウス=ザイデル法
x(i)y(i)z(i)x(i)y(i)z(i)
0000000
11.5-0.333333-0.7777781.5-0.583333-1.04630
21.19444-0.194444-1.074071.122690.002701-1.02756
31.011570.004630-1.021600.9855430.016191-0.998586
40.9880400.008873-1.003090.996659-0.000150-0.999241
50.9962380.003537-0.9983281.00042-0.000449-1.00004
60.999952-0.000209-0.9995571.000090.000006-1.00002
71.00027-0.000213-0.9999660.9999880.000012-0.999999
81.00007-0.000063-1.000040.9999980.000000-0.999999
90.9999970.000007-1.0000110.000000-1
100.9999940.000005-110.000000-1
110.9999990.000001-0.99999910.000000-1
1210.000000-110.000000-1

この結果を見る限りは、ガウス=ザイデル法の方がヤコビ法よりも早く収束しています。

ガウス=ザイデル法のサンプル・プログラムを以下に示します。

/**
   LES_GaussSeidel : ガウス=ザイデル法による連立方程式の計算
**/
template< class T, template< class T > class Test = ConvergenceTest_BySum >
class LES_GaussSeidel
{
  Test< T > test_; // 収束判定
  unsigned int maxCnt_; // 収束しなかった場合の最大処理回数

 public:

  /*
    しきい値 e と最大処理回数 maxCnt を指定して構築
  */
  LES_GaussSeidel( T e, unsigned int maxCnt )
  : test_( e ), maxCnt_( maxCnt ) {}

  /*
    連立方程式 les を解く

    計算回数が最大処理回数 maxCnt を超えた場合は例外 ExceptionExcessLimit を投げる。
    但し、それまでに得られた近似解は連立方程式 les の解ベクトルに返される。
  */
  void operator()( LinearEquationSystem< T >& les );
};

/*
  LES_GaussSeidel< T, Test >::operator() : ガウス=ザイデル法による連立方程式 les の計算

  計算回数が maxCnt_ を超えた場合は例外 ExceptionExcessLimit を投げる。
*/
template< class T, template< class T > class Test >
  void LES_GaussSeidel< T, Test >::operator()( LinearEquationSystem< T >& les )
{
  typedef typename LinearEquationSystem< T >::size_type size_type;

  assert( &les != 0 );

  std::vector< T > ans( les.size(), T() ); // 仮の解
  unsigned int cnt = 1;                    // 処理回数

  // 対角成分が 0 なら処理失敗
  for ( size_type i = 0 ; i < les.size() ; ++i )
    if ( les[i][i] == 0 ) throw ExceptionUnsolvedEquation();

  while ( cnt - 1 < maxCnt_ ) {
    test_.init();
    for ( size_type i = 0 ; i < les.size() ; ++i ) {
      T preAns = ans[i]; // 前回の解を保持しておく

      // 近似解の計算
      ans[i] = les.ans( i );
      for ( size_type j = 0 ; j < les.size() ; ++j ) {
        if ( i == j ) continue;
        ans[i] -= les[i][j] * ans[j];
      }
      ans[i] /= les[i][i];

      // 前回の解との差を求める
      test_( preAns, ans[i] );
    }
    if ( test_.isConvergent() )
      break;
    ++cnt;
  }

  // 求めた解を代入
  for ( size_type i = 0 ; i < ans.size() ; ++i )
    les.ans( i ) = ans[i];

  if ( cnt - 1 >= maxCnt_ )
    throw ExceptionExcessLimit( maxCnt_ );
}

計算方法はヤコビ法と変わりません。大きな違いは、前回求めた解を保持するための領域 ( ヤコビ法での buff ) が不要になる部分です。求めた解はそのまま順次計算に利用していくので、前回求めた解を計算のために残しておく必要はなくなります。
収束が速い上に領域も半分で済むため、通常はガウス=ザイデル法を利用した方が有利です。しかし、それぞれの計算処理を並列に行うような場合には利用できないので、そのような場合はヤコビ法が用いられるようです。


3) SOR法 ( Successive Over-Relaxation )

ヤコビ法の漸化式をもう一度以下に示します。

xr(i+1) = { br - Σk{k≠r}( arkxk(i) ) } / arr

前回求めた仮の解との差を er(i) で表したとき、関係式 xr(i+1) = xr(i) + er(i+1) が得られますが、er(i) に対して適当な値を掛けると、仮の解は真の解に対してより収束する可能性があります ( 逆に収束が遅くなったり発散する可能性もあります )。このように、差 er(i) を適当な値で補正して収束を速くする方法を「逐次過緩和法 ( Successive Over-Relaxation Method ; SOR 法 )」と言います。

具体的には、ヤコビ法を使って求めた解を前回求めた解の値から減算して差を求め、その値に補正値を掛けてから前の解に加算することによって次の解を求めます。式に表すと、次のようになります。

er(i+1) = { br - Σk{k≠r}( arkxk(i) ) } / arr - xr(i)

xr(i+1) = xr(i) + ωer(i+1)

上式で、ω が補正値 ( 緩和係数と呼ばれます ) を表しています。問題は、ω としてどのような値が最適なのかということですが、ω = 1 の場合はヤコビ法と等しく、1 より小さくなると収束が逆に遅くなることは明らかです。しかし、あまり大きな値にすると今度は発散してしまいます ( 理論上、ω > 2 のときは発散するそうです )。方程式の数などによっても、ω の最適値は変化するようです。

SOR 法のサンプル・プログラムを以下に示します。

/**
   LES_Jacobi_SOR : ヤコビ法に SOR 法を適用した連立方程式の計算
**/
template< class T, template< class T > class Test = ConvergenceTest_BySum >
class LES_Jacobi_SOR
{
  T w_;                 // 緩和係数
  Test< T > test_;      // 収束判定
  unsigned int maxCnt_; // 収束しなかった場合の最大処理回数

 public:

  /*
    緩和係数 w、しきい値 e、最大処理回数 maxCnt を指定して構築
  */
  LES_Jacobi_SOR( T w, T e, unsigned int maxCnt )
  : w_( w ), test_( e ), maxCnt_( maxCnt ) {}

  /*
    連立方程式 les を解く

    ・計算回数が最大処理回数 maxCnt を超えた場合は例外 ExceptionExcessLimit を投げる。
      但し、それまでに得られた近似解は連立方程式 les の解ベクトルに返される。
  */
  void operator()( LinearEquationSystem< T >& les );
};

/*
  LES_Jacobi_SOR< T, Test >::operator() : Jacobi法にSOR法を適用

  les : 対象の連立方程式

  計算回数が maxCnt_ を超えた場合は例外 ExceptionExcessLimit を投げる。
*/
template< class T, template< class T > class Test >
  void LES_Jacobi_SOR< T, Test >::operator()( LinearEquationSystem< T >& les )
{
  typedef typename LinearEquationSystem< T >::size_type size_type;

  assert( &les != 0 );

  std::vector< T > ans( les.size(), T() ); // 仮の解
  unsigned int cnt = 1;                    // 処理回数

  // 対角成分が 0 なら処理失敗
  for ( size_type i = 0 ; i < les.size() ; ++i )
    if ( les[i][i] == 0 ) throw ExceptionUnsolvedEquation();

  while ( cnt - 1 < maxCnt_ ) {
    std::vector< T > buff( les.size() ); // 次に計算する近似解

    // 近似解の計算
    for ( size_type i = 0 ; i < les.size() ; ++i ) {
      T t = les.ans( i );
      for ( size_type j = 0 ; j < les.size() ; ++j ) {
        if ( i == j ) continue;
        t -= les[i][j] * ans[j];
      }
      t = t * w_ / les[i][i];
      buff[i] = ans[i] * ( 1.0 - w_ ) + t;
    }

    // 収束判定
    test_.init();
    for ( typename std::vector< T >::size_type i = 0 ; i < buff.size() ; ++i ) {
      test_( ans[i], buff[i] );
      ans[i] = buff[i];
    }
    if ( test_.isConvergent() )
      break;
    ++cnt;
  }

  // 求めた解を代入
  for ( size_type i = 0 ; i < ans.size() ; ++i )
    les.ans( i ) = ans[i];

  if ( cnt - 1 >= maxCnt_ )
    throw ExceptionExcessLimit( maxCnt_ );
}

/**
   LES_GaussSeidel_SOR : ガウス=ザイデル法に SOR 法を適用した連立方程式の計算
**/
template< class T, template< class T > class Test = ConvergenceTest_BySum >
class LES_GaussSeidel_SOR
{
  T w_;                 // 緩和係数
  Test< T > test_;      // 収束判定
  unsigned int maxCnt_; // 収束しなかった場合の最大処理回数

 public:

  /*
    緩和係数 w、しきい値 e、最大処理回数 maxCnt を指定して構築
  */
  LES_GaussSeidel_SOR( T w, T e, unsigned int maxCnt )
    : w_( w ), test_( e ), maxCnt_( maxCnt ) {}

  /*
    連立方程式 les を解く

    ・計算回数が最大処理回数 maxCnt を超えた場合は例外 ExceptionExcessLimit を投げる。
      但し、それまでに得られた近似解は連立方程式 les の解ベクトルに返される。
  */
  void operator()( LinearEquationSystem< T >& les );
};

/*
  LES_GaussSeidel_SOR< T, Test >::operator() : ガウス=ザイデル法にSOR法を適用

  les : 対象の連立方程式

  計算回数が maxCnt_ を超えた場合は例外 ExceptionExcessLimit を投げる。
*/
template< class T, template< class T > class Test >
  void LES_GaussSeidel_SOR< T, Test >::operator()( LinearEquationSystem< T >& les )
{
  typedef typename LinearEquationSystem< T >::size_type size_type;

  assert( &les != 0 );

  std::vector< T > ans( les.size(), T() ); // 仮の解
  unsigned int cnt = 1;                    // 処理回数

  // 対角成分が 0 なら処理失敗
  for ( size_type i = 0 ; i < les.size() ; ++i )
    if ( les[i][i] == 0 ) throw ExceptionUnsolvedEquation();

  while ( cnt - 1 < maxCnt_ ) {
    test_.init();
    for ( size_type i = 0 ; i < les.size() ; ++i ) {
      T preAns = ans[i]; // 前回の解を保持しておく

      // 近似解の計算
      T t = les.ans( i );
      for ( size_type j = 0 ; j < les.size() ; ++j ) {
        if ( i == j ) continue;
        t -= les[i][j] * ans[j];
      }
      t = t * w_ / les[i][i];
      ans[i] = ans[i] * ( 1.0 - w_ ) + t;

      // 前回の解との差を求める
      test_( preAns, ans[i] );
    }
    if ( test_.isConvergent() )
      break;
    ++cnt;
  }

  // 求めた解を代入
  for ( size_type i = 0 ; i < ans.size() ; ++i )
    les.ans( i ) = ans[i];

  if ( cnt - 1 >= maxCnt_ )
    throw ExceptionExcessLimit( maxCnt_ );
}

SOR 法は、ヤコビ法だけでなくガウス=ザイデル法にも適用可能です。計算処理はそのままで、最後に前の解との差分を計算してから緩和係数を掛け、それを前の解に加算する部分だけが追加されているだけです。精度を向上させるため、先ほど説明した計算式はそのままの形で適用していないことに注意してください。

収束が最も早くなるような最適な緩和係数は、連立方程式の数だけでなく要素によっても変化します。上記プログラムでテストした結果、緩和係数 w の値として 1 より小さな値である場合に収束が早くなる場合もありました ( 通常は緩和係数が 1 より大きい場合に収束が早くなるようです)。式の数が少ない場合は特にその傾向が大きいようです。式の数が 100 程度になると、緩和係数が 1 より大きい方が収束は早くなる傾向にあります。また、緩和係数を小さくすると、通常では発散してしまうような連立方程式を解くことができるようになる場合があります。


4) 反復法の収束条件

ヤコビ法では、連立方程式を Dx = b - ( L + U )x の形に変形して漸化式を作り出し、これを利用して解を求めていました。もっと一般的に、連立方程式 Ax = b は、A を任意の二つの行列 S, T の差 S - T に分解して、次のように表すことができます。

Sx = Tx + b ... (1)

右辺に仮の解 x' を代入して左辺の解 x を求める操作を繰り返し、x' = x を満たせばその時の解 x が真の解になるというのが反復法の考え方でした。
k 回目に得られた仮の解を x(k) としたとき、x(k)x(k+1) の間には次のような関係式が成り立ちます。

Sx(k+1) = Tx(k) + b ... (2)

(1) 式の x は真の解を表しているとします ( その時は当然、(1)式の等号は成り立ちます )。(1) - (2) を求めると、

S( x - x(k+1) ) = T( x - x(k) )

ここで、x - x(k) は k 番目の仮の解と真の解との誤差を表すので、これを e(k) とすると

Se(k+1) = Te(k)

両辺に左側から S の逆行列を掛けると

e(k+1) = S-1Te(k) ... (3)

になります。右辺にある e(k) は、e(k-1) の間に (3) 式と同じ関係が成り立つので、次のように式を変形することができます

e(k+1) = S-1T( S-1Te(k-1) ) = (S-1T)2e(k-1)

この操作を繰り返すと、右辺の誤差項は最終的に e(0) ( つまり最初に決めた仮の解と真の解との誤差 ) になり、以下のような関係式が得られることになります。

e(k+1) = (S-1T)k+1e(0)

上式から、反復法とは、e(0)S-1T のべき乗を掛け合わせることによって、全ての成分が 0 になるようにする手法であることが分かります。そのためには、S-1T のべき乗の全成分が 0 に収束しなければなりません。まずはこれが、反復法で解が収束するために必要な条件となります。

以下、対角行列・下三角行列・上三角行列をそれぞれ D, L, U で表します。但し、LU は対角成分もゼロであるとします。任意の行列 A は、D, L, U の三つの行列の和で表すことができます。

D=|a11,0,...0|  L=|0,0,...0|  U=|0,a12,...a1n|
|0,a22,...0||a21,0,...0||0,0,...a2n|
|::...:||::...:||::...:|
|0,0,...ann||an1,an2,...0||0,0,...0|

ヤコビ法の場合、S = D, T = -L - U なので、S-1T = -D-1( L + U ) です。ガウス=ザイデル法では L がすでに求められた仮の解になるので、次のような関係式になります。

Dx(k+1) = b - Lx(k+1) - Ux(k)

右辺にある x(k+1) の項を移項して、

( D + L )x(k+1) = b - Ux(k)

よって、S-1T = -( D + L )-1U が得られます。

反復法に SOR 法を適用した場合、

x(k+1)=x(k) + ω{ S-1( Tx(k) + b ) - x(k) }
={ ( 1 - ω )E + ωS-1T }x(k) + ωS-1b

ヤコビ法に上式を適用すると、

x(k+1) = [ ( 1 - ω )E - ωD-1( L + U ) ]x(k) + ωD-1b

両辺に左側から D / ω を掛けて整理すると

(D/ω)x(k+1)={ ( 1/ω - 1 )D - ( L + U ) }x(k) + b
={ D/ω - ( D + L + U ) }x(k) + b
=( D/ω - A )x(k) + b

よって、ヤコビ法に SOR 法を適用した場合は S-1T = ωD-1( D/ω - A ) = E - ωD-1A になります。ω = 1 のときは、

S-1T=D-1( D - A )
=-D-1( L + U )

となるので、ヤコビ法での式と一致します。

ガウス=ザイデル法に SOR 法を適用した場合、

x(k+1)=x(k) + ω{ D-1( b - Lx(k+1) - Ux(k) ) - x(k) }
=ωD-1b + { ( 1 - ω )E - ωD-1U }x(k) - ωD-1Lx(k+1)

右辺にある x(k+1) の項を移項して、

( E + ωD-1L )x(k+1) = { ( 1 - ω )E - ωD-1U }x(k) + ωD-1b

両辺に左側から D / ω を掛けると

( D/ω + L )x(k+1) = { ( 1/ω - 1 )D - U }x(k) + b

ここで、

G(ω) = D/ω + L

としたとき、

( 1/ω - 1 )D - U=D/ω - ( D + U )
=( D/ω + L ) - ( D + U + L )
=G(ω) - A

よって、

G(ω)x(k+1) = { G(ω) - A }x(k) + b

より、ガウス=ザイデル法に SOR 法を適用した場合は S-1T = G(ω)-1{ G(ω) - A } = E - G(ω)-1A になります。ω = 1 のときは、G(1) = D + L より

S-1T=( D + L )-1( D + L - A )
=( D + L )-1( -U )

となるので、ガウス=ザイデル法での式と一致します。

以上をまとめると、各反復法が収束するためには以下の行列の積が零行列に収束することが必要になります。

反復法S-1T
ヤコビ法-D-1(L + U)
ガウス=ザイデル法-( D + L )-1U
ヤコビ + SOR法E - ωD-1A
ガウス=ザイデル + SOR法E - G(ω)-1A ( 但し、G(ω) = D/ω + L )

5) 固有値と固有ベクトル

連立方程式 Ax = b は、行列 A によってベクトル xb に変換されたことを表しているものと見ることができます。このような変換を「一次変換(Linear Transformation)」といいます。行列 A の各列をベクトル ai = ( a1i, a2i, ... ani )T で表すと、一次変換は次のように表すこともできます。

x1a1 + x2a2 + ... + xnan = b

左辺のように、ベクトル ai と数値 ( スカラー ) xi を使って和の形で表したものを、xi を係数とする「線形結合 ( Linear Combination )」と言います。また、ベクトルの中のいずれかが、残る他のベクトルの線形結合で表せるとき「線形従属 ( Linearly Dependent )」、そうでないときは「線形独立 ( Linearly Independent )」であると言います。
ベクトルが線形独立であるためには、x1a1 + x2a2 + ... + xnan = 0 となるのが、全ての係数 xi が 0 である場合に限ることが必要十分条件となります。もし、ある係数 xi が 0 でないとすると、

ai = -( x1a1 + ... + xi-1ai-1 + xi+1ai+1 + ... + xnan ) / xi

のように ai が他のベクトルの線形結合で表せてしまい、逆に ai がこのように表せる場合、与式に対して係数 xi は 0 でなくなります ( このことから、線形独立なベクトルは 0 でないことも分かります)。

ベクトル ai が線形独立であるならば、任意のベクトル b は一通りの線形結合で表せます。もし、二通りの線形結合で次のように表せたとすると、

x1a1 + x2a2 + ... + xnan = b
y1a1 + y2a2 + ... + ynan = b

二式を引くと

( x1 - y1 )a1 + ( x2 - y2 )a2 + ... + ( xn - yn )an = 0

ベクトル ai は線形独立なので、各係数は 0 でなくてはならず、xi = yi となります。このことは、ベクトル ai が線形独立ならば、連立方程式の解が一意に決まることを意味しています。行列 A の各列ベクトル ai の中で線形独立なものの個数をその行列の「ランク」と言うため、n 行 n 列の行列 A で表された連立方程式 Ax = b が唯一の解を持つための条件は、次の条件と同値であるということになります。

  1. A のランクが n である ( 各列ベクトルが全て線形独立である )。
  2. 行列式 |A| が 0 でない。
  3. 逆行列 A-1 が存在する。

ここで、一次変換 Ax = b によって、bx ( ≠ 0 ) の定数倍 λ で表されるような場合を考えます。式に書くと、次のようになります。

Ax = λx

この時のスカラー λ を行列 A の「固有値 ( Eigenvalue )」、ベクトル x を「固有ベクトル ( Eigenvector )」と呼びます。上式は次のように変形できます。

λx - Ax=0
λEx - Ax=0
( λE - A )x=0

ここで、E は単位行列を表しています。これも連立方程式を表しており、解は明らかに x = 0 です。しかし、x0 としているため、それ以外に解を持たなければなりません。ここで、連立方程式が唯一の解を持つのは係数の行列の行列式が 0 でないときなので、x = 0 以外に解を持つためには

| λE - A | = 0

でなければならないことになります。この式の左辺は、λ を変数とする多項式となります。これを、行列 A の「固有方程式 ( Characteristic Polynomial )」と言います。行列 A の行列数が n のとき、固有方程式の解は最大で n 個の解を持ち、それぞれの固有値に対して固有ベクトルを求めることができます。

例として、次の行列の固有値と固有ベクトルを求めてみます。

| 6, 4 |
| 1, 3 |

固有方程式は次のようになります。

| λ - 6     -4 | = 0
|     -1 λ - 3 |     より
( λ - 6 )( λ - 3 ) - ( -4 )( -1 )=λ2 - 9λ + 14
=( λ - 2 )( λ - 7 )
=0

よって、固有値 λ = 2, 7 になります。これらを連立方程式に代入すると、

λ = 2 のとき
|-4,-4||x1| = |0|
|-1,-1||x2||0|

この連立方程式は不定になりますが、解の一つとして ( x1, x2 ) = ( 1, -1 )があります。

λ = 7 のとき
|1,-4||x1| = |0|
|-1,4||x2||0|

こちらは解の一つとして ( x1, x2 ) = ( 4, 1 ) があります。

従って、上記行列の固有値は 2 と 7 で、それぞれの固有ベクトルは ( 1, -1 ) と( 4, 1 ) になります。

前述の通り、一次変換は任意のベクトルを別の成分のベクトルに変換します。n 行 n 列の行列は、n 次元空間にあるベクトルの向きと大きさを変化させますが、固有ベクトルに対しては大きさを変化させるだけで、向きを変えることはありません。

図 5-1. 行列によるベクトルの一次変換の様子
行列による図形の変換

上の図は、固有ベクトル上に頂点を持つ四角形が行列によって変換された時の様子を示したものです。中央の塗りつぶされた長方形が平行四辺形に変換されていますが、頂点は固有ベクトル上から外れていません。次元が多くなってもこの性質は保たれ、各固有ベクトル上の座標は向きを変えることなく大きさだけ固有値倍されます。

次に、行列 A の固有ベクトル x1, x2, ... xn を列ベクトルとする行列を X = ( x1 x2 ... xn ) としたとき、行列の積 AX がどうなるのかを見てみます。

AX=A( x1 x2 ... xn )
=( Ax1 Ax2 ... Axn )
=( λ1x1 λ2x2 ... λnxn )
=( x1 x2 ... xn )|λ1,0,...0|
|0,λ2,...0|
|::...:|
|0,0,...λn|

前の例で、実際に計算してみます。

X=|1,4|
|-1,1| より
AX=|6,4||1,4|
|1,3||-1,1|
=|2,28|
|-2,7|
=|1,4||2,0|
|-1,1||0,7|
= X|2,0|
|0,7|

固有値を対角成分とする行列を D としたとき、AX = XD が成り立つので、両辺に右側から X-1 を掛けると A = XDX-1 となります。これを「行列の対角化」といいます。但し、行列の対角化が可能になるには X の逆行列が存在する必要があるため、各固有ベクトルは線形独立の関係でなくてはなりません。
固有値が全て異なれば、固有ベクトルは必ず線形独立になります。これは、以下のように証明することができます。

行列 A の固有ベクトルの一つ xk が、線形独立な他の固有ベクトルによって線形結合の形に表せると仮定する。すなわち

xk = Σi( cixi )

左辺に左側から行列 A を掛けると、

Axk = Σi( ciAxi )

固有ベクトル xi に対する固有値を λi とすると、上式は

Axk = Σi( ciλixi )

しかし、Axk = λkxk でもあるので、

Axk = λkxk = Σi( ciλkxi )

二式の差を取ると

Σi( ci( λi - λk )xi ) = 0

xi は線形独立なので、各固有値 λi が互いに異なる場合、上式が成り立つためには ci = 0 でなければならない。しかし、この時 xk = 0 となり、xk は固有ベクトルとならない。

重解が存在する場合、固有値は行列の行列数よりも少なくなります。このときでも、線形独立な固有ベクトルを行列数と等しい数だけ得ることができれば、X は逆行列を持つことができます。しかし、例えば 2 行 2 列の行列に対して固有値が一つしかなければ、直線の方程式 ax + by = 0 一つから固有ベクトルを求めなければならず、固有ベクトルはただ一つのみになります。3 行 3 列の場合でも、固有方程式が三重解になって固有値が一つしかない場合は、平面を表す方程式 ax + by + cz = 0 一つから固有ベクトルを求めることになるため、線形独立な固有ベクトルは二つまでしか求めることができません。このように、重解がある場合、対角化が常に可能であるとは限らないことになります。

対角化が可能な場合、A = XDX-1 から、A のべき乗 Ak は次のようになります。

Ak=( XDX-1 )k
=XDX-1XDX-1 ... XDX-1
=XDkX-1

Dk は、対角成分が k 乗されるだけで他の成分は 0 のままです。固有値を k 乗するだけで求めることができるので、X の逆行列さえ分かればべき乗は簡単に計算することができます。さらに、もし固有値の絶対値の最大値が 1 より小さい場合、べき乗は 0 に近づくことになります。反復法における収束条件の式は、以下のようなものでした。

e(k+1) = (S-1T)k+1e(0)

S-1T の固有値からなる対角行列を D、固有ベクトルを列成分とする行列を X としたとき、( S-1T )k+1 = XDk+1X-1 なので、S-1T の固有値 λ の全てに対して |λ| < 1 が成り立てば、反復法において解は収束します。絶対値が最大になる固有値のことを「スペクトル半径 ( Spectral Radius )」といいます。

残念ながら、固有値を求めるためには行列式を求める必要があります。しかも、行列式の中には未知数 λ が含まれており、代数方程式を解かなければ固有値を得ることができません。反復法によって誤差が収束する条件を満たす行列をあらかじめ探すことは難しく、処理時間にも影響を及ぼすので、試しに計算してみて発散するようであれば他の方法に切り替えるようにするのが通常のやり方のようです。しかしヤコビ法の場合、行列が対角優位であるならば必ず収束することを示すことができます。

A=|a11,a12,...a1n|
|a21,a22,...a2n|
|::...:|
|an1,an2,...ann| のとき、
 
S=|a11,0,...0|  T=|0,-a12,...-a1n|
|0,a22,...0||-a21,0,...-a2n|
|::...:||::...:|
|0,0,...ann||-an1,-an2,...0|

対角行列の逆行列も対角行列となり、その要素は逆数となるので ( 補足 1 )、

S-1T=|0,-a12/a11,...-a1n/a11|
|-a21/a22,0,...-a2n/a22|
|::...:|
|-an1/ann,-an2/ann,...0|

Gershgorin の円定理 ( 補足 2 )」から、行列の任意の固有値 λ に対して | mrr - λ | ≤ Rr ≡ Σc{ 1→n ; c≠r }( |mrc| ) ( 但し、mrc は n 行 n 列の正方行列の r 行 c 列の要素 ) を満たす Rr が存在するので、上記の行列の場合、対角成分が 0 であることから、各固有値に対して

| 0 - λ | = |λ| ≤ Σc( |arc| / |arr| ) ( r = 1, 2, ... , n )

を満たす r が必ず存在する。従って、全ての行 ( r = 1, 2, ... , n ) に対し、Σc( |arc| / |arr| ) < 1 であれば、任意の固有値 λ に対しても |λ| < 1 が成り立つ。

つまり、Σc( |arc| / |arr| ) の最大値が 1 より小さければ、ヤコビ法による処理によって解は必ず収束します。このような行列を「狭義優対角行列」といいます。ちなみに、「(広義)優対角行列」は Σj( |aij| / |aii| ) の最大値が 1 以下となる行列のことを指します。


6) 性能評価

直接法と反復法のそれぞれのアルゴリズムを使ったときの処理時間を計測した結果を以下に示します。表の N 及びグラフの横軸は未知数と式の数を、グラフの縦軸は一回の処理にかかった時間(sec.)を表しています。計測は Windows7 上の VirtualBox でゲスト OS を Xubuntu 16.04 LTS として行っています。CPU は Intel Core i7-2600 3.40GHz で、ゲスト OS 上ではコア数を 1 としています。

表 6-1. 各解法の性能評価結果
N直接法反復法
ガウスの消去法ガウス・ジョルダン法ヤコビ法ガウス=ザイデル法
処理時間処理時間処理時間処理回数処理時間処理回数
1000.0010620.0015850.000241100.0001058
2000.0087920.0294430.001122110.0003888
4000.2533480.4105460.030131280.02836319
6001.916644.224720.106021260.06559315
8005.6249913.39400.166268180.12168210
100012.524829.01570.376266220.22230915
 
図 6-1. 処理性能評価グラフ
直接法反復法
処理性能評価グラフ(直接法) 処理性能評価グラフ(反復法)

ガウスの消去法では、未知数の数が N の時、一つの方程式にある N 個の係数に対する加減乗除の操作を約 N( N + 1 ) / 2 回繰り返すことになるので、処理量はだいたい N3 に比例することになります。上記結果を見ると、N が小さいうちは未知数が 2 倍になると処理量は約 8 倍に増えていて理論値に近い値を示していますが、N が大きくなると増加率はある程度抑えられています。CPU のキャッシュ等が有効に働いたのか、それとも他の理由によるものかは不明です。また、ガウス・ジョルダン法は、一つの式に対する加減乗除の操作が約 N2 回とガウスの消去法よりも多くなる分、処理速度は遅くなります。
直接法に比べると、反復法は圧倒的に早い速度で処理ができます。今回のテストでは、対角成分が他の成分よりもかなり大きく ( 103 倍程度に ) なるように係数に細工をしてあるため、多くても数十回程度のループで収束します。しかし、対角成分に対して変更を行わない場合、N = 100 の場合でも収束しなくなりました。少なくともこのサンプル・プログラムに関しては、反復法は用途がかなり限定され、汎用的に利用することはできないようです。それでも、対角成分が他の成分に比べて大きいことが既知であれば反復法の方がパフォーマンス的に有利であることは間違いないようです。直接法で実用に耐えうるのは、せいぜい 2000程度ではないでしょうか。この場合、ガウスの消去法でも、約 2 分程度かかる計算になります。

以下は、SOR 法を利用した場合の処理時間と処理回数を表にまとめたものになります。一つの条件に対し、処理は N = 1000 の連立方程式をランダムに変えながら 10 回行い、その平均値を示しています。

表 6-2. SOR 法における処理時間と処理回数の推移
反復法ヤコビ法ガウス=ザイデル法
緩和係数処理時間処理回数処理時間処理回数
0.850.08410131.10.07321021.8
0.860.11285130.70.07134121.3
0.870.10801130.20.06384420.8
0.880.12347230.00.07975420.3
0.890.13592129.70.07102019.8
0.900.09775829.30.06447519.3
0.910.09564629.20.05646518.9
0.920.09907128.80.05685818.5
0.930.09800528.60.05686218.0
0.940.09394028.30.05215317.7
0.950.09339028.10.05119217.2
0.960.09891928.00.05031117.0
0.970.09412727.70.04854516.7
0.980.09357427.70.04272016.7
0.990.10123627.90.04805716.2
1.000.12220627.70.05614516.2
1.010.11426828.00.05852416.5
1.020.08904228.30.04789116.9
1.030.09191128.40.04182417.2
1.040.09056029.00.05369317.6
1.050.09354029.20.05993818.1
 
図 6-2. 緩和係数と処理性能の相関グラフ
処理時間処理回数
緩和係数と処理時間の相関図 緩和係数と処理回数の相関図

処理時間はバラツキが大きいため ( さらに少しのソース変更でも処理時間が変化してしまいます )、処理回数の方が傾向がわかりやすいと思います。緩和係数の小さい側から 1 に近づくに従って処理回数は減少する傾向にあり、0.97 〜 0.98 あたりで最小値となりました。また、その後は増大する傾向にあります。ヤコビ法と比較してガウス=ザイデル法の方が処理時間・処理回数ともに小さくなります。処理の効率がよく、収束もしやすいことをこの結果は表しています。


反復法は、「正則化最小二乗法」の「LASSO」でも応用されています。特定の条件が必要になるものの、うまく活用すれば直接法よりも高速に処理を行うことができます。欠点としては収束しない場合の対応で、対角成分が大きい場合でも収束しないときがあるので他の手法に切り替えるなどの対処が必要になります。


補足 1) 対角行列の逆行列

n 行 n 列の対角行列 D の逆行列 D-1 に対し、DD-1 = E が成り立ちます。D の r 行めのベクトル成分を dr = ( 0, 0, ... , drr, ... , 0 )、D-1 の c 列めのベクトル成分を d'c = ( d'1c, d'2c, ... d'nc )T とした時、

DD-1 = ( d1, d2, ... dn )T( d'1, d'2, ... d'n )

よって、drd'c = drrd'rc と計算することができます。右辺は単位行列であり、対角成分のみが 1 で残りの成分は全て 0 なので、結局 d'rc は r ≠ c の場合 0 になり、r = c のときは drrd'rc = 1 より

d'rc = 1 / drr

となります。よって、逆行列も対角行列となり、その成分は元の対角行列の要素の逆数となります。


補足 2) Gershgorinの円定理 ( Gershgorin Circle Theorem )

Gershgorinの円定理」は、次のような内容になります。

行列 A の r 行 c 列の要素を arc とする。行列の r 行めのベクトル成分から対角要素を除いた値の絶対値の和を取り Rr とした時、行列 A の任意の固有値 λ に対し、| arr - λ | ≤ Rr を満たす Rr が必ず存在する。

証明は次のようになります。

行列 A の固有値を λ、その時の固有ベクトルを x とし、x の成分の中で絶対値が最大のものを xr とする。A の r 行 c 列の要素を arcx = ( x1, x2, ... )T とした時、Ax = λx より、

Σc( arcxc ) = λxr

x0 から |xr| > 0 なので、両辺を xr で割って、

Σc( arcxc / xr ) = λ

左辺から、対角成分 arc を含む項を右辺に移項すると、

Σc{r≠c}( arcxc / xr ) = λ - arrxr / xr = λ - arr

xrx の成分の中で絶対値が最大となるため、r ≠ c のとき | xc/xr | ≤ 1。よって、左辺の絶対値は必ず | Σc{r≠c}( arc ) |以下になる。従って、

| λ - arr |=| Σc{r≠c}( arcxc / xr ) |
| Σc{r≠c}( arc ) |
Σc{r≠c}( | arc | ) = Rr

複素平面上において、対角要素 arr を中心とする半径 Rr の円を描いたとき、固有値はいずれかの円の中に含まれ、その外部に存在することはないことを、この定理は示しています。この円のことを「Gershgorin の円 ( Gershgorin Disc )」といいます。この円は、非対角成分の値に応じて大きさが変化し、対角行列に対しては全ての円の半径は 0 になります。このことは、固有値が対角成分と等しくなることを示しています。優対角行列では、各行において、非対角成分の和が対角成分以下になります。Gershgorin の円定理から、優対角行列の場合は対角成分と固有値の差が 1 以下であることになり、対角成分を固有値と見なしてもそれほど大きなズレはないということになります。


<参考文献>
  1. 「これなら分かる応用数学教室」 金谷健一著 (共立出版)
  2. 川上一郎「数値計算の基礎」(http://www7.ocn.ne.jp/~kawa1/) : ジョルダン標準系(PDF形式)
  3. 「Yamamoto's Laboratory」 「連立一次方程式(反復法)」
  4. Wikipedia ... 「Gershgorinの円定理」は Wikipediaを参考にしました。しかし、日本語のページはないようです...

[Go Back]前に戻る [Back to HOME]タイトルに戻る
inserted by FC2 system