Limited Memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS)

lbfgs.cc

#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 lbfgs(F f, G df, int n) {
  int m = min(10, n);

  vector<double> x(n);
  deque< vector<double> > y, s;

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

    // compute direction by lbfgs
    vector<double> q = g;
    double alpha[m];
    for (int i = 0; i < s.size(); ++i) {
      alpha[i] = dot(s[i], q) / dot(y[i], s[i]);
      q = q - alpha[i] * y[i];
    }
    for (int i = s.size()-1; i >= 0; --i) {
      double beta = dot(y[i], q) / dot(y[i], s[i]);
      q = q + (alpha[i] - beta) * s[i];
    }

    // optimize step size by golden section search
    vector<double> px = x;
    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*q), fc = f(x - c*q);
    while (d -a > eps*(c-d)) {
      if (fb < fc) {
        d = c; c = b; fc = fb;
        b = a + (d-c); fb = f(x-b*q);
      } else {
        a = b; b = c; fb = fc;
        c = d - (b-a); fc = f(x-c*q);
      }
    }
    x = x - c * q;

    // update parameters
    vector<double> pg = g;
    g = df(x);
    s.push_front(x - px);
    y.push_front(g - pg);
    if (s.size() > m) {
      y.pop_back();
      s.pop_back();
    }
  }
}


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