数値演算法

(7) 連立方程式を解く -1-

連立方程式の解き方は、中学校で習ったという方が大半ではないでしょうか。といっても、そこで教えてもらうのは二元一次方程式であり、求める未知数も式も二つまででした。いわゆる鶴亀算を解いたり、二直線の交点を求めたりするのに用いられ、その解き方も、代入法や加減法といったパターンがあって、計算ミスさえなければそれほど難しくはありません。
もっと一般的な連立一次方程式は、データ解析や固有値問題など多くの分野で利用されるため、様々な解の求め方が誕生してきました。解くことは難しくないものの、中学校で習った方法を利用した場合、未知数の数が多くなるに従い解くのに非常に時間がかかるので、少しでも短時間で解を求めるための手法が開発されてきたわけです。ここでは、連立方程式を解くための手法について紹介したいと思います。

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

1) ガウスの消去法 ( Gaussian elimination )

まずは中学校で習った「消去法」を復習しましょう。次の連立方程式を解くとします。

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

まずは (1) から (2) × 3 を引きます。

   3x + 4y = 7
-) 3x - 9y = -6
----------------
       13y = 13
         y = 1

これで y の値が求められました。次に、(1) に y = 1 を代入します。

3x + 4 = 7
    3x = 3
     x = 1

よって ( x, y ) = ( 1, 1 ) になります。

消去法は、二式の加減算によって未知数の片方を消去することで解を求めていました。この方法は、未知数が三つ以上になっても利用することができます。

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

(1) × 2 から (2) を引く

   4x + 6y -  8z = 12
-) 4x -  y +  7z = -3
---------------------
        7y - 15z = 15 ... (2')

(1) から (3) × 2 を引く

   2x + 3y - 4z = 6
-) 2x + 2y - 6z = 8
---------------------
         y + 2z = -2 ... (3')

(2') から (3') × 7 を引く

   7y - 15z = 15
-) 7y + 14z = -14
---------------------
       -29z = 29
          z = -1

(2') に z = -1 を代入する

7y + 15 = 15 より y = 0

(1) に y = 0, z = -1 を代入する

2x + 4 = 6 より x = 1

よって ( x, y, z ) = ( 1, 0, -1 )

消去法では、異なった式同士を加減算することによって、特定の未知数を消去した式を新たに作成する操作を繰り返し、最終的には、式と未知数の数を一つ減らします。これを繰り返していくと、最後に未知数と式がどちらも一つだけ残ることになり、これで解の一つが求められます。今までの処理で、未知数が消去された方程式がいくつも作成されているので、これらに求めた解を代入して他の解を求める操作を続けていくと、最終的に全ての解が求められます。どの式を組み合わせてどのように未知数を消去するか、またどの式に求めた解を代入して他の解を求めていくかは、経験と勘に頼っている部分が大きいと思いますが、これを機械的に処理できるようにしたのが「ガウスの消去法 ( Gaussian elimination )」です。

ガウスの消去法は、「前進消去 ( Forward Elimination )」と「後退代入 ( Back Substitution )」の二つのステップから成ります。前の説明にある、式と未知数の数を一つずつ減らしていく操作が前進消去、解を一つ得てから他の解を求めていく操作が後退代入にあたります。
以下の説明において、n 行目にある式は「式(n)」、n 番目の未知数の係数は「係数[n]」と表すようにします。まずはじめに、式(1)にある係数[1]の値が 1 になるように、式(1)全体を係数[1]で割ります。次に、二行目以降の各式について、その式の係数[1]を式(1)に掛けた結果で減算する処理を行うと、それぞれの式の一つ目の未知数は消去されます。以下、式(2)を使って以降の式から二番目の未知数を、式(3)を使って以降の式から三番目の未知数を、という形で次々と係数を消去すると、最後の式にまで達したときには、未知数が一つだけ残る形になります。ちょうど、式を並べたときに三角形のイメージになると考えると分かりやすいと思います。
最後の式は未知数が一つだけなので、その解はすぐに計算できます。この結果を一つ上の式に代入すると、残った未知数がやはり一つとなって、これも簡単に計算できます。このように、順次求めた解を全て上側の式に代入すれば、未知数が一つだけ残り、解を求めることができます。一番上にある式の未知数が解けたら、全ての結果が得られたことになります。

先程と同じ連立方程式を、ガウスの消去法を使って解いてみます。まずは、前進消去の処理です。

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

(1) を 2 で割って、x の係数が 1 になるようにする。

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

(1) × 4 を (2) から引く。

 x + (3/2)y -  2z = 3   ... (1)
        -7y + 15z = -15 ... (2)
 x +      y -  3z = 4   ... (3)

(1) を (3) から引く。

 x + (3/2)y -  2z = 3   ... (1)
        -7y + 15z = -15 ... (2)
    -(1/2)y -   z = 1   ... (3)

(2) を -7 で割って、y の係数が 1 になるようにする。

 x + (3/2)y -      2z = 3    ... (1)
          y - (15/7)z = 15/7 ... (2)
    -(1/2)y -       z = 1    ... (3)

(2) × (-1/2) を (3) から引く。

 x + (3/2)y -       2z = 3     ... (1)
          y -  (15/7)z = 15/7  ... (2)
            - (29/14)z = 29/14 ... (3)

前進消去の処理を行った結果、最下位にある式は未知数が一つのみで、上位側に進むに従い、一つずつ未知数が増えていくような形になります。この状態になったら、次に後退代入の処理を行います。

(3) から z の値を求める。

- (29/14)z = 29/14 より z = -1

(2) に z = -1 を代入して y の値を求める。

y - (15/7)×(-1) = 15/7 より y = 0

(1) に y = 0, z = -1 を代入して x の値を求める。

x + (3/2)×0 - 2×(-1) = 3 より x = 1

以上で、連立方程式の解を求めることできました。

さて、連立方程式を解く場合に気をつけるべきこととして、解を求めることができない場合(不能)と解が定まらない場合(不定)の存在があります。例えば、次の連立方程式

6x + 2y = 4 ... (1)
3x +  y = 1 ... (2)

を解く場合、(2) を 2 倍して二式を引くと、左辺が 0、右辺が 2 となって等号が成り立たなくなります。この場合、( x, y ) としてどのような値を組み合わせても二式を満たすことができないので、解が存在しないことになり不能になります。また、次の連立方程式

3x +  4y = 1 ... (1)
9x + 12y = 3 ... (2)

では (1) を 3 倍して二式を引くと、両辺とも 0 になります。この場合、( x, y ) としてどのような値を組み合わせても常に二式を満たすことになるので、解が一意に決められず不定となります。前進消去によって未知数を一つ消去したときに、消去を行った全ての式に対して他の未知数も同時に消去されてしまった場合、解は不定または不能の状態になります。そのまま処理を進めると、係数で式を徐算する際にゼロ徐算を行うことになるので、係数をチェックすることで不定・不能の状態であるかが分かります。以下に、不能となる連立方程式をガウスの消去法で処理した場合の例を示します。

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

(1) を 2 で割って、x の係数が 1 になるようにする。

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

(1) を (2) から、(1) × 4 を (3) から引く。

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

(2) × 0 を (3) から引く。

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

(3) 式を満たす z の値が存在しないため、不能となる。

ところで、ゼロ徐算が発生する場合は、不能・不定の場合のみであるとは限りません。例えば、次の連立方程式は解が存在するにも係わらず、ガウスの消去法を利用するとゼロ徐算が発生してしまいます。

   4y = 8 ... (1)
x + y = 6 ... (2)

消去する対象の未知数が持つ係数が 0 であった場合、下側にある式と位置を交換することで処理ができる状態になるように調整を行います。上の例では、(1) と (2) を交換すれば処理を進めることができます。このような処理を「ピボット選択 ( Pivot )」と言います。ピボット選択には「部分選択法」と「完全選択法」の二種類あり、式を交換するやり方は「部分選択法」にあたります。式だけではなく、その中の各項(係数)も交換する方法が「完全選択法」になりますが、式を交換しても結果に変化がないのに対し、項を交換した場合は求める解の位置が変わるので、あとで補正するために交換した位置を記憶しておく必要があります。

ピボット選択は、ゼロ徐算を防ぐだけではなく、解の誤差を小さくする目的にも利用されます。

εx1 + a2x2 + a3x3 + ... + anxn = b

上式で、ε は他の係数に比べて非常に小さいとします。x1 の係数を 1 にするために、両辺を ε で割ると、

x1 + (a2/ε)x2 + (a3/ε)x3 + ... + (an/ε)xn = b/ε

ここで、x2 以降の係数は他式のそれに対して非常に大きくなる可能性があります。すると、他式との差を計算した時に、値の差が大きくなりすぎて丸め誤差が発生することになり、この式一つだけで他の式の係数が丸められてしまうことによって、実際の解から大きく外れてしまうことになります。このような現象は、ピボット選択によって係数がなるべく大きいものを選択するすることで回避することができます。たとえこの式の中の係数が無視できるほど小さくなったとしても、他の式の係数は影響を受けないことになるため誤差を最小限に抑えることができるわけです。
しかし、単純に係数の大小を比較するだけでは不十分で、徐数となる係数は大きいが、他の項の係数はさらに大きいというような場合は、徐算した後も係数は大きいままになります。よって、式の中で最も大きな値を抽出して、その値で徐算した上で大小比較すると、さらに誤差を抑えることができます。


以下に、サンプル・プログラムを示します。まずは、連立方程式を表現するクラスから示します。

/**
   連立一次方程式

   連立方程式 Ax = y を定義するためのクラス。A は係数行列、x は未知数、y は解をそれぞれ表す。
**/
template< typename T >
class LinearEquationSystem
{
 public:

  // 係数や解の型
  typedef T value_type;
  // 係数行列のサイズに対する型
  typedef typename SquareMatrix< T >::size_type size_type;

  // 係数行列への反復子の型(行方向のみ利用可能)
  typedef typename SquareMatrix< T >::iterator iterator;
  // 係数行列への定数反復子の型(行方向のみ利用可能)
  typedef typename SquareMatrix< T >::const_iterator const_iterator;

 private:

  SquareMatrix< T > coef_; // 係数行列 A (左辺)
  std::valarray< T > ans_; // 解ベクトル y (右辺)

 public:

  /*
    デフォルト・コンストラクタ

    sz : 方程式(未知数)の数
  */
  explicit LinearEquationSystem( size_type sz = 0 )
    : coef_( sz ), ans_( sz ) {}

  /*
    係数行列 coef と解ベクトル [ s, e ) からのコピー
  */
  template< class In >
    LinearEquationSystem( const SquareMatrix& coef, In s, In e )
  { init( coef, s, e ); }

  /*
    係数行列 coef からのコピー

    解ベクトルの内容は全て T() に初期化される。
  */
  explicit LinearEquationSystem( const SquareMatrix< T >& coef )
  {
    assert( &coef != 0 );

    coef_ = coef;
    ans_.resize( coef_.size() );
  }

  /*
    コピー・コンストラクタ
  */
  LinearEquationSystem( const LinearEquationSystem< T >& les )
  {
    coef_ = les.coef_;
    ans_.resize( les.size() );
    ans_ = les.ans_;
  }

  /*
    係数行列 coef と解ベクトル [ s, e ) による初期化
  */
  template< class In >
    void init( const SquareMatrix& coef, In s, In e )
  {
    assert( &coef != 0 );

    ErrLib::CheckLinearModel< SquareMatrix >( coef, s, e );

    coef.clone( &coef_ );
    ans_.resize( coef.size() );
    for ( size_type i = 0 ; i < ans_.size() ; ++i  )
      ans_[i] = *s++;
  }

  /*
    未知数の数(=方程式の数)を返す
  */
  size_type size() const
  { return( coef_.size() ); }

  /*
    未知数の数(=方程式の数)を変更する

    sz : 未知数の数(=方程式の数)
    val : 係数と解の初期値
  */
  void resize( size_type sz, value_type val = value_type() )
  {
    coef_.assign( sz, val );
    ans_.resize( sz, val );
  }

  /*
    方程式の入れ換え(係数行列の行と解の入れ替え)

    式の入れ替えだけなので、連立方程式としては変化しない。

    e1 : 入れ替える対象の方程式 1
    e2 : 入れ替える対象の方程式 2
  */
  void swapEquation( size_type e1, size_type e2 )
  {
    if ( e1 == e2 ) return;

    iterator row1 = (*this)[e1];
    iterator row2 = (*this)[e2];
    while ( row1 != row1.end() ) {
      value_type val = *row1;
      *row1++ = *row2;
      *row2++ = val;
    }
    std::swap( ans_[e1], ans_[e2] );
  }

  /*
    係数の入れ換え(係数行列の列の入れ替え)

    係数の入れ替えなので未知数の位置が変化するが、その管理は行わない。
    そのため、入れ替え時に外部で管理する必要がある。

    c1 : 入れ替える対象の係数 1
    c2 : 入れ替える対象の係数 2
  */
  void swapCoef( size_type c1, size_type c2 )
  {
    if ( c1 == c2 ) return;

    iterator col1 = coef_.column( c1 );
    iterator col2 = coef_.column( c2 );
    while ( col1 != col1.end() ) {
      value_type val = *col1;
      *col1++ = *col2;
      *col2++ = val;
    }
  }

  /*
    指定した方程式 src を val で定数倍して他の方程式 dest から引く

    連立方程式としては変化しない。

    src : 対象の方程式の添字
    val : 定数
    dest : 被演算対象の方程式の添字
  */
  void sub( size_type src, T val, size_type dest )
  {
    const_iterator sRow = (*this)[src];
    iterator dRow = (*this)[dest];

    while ( sRow != sRow.end() )
      *dRow++ -= *sRow++ * val;

    ans_[dest] -= ans_[src] * val;
  }

  /*
    指定した方程式 r を定数 val で割る

    連立方程式としては変化しない。

    r : 対象の方程式の添字
    val : 定数
  */
  void div( size_type r, T val )
  {
    for ( iterator row = (*this)[r] ; row != row.end() ; ++row )
      *row /= val;
    ans_[r] /= val;
  }

  /* 要素へのアクセス */

  /*
    指定した行番号 r の係数行列への反復子を返す

    反復子を使って係数を変更した場合、変更前の連立方程式と同一である保証はなくなる。
  */
  iterator operator[]( size_type r )
  { return( coef_[r] ); }

  /*
    指定した行番号 r の係数行列への定数反復子を返す
  */
  const_iterator operator[]( size_type r ) const
  { return( coef_[r] ); }

  /*
    係数の値への参照を返す

    係数を変更した場合、変更前の連立方程式と同一である保証はなくなる。

    r : 方程式の番号(行番号)
    c : 未知数の番号(列番号)
  */
  value_type& coef( size_type r, size_type c )
  { return( coef_[r][c] ); }

  /*
    係数の値を返す

    r : 方程式の番号(行番号)
    c : 未知数の番号(列番号)
  */
  value_type coef( size_type r, size_type c ) const
  { return( coef_[r][c] ); }

  /*
    解の値への参照を返す

    解を変更した場合、変更前の連立方程式と同一である保証はなくなる。

    i : 対象の解に対する添字
  */
  value_type& ans( size_type i )
  {
    if ( i >= size() )
      throw std::out_of_range( "Exception : out of range @LinearEquationSystem::ans" );

    return( ans_[i] );
  }

  /*
    解の値を返す

    i : 対象の解に対する添字
  */
  value_type ans( size_type i ) const
  {
    if ( i >= size() )
      throw std::out_of_range( "Exception : out of range @LinearEquationSystem::ans" );

    return( ans_[i] );
  }

  /*
    代入演算子の多重定義
  */
  LinearEquationSystem& operator=( const LinearEquationSystem< T >& les )
  {
    if ( this == &les ) return( *this );

    coef_ = les.coef_;
    ans_.resize( les.size() );
    ans_ = les.ans_;

    return( *this );
  }
};

連立方程式は、クラス LinearEquationSystem の中で係数行列 coef_ と解ベクトル ans_ を使って表現されています。連立方程式が

a11x1 + a12x2 + ... + a1nxn = b1
a21x1 + a22x2 + ... + a2nxn = b2
:
an1x1 + an2x2 + ... + annxn = bn

で表されるとき、n 行 n 列の正方行列 A 及び n 次のベクトル x, b を以下のように定義すると、Ax = b と表すことができます。coef_A に、ans_b に対応し、未知数 x は処理によって ans_ に出力されるのでメンバ変数としては不要になります。

A = | a11, a12, ... a1n |
| a21, a22, ... a2n |
|:: ...:|
| an1, an2, ... ann |

x = ( x1, x2, ... xn )T

b = ( b1, b2, ... bn )T

LinearEquationSystem では、係数行列や解ベクトルを初期化するためのメンバ関数 init、ガウスの消去法での前進消去や後退代入を行うために必要なメンバ関数 subdiv、ピボットに必要なメンバ関数 swapEquation, swapCoef が用意されています。また、係数や解へのアクセスを行うためのメンバ関数として operator[]coefans があります。
係数行列の型である SquareMatrix は、任意の型を要素とする正方行列を表したクラスです。定義や実装はサンプル・プログラム内では省略していますが、任意の行や列への反復子 ( iterator ) を返すメンバ関数 operator[], column を持ち、各要素への反復処理ができるようになっています。また、coef_[r][c] によって r 行 c 列の要素へのアクセスができます。
解ベクトルには、STL ( Standard Template Library ) で提供されている valarray を使用しています。vector と同様に配列 ( ベクトル ) を表現したコンテナ・クラスですが、vector のように任意の要素数のインスタンスをコピーすることはできず、あらかじめ同じ要素数にリサイズした上でコピーする必要があります。そのため、コピー・コンストラクタと代入演算子の多重定義はデフォルトのものが利用できず、改めて実装を行っています。vector を使えばその必要はありませんが、数値演算に対しては valarray の方が高速化が期待できるため、今回は valarray を採用しました。


これらの道具が一通り揃えば、ガウスの消去法を処理することが可能になります。

/*
  ガウスの消去法により連立方程式 sysLE の解を求める

  解が不定または不能の場合は false を返す
*/
template< typename T >
bool GaussianElimination( LinearEquationSystem< T >& sysLE )
{
  typedef typename LinearEquationSystem< T >::size_type size_type;

  assert( &sysLE != 0 );

  // 前進消去
  for ( size_type i = 0 ; i < sysLE.size() ; ++i ) {
    Pivot( &sysLE, i );
    T k = sysLE[i][i];
    if ( k == 0 ) return( false ); // 最大係数が 0 の場合は不定または不能となるため false
    sysLE.div( i, k );
    for ( size_type r = i + 1 ; r < sysLE.size() ; ++r )
      sysLE.sub( i, sysLE[r][i], r );
  }

  // 後退代入
  for ( size_type i = sysLE.size() ; i > 0 ; --i ) {
    for ( size_type j = i ; j < sysLE.size() ; ++j )
      sysLE.ans( i - 1 ) -= sysLE.ans( j ) * sysLE[i - 1][j];
    sysLE.ans( i - 1 ) /= sysLE[i - 1][i - 1];
  }

  return( true );
}

処理の内容は、前に説明したものと全く同じです。前進消去では、消去すべき項の係数で式全体を徐算した後、それより下側にある式に対して係数を消去するための減算処理を順に行っています。係数での徐算前にその値が 0 でないかをチェックして、0 であれば解が不定または不能となると見なして処理を中断しています。
後退代入では、それまでに求められた解を代入することで、未知数となっている変数の値を求めています。

前進消去での各式に対する処理に先立って、Pivotを呼び出してピボット選択を行っています。そのプログラムは次のようになります。

/*
  Normalize : 正規化した係数を返す

  正規化する ( r, c ) の位置の係数から行方向に末尾までたどって最大の値(スケール因子;scaled factor)を探し、
  ( r, c ) の位置の係数をその値で除算した結果を絶対値で返す。
*/
template< typename T >
T Normalize( const LinearEquationSystem< T >& sysLE,
             typename LinearEquationSystem< T >::size_type r, typename LinearEquationSystem< T >::size_type c )
{
  typedef typename LinearEquationSystem< T >::size_type size_type;

  T coef = std::abs( sysLE[r][c] ); // 正規化対象の係数

  T max = coef; // 現在の最大値を正規化対象の係数で初期化
  for ( size_type i = c + 1 ; i < sysLE.size() ; ++i ) {
    T t = std::abs( sysLE[r][i] );
    if ( max < t ) max = t;
  }

  return( ( max > 0 ) ? coef / max : T() );
}

/*
  Pivot : 部分ピボットを行う(正規化処理付き)

  i 行目の対角成分より列方向下側の係数の中で最大値を求め、最大値を持つ行と入れ替える。

  戻り値 : i行目と入れ替えした行の番号
*/
template< typename T >
typename LinearEquationSystem< T >::size_type
Pivot( LinearEquationSystem< T >* sysLE, typename LinearEquationSystem< T >::size_type i )
{
  typedef typename LinearEquationSystem< T >::size_type size_type;

  size_type maxIdx = i; // 最大値を持つ行をピボット対象の行番号で初期化
  T max = Normalize( *sysLE, maxIdx, i ); // 現在の最大値

  for ( size_type r = i + 1 ; r < sysLE->size() ; ++r ) {
    T t = Normalize( *sysLE, r, i );
    if ( max < t ) {
      maxIdx = r;
      max = t;
    }
  }
  if ( maxIdx != i ) sysLE->swapEquation( maxIdx, i );

  return( maxIdx );
}

ここでは「部分選択法」が使われています。i 番目の式が対象であるとした場合、i 番目の未知数が持つ係数を比較することになるので、それ以降の式に対して同じく i 番目の係数をチェックして、その中で最大の値を持つものを決定し、それと i 番目の式を交換します。係数を比較する際、その式の係数の最大値であらかじめ徐算した上で行った方がよいことは前に説明しましたが、それを行っているのが Normalizeです。ここでは係数の最大値を求めてから、その値で対象の係数を徐算した結果を返す処理を行っています。なお、対象の係数を持つ項より前の項は、消去法によってすでに消去済みなので、除外することができます。


2) ガウス・ジョルダン法 ( Gauss-Jordan Elimination )

ガウスの消去法では、式同士の減算はより下側にある式だけを対象に行っていました。これを、すでに処理された上側の式にも適用することを検討してみます。

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

      ↓

x + (3/2)y -  2z = 3   ... (1)
       -7y + 15z = -15 ... (2)
   -(1/2)y -   z = 1   ... (3)

上の例は、ガウスの消去法を使って x の項を消去した後の状態です。次に (2) を対象に処理を続けます。

x + (3/2)y -      2z = 3    ... (1)
         y - (15/7)z = 15/7 ... (2)
   -(1/2)y -       z = 1    ... (3)

ここで、(3) だけでなく (1) にある y の項も同様に処理します。(2) × (3/2) を (1) から、(2) × (-1/2)を (3) から減算すると、次のようになります。

x   - (17/14)z = -3/14 ... (1)
  y -  (15/7)z = 15/7  ... (2)
    - (29/14)z = 29/14 ... (3)

最後に、(3) を -29/14 で徐算してから、(3) × (-15/7) - (2)、(3) × (-17/14) - (1) を計算して両式の z の項を消去します。

x     = 1  ... (1)
  y   = 0  ... (2)
    z = -1 ... (3)

前進消去処理が完了した段階で、各式にはそれぞれ異なる項が一つだけ残ることになり、しかもその係数は全て 1 になるので、後退代入をしなくてもすでに解が求められた状態になっています。この手法を「ガウス・ジョルダン法」と言います。ガウス・ジョルダン法は、ガウスの消去法に比べて演算回数は増えますが、後退代入が不要な分、実装がシンプルになるという利点があります。


ガウス・ジョルダン法のサンプル・プログラムを以下に示します。

/*
  ガウス・ジョルダン法により連立方程式 sysLE の解を求める

  解が不定または不能の場合は false を返す
*/
template< typename T >
bool GaussJordanElimination( LinearEquationSystem< T >& sysLE )
{
  typedef typename LinearEquationSystem< T >::size_type size_type;

  assert( &sysLE != 0 );

  for ( size_type i = 0 ; i < sysLE.size() ; ++i ) {
    Pivot( &sysLE, i );
    T k = sysLE[i][i];
    if ( k == 0 ) return( false ); // 最大係数が 0の場合は不定または不能となるためfalse
    sysLE.div( i, k );
    for ( size_type r = 0 ; r < sysLE.size() ; ++r ) {
      if ( r == i ) continue;
      sysLE.sub( i, sysLE[r][i], r );
    }
  }

  return( true );
}

内容は、ガウスの消去法における前進消去の部分とほとんど変わりません。式同士の減算処理の範囲が変わる以外は全く同じ内容です。


3) 連立方程式による逆行列の計算

数値演算法 (4) 高速フーリエ変換」の章で、逆行列の計算による連立方程式の解法について説明をしました。この時は逆行列の計算に余因子行列を利用したのですが、あまりに処理が重くて使用に耐えないような代物でした。そもそも通常は、逆行列の計算自体に連立方程式が利用されているため、全く逆のことをしていたわけです。

前述の通り、連立方程式は n 行 n 列の正方行列 A 及び n 次のベクトル x, b によって Ax = b と表すことができます。x の成分を未知数としたとき、ガウスの消去法やガウス・ジョルダン法によって任意の b に対して x を求めることができます。b の中の一つの要素を 1 とおき、残りを 0 にした上で連立方程式を解く操作を、各要素が 1 となるように n 回繰り返したとき、b の第 i 成分が 1 であるときの解を ci = ( ci1, ci2, ... cin )T とすると、A( c1, c2, ... cn ) = E が成り立つことになります。ここで、E は単位行列 ( Identity Matrix ) を表します。各成分を使って書くと、次のようになります。

| a11, a12, ... a1n || c11, c12, ... c1n | = |1, 0, ... 0|
| a21, a22, ... a2n || c21, c22, ... c2n ||0, 1, ... 0|
|:: ...:||:: ...:||:: ... :|
| an1, an2, ... ann || cn1, cn2, ... cnn ||0, 0, ... 1|

上式から、行列 C = ( c1, c2, ... cn ) は行列 A の逆行列を表していることになります。

ガウスの消去法を利用して逆行列を求める処理のサンプル・プログラムを以下に示します。

/*
  連立方程式 les の係数行列の逆行列を求める

  les : 対象の行列を係数行列とする LinearEquationSystem クラスのインスタンス
  inv 求めた逆行列を書き込む正方行列へのポインタ

  lesの係数行列が正則行列ではない(連立方程式の解が不定または不能の)場合は false を返す
*/
template< typename T >
bool InverseMatrix( LinearEquationSystem< T >& les, SquareMatrix< T >* inv )
{
  typedef typename LinearEquationSystem< T >::size_type size_type;

  size_type sz = les.size();
  inv->assign( sz );

  for ( size_type i = 0 ; i < sz ; ++i ) {
    LinearEquationSystem< T > buff( les );
    for ( size_type r = 0 ; r < sz ; ++r )
      buff.ans( r ) = ( r == i ) ? 1 : 0;
    if ( ! GaussianElimination( buff ) ) return( false );
    for ( size_type r = 0 ; r < sz ; ++r )
      (*inv)[r][i] = buff.ans( r );
  }

  return( true );
}

ループ内では、行列数と等しいサイズの連立方程式を、buff にコピーして、単位行列の第一列から順番に列の要素を解として未知数を解いています。連立方程式が解けない ( 解が不定あるいは不能になった ) 場合、逆行列は存在しないことになるため処理を中断します。なお、逆行列が存在する行列は「正則行列 ( Regular Matrix )」、それに対して逆行列の存在しない行列は「特異行列 ( Singular Matrix )」と呼ばれます。


4) LU分解 ( LU Decomposition )

連立方程式を解くことによって逆行列を計算する場合、各項の係数はそのままで右辺の値だけが変化するため、同じ係数に対して前進消去を繰り返していることになります。あらかじめ前進消去を済ませた状態で、右辺の値を変化させながら後退代入だけを行うようにすれば、処理量を減らすことができます。しかし、前進消去によって右辺の値も変化するため、ガウスの消去法ではそのようなことは簡単にはできません。そこで、別の手段を検討します。

正方行列に対し、行番号を r、列番号を c、その要素を arc としたとき、

r > c となる全ての要素に対して arc = 0 を満たす行列を上三角行列

r < c となる全ての要素に対して arc = 0 を満たす行列を下三角行列

といいます。各成分を使って書くと、次のようになります。

上三角行列 U = | u11, u12, ... u1(n-1), u1n |
| 0, u22, ... u2(n-1), u2n |
|:: ...::|
| 0, 0, ... u(n-1)(n-1), u(n-1)n |
| 0, 0, ... 0, unn |

下三角行列 L = | l11, 0, ... 0, 0 |
| l21, l22, ... 0, 0 |
|:: ...::|
| l(n-1)1, l(n-1)2, ... l(n-1)(n-1), 0 |
| ln1, ln2, ... ln(n-1), lnn |

上三角行列の形はちょうど、ガウスの消去法によって前進消去を行った後の状態を表しています。下三角行列は行の並びが逆転しているだけなので、どちらの三角行列も後退代入だけで解を得ることができます。行列がそれぞれの三角行列の積の形 A = LU で表せたと仮定すると、任意の連立方程式 Ax = bLUx = b になります。連立方程式 Ux = y は後退代入だけで解くことができます。この結果を元の連立方程式に代入すると Ly = b となり、これも後退代入だけで解けるので、前進消去を行うことなく解を得ることができることになります。

次に問題になるのは、どうやって行列を上三角行列と下三角行列に分解するのかということになります。行列のサイズが N のとき、分解後の行列において、0 以外の要素の数は Σk{1→N}( k ) = N( N + 1 ) / 2 なので、解くべき未知数はその倍の N( N + 1 ) = N2 + N になります。元の行列の要素数が N2 で、それぞれの要素に対して方程式が作られるので、N2 個の方程式を解いて N2 + N 個の未知数を解くことになり、一意には値を決めることができません。逆に言えば、計算しやすい形で処理できることになり、解法として「ドゥーリトル法 ( Doolittle Algorithm )」と「クラウト法 ( Crout Algorithm )」の二つが一般的に知られています。

ドゥーリトル法では、下三角行列 L の対角成分を 1 として、残りの要素を求めていきます。具体的な例として、5 次の正方行列を LU 分解してみましょう。

A = | a00, a01, a02, a03, a04 |
| a10, a11, a12, a13, a14 |
| a20, a21, a22, a23, a24 |
| a30, a31, a32, a33, a34 |
| a40, a41, a42, a43, a44 |

上三角行列を U、下三角行列を L としたとき、それぞれの要素を以下のように定義します。

L = | 1, 0, 0, 0, 0 |
| l10, 1, 0, 0, 0 |
| l20, l21, 1, 0, 0 |
| l30, l31, l32, 1, 0 |
| l40, l41, l42, l43, 1 |

U = | u00, u01, u02, u03, u04 |
| 0, u11, u12, u13, u14 |
| 0, 0, u22, u23, u24 |
| 0, 0, 0, u33, u34 |
| 0, 0, 0, 0, u44 |

LU = A から、

| 1, 0, 0, 0, 0 || u00, u01, u02, u03, u04 | = | a00, a01, a02, a03, a04 |
| l10, 1, 0, 0, 0 || 0, u11, u12, u13, u14 || a10, a11, a12, a13, a14 |
| l20, l21, 1, 0, 0 || 0, 0, u22, u23, u24 || a20, a21, a22, a23, a24 |
| l30, l31, l32, 1, 0 || 0, 0, 0, u33, u34 || a30, a31, a32, a33, a34 |
| l40, l41, l42, l43, 1 || 0, 0, 0, 0, u44 || a40, a41, a42, a43, a44 |

この中には 25 組の方程式があります。まず、L の 0 行目の要素と U の各列の要素との積から

u00 = a00, u01 = a01, u02 = a02, u03 = a03, u04 = a04

これで、u00 〜 u04 の各値が決まります。次に L の 1 行目から 4 行目までの要素と U の 0 列目の要素との積から

l10u00 = a10, l20u00 = a20, l30u00 = a30, l40u00 = a40

となり、u00 の値がすでに判明しているので l10 〜 l40 の各値を求めることができます。以下、次のように積を求めていきます。

L の 1 行目の要素と U の 1 列目から 4 列目までの要素との積から

l10u01 + [u11] = a11, l10u02 + [u12] = a12, l10u03 + [u13] = a13, l10u04 + [u14] = a14

L の 2 行目から 4 行目までの要素と U の 1 列目の要素との積から

l20u01 + [l21]u11 = a21, l30u01 + [l31]u11 = a31, l40u01 + [l41]u11 = a41

L の 2 行目の要素と U の 2 列目から 4 列目までの要素との積から

l20u02 + l21u12 + [u22] = a22, l20u03 + l21u13 + [u23] = a23, l20u04 + l21u14 + [u24] = a24

L の 3 行目から 4 行目までの要素と U の 2 列目の要素との積から

l30u02 + l31u12 + [l32]u22 = a32, l40u02 + l41u12 + [l42]u22 = a42

L の 3 行目の要素と U の 3 列目から 4 列目までの要素との積から

l30u03 + l31u13 + l32u23 + [u33] = a33, l30u04 + l31u14 + l32u24 + [u34] = a34

L の 4 行目の要素と U の 3 列目の要素との積から

l40u03 + l41u13 + l42u23 + [l43]u33 = a43

L の 4 行目の要素と U の 4 列目の要素との積から

l40u04 + l41u14 + l42u24 + l43u34 + [u44] = a44

それぞれの式で求めることのできる未知数は "[ ]" で囲ってあります。この例から、n 次の行列に対する一般式は次のようになることが分かります。

L の r 行目の要素と U の c 列目 ( c = r, r + 1, ... n - 1 ) の要素との積から

urc = arc - Σk{0→r-1}( lrkukc )

L の r 行目 ( r = c + 1, c + 2, ... n - 1 ) の要素と U の c 列目の要素との積から

lrc = [ arc - Σk{0→c-1}( lrkukc ) ] / ucc

但し、和の計算において、r - 1 < 0, c - 1 < 0 であれば 0 と見なします。

クラウト法の場合は、上三角行列 U の対角成分を 1 として残りの要素を求める形になり、考え方はドゥーリトル法と大差ありません。n 次の行列に対する一般式は次のようになります。

L の r 行目 ( r = c, c + 1, ... n - 1 ) の要素と U の c 列目の要素との積から

lrc = arc - Σk{0→c-1}( lrkukc )

L の r 行目の要素と U の c 列目 ( r = c + 1, c + 2, ... n - 1 ) の要素との積から

urc = [ arc - Σk{0→r-1}( lrkukc ) ] / lrr


ドゥーリトル法とクラウト法による LU 分解のサンプル・プログラムを以下に示します。

/*
  PivotRow : 行 i に対して部分ピボットを行う

  戻り値 : i 行目と入れ替えした行の番号
*/
template< typename T >
typename Matrix< T >::size_type
PivotRow( LinearEquationSystem< T >* sysLE, typename LinearEquationSystem< T >::size_type i )
{
  typedef typename LinearEquationSystem< T >::size_type size_type;

  size_type maxIdx = i;
  T max = std::abs( (*sysLE)[maxIdx][i] );

  for ( size_type r = i + 1 ; r < sysLE->size() ; ++r ) {
    T t = std::abs( (*sysLE)[r][i] );
    if ( max < t ) {
      maxIdx = r;
      max = t;
    }
  }
  if ( maxIdx != i ) sysLE->swapEquation( maxIdx, i );

  return( maxIdx );
}

/*
  PivotCol : 列 i に対して部分ピボットを行う

  戻り値 : i 列目と入れ替えした列の番号
*/
template< typename T >
typename Matrix< T >::size_type
PivotCol( LinearEquationSystem< T >* sysLE, typename LinearEquationSystem< T >::size_type i )
{
  typedef typename LinearEquationSystem< T >::size_type size_type;

  size_type maxIdx = i;
  T max = std::abs( (*sysLE)[i][maxIdx] );

  for ( size_type c = i + 1 ; c < sysLE->size() ; ++c ) {
    T t = std::abs( (*sysLE)[i][c] );
    if ( max < t ) {
      maxIdx = c;
      max = t;
    }
  }
  if ( maxIdx != i ) sysLE->swapCoef( maxIdx, i );

  return( maxIdx );
}

/*
  Doolittle : LU 分解 ドゥーリトル法

  ・LU 分解に失敗した場合 ( false が返された場合 ) の l と u の内容は不定である。
  ・精度向上のため未知数の入れ替え ( ピボット選択 ) を行なっている。よって、行列の LU 分解の用途には利用できない。

  a : LU 分解する対象の行列
  l : 分解した下三角行列 L
  u : 分解した上三角行列 U
  index : 処理前に対する処理後の未知数の順番 ( NULL ならば無視 )

  行列 a が LU 分解できない場合 ( 処理途中で対角成分にゼロが見つかった場合 ) は false を返す
*/
template< typename T >
bool Doolittle( LinearEquationSystem< T > a,
                LinearEquationSystem< T >* l, LinearEquationSystem< T >* u,
                std::vector< typename LinearEquationSystem< T >::size_type >* index = 0 )
{
  assert( l != 0 && u != 0 );

  typedef typename LinearEquationSystem< T >::size_type size_type;

  size_type sz = a.size();
  l->resize( sz );
  u->resize( sz );
  if ( index != 0 ) index->resize( sz );

  for ( size_type i = 0 ; i < sz ; ++i ) {
    (*l)[i][i] = 1;         // Lの対角成分
    (*u)[0][i] = a[0][i];   // Uの一行目
    l->ans( i ) = u->ans( i ) = a.ans( i );
    if ( index != 0 ) (*index)[i] = i;
  }

  for ( size_type i = 1 ; i < sz ; ++i ) {
    size_type swapped = PivotCol( u, i - 1 );
    if ( (*u)[i - 1][i - 1] == 0 ) return( false );
    if ( swapped != i - 1 ) {
      a.swapCoef( swapped, i - 1 );
      if ( index != 0 )
        std::swap( (*index)[swapped], (*index)[i - 1] );
    }
    for ( size_type j = i ; j < sz ; ++j ) {
      // L成分の算出
      (*l)[j][i - 1] = a[j][i - 1];
      for ( size_type k = 0 ; k < i - 1 ; ++k )
        (*l)[j][i - 1] -= (*l)[j][k] * (*u)[k][i - 1];
      (*l)[j][i - 1] /= (*u)[i - 1][i - 1];

      // U成分の算出
      (*u)[i][j] = a[i][j];
      for ( size_type k = 0 ; k < i ; ++k )
        (*u)[i][j] -= (*l)[i][k] * (*u)[k][j];
    }
  }

  return( true );
}

/*
  Crout : LU 分解 クラウト法

  ・LU 分解に失敗した場合 ( false が返された場合 ) の l と u の内容は不定である。
  ・精度向上のため式の入れ替え ( ピボット選択 ) を行なっている。よって、行列の LU 分解の用途には利用できない。

  a : LU 分解する対象の行列
  l : 分解した下三角行列 L
  u : 分解した上三角行列 U
  index : 処理前に対する処理後の式の位置 ( NULL ならば無視 )

  行列 a が LU 分解できない場合 ( 処理途中で対角成分にゼロが見つかった場合 ) は false を返す
*/
template< typename T >
bool Crout( LinearEquationSystem< T > a,
            LinearEquationSystem< T >* l, LinearEquationSystem< T >* u,
            std::vector< typename LinearEquationSystem< T >::size_type >* index = 0 )
{
  assert( l != 0 && u != 0 );

  typedef typename LinearEquationSystem< T >::size_type size_type;

  size_type sz = a.size();
  l->resize( sz );
  u->resize( sz );
  if ( index != 0 ) index->resize( sz );

  for ( size_type i = 0 ; i < sz ; ++i ) {
    (*u)[i][i] = 1;         // Uの対角成分
    (*l)[i][0] = a[i][0];   // Lの一列目
    l->ans( i ) = u->ans( i ) = a.ans( i );
    if ( index != 0 ) (*index)[i] = i;
  }

  for ( size_type i = 1 ; i < sz ; ++i ) {
    size_type swapped = PivotRow( l, i - 1 );
    if ( (*l)[i - 1][i - 1] == 0 ) return( false );
    if ( swapped != i - 1 ) {
      a.swapEquation( swapped, i - 1 );
      if ( index != 0 )
        std::swap( (*index)[swapped], (*index)[i - 1] );
    }
    for ( size_type j = i ; j < sz ; ++j ) {
      // U成分の算出
      (*u)[i - 1][j] = a[i - 1][j];
      for ( size_type k = 0 ; k < i - 1 ; ++k )
        (*u)[i - 1][j] -= (*l)[i - 1][k] * (*u)[k][j];
      (*u)[i - 1][j] /= (*l)[i - 1][i - 1];

      // L成分の算出
      (*l)[j][i] = a[j][i];
      for ( size_type k = 0 ; k < i ; ++k )
        (*l)[j][i] -= (*l)[j][k] * (*u)[k][i];
    }
  }

  return( true );
}

LU 法でもピボット選択は有効になります。ドゥーリトル法の場合、上三角行列 U の対角成分を使って徐算をしているので、徐算に使う対角成分を含む行にある 0 以外の成分 ( 対角成分より右側 ) の中から最大のものを利用するのが最適となります。よって、列の交換が必要になり、そのために関数 PivotCol を利用しています。また、クラウト法は、下三角行列 L の対角成分を使うため、対角成分を含む列にある要素を使うことになります。行の交換が必要になるので、先に用意した Pivot がそのまま利用できそうですが、Pivot では対象の式の係数の最大値を使って除算した上で比較するようになっていました。今回は対角成分で除算するだけなので、単純に最大値を選択すればよく正規化は不要です。そのため、正規化処理を除いた PivotRow を新たに用意しています。
ドゥーリトル法の場合、列の交換をするので未知数自体が交換されます。よって、どの列を交換したのかを保持して、解を求めた後に位置を修正する必要があります。クラウト法では行の交換になるので、右辺があらかじめ代入されていれば後で解を入れ替える必要はありません。しかし、逆行列を求めるような場合は、LU 分解をしてから右辺に値を代入して連続的に解く形になるので、やはり入れ替え位置を記憶しておいて、右辺に値を代入するときに補正をする必要が生じます。もし、左辺の式の位置だけ変更して右辺の値をそのままにしてしまうと、元の連立方程式と一致しなくなります。


LU 分解ができれば、あとは解を求める処理のみです。

/*
  ForwardSub : 上側から解を求める ( 下三角行列用 )

  l : 下三角行列
  u : 上三角行列

  u は、求めた解を代入するだけのために渡される。
*/
template< class T >
void ForwardSub( LinearEquationSystem< T >* l, LinearEquationSystem< T >* u )
{
  typedef typename LinearEquationSystem< T >::size_type size_type;

  for ( size_type r = 0 ; r < l->size() ; ++r ) {
    T t = l->ans( r );
    for ( size_type i = 0 ; i < r ; ++i )
      t -= (*l)[r][i] * l->ans( i );
    l->ans( r ) = u->ans( r ) = t / (*l)[r][r];
  }
}

/*
  BackSub : 下側から解を求める ( 上三角行列用 )

  u : 上三角行列
*/
template< class T >
void BackSub( LinearEquationSystem< T >* u )
{
  typedef typename LinearEquationSystem< T >::size_type size_type;

  for ( size_type r = u->size() ; r > 0 ; --r ) {
    T t = u->ans( r - 1 );
    for ( size_type i = u->size() ; i > r ; --i )
      t -= (*u)[r - 1][i - 1] * u->ans( i - 1 );
    u->ans( r - 1 ) = t / (*u)[r - 1][r - 1];
  }
}

上の二つの関数は、ガウスの消去法の中の後退代入を、上側または下側から行って解を求めているだけです。先に ForwardSub へ上三角行列 L と下三角行列 U を渡すと、L を使って解を求め、その結果を U の右辺に代入します。次に BackSubU を使って解を求めます。


LU 分解と後退代入を組み合わせて連立方程式を解くサンプル・プログラムです。

/*
  LUdecomposition_Doolittle : LU 分解 ( ドゥーリトル法 ) を利用した連立方程式の計算

  les : 対象の連立方程式

  係数行列が LU 分解できない場合は false を返す
*/
template< class T >
bool LUdecomposition_Doolittle( LinearEquationSystem< T >& les )
{
  assert( &les != 0 );

  typedef typename LinearEquationSystem< T >::size_type size_type;

  LinearEquationSystem< T > l, u;
  std::vector< size_type > index;
  if ( ! Doolittle( les, &l, &u, &index ) )
    return( false );
  ForwardSub( &l, &u );
  BackSub( &u );

  for ( size_type i = 0 ; i < u.size() ; ++i )
    les.ans( index[i] ) = u.ans( i );

  return( true );
}

/*
  LUdecomposition_Crout : LU 分解 ( クラウト法 ) を利用した連立方程式の計算

  les : 対象の連立方程式

  係数行列がLU分解できない場合は false を返す
*/
template< class T >
bool LUdecomposition_Crout( LinearEquationSystem< T >& les )
{
  assert( &les != 0 );

  typedef typename LinearEquationSystem< T >::size_type size_type;

  LinearEquationSystem< T > l, u;
  if ( ! Crout( les, &l, &u ) )
    return( false );
  ForwardSub( &l, &u );
  BackSub( &u );

  for ( size_type i = 0 ; i < u.size() ; ++i )
    les.ans( i ) = u.ans( i );

  return( true );
}

前述の通り、ドゥーリトル法を利用した場合は未知数の入れ替えが発生するため、元の位置に戻す処理が必要になります。クラウト法ではその必要はありませんが、LU 分解をして複数の解を使って一気に後退代入を行う場合は ( 後で右辺を代入して処理することになるので ) やはり元に戻す処理が必要となることに注意してください。


今回紹介したアルゴリズムは、連立方程式を手作業で解くやり方をそのまま利用しています。このような処理方法は「直接法」と呼ばれます。それとは別に、「反復解法」と呼ばれるアルゴリズムがあるので、次はこの「反復解法」を紹介したいと思います。


少し、宣伝です。

趣味で小説を書いています。以前から、書いてみたいという思いがあり、アマチュアの小説を扱うサイトを見つけたことから実際に始めてみました。

最初は、小説を書ければ満足で、アクセス数が少なくてもいいかなと思っていましたが、やはりせっかく書いたからには読んでもらいたいし、感想をいただきたいという欲も出てきました。
もし、よろしければ、読んでもらえると嬉しいですし、今後の励みにもなります。

小説は、以下のリンク先から読むことができます。

「鬼退治」
「瑠璃色の海と空」

<参考文献>
  1. 「これなら分かる応用数学教室」 金谷健一著 (共立出版)
  2. 「Yamamoto's Laboratory」 「4 ガウス・ジョルダン法」 ... 他にも様々な内容が「講義ノート」として公開されています。
  3. Wikipedia ... 今回はほとんどの内容が Wikipediaから得られました。

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