数値演算法

(10) 極値計算法

関数の最大・最小値を求める問題は様々な場面で必要になります。このとき、全定義域に関して最大・最小値を求める代わりに局所的な最大・最小値、すなわち極値を計算する手法がよく用いられます。今回は、極値の計算を行うための方法について紹介します。

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

1) 勾配法 ( Gradient Method )

一変数の関数 f(x) に対し、その極大値・極小値を求めるためには、f(x) の導関数 f'(x) を計算して f'(x) = 0 となる x を求め、その点での f(x) から得るのが一般的な方法です。例えば

f(x) = x( x - 3 )2

の導関数は

f'(x) = ( x - 3 )2 + 2x( x - 3 ) = 3( x - 1 )( x - 3 )

なので、極値は x = 1, 3 のときで、f(1) = 4, f(3) = 0 となります。極小値か極大値かを判別したければ、f'(x) の正負から関数の増減を調べたり、場合によっては二階の導関数 f''(x) を求めて利用することもあります。

しかし、導関数 f'(x) が複雑で極値をとるときの x の値を求めることができない場合、この方法を使うことはできません。このとき利用される数値計算法に「勾配法 ( Gradient Method )」があります。

まず、極値に近いと思われる点 x0 を初期値とします。f'(x0) ≅ 0 ならばその点で極値をとるのでそれで処理は終了です。極大値を求めたい場合、f'(x0) > 0 なら x 軸上を正の方向へ、f'(x0) < 0 なら負の方向へ進みます。極小値を求めたい場合はその逆になります。このとき、どの程度進めばよいか、その距離が問題となりますが、一つの方法として正の定数値 h を決めて進んだ点の導関数値 f'( x0 + h ) ( 負の方向へ進むのなら f'( x0 - h ) ) を調べ、極大値を求める場合 f'( x0 + h ) > 0 ならばさらに 2h 進め、それでも符号が反転しなければさらに 4h ... という具合に導関数値が正である限り h を倍々していき、符号が反転したら直前で停止して、その位置を次の位置 x1、同じく直前の h を次の h とします。逆にもし f'( x0 + h ) < 0 なら進み過ぎということになるので半分の x0 + h / 2 を調べ、それでも符号が反転しなければ x0 + h / 4 ... という具合に h を半分ずつにしていき、符号が反転したらその位置を次の位置 x1、直前の h を次の h とします。これを繰り返し、数列 x0, x1, ... を計算していくと、やがて f'(xn) ≅ 0 となる xn が求められます。実際には完全に一致させることはできないので、収束判定のためのしきい値 ε を定めてそれより小さくなったら処理を終了させます。

図 1-1. 単純な勾配法
単純な勾配法

極値の求め方は他にも、以前紹介した「二分法 ( Bisection Method )」や「ニュートン-ラフソン法 ( Newton-Raphson Method )」などがあります。例えば「ニュートン-ラフソン法」では f'(x) と f''(x) を使って次の漸化式が成り立ちます。

xn+1 = xn - f'(xn) / f''(xn)

しかし、この場合は二階導関数 f''(x) が必要になるため、二階導関数に近似式

f''(xn) ≅ [ f'(xn) - f'(xn-1) ] / ( xn - xn-1 )

を利用すると、以下のように式を変形することができます。

xn+1=xn - f'(xn) / f''(xn)
=xn - f'(xn)( xn - xn-1 ) / [ f'(xn) - f'(xn-1) ]
={ [ f'(xn) - f'(xn-1) ]xn - f'(xn)( xn - xn-1 ) } / [ f'(xn) - f'(xn-1) ]
={ f'(xn)xn-1 - f'(xn-1)xn } / [ f'(xn) - f'(xn-1) ]

こうすると計算に xn と xn-1 の二つが必要になりますが、二階導関数なしで処理することが可能になります。この手法を「セカント法 ( Secant Method )」といいます。

図 1-2. セカント法
セカント法

二分法を改良した手法に「はさみうち法 ( False Position Method )」というものがあります。二分法の場合は f(x) = 0 となる点 x を含む範囲を半分ずつ小さくしていく方式でした。はさみうち法では範囲を示す二点 ( x0, f'(x0) ) と ( x1, f'(x1 ) ) を結んだ直線と x 軸との交点の x 座標 xc を次の範囲の始点または終点とします。f'(xc) の符号が始点と等しい時、f'(xc) = 0 となる点は終点側にあるので xc を新たな始点とし、そうでなければ新たな終点とします。

図 1-3. はさみうち法
はさみうち法

「セカント法」や「はさみうち法」はどちらも直線に近似して x 軸との交点を求める方法と考えることができます。直線ではなく二次関数に近似して交点を求める方法に「逆二次補間法 ( Inverse Quadratic Interpolation )」というものがあります。しかし、二次関数に近似するためには二点だけでは足りないためもう一点追加する必要があります。また、三点を定めたとしてもそれを通る二次関数を求めるには連立方程式を解く必要があるため、それを避けるために「ラグランジュの多項式 ( Lagrange Polynomial )」を利用します (*1-1)。

三次のラグランジュ多項式は次のようになります。

f(x) = ( x - x1 )( x - x2 )y0 / ( x0 - x1 )( x0 - x2 ) + ( x - x2 )( x - x0 )y1 / ( x1 - x2 )( x1 - x0 ) + ( x - x0 )( x - x1 )y2 / ( x2 - x0 )( x2 - x1 )

上式において、f(x0) = y0, f(x1) = y1, f(x2) = y2 が成り立つので、f(x) は三点 ( x0, y0 ), ( x1, y1 ), ( x2, y2 ) を通る二次関数になります。しかし、この式では x = 0 のときの y の値 f(0) は求められても y = 0 のときの x の値は求められないため、x と y を反転した f(y) を使うことにします。

f(y) = ( y - y1 )( y - y2 )x0 / ( y0 - y1 )( y0 - y2 ) + ( y - y2 )( y - y0 )x1 / ( y1 - y2 )( y1 - y0 ) + ( y - y0 )( y - y1 )x2 / ( y2 - y0 )( y2 - y1 )

この式に y = 0 を代入した式 y(0) を計算すると x 軸との交点を求めることができるので、セカント法のようにそれを次の点として更新すれば解に収束していきます。

図 1-4. 逆二次補間法
逆二次補間法

いずれの方法においても極値を求めるため f'(x) を対象に説明していますが、f(x) = 0 となる x の値を求めるときにも当然利用可能です。特に、f(x) = x2 - C とすれば、√C を数値計算として求めるアルゴリズムとなります。


それぞれの手法について、サンプル・プログラムを以下に示します。

/*
  GradientMethod : 勾配法 ( Gradient Method )

  テンプレート引数 F は関数の型を表す。
  収束判定値 e は f(x) がゼロに収束したか、ステップ幅がゼロに収束したかのいずれにも利用される。

  f : 関数 f(x)
  x0 : x の初期値
  h : ステップ幅
  e : 収束判定値
*/
template< class F >
  double GradientMethod( F f, double x0, double h, double e )
{
  while ( std::isfinite( x0 ) && std::abs( f( x0 ) ) > e && std::abs( h ) > e ) {
    double x1 = x0 + h;
    // f(x0) と f(x1) が同符号なら h を二倍にしてさらに進める
    if ( f( x0 ) * f( x1 ) > 0 ) {
      do {
        h *= 2;
        x0 = x1;
        x1 += h;
      } while ( std::isfinite( x0 ) && f( x0 ) * f( x1 ) > 0 );
      h /= 2;
    // f(x0) と f(x1) が異符号なら h を半分にして戻す
    } else {
      do {
        h /= 2;
        x1 -= h;
      } while ( std::isfinite( x0 ) && f( x0 ) * f( x1 ) < 0 );
      x0 = x1;
      h *= 2;
    }
  }

  return( x0 );
}

/*
  BisectionMethod : 二分法 ( Bisection Method )

  テンプレート引数 F は関数の型を表す。
  収束判定値 e は f(x) がゼロに収束したか、x0, x1 が充分近くなったかのいずれにも利用される。
  x0, x1 の間に解があることを想定している。

  f : 関数 f(x)
  x0, x1 : 解のある範囲
  e : 収束判定値
*/
template< class F >
  double BisectionMethod( F f, double x0, double x1, double e )
{
  if ( f( x0 ) > f( x1 ) ) std::swap( x0, x1 );
  double xm = ( x0 + x1 ) / 2; // x0, x1 の中点

  while ( std::abs( f( xm ) ) > e && std::abs( x1 - x0 ) > e ) {
    // 中点が負の場合は解は x1 側にあるので範囲を ( xm, x1 ) にする
    if ( f( xm ) < 0 )
      x0 = xm;
    // 中点が正の場合は解は x0 側にあるので範囲を ( x0, xm ) にする
    else
      x1 = xm;
    if ( f( x0 ) > f( x1 ) ) std::swap( x0, x1 );
    xm = ( x0 + x1 ) / 2;
  }

  return( xm );
}

/*
  SecantMethod : セカント法 ( Secant Method )

  テンプレート引数 F は関数の型を表す。
  収束判定値 e は f(x) がゼロに収束したかに利用される。

  f : 関数 f(x)
  x0 : x の初期値
  x1 : 直線を求めるためのもう一つの点
  e : 収束判定値
*/
template< class F >
  double SecantMethod( F f, double x0, double x1, double e )
{
  while ( std::abs( f( x1 ) ) > e ) {
    // 直線 ( x0, x1 ) と x 軸との交点を求める
    double x2 = ( f( x1 ) * x0 - f( x0 ) * x1 ) / ( f( x1 ) - f( x0 ) );
    x0 = x1;
    x1 = x2;
  }

  return( x1 );
}

/*
  FalsePositionMethod : はさみうち法 ( False Position Method )

  テンプレート引数 F は関数の型を表す。
  収束判定値 e は f(x) がゼロに収束したか、x0, x1 が充分近くなったかのいずれにも利用される。
  x0, x1 の間に解があることを想定している。

  f : 関数 f(x)
  x0, x1 : 解のある範囲
  e : 収束判定値
*/
template< class F >
  double FalsePositionMethod( F f, double x0, double x1, double e )
{
  if ( f( x0 ) > f( x1 ) ) std::swap( x0, x1 );

  // 直線 ( x0, x1 ) と x 軸との交点
  double xm = ( f( x1 ) * x0 - f( x0 ) * x1 ) / ( f( x1 ) - f( x0 ) );

  while ( std::abs( f( xm ) ) > e && std::abs( x1 - x0 ) > e ) {
    // 交点が負の場合は解は x1 側にあるので範囲を ( xm, x1 ) にする
    if ( f( xm ) < 0 )
      x0 = xm;
    // 交点が正の場合は解は x0 側にあるので範囲を ( x0, xm ) にする
    else
      x1 = xm;
    if ( f( x0 ) > f( x1 ) ) std::swap( x0, x1 );
    xm = ( f( x1 ) * x0 - f( x0 ) * x1 ) / ( f( x1 ) - f( x0 ) );
  }

  return( xm );
}

/*
  InverseQuadraticInterpolation : 逆二次補間法 ( Inverse Quadratic Interpolation )

  テンプレート引数 F は関数の型を表す。
  収束判定値 e は f(x) がゼロに収束したかに利用される。

  f : 関数 f(x)
  x0, x2 : 二次関数を求めるときに利用する二点
  e : 収束判定値
*/
template< class F >
  double InverseQuadraticInterpolation( F f, double x0, double x2, double e )
{
  double x1 = ( x0 + x2 ) / 2; // 二次関数を求めるためのもう一点

  while ( std::abs( f( x2 ) ) > e ) {
    double f0 = f( x0 );
    double f1 = f( x1 );
    double f2 = f( x2 );
    double x3 = ( f0 != f1 && f1 != f2 ) ?
      -( x0 * f1 * f2 * ( f1 - f2 ) + x1 * f2 * f0 * ( f2 - f0 ) + x2 * f0 * f1 * ( f0 - f1 ) )
      / ( ( f0 - f1 ) * ( f1 - f2 ) * ( f2 - f0 ) ) :
      ( f2 * x0 - f0 * x2 ) / ( f2 - f0 );
    x0 = x1;
    x1 = x2;
    x2 = x3;
  }

  return( x2 );
}

処理の流れに関しては前述した通りですが、いくつか注意する必要のあるところがあります。

まず、勾配法 ( GradientMethod ) に対しては進む方向は変数 h のステップ幅で決まります。正しい方向に進んでいるかどうかはチェックしていないので、解のない方向へ進めてしまった場合はいつまでも処理が戻らないことになります。サンプル・プログラムの中では x の値が発散してしまったら処理を終了するようにしていますが、このときは当然、正しい結果は得られないことになります。なお、f(x) にだけでなく導関数 f'(x) も引数として渡せるなら、初期値に対して f'( x0 ) と f( x0 ) が異符号なら正の方向へ、同符号なら負の方向へ進むように h を補正すれば正しい方向へ進むことができます。

二分法 ( BisectionMethod ) の場合、引数として渡す二つの x 座標 x0, x1 の間に解が存在することを前提としています。そうではない場合、正しい結果を得ることはできません。セカント法 ( SecantMethod ) とはさみうち法 ( FalsePositionMethod )、逆二次補間法 ( InverseQuadraticInterpolation ) においても、初期値の与え方によっては解が収束せず発散してしまう可能性はあります。

セカント法とはさみうち法はどちらも直線と x 軸の交点を求めるため式が全く同じですが、前者が x0, x1 を順次更新していくのに対して後者は二分法と同様に二座標間の点を次の新たな点としているところが異なります。また、逆二次補間法は二次関数を使ったセカント法と考えることができます。


それぞれの手法を用いて 2 の平方根を求めてみます。使用する関数は以下のようになります。

double sq2( double x )
{
  return( x * x - 2.0 );
}

下表は、各手法での収束の様子をまとめたものです。どの手法も初期値は 0 としています。

表 1-1. 各手法の収束の様子
勾配法 ( x0 = 0, h = 0.5 )二分法 ( x0 = 0, x1 = 2 )セカント法 ( x0 = 0, x1 = 2 )はさみうち法 ( x0 = 0, x1 = 2 )逆二次補間法 ( x0 = 0, x1 = 2 )
xf(x)xf(x)xf(x)xf(x)xf(x)
10.5-1.751-11-11-11.6666666670.7777777778
21-11.50.251.333333333-0.22222222221.333333333-0.22222222221.401515152-0.03575528007
31.25-0.43751.25-0.43751.4285714290.040816326531.4-0.041.41389545-0.0008996552786
41.375-0.1093751.375-0.1093751.413793103-0.0011890606421.411764706-0.0069204152251.4142137886.391182188E-07
51.40625-0.02246093751.43750.066406251.414211438-6.007286839E-061.413793103-0.0011890606421.4142135626.49702514E-13
61.4140625-0.00042724609381.40625-0.02246093751.4142135638.931460016E-101.414141414-0.0002040608101
71.41418457-8.200109005E-051.4218750.021728515621.414213562-4.440892099E-161.414201183-3.501277966E-05
81.414199829-3.884336911E-051.4140625-0.00042724609381.414211438-6.007286839E-06
91.414207458-1.726433402E-051.417968750.010635375981.414213198-1.030688757E-06
101.414211273-6.474772817E-061.4160156250.0051002502441.4142135-1.768382716E-07
111.414213181-1.079981303E-061.4150390620.0023355484011.414213552-3.034065155E-08
121.414213419-4.056318517E-071.4145507810.0009539127351.414213561-5.205633036E-09
131.414213538-6.845708356E-081.4143066410.0002632737161.414213562-8.931462236E-10
141.414213553-2.631023555E-081.41418457-8.200109005E-051.414213562-1.532394212E-10
151.414213561-5.236811429E-091.4142456059.063258767E-051.414213562-2.629119145E-11
161.414213561-2.602633442E-091.4142150884.314817488E-06
171.414213562-1.285544338E-091.414199829-3.884336911E-05
181.414213562-6.269997854E-101.414207458-1.726433402E-05
191.414213562-2.977276203E-101.414211273-6.474772817E-06
201.414213562-1.330915378E-101.414213181-1.079981303E-06
211.414213562-5.07733855E-111.4142141341.617417183E-06
221.4142136572.68717713E-07
231.414213419-4.056318517E-07
241.414213538-6.845708356E-08
251.4142135981.001303112E-07
261.4142135681.583661291E-08
271.414213553-2.631023555E-08
281.414213561-5.236811429E-09
291.4142135645.299900963E-09
301.4142135623.154454475E-11

セカント法と逆二次補間法の収束の速さが際立っています。この二つの手法は初期値が適切であれば収束が速いという特徴があります。しかし前述の通り、初期値の与え方によっては値が発散してしまうこともあるので注意が必要です。二分法は収束が遅く、指定した二つの値の間に解がある必要がありますが、それが保証されている限りほぼ確実に正しい解へ収束します。


*1-1) 「ラグランジュの多項式 ( Lagrange Polynomial )」については「グラフィック・パターンの扱い (5) サンプル補間」の中の「6) 多項式補間(Polynomial Interpolation)」で詳しく紹介しています。


2) 多変量の勾配法

まずは二変数の場合を考えてみます。例えば、

f( x, y ) = ( x - 2 )2 + ( y - 1 )2

に対してその極値をとる点は ( 2, 1 ) であり、その値が 0 になることはすぐにわかります。これを計算するためには、偏微分 ∂f / ∂x, ∂f / ∂y を求めて

∂f / ∂x = 2( x - 2 )

∂f / ∂y = 2( y - 1 )

より ∂f / ∂x = 0, ∂f / ∂y = 0 から得ることができます。しかし、式が複雑になると簡単に求めることはできなくなります。そこで使われる手法に「山登り法 ( Hill Climbing )」があります。

まず、平面上に初期値 ( x0, y0 ) を定めます。その点からある方向に向けて直線をひいてその断面を見ると、それは一変数関数と見ることができます。そこで、ある方向へ関数の値を調べながらたどり、極大値を求めたいのなら値の最大となる点へ登って行き、逆に極小値の場合なら値の最小となる点へ下っていきます。極大・極小点が見つかったら、また同じように今度は違う方向へ走査していきます。一変数関数の極大・極小値なら、今まで紹介した方法をそのまま利用することができます。

ある点から走査する方向ベクトルを a = ( ax, ay ) としたとき、その走査線上にある関数の値は、任意の値 t を使って f( x0 + axt, y0 + att ) で求めることができます。これを t の一変数関数 F(t) とおいて計算を行います。この時、dF / dt は合成関数の微分の公式から

dF / dt = ( ∂f / ∂x )( dx / dt ) + ( ∂f / ∂y )( dy / dt ) = ax( ∂f / ∂x ) + ay( ∂f / ∂y )

なので、dF / dt = 0 となる t を求め、それを使って次の点 ( x0 + axt, y0 + att ) が得られます。

図 2-1. 線形探索の様子
線形探索

問題はどの方向に走査すればよいかということになります。ランダムに走査しても極点に到達できないので、( x0, y0 ) において最も傾きの大きな方向へ進むことにします。f( x, y ) の ( x0, y0 ) でのテイラー-マクローリン展開は

f( x, y ) = f( x0, y0 ) + ( ∂f0 / ∂x )( x - x0 ) + ( ∂f0 / ∂y )( y - y0 ) + ...

となります。但し、∂f0 / ∂x, ∂f0 / ∂y はそれぞれ ∂f / ∂x, ∂f / ∂y の ( x0, y0 ) での値を表し、... 以降は ( x - x0 )2, ( x - x0 )( y - y0 ), ( y - y0 )2 ... と二次以上の高次項が並ぶのを省略して書いています。これらの項は、( x - x0 ), ( y - y0 ) が非常に小さいことから無視できるとすれば、上式は平面の方程式に近似することができます。

点 ( x0, y0 ) から x 軸、y 軸方向へそれぞれ Δx, Δy だけ移動したとき、f( x, y ) の大きさは

f( x0 + Δx, y0 + Δy ) - f( x0, y0 )={ f( x0, y0 ) + ( ∂f0 / ∂x )[ ( x0 + Δx ) - x0 ] + ( ∂f0 / ∂y )[ ( y0 + Δy ) - y0 ] } - f( x0, y0 )
=( ∂f0 / ∂x )Δx + ( ∂f0 / ∂y )Δy

だけ変化します。∇f0 = ( ∂f0 / ∂x, ∂f0 / ∂y )、Δ = ( Δx, Δy ) とすれば、上式は二つのベクトルの内積 ( ∇f0, Δ ) で表されます。二つのベクトルのなす角を θ としたとき、

( ∇f0, Δ ) = ||∇f0||・||Δ||cosθ

なので、f( x, y ) の変化量は θ = 0 のとき、すなわち Δ∇f0 の方向へ進むときに最大となり、最も傾きが大きくなります。これを「勾配 (Gradient)」といいます。逆に、f( x, y ) が変化しない Δ は「等高線 (Contour)」と呼ばれます。このとき θ = π / 2 であり、勾配 ∇f0 と等高線は直交することになります。

図 2-2. 平面上の増加量は ∇f0 へ進むときに最大
平面での最大増加量

∇f0 の方向に進むとき、dF/ dt は

dF / dt = ( ∂f / ∂x )( ∂f0 / ∂x ) + ( ∂f / ∂y )( ∂f0 / ∂y ) = ( ∇f, ∇f0 )

と表されますが、dF / dt = 0 となることは ∇f∇f0 が直交することを意味します。走査する直線は勾配 ∇f0 方向なので、dF / dt = 0 となる位置で等高線に接し、次はその接点から直交する方向へ ( その方向は ∇f で決まります ) 走査すればよいことになります。この様子を下図に示します。

図 2-3. 二変量での勾配法
二変量での勾配法

変数が増えた場合も考え方は同様で、∇fk = ( ∂f(k) / ∂x1, ∂f(k) / ∂x2, ... ∂f(k) / ∂xn ) ( 但し、∂f(k) / ∂xi は点 x(k) = ( x1(k), x2(k), ... xn(k) ) での ∂f / ∂xi の値 ) 方向へ線形探索を行い、極点が見つかったら ( ∇fk+1, ∇fk ) = 0 となる ∇fk+1 を次の探索方向として走査を繰り返します。走査の終了は ∇fk+1 = 0 となったときですが、厳密にゼロになることはないのであるしきい値より小さくなったらそれで処理を終了します。これが、多変量での勾配法のアルゴリズムとなります。


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

/**
   HC_Operator : 山登り法 ( Hill Climbing ) 用多変数関数オブジェクト
**/
class HC_Operator
{
  std::vector< double > x0_; // 初期値

  // fx : f(x) の値を求める
  virtual double fx( const std::vector< double >& x ) const = 0;

public:

  // 仮想デストラクタ(何もしない)
  virtual ~HC_Operator() {}

  /*
    nabla : ∇f の値を求める

    x : x の値(多変量)
    i : x の何番目のパラメータに対する偏微分を返すかを指定する添字
  */
  virtual double nabla( const std::vector< double >& x, size_t i ) const = 0;

  // assign : 初期値 x0 を代入する
  void assign( const std::vector< double >& x0 )
  { x0_.assign( x0.begin(), x0.end() ); }

  // ft : F(t) = f(x0 + t∇f(x0)) の値を求める
  double ft( double t ) const
  {
    std::vector< double > x( x0_ );
    for ( std::vector< double >::size_type i = 0 ; i < x0_.size() ; ++i )
      x[i] += t * nabla( x0_, i );

    return( fx( x ) );
  }

  // operator() : dF(t)/dt = ( ∇f(x0 + t∇f0), ∇f(x0) ) の値を求める
  double operator()( double t ) const
  {
    std::vector< double > x( x0_ );
    for ( std::vector< double >::size_type i = 0 ; i < x0_.size() ; ++i )
      x[i] += t * nabla( x0_, i );

    double ans = 0;
    for ( std::vector< double >::size_type i = 0 ; i < x0_.size() ; ++i )
      ans += nabla( x0_, i ) * nabla( x, i );

    return( ans );
  }
};

/*
  HillClimbing : 山登り法 ( Hill Climbing ) による多変量関数の極値計算

  テンプレート引数の Op は多変量関数オブジェクト HC_Operator の派生型であることを想定

  op : 関数 f(x)
  x0 : 初期値へのポインタ (結果をそのまま返す)
  h : ステップ幅
  e : 収束判定値
*/
template< class Op >
  void HillClimbing( Op& op, std::vector< double >* x0, double h, double e )
{
  std::vector< double > dx( x0->size() );

  for ( ; ; ) {
    op.assign( *x0 );
    h = ( ( op( 0.0 ) * op.ft( 0.0 ) < 0 ) ? 1 : -1 ) * std::abs( h );
    double t = GradientMethod( op, 0.0, h, e );

    for ( std::vector< double >::size_type i = 0 ; i < x0->size() ; ++i )
      dx[i] = t * op.nabla( *x0, i );

    bool isConvergence = true;
    for ( std::vector< double >::size_type i = 0 ; i < x0->size() ; ++i ) {
      (*x0)[i] += dx[i];
      if ( dx[i] > e ) isConvergence = false;
    }

    if ( isConvergence ) break;
  }
}

一変量の勾配法をそのまま利用するためかなり強引な実装になっています。HC_Operator は多変量関数用の基底クラスです。可変長配列の x0_ がメンバ変数で、これは x の初期値を表します。fx が f(x) の値を、nabla∇f を求めるためのメンバ関数で、どちらも純粋仮想関数として定義されています。∇f = ( ∂f / ∂x1, ∂f / ∂x2, ... ∂f / ∂xi, ... ) より nablax の値の他に何番目の要素かを表す添字 i を引数として渡します。また、ft は引数 t に対して F(t) = f( x0 + t∇f(x0) ) を求めるための関数、operator() は dF / dt = ( ∇f(x0 + t∇f0), ∇f(x0) ) を求めるための関数で、fx, nabla を利用して計算しています。最後に、assign で初期値 x0_ を代入することができます。

HillClimbing は山登り法を行うための関数です。テンプレート引数の Op は、先ほど定義したクラス HC_Operator から派生したインスタンスであることを前提としています。そうであれば最初から HC_Operator への参照かポインタで渡せば済むように見えますが、一変数関数用の勾配法 GradientMethod がテンプレート引数を持つためそれはできません。そうするためには、一変数の勾配法からテンプレートではなく基底クラスへの参照かポインタを介して関数を渡すように再定義する必要があります。
関数の中で行っていることは、xk を初期値として一変数用の勾配法 GradientMethod を呼び出して dF / dt = 0 になる t を求め、xk+1 = xk + t∇f(xk) として初期値を更新する処理を、t∇f(xk) が収束するまで繰り返しているだけです。注意点として、勾配法で利用するステップ幅 h を、F'(t0) と F(t0) が異符号なら正に、同符号なら負に変換しています。GradientMethod そのものは h の符号に関してチェックしていないので、外側でチェックをしているわけです。なお、t の初期値 t0 はゼロ固定としています。


HC_Operator からの派生クラスの例として、楕円体の関数オブジェクトのサンプル・プログラムを以下に示します。

/**
   HC_Oval : 楕円体の関数オブジェクト
**/
class HC_Oval : public HC_Operator
{
  double a_, b_;   // 扁平率
  double x0_, y0_; // 中心の座標

  /*
    fx : f(x) = a( x - x0 )^2 + b( y - ya )^2 を求める

    ( x, y ) = ( x[0], x[1] ) として計算する
  */
  virtual double fx( const std::vector< double >& x ) const
  {
    assert( x.size() >= 2 );

    return( a_ * pow( x[0] - x0_, 2 ) + b_ * pow( x[1] - y0_, 2 ) );
  }

public:

  // 扁平率 a, b と中心座標 ( x0, y0 ) を指定して構築
  HC_Oval( double a, double b, double x0, double y0 )
    : a_( a ), b_( b ), x0_( x0 ), y0_( y0 ) {}

  /*
    ∇f : ( ∂f/∂x, ∂f/∂y ) を求める

    ( x, y ) = ( x[0], x[1] ) として計算する
    i = 0 で ∂f/∂x、i = 1 で ∂f/∂y を返す
  */
  virtual double nabla( const std::vector< double >& x, size_t i ) const
  {
    assert( x.size() >= 2 );
    assert( i < 2 );

    if ( i == 0 )
      return( 2 * a_ * ( x[0] - x0_ ) );
    else
      return( 2 * b_ * ( x[1] - y0_ ) );
  }
};

楕円体は以下に示すような式で表されます。

f( x, y ) = a( x - x0 )2 + b( y - y0 )2

その偏微分は次のとおりです。

∂f / ∂x = 2a( x - x0 )

∂f / ∂y = 2b( y - y0 )

但し、a と b の符号が等しくなければ楕円体にはなりません。異符号の場合は双曲型、どちらかがゼロの場合は放物型となります。

図 2-4. a( x - x0 )2 + b( y - y0 )2 の形
楕円型放物型双曲型
楕円型 放物型 双曲型

楕円体ならば、極値は明らかに ( x0, y0 ) となります。下記のような例での出力結果は次のようになります。

■ サンプル・コード (1)
HC_Oval oval( 2, 1, 1.5, 2.5 );
std::vector< double > x0( 2 );
HillClimbing( oval, &x0, 0.5, 1E-10 );
cout << "( " << x0[0] << ", " << x0[1] << " )" << endl;
■ 出力結果 (1)
( 1.5, 2.499997856 )

ところが、双曲型になるように引数を変更すると、値は発散してしまいます。双曲型は極値を持たないためです。

■ サンプル・コード (2)
HC_Oval oval( 2, -1, 1.5, 2.5 );
std::vector< double > x0( 2 );
HillClimbing( oval, &x0, 0.5, 1E-10 );
cout << "( " << x0[0] << ", " << x0[1] << " )" << endl;
■ 出力結果 (2)
( -3, 3.678732859e+297 )

3) 共役勾配法 ( Conjugate Gradient Method )

線形重回帰モデル」における正規方程式は

XTXa = XTy

で表されます。ここで X は独立変数ベクトルからなる行列で「デザイン行列」と呼ばれ、y は従属変数からなるベクトル、a は求めたい係数です。この式は連立方程式の形をしていますが、a に対する係数となる XTX は「半正値対称行列」であり、固有値は必ずゼロ以上となります ( 補足 1 )。上式をもっと一般化して、任意の半正値対称行列 A を係数とする連立方程式

Ax = b

の解を求めることを考えると、これは

f(x) = (1/2)( x, Ax ) - ( x, b )

を最小にすることと同値となります。そこで、連立方程式を解く代わりに f(x) を最小にすることで解を得られないか検討してみます。

f(x) の偏微分からなるベクトル ∇f = ( ∂f / ∂x1, ∂f / ∂x2, ... ∂f / ∂xn ) は、

∇f = Ax - b

であり、∇f = 0 のとき f(x) は最小値を持ちます ( 補足 2 )。しかし、これは結局、元の連立方程式を解くことと変わりません。そこで、Ax - b = 0 の真の解 x が近似解 x(k) に対して t(k)m(k) だけずれている、すなわち

x = x(k) + t(k)m(k)

が成り立つとして、この式を代入すると

A( x(k) + t(k)m(k) ) - b = 0 より

( Ax(k) - b ) = -t(k)Am(k) --- (1)

となります。Ax(k) - b∇f(k) としたとき、∇f(k) は点 x(k) での勾配を意味します。上式 (1) は、Am(k)∇f(k) と平行であること、すなわち

Am(k)∇f(k) --- (2)

であることを表しています。

x(k) での等高線の接線ベクトルを u(k) とすると、u(k)∇f(k) は直交するので、

( u(k), Am(k) ) = 0 --- (3)

が成り立ちます。u(k)m(k) は、m(k)A を作用させることによって直交します。このようなとき、u(k)m(k)A に関して互いに「共役 (Conjugation)」であるといいます。勾配法を用いて真の解 x に近づきたい場合、u(k) に直交した勾配 ∇f(k) の方向ではなく、共役の関係にある m(k) の方向へ進まなければならないことになります。

(2) より、Am(k)∇f(k) と向きが等しいことは先程述べた通りです。A = E ( 単位行列 ) ならば、m(k)∇f(k) となりますが実際にはそうはなりません。そのズレ量を α(k) で表して

m(k) = ∇f(k) + α(k)u(k) --- (4)

として、これを (3) に代入すると

( u(k), A( ∇f(k) + α(k)u(k) ) ) = ( u(k), A∇f(k) ) + α(k)( u(k), Au(k) ) = 0 より

α(k) = - ( u(k), A∇f(k) ) / ( u(k), Au(k) ) --- (5)

という結果が得られます。

勾配法では、前にたどった勾配方向が次の勾配と直交するのでした。すなわち、u(k) = m(k-1) ということになるので、(4), (5) 式は

m(k)=∇f(k) + α(k)m(k-1)
=( Ax(k) - b ) + α(k)m(k-1)
α(k)=( m(k-1), A( Ax(k) - b ) ) / ( m(k-1), Am(k-1) )

また、(1) より両辺に対して m(k) との内積をとって整理すると

t(k) = -( m(k), Ax(k) - b ) / ( m(k), Am(k) )

となります。以上をまとめて

  1. α(k) = -( m(k-1), A( Ax(k) - b ) ) / ( m(k-1), Am(k-1) )
  2. m(k) = ( Ax(k) - b ) + α(k)m(k-1)
  3. t(k) = -( m(k), Ax(k) - b ) / ( m(k), Am(k) )
  4. x(k+1) = x(k) + t(k)m(k)

の 4 つの反復公式を繰り返し適用することで解を得ることができます。但し、初期値は α(0) = 0, x(0) = 0 として 2 から処理を開始します。

図 3-1. 共役勾配法での各ベクトルの関係
共役勾配法での各ベクトルの関係

前述の通り、極値のある方向は勾配 ∇f ではなく m なので、通常の勾配法の場合ではなかなか極値にたどり着けないときがあります。二変数関数の場合、共役勾配法であれば 2 回の反復処理で極値へ到達してしまいます。一般に、n 変数関数に対しては n 回の反復処理で済むことから、n 個の連立方程式が線形の速さで処理できることを意味します。これは、今まで紹介した直接法や反復法と比較しても非常に有利であることになります。

図 3-2. 勾配法と共役勾配法の進み方の違い
勾配法と共役勾配法の進み方

上述の方法は、係数が半正値対称行列である場合に限られていましたが、これは一般の連立方程式にも応用できます。一般の連立方程式 Ax = b の解を求めるためには

f(x) = (1/2)|| Ax - b ||2

の最小値を計算すればよいので、これを展開すると

f(x)=(1/2)( Ax - b, Ax - b )
=(1/2)( Ax, Ax ) - (1/2)( Ax, b ) - (1/2)( b, Ax ) + (1/2)( b, b )
=(1/2)( x, ATAx ) - ( x, ATb ) + (1/2)|| b ||2

となります。但し、ここで ( Ax, y ) = ( x, ATy ) であることを利用しました。勾配 ∇f を計算すると

∇f = ATAx - ATb

となるので、前述と同様の計算を行うと各反復公式は次のようになります。

  1. α(k) = -( m(k-1), ATAAT( Ax(k) - b ) ) / ( m(k-1), ATAm(k-1) )
  2. m(k) = AT( Ax(k) - b ) + α(k)m(k-1)
  3. t(k) = -( m(k), AT( Ax(k) - b ) ) / ( m(k), ATAm(k) )
  4. x(k+1) = x(k) + t(k)m(k)

共役勾配法による連立方程式の計算用サンプル・プログラムを以下に示します。

/*
  CGM_CalcAlpha : 共役勾配法による連立方程式計算 ( αを求める )

  α(k) = -( m(k-1), aTaaT・ax_b(k) ) / ( m(k-1), aTa・m(k-1) )
*/
template< class T >
  T CGM_CalcAlpha( const std::vector< T >& m, const Matrix< T >& aTaaT, const std::vector< T >& ax_b, const Matrix< T >& aTa )
{
  // 分子の計算
  std::vector< T > aTaaTax_b( m.size() );
  MatrixLib::mult( aTaaT, ax_b.begin(), ax_b.end(), aTaaTax_b.begin(), aTaaTax_b.end() );
  T n = std::inner_product( m.begin(), m.end(), aTaaTax_b.begin(), T() );

  // 分母の計算
  std::vector< T > aTam( m.size() );
  MatrixLib::mult( aTa, m.begin(), m.end(), aTam.begin(), aTam.end() );
  T d = std::inner_product( m.begin(), m.end(), aTam.begin(), T() );

  return( ( d == 0 ) ? 0 : -n / d );
}

/*
  CGM_CalcM : 共役勾配法による連立方程式計算 ( mを求める )

  m(k) = aT・ax_b(k) + alpha(k)・m(k-1)
*/
template< class T >
  void CGM_CalcM( const SquareMatrix< T >& aT, const std::vector< T >& ax_b, T alpha, std::vector< T >* m )
{
  // alpha・m(k-1) の計算
  std::transform( m->begin(), m->end(), m->begin(), std::bind2nd( std::multiplies< T >(), alpha ) );

  // aT・ax_b(k) の計算
  std::vector< T > aTax_b( m->size() );
  MatrixLib::mult( aT, ax_b.begin(), ax_b.end(), aTax_b.begin(), aTax_b.end() );

  // 両者の和を求める
  std::transform( aTax_b.begin(), aTax_b.end(), m->begin(), m->begin(), std::plus< T >() );
}

/*
  CGM_CalcT : 共役勾配法による連立方程式計算 ( tを求める )

  t(k) = -( m(k), aT・ax_b(k) ) / ( m(k), aTa・m(k) )
*/
template< class T >
  T CGM_CalcT( const std::vector< T >& m, const SquareMatrix< T >& aT, const std::vector< T >& ax_b, const Matrix< T >& aTa )
{
  // 分子の計算
  std::vector< T > aTax_b( m.size() );
  MatrixLib::mult( aT, ax_b.begin(), ax_b.end(), aTax_b.begin(), aTax_b.end() );
  T n = std::inner_product( m.begin(), m.end(), aTax_b.begin(), T() );

  // 分母の計算
  std::vector< T > aTam( m.size() );
  MatrixLib::mult( aTa, m.begin(), m.end(), aTam.begin(), aTam.end() );
  T d = std::inner_product( m.begin(), m.end(), aTam.begin(), T() );

  return( ( d == 0 ) ? 0 : -n / d );
}

/*
  CGM_CalcX : 共役勾配法による連立方程式計算 ( xを求める )

  x が収束したかを同時に判定する

  x(k+1) = x(k) + t(k)・ m(k)
*/
template< class T >
  bool CGM_CalcX( T t, const std::vector< T >& m, std::vector< T >* x, T e )
{
  // t(k)・m(k) の計算
  std::vector< T > tm( m );
  std::transform( tm.begin(), tm.end(), tm.begin(), std::bind2nd( std::multiplies< T >(), t ) );

  // x(k+1) = x(k) + t(k)・ m(k) の計算
  std::transform( x->begin(), x->end(), tm.begin(), x->begin(), std::plus< T >() );

  // 収束したかを t(k)・m(k) の値で判断
  for ( typename std::vector< T >::const_iterator i = tm.begin() ; i != tm.end() ; ++i )
    if ( *i > e ) return( false );

  return( true );
}

/*
  ConjugateGradientMethod : 共役勾配法による連立方程式計算

  a : 係数行列 A
  b : 解ベクトル b
  x : 未知数 x を保持するための配列
  e : 収束判定値
*/
template< class T >
  void ConjugateGradientMethod( const SquareMatrix< T >& a, const std::vector< T >& b, std::vector< T >* x, T e )
{
  ErrLib::CheckLinearModel( a, b.begin(), b.end() );

  T alpha = 0;
  SquareMatrix< T > aT( a ); // At
  aT.transpose();
  Matrix< T > aTa; // At・A
  MatrixLib::mult( aT, a, &aTa );
  Matrix< T > aTaaT; // At・A・At
  MatrixLib::mult( aTa, aT, &aTaaT );

  std::vector< T > m( b.size() ); // m
  x->assign( b.size(), T() );

  for ( ; ; ) {
    std::vector< T > ax_b( b.size() ); // Ax - b
    MatrixLib::mult( a, x->begin(), x->end(), ax_b.begin(), ax_b.end() );
    std::transform( ax_b.begin(), ax_b.end(), b.begin(), ax_b.begin(), std::minus< T >() );

    alpha = CGM_CalcAlpha( m, aTaaT, ax_b, aTa ); // α(k) = -( m(k-1), AtAAt( Ax(k) - b ) ) / ( m(k-1), AtA・m(k-1) )
    CGM_CalcM( aT, ax_b, alpha, &m );             // m(k) = At( Ax(k) - b ) + α(k)m(k-1)
    T t = CGM_CalcT( m, aT, ax_b, aTa );          // t(k) = -( m(k), At( Ax(k) - b ) ) / ( m(k), AtA・m(k) )
    if ( CGM_CalcX( t, m, x, e ) ) break;         // x(k+1) = x(k) + t(k)・m(k)
  }
}

内容は前述の反復公式をそのまま当てはめただけなので難しいところはないと思います。補足として、Matrix は任意の型を要素とする行列、SquareMatrix は任意の型を要素とする正方行列を表したクラスです。定義や実装はサンプル・プログラム内では省略していますが、以下の関数やメンバ関数を用いています。

コピー・コンストラクタ
引数の行列の内容をそのままコピーして構築。
void transpose()
行列の転置を行うメンバ関数。
template< class MATRIX1, class MATRIX2, typename Res >
  void MatrixLib::mult( const MATRIX1& mat1, const MATRIX2& mat2, Matrix< Res >* res )
行列どうしの乗算を行う関数。mat1mat2 の乗算結果を res に返す。
template< class MATRIX, class In, class Out >
  void MatrixLib::mult( const MATRIX& mat, In is, In ie, Out os, Out oe )
行列と列ベクトルの乗算を行う関数。行列 mat と列ベクトル [ is, ie ) の乗算結果を配列 [ os, oe ) に返す。
template< class MATRIX, class In >
  void ErrLib::CheckLinearModel( const MATRIX& x, In ys, In ye )
線形モデルとして適切か ( 連立方程式の形になっているか ) をチェックする。そうでない場合は例外を投げる。

その他に、STL ( Standard Template Library ) で提供されている関数 transforminner_product を多用しています。前者は、配列に定数を作用させたり、他の配列の要素どうしと演算させるために利用しています。また、後者は二つのベクトルの内積を計算するための関数です。特に transform の使い方は複数あって分かりづらいかもしれませんが、例えば

std::transform( m->begin(), m->end(), m->begin(), std::bind2nd( std::multiplies< T >(), alpha ) );

の場合、第一・第二引数が演算対象の要素範囲、第三引数が結果を代入する位置の先頭 ( この場合、演算対象と結果の格納位置が同じであることになります )、最後の引数が演算の内容になります。multiplies は二つの引数どうしを乗算する単純な関数ですが、この関数の第二引数を alpha 固定とするためにバインダ生成用関数 bind2nd を利用しています。これにより、各要素に定数 alpha を乗算することができます。もう一つの使い方の例として

std::transform( aTax_b.begin(), aTax_b.end(), m->begin(), m->begin(), std::plus< T >() );

の場合、第一・第二引数が一つめの配列の要素の範囲、第三引数が二つめの配列の先頭 ( 末尾の位置は一つめの配列の範囲で決まるため省略されています )、第四引数が結果を代入する位置の先頭、最後の引数が演算の内容 ( plus は加算を行う関数 ) を表します。よって、二つの配列を加算する処理を行っていることになります。


共役勾配法は、一般の n 変量関数に対する極値計算にも応用することができます。多変量のテーラー・マクローリン展開から、n 変数の関数 f(x) は二次近似 fII(x) として以下の式で表すことができます。

fII(x) = f(x(k)) + ( ∇f(k), x - x(k) ) + (1/2)( x - x(k), H(k)( x - x(k) ) )

ここで、x は真の解を表し、x(k) はそれに対する近似解を意味します。∇f は勾配、H は二階導関数 ∂2f / ∂xi∂xj ( i, j = 1, 2, ... n ) を i 行 j 列の要素とする行列で「ヘッセ行列 (Hessian Matrix)」と呼ばれます。また、添字の (k) は xx(k) を代入した値であることを表しています。勾配 ∇fII を計算すると、

∇fII = ∇f(k) + H(k)( x - x(k) ) = 0

となるので、先ほどと同様に x = x(k) + t(k)m(k) とおいて代入すると

∇f(k) = -t(k)H(k)m(k) --- (1)

となり、

H(k)m(k)∇f(k) --- (2)

が成り立ちます。同様に (3), (4) 式は

( u(k), H(k)m(k) ) = 0 --- (3)

m(k) = ∇f(k) + α(k)u(k) --- (4)

なので、

α(k) = -( u(k), H(k)∇f(k) ) / ( u(k), H(k)u(k) ) --- (5)

という結果が得られます。もし、f(x) が二次式なら、fII(x) と f(x) は等しいことになるので、(1) より

t(k) = -( m(k), ∇f(k) ) / ( m(k), H(k)m(k) )

として 4 つの反復公式で計算すればよいのですが、そうでない場合は t(k) を求めるために元の関数 f(x) を利用した方がよく、これには勾配法がそのまま適用できます。以上から、次のような反復処理を得ることができます。

  1. α(k) = -( m(k-1), H(k)∇f(k) ) / ( m(k-1), H(k)m(k-1) )
  2. m(k) = ∇f(k) + α(k)m(k-1)
  3. x(k) を初期値として、f(x) に対して m(k) の方向へ勾配法により線形探索を行い t(k) を得る
  4. x(k+1) = x(k) + t(k)m(k)

共役勾配法による極値計算用のサンプル・プログラムを以下に示します。

/**
   CGM_Operator : 共役勾配法 ( Conjugate Gradient Method ) 用多変数関数オブジェクト
**/
class CGM_Operator
{
  std::vector< double > x0_; // 初期値
  std::vector< double > m_;  // 探索方向ベクトル

  // fx : f(x) の値を求める
  virtual double fx( const std::vector< double >& x ) const = 0;

public:

  // 仮想デストラクタ(何もしない)
  virtual ~CGM_Operator() {}

  /*
    nabla : ∇f の値を求める

    x : x の値(多変量)
    i : x の何番目のパラメータに対する偏微分を返すかを指定する添字
  */
  virtual double nabla( const std::vector< double >& x, size_t i ) const = 0;

  /*
    hessian : ヘッセ行列 H の値を求める

    x : x の値(多変量)
    i,j : H の何行・何列目の二階偏微分を返すかを指定する添字
  */
  virtual double hessian( const std::vector< double >& x, size_t i, size_t j ) const = 0;

  // assign : 初期値 x0 と探索方向ベクトル m を代入する
  void assign( const std::vector< double >& x0, const std::vector< double >& m )
  {
    x0_.assign( x0.begin(), x0.end() );
    m_.assign( m.begin(), m.end() );
  }

  // ft : F(t) = f(x0 + t∇f(x0)) の値を求める
  double ft( double t ) const
  {
    std::vector< double > x( x0_ );
    for ( std::vector< double >::size_type i = 0 ; i < x0_.size() ; ++i )
      x[i] += t * m_[i];

    return( fx( x ) );
  }

  // operator() : dF(t)/dt = ( ∇f(x0 + t∇f0), ∇f(x0) ) を返す
  double operator()( double t ) const
  {
    std::vector< double > x( x0_ );
    for ( std::vector< double >::size_type i = 0 ; i < x0_.size() ; ++i )
      x[i] += t * m_[i];

    double ans = 0;
    for ( std::vector< double >::size_type i = 0 ; i < x0_.size() ; ++i )
      ans += m_[i] * nabla( x, i );

    return( ans );
  }
};

/*
  CGM_CalcAlpha : 共役勾配法による連立方程式計算 ( αを求める )

  α(k) = -( m(k-1), h(k)・nabla(k) ) / ( m(k-1), h(k)・m(k-1) )
*/
double CGM_CalcAlpha( const std::vector< double >& m, const SquareMatrix< double >& h, const std::vector< double >& nabla )
{
  // 分子の計算
  std::vector< double > hNabla( m.size() );
  MatrixLib::mult( h, nabla.begin(), nabla.end(), hNabla.begin(), hNabla.end() );
  double n = std::inner_product( m.begin(), m.end(), hNabla.begin(), double() );

  // 分母の計算
  std::vector< double > hm( m.size() );
  MatrixLib::mult( h, m.begin(), m.end(), hm.begin(), hm.end() );
  double d = std::inner_product( m.begin(), m.end(), hm.begin(), double() );

  return( ( d == 0 ) ? 0 : -n / d );
}

/*
  CGM_CalcM : 共役勾配法による連立方程式計算 ( mを求める )

  m(k) = nabla(k) + alpha(k)・m(k-1)
*/
void CGM_CalcM( const std::vector< double >& nabla, double alpha, std::vector< double >* m )
{
  // alpha・m(k-1) の計算
  std::transform( m->begin(), m->end(), m->begin(), std::bind2nd( std::multiplies< double >(), alpha ) );

  // 両者の和を求める
  std::transform( nabla.begin(), nabla.end(), m->begin(), m->begin(), std::plus< double >() );
}

/*
  ConjugateGradientMethod : 共役勾配法による多変量関数の極値計算

  テンプレート引数の Op は多変量関数オブジェクト CGM_Operator の派生型であることを想定

  op : 関数 f(x)
  x0 : 初期値へのポインタ (結果をそのまま返す)
  h : ステップ幅
  e : 収束判定値
*/
template< class Op >
  void ConjugateGradientMethod( Op& op, std::vector< double >* x0, double h, double e )
{
  double alpha = 0;                             // α
  SquareMatrix< double > hessian( x0->size() ); // H
  std::vector< double > nabla( x0->size() );    // ∇f
  std::vector< double > m( x0->size() );        // m

  for ( ; ; ) {
    // ∇f(k) の計算
    for ( std::vector< double >::size_type i = 0 ; i < x0->size() ; ++i )
      nabla[i] = op.nabla( *x0, i );
    // H(k) の計算
    for ( SquareMatrix< double >::size_type r = 0 ; r < hessian.rows() ; ++r )
      for ( SquareMatrix< double >::size_type c = 0 ; c < hessian.cols() ; ++c )
        hessian[r][c] = op.hessian( *x0, r, c );

    alpha = CGM_CalcAlpha( m, hessian, nabla ); // α(k) = -( m(k-1), H(k)・∇f(k) ) / ( m(k-1), H(k)・m(k-1) )
    CGM_CalcM( nabla, alpha, &m );              // m(k) = ∇f(k) + α(k)m(k-1)

    op.assign( *x0, m );
    h = ( ( op( 0.0 ) * op.ft( 0.0 ) < 0 ) ? 1 : -1 ) * std::abs( h );
    double t = GradientMethod( op, 0.0, h, e );

    if ( CGM_CalcX( t, m, x0, e ) ) break;      // x(k+1) = x(k) + t(k)・m(k)
  }
}

CGM_Operator は、山登り法のところでも登場した多変量関数用の基底クラスですが、違いとしてはヘッセ行列の要素を返す hessian という純粋仮想関数が一つ追加されたところです。ConjugateGradientMethod が極値計算のためのメイン関数となりますが、これも前述の連立方程式計算用サンプル・プログラムとほぼ同じ構成です。違いとしては、t の計算に勾配法を利用している部分です。


山登り法の時と同様に、CGM_Operator からの派生クラスである楕円体の関数オブジェクトを作成します。

/**
   CGM_Oval : 楕円体の関数オブジェクト
**/
class CGM_Oval : public CGM_Operator
{
  double a_, b_;   // 扁平率
  double x0_, y0_; // 中心の座標

  /*
    fx : f(x) = a( x - x0 )^2 + b( y - ya )^2 を求める

    ( x, y ) = ( x[0], x[1] ) として計算する
  */
  virtual double fx( const std::vector< double >& x ) const
  {
    assert( x.size() >= 2 );

    return( a_ * pow( x[0] - x0_, 2 ) + b_ * pow( x[1] - y0_, 2 ) );
  }

public:

  // 扁平率 a, b と中心座標 ( x0, y0 ) を指定して構築
  CGM_Oval( double a, double b, double x0, double y0 )
    : a_( a ), b_( b ), x0_( x0 ), y0_( y0 ) {}

  /*
    ∇f : ( ∂f/∂x, ∂f/∂y ) を求める

    ( x, y ) = ( x[0], x[1] ) として計算する
    i = 0 で ∂f/∂x、i = 1 で ∂f/∂y を返す
  */
  virtual double nabla( const std::vector< double >& x, size_t i ) const
  {
    assert( x.size() >= 2 );
    assert( i < 2 );

    if ( i == 0 )
      return( 2 * a_ * ( x[0] - x0_ ) );
    else
      return( 2 * b_ * ( x[1] - y0_ ) );
  }

  /*
    H : | ∂^2f/∂x^2,  ∂^2f/∂x∂y | を求める
        | ∂^2f/∂x∂y, ∂^2f/∂y^2  |

    ( x, y ) = ( x[0], x[1] ) として計算する
    ( i, j ) = ( 0, 0 ) で ∂^2f/∂x^2、( i, j ) = ( 1, 1 ) で ∂^2f/∂y^2
    ( i, j ) = ( 1, 0 ), ( 0, 1 ) で ∂^2f/∂x∂y を返す
  */
  virtual double hessian( const std::vector< double >& x, size_t i, size_t j ) const
  {
    assert( x.size() >= 2 );
    assert( i < 2 && j < 2 );

    if ( i == 0 && j == 0 )
      return( 2 * a_ );
    else if ( i == 1 && j == 1 )
      return( 2 * b_ );
    else
      return( 0 );
  }
};

山登り法と共役勾配法を使い、楕円体の極値を求めたときの収束の様子を以下に示します。楕円体の方程式は f( x, y ) = a( x - 5 )2 + b( y - 5 )2 で、a, b をいくつかの値に切り替えて計算しています。いずれの場合も初期値は ( x, y ) = ( 0, 0 ) で、収束判定のためのしきい値は 1E-10 としています。

表 3-1. 楕円体の極値計算の収束
a = 1, b = 1a = 2, b = 1a = 3, b = 1
山登り法共役勾配法山登り法共役勾配法山登り法共役勾配法
xyxyxyxyxyxy
000000000000
55555.5555555552.7777777785.5555555552.7777777785.3571428571.7857142865.3571428571.785714286
4.629629634.62962963554.4642857154.4642857164.9999999995
5.0411522634.8353909475.0382653064.655612246
4.9725651584.9725651584.9426020414.942602041
5.0030483164.9878067375.0040998544.963101312
4.997967794.9979677884.9938502194.993850219
5.00022584.9990967945.000439274.996046569
4.9998494694.9998494614.9993410964.999341094
5.0000167224.9999330925.0000470624.999576416
4.9999888964.9999887624.9999294194.999929368
5.0000010414.9999949085.0000050234.999954588
54.9999974544.9999925444.999992195
5.0000002334.999994878
4.9999998834.999997439

a = 1, b = 1 の場合は球面であり、∇fm は同じ方向となるため、どちらの場合も 2 回で収束していますが、a の値が増えるに従い楕円体は x の方向につぶれた形になるので、山登り法では ∇fm のズレが大きくなりなかなか収束しなくなります。それに対して共役勾配法は常に 2 回で収束します。


今回は極値計算手法について紹介しました。地味ですが、様々な場面で利用される非常に重要なテーマだと思います。文献なんかを読むと、極値計算としてよく紹介されているのが「最急降下法 ( Steepest Descent Method )」という「山登り法」とよく似た手法ですが、ステップ幅の値はたいてい定数となっていて、その更新方法などについては触れていないのがほとんどです。
今回紹介していない手法として「準ニュートン法 ( Quasi-Newton Method )」があります。「ニュートン - ラフソン法」がヘッセ行列の逆行列を必要とするのに対し、準ニュートン法では別の方法で置き換えることで計算の負荷を軽減しています。


補足 1) XTX が半正値対称行列であることの証明

まず、「固有値」や「対称行列」については「固有値問題 (1) 対称行列の固有値」にて、「半正値対称行列」については「固有値問題 (2) カルーネン・レーベ展開」にて詳しく紹介していますので、そちらを参照して下さい。

XTX が半正値対称行列であることを証明するために、まず以下の命題を証明します。

対称行列 A が半正値であるときに限って、任意のベクトル y に対して ( y, Ay ) ≥ 0 が成り立つ。

A は対称行列なので、固有値分解により QDQT の形に変換することができます。但し、Q は固有ベクトルからなる直交行列、D は固有値を対角成分とする対角行列です。これを ( y, Ay ) に代入すると

( y, Ay ) = ( y, QDQTy ) = ( QTy, DQTy )

となります。但し、( x, Ax ) = ( ATx, x ) であることを利用しています。QTyy' = ( y1, y2, ... yn )、i 番目の固有値を λi とすると、

( y, Ay ) = λ1y12 + λ2y22 + ... + λnyn2

なので、全ての固有値について λi ≥ 0 ならば任意の y について ( y, Ay ) ≥ 0 であり、またその逆も成り立ちます。よって、命題が証明されました。

ここまで証明できればあとは簡単です。XTX

( XTX )T = XTX

なので明らかに対称行列です。但し、( AB )T = BTAT であることを利用しています。また、任意のベクトル y について

( y, XTXy ) = ( Xy, Xy ) = ||Xy||2 ≥ 0

が成り立つので、先ほどの命題から XTX は半正値対称行列であることになります。


補足 2) 内積の微分

x, b を i 番目の要素をそれぞれ xi, bi とするベクトル、A を r 行 c 列目の要素を arc とする行列として、( x, b ), ( x, Ax ) の xi による偏微分を求めてみます。まず、( x, b ) は

( x, b ) = Σk( bkxk )

なのでその xi による偏微分は bi であることは簡単にわかります。次に、( x, Ax ) の場合、Ax の r 番目の成分 { Ax }r

{ Ax }r = Σc( arcxc )

なので、

( x, Ax ) = Σr( xr・{ Ax }r ) = Σr( Σc( arcxrxc ) )

となります。従って、xi による偏微分は

∂( x, Ax ) / ∂xi = Σr( arixr ) + Σc( aicxc )

と求めることができます。Σc( aicxc ) = { Ax }i、Σr( arixr ) = { ATx }i なので、

∂( x, Ax ) / ∂xi = { Ax }i + { ATx }i

と書くこともできます。∂( x, b ) / ∂xi を i 番目の要素とするベクトルを ∇( x, b )、∂( x, Ax ) / ∂xi を i 番目の要素とするベクトルを ∇( x, Ax ) とすれば、上記結果から

∇( x, b ) = b

∇( x, Ax ) = Ax + ATx

であり、特に A が対称行列ならば

∇( x, Ax ) = 2Ax

ということになります。


<参考文献>
  1. 「これなら分かる 最適化数学」 金谷健一 著 (共立出版)
  2. 秋田県立大学 シミュレーション工学 (http://www.pna.eis.akita-pu.ac.jp/) - 「数値解析 ( Numerical Analysis )」 - 「数値解析 ( 第二章 ) 非線形方程式を解く
  3. Wikipedia

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