Newton-Krylov method for function minimization.

newtonkrylov.cc

// 
// Newton-Krylov method
//
// The Newton method is:
//   Step 1. solve H(x) u = g(x), where H is a Hessian and g is a gradient.
//   Step 2. update x = x - u.
// This is a second order method, which requires a Hessian information.
// On the Newton-Krylov method, Step 1 is implemented as
//   (a) solve this equation by iterative (Krylov subspace) method.
//   (b) use numerical difference to LHS evaluation, i.e.,
//     (g(x + h u) - g(x)) / h -> H(x) u, where h is a small number.
// This is a first order method, which requires only gradient information.
// 
// Complexity: Typically O(n log n) per iteration.
//
#include <iostream>
#include <vector>
#include <cstdio>
#include <cstdlib>
#include <map>
#include <cmath>
#include <cstring>
#include <functional>
#include <algorithm>

using namespace std;

#define ALL(c) c.begin(), c.end()
#define FOR(i,c) for(typeof(c.begin())i=c.begin();i!=c.end();++i)
#define REP(i,n) for(int i=0;i<n;++i)
#define REPD(i,n) for(int i=n;i>=0;--i)
#define REPS(i,s,n) for(int i=s;i<n;++i)
#define fst first
#define snd second


// Lv.1 BLAS operations
typedef vector<double> Vector;
Vector axpy(double a, Vector x, Vector y = Vector(0)) {
  REP(i, x.size()) x[i] *= a;
  REP(i, y.size()) x[i] += y[i];
  return x;
}
Vector operator+(Vector x, Vector y) { return axpy(1, x, y); }
Vector operator-(Vector x, Vector y) { return axpy(-1, y, x); }
Vector operator*(double a, Vector x) { return axpy(a, x); }
Vector operator/(Vector x, double a) { return axpy(1/a, x); }
double dot(Vector x, Vector y) {
  double d = 0;
  REP(i, x.size()) d += x[i] * y[i];
  return d;
}
double norm(Vector x) { return sqrt(dot(x,x)); }


// test function
double f(Vector x) { // f(x) = (x-1)^2 + 2 (y-2)^4
  return pow(x[0]-1, 2) + 2 * pow(x[1]-2, 4);
}
Vector g(Vector x) { // g(x) = grad f(x) = (2 (x-1), 8 (y-2)^3)
  x[0] = 2 * pow(x[0]-1, 1);
  x[1] = 8 * pow(x[1]-2, 3);
  return x;
}

// Newton-Krylov Optimization
//
// solve H(x) u == g(x) by iterative method
//   LHS evaluation: numerical difference H(x) u = (g(x + h u) - g(x)) / h
//   linear solver:  simpler GMRES
//
bool newtonDirection(Vector x, Vector &dx) {
  const double EPS = 1e-6, h = EPS / 100;

  vector<Vector> V, T;
  Vector gx = g(x), b = gx, r, c;
  double b0 = norm(b);
  if (b0 < EPS) return false;
  r = b / b0;

  REP(k, x.size()) {
    Vector &u = k == 0 ? r : V[k-1];
    V.push_back((g(x+h*u)-gx)/h); // LHS evaluation

    T.push_back( Vector(k+1) );
    REP(i, k) {
      T[k][i] = dot(V[i], V[k]);
      V[k] = V[k] - T[k][i] * V[i];
    }
    T[k][k] = norm(V[k]);
    V[k] = V[k] / T[k][k];
    c.push_back(dot(r, V[k]));
    r = r - c[k] * V[k];

    if (b0 * norm(r) < EPS || k == x.size()-1) {
      REPD(i, k) {
        REPS(j, i+1, k+1) c[i] -= T[j][i] * c[j];
        c[i] /= T[i][i];
      }
      dx = c[0] * b;
      REP(i, k) dx = dx + b0 * c[i+1] * V[i];
      return true;
    }
  }
}
double minimize() {
  int n = 2;
  Vector x(n);
  REP(iter, 100) {
    Vector dx;
    if (!newtonDirection(x, dx)) break;
    x = x - dx;
  }
  REP(i, n) cout << x[i] << " "; cout << endl;
}

int main() { minimize(); }