Nonlinear Conjugate Gradient Descent (Polak-Ribiere)

nlcg.cc

// nonlinear conjugate gradient

#include <iostream>
#include <vector>
#include <deque>
#include <cstdio>
#include <cstdlib>
#include <map>
#include <cmath>
#include <cstring>
#include <functional>
#include <algorithm>
#include <unordered_map>

using namespace std;

double dot(vector<double> a, vector<double> b) {
  double c = 0;
  for (int i = 0; i < a.size(); ++i)
    c += a[i] * b[i];
  return c;
}
vector<double> operator*(double a, vector<double> x) {
  for (int i = 0; i < x.size(); ++i) x[i] *= a;
  return x;
}
vector<double> operator+(vector<double> x, vector<double> y) {
  for (int i = 0; i < x.size(); ++i) x[i] += y[i];
  return x;
}
vector<double> operator-(vector<double> x, vector<double> y) {
  return x + ((-1) * y);
}


double f(vector<double> x) {
  double score = 0;
  for (int i = 0; i < x.size(); ++i)
    score += (i + 1) * (x[i] - i)*(x[i] - i) / 2;
  return score;
}
vector<double> df(vector<double> x) {
  vector<double> g(x.size());
  for (int i = 0; i < x.size(); ++i)
    g[i] = (i + 1) * (x[i] - i);
  return g;
}


template <class F, class G>
double nlcg(F f, G df, int n) {
  vector<double> x(n);

  vector<double> pg;
  for (int k = 0; k < 20; ++k) {
    cout << f(x) << endl;

    // conjugate gradient direction (Polak-Ribiere)
    vector<double> g = df(x);
    if (!pg.empty()) {
      double beta = max(0.0, (dot(g, g)-dot(g,pg))/(dot(pg, pg)));
      g = g + beta * pg;
    }
    pg = g;

    // optimize step size by golden section search
    double r = 2/(3+sqrt(5)), eps = 1e-5;
    double a = 0, d = 1, b = a + (d-a)*r, c = d - (b-a);
    double fb = f(x - b*g), fc = f(x - c*g);
    while (d -a > eps*(c-d)) {
      if (fb < fc) {
        d = c; c = b; fc = fb;
        b = a + (d-c); fb = f(x-b*g);
      } else {
        a = b; b = c; fb = fc;
        c = d - (b-a); fc = f(x-c*g);
      }
    }

    // update parameters
    x = x - c * g;
  }
}


int main() {
  int n = 50;
  nlcg(f, df, n);
}