無制約非線型最適化 (準ニュートン法)

説明

準ニュートン法は，二回連続微分可能な R^n 上の凸関数 f(x) の R^n 上の最小値を求めるアルゴリズムである（凸性は，大域収束性の保証と，ヘッセ行列が正定値になることに要請される）．

f(x) が二回連続微分可能と仮定すると，f の Taylor 展開は高次の項を無視して次のようになる：

f(x+dx) 〜 f(x) + ∇f(x)・dx + 1/2 dx^T ∇^2 f(x)・dx・dx

ここに停留条件 f(x+dx) 〜f(x) を代入して整理すれば

dx 〜 - [ 1/2 ∇^2 f(x) ]^{-1} ∇f(x)

となり，これで逐次更新してするというアルゴリズムが考えられる．これをニュートン法という．ニュートン法はヘッセ行列の逆行列を毎ステップごとに陽に計算する必要があるため，ヘッセ行列を求めにくい場合に適用しづらく，計算量的にもうれしくない．そこで，ここもヘッセ行列の逆行列に収束するような行列を逐次更新するすることが考えられる．

逆ヘッセ行列の逐次更新はいくつか知られているが，現在広く用いられているものとして Broyden-Fletcher-Goldfarb-Shanno 法がある．これは，逆ヘッセ行列を次で更新する：

dB = (y y^T)/(y^T s) - (B s)(B s)^T/(s^T B s)

さらに，ヘッセ行列の不正確を補うため，x の更新則も緩和して

dx = - t B ∇f(x)

とスカラー倍 t を含める．t は更新方向 v = B ∇f(x) を固定したとき，f(x + t v) を最小化するような値で定める．これは黄金探索法で求めることができる（厳密に最小化しなくとも，それなりに小さくしていれば十分動作する）．

計算量

O(M n^2) 回の関数 f と df の評価．ただし M は反復回数．

ソースコード

array quasi_newton(array x) {
  const int n = x.size();
  matrix B(n, array(n));
  for (int i = 0; i < n; ++i)
    B[i][i] = 1;
  double fx =  f(x), fz;
  array   g = df(x);
  do {
    fz = fx;
    array s(n);
    for (int i = 0; i < n; ++i)
      for (int j = 0; j < n; ++j)
        s[i] -= B[i][j] * g[j];

    array p = x; // goldsection-minimum
    const double r = 2 / (3 + sqrt(5));
    double a = 0, b = 1, t = r * (b - a), c = a + t, d = b - t;
    double fc, fd;
    for (int i = 0; i < n; ++i)
      x[i] = p[i] + c * s[i];
    fc = f(x);
    for (int i = 0; i < n; ++i)
      x[i] = p[i] + d * s[i];
    fd = f(x);
    while (d - c > EPS) {
      if (fc > fd) {
        a = c; c = d; d = b - r * (b - a);
        for (int i = 0; i < n; ++i)
          x[i] = p[i] + d * s[i];
        fc = fd; fd = f(x);
      } else {
        b = d; d = c; c = a + r * (b - a);
        for (int i = 0; i < n; ++i)
          x[i] = p[i] + c * s[i];
        fd = fc; fc = f(x);
      }
    }
    fx = f(x);
    array y = g; g = df(x);
    for (int i = 0; i < n; ++i)
      y[i] = g[i] - y[i];
    array By(n);
    for (int i = 0; i < n; ++i)
      for (int j = 0; j < n; ++j)
        By[i] += B[i][j] * y[j];
    double sy = inner_product(s, y) + EPS;
    double yBy = inner_product(y, By) + EPS;
    array u(n);
    for (int i = 0; i < n; ++i)
      u[i] = s[i]/sy - By[i]/yBy;
    for (int i = 0; i < n; ++i)
      for (int j = 0; j < n; ++j)
        B[i][j] += s[i]*s[j]/sy - By[i]*By[j]/yBy + u[i]*u[j]*yBy;
  } while ( abs(fx - fz) > EPS );
  return x;
}

Verified