Merge branch 'const_reorder_no_state' into const_reorder_2

1 files changed, 193 insertions, 0 deletions

diff --git a/utils/sgd.cpp b/utils/sgd.cpp
new file mode 100644
index 00000000..8613edca
--- /dev/null
+++ b/utils/sgd.cpp
@@ -0,0 +1,193 @@
+#include "maxent.h"
+#include <cmath>
+#include <stdio.h>
+
+using namespace std;
+
+// const double SGD_ETA0 = 1;
+// const double SGD_ITER = 30;
+// const double SGD_ALPHA = 0.85;
+
+//#define FOLOS_NAIVE
+//#define FOLOS_LAZY
+#define SGD_CP
+
+inline void apply_l1_penalty(const int i, const double u, vector<double>& _vl,
+                             vector<double>& q) {
+  double& w = _vl[i];
+  const double z = w;
+  double& qi = q[i];
+  if (w > 0) {
+    w = max(0.0, w - (u + qi));
+  } else if (w < 0) {
+    w = min(0.0, w + (u - qi));
+  }
+  qi += w - z;
+}
+
+static double l1norm(const vector<double>& v) {
+  double sum = 0;
+  for (size_t i = 0; i < v.size(); i++) sum += abs(v[i]);
+  return sum;
+}
+
+inline void update_folos_lazy(const int iter_sample, const int k,
+                              vector<double>& _vl,
+                              const vector<double>& sum_eta,
+                              vector<int>& last_updated) {
+  const double penalty = sum_eta[iter_sample] - sum_eta[last_updated[k]];
+  double& x = _vl[k];
+  if (x > 0)
+    x = max(0.0, x - penalty);
+  else
+    x = min(0.0, x + penalty);
+  last_updated[k] = iter_sample;
+}
+
+int ME_Model::perform_SGD() {
+  if (_l2reg > 0) {
+    cerr << "error: L2 regularization is currently not supported in SGD mode."
+         << endl;
+    exit(1);
+  }
+
+  cerr << "performing SGD" << endl;
+
+  const double l1param = _l1reg;
+
+  const int d = _fb.Size();
+
+  vector<int> ri(_vs.size());
+  for (size_t i = 0; i < ri.size(); i++) ri[i] = i;
+
+  vector<double> grad(d);
+  int iter_sample = 0;
+  const double eta0 = SGD_ETA0;
+
+  //  cerr << "l1param = " << l1param << endl;
+  cerr << "eta0 = " << eta0 << " alpha = " << SGD_ALPHA << endl;
+
+  double u = 0;
+  vector<double> q(d, 0);
+  vector<int> last_updated(d, 0);
+  vector<double> sum_eta;
+  sum_eta.push_back(0);
+
+  for (int iter = 0; iter < SGD_ITER; iter++) {
+
+    random_shuffle(ri.begin(), ri.end());
+
+    double logl = 0;
+    int ncorrect = 0, ntotal = 0;
+    for (size_t i = 0; i < _vs.size(); i++, ntotal++, iter_sample++) {
+      const Sample& s = _vs[ri[i]];
+
+#ifdef FOLOS_LAZY
+      for (vector<int>::const_iterator j = s.positive_features.begin();
+           j != s.positive_features.end(); j++) {
+        for (vector<int>::const_iterator k = _feature2mef[*j].begin();
+             k != _feature2mef[*j].end(); k++) {
+          update_folos_lazy(iter_sample, *k, _vl, sum_eta, last_updated);
+        }
+      }
+#endif
+
+      vector<double> membp(_num_classes);
+      const int max_label = conditional_probability(s, membp);
+
+      const double eta =
+          eta0 * pow(SGD_ALPHA,
+                     (double)iter_sample / _vs.size());  // exponential decay
+      //      const double eta = eta0 / (1.0 + (double)iter_sample /
+      // _vs.size());
+
+      //      if (iter_sample % _vs.size() == 0) cerr << "eta = " << eta <<
+      // endl;
+      u += eta * l1param;
+
+      sum_eta.push_back(sum_eta.back() + eta * l1param);
+
+      logl += log(membp[s.label]);
+      if (max_label == s.label) ncorrect++;
+
+      // binary features
+      for (vector<int>::const_iterator j = s.positive_features.begin();
+           j != s.positive_features.end(); j++) {
+        for (vector<int>::const_iterator k = _feature2mef[*j].begin();
+             k != _feature2mef[*j].end(); k++) {
+          const double me = membp[_fb.Feature(*k).label()];
+          const double ee = (_fb.Feature(*k).label() == s.label ? 1.0 : 0);
+          const double grad = (me - ee);
+          _vl[*k] -= eta * grad;
+#ifdef SGD_CP
+          apply_l1_penalty(*k, u, _vl, q);
+#endif
+        }
+      }
+      // real-valued features
+      for (vector<pair<int, double> >::const_iterator j = s.rvfeatures.begin();
+           j != s.rvfeatures.end(); j++) {
+        for (vector<int>::const_iterator k = _feature2mef[j->first].begin();
+             k != _feature2mef[j->first].end(); k++) {
+          const double me = membp[_fb.Feature(*k).label()];
+          const double ee = (_fb.Feature(*k).label() == s.label ? 1.0 : 0);
+          const double grad = (me - ee) * j->second;
+          _vl[*k] -= eta * grad;
+#ifdef SGD_CP
+          apply_l1_penalty(*k, u, _vl, q);
+#endif
+        }
+      }
+
+#ifdef FOLOS_NAIVE
+      for (size_t j = 0; j < d; j++) {
+        double& x = _vl[j];
+        if (x > 0)
+          x = max(0.0, x - eta * l1param);
+        else
+          x = min(0.0, x + eta * l1param);
+      }
+#endif
+    }
+    logl /= _vs.size();
+//    fprintf(stderr, "%4d logl = %8.3f acc = %6.4f ", iter, logl,
+// (double)ncorrect / ntotal);
+
+#ifdef FOLOS_LAZY
+    if (l1param > 0) {
+      for (size_t j = 0; j < d; j++)
+        update_folos_lazy(iter_sample, j, _vl, sum_eta, last_updated);
+    }
+#endif
+
+    double f = logl;
+    if (l1param > 0) {
+      const double l1 =
+          l1norm(_vl);  // this is not accurate when lazy update is used
+      //      cerr << "f0 = " <<  update_model_expectation() - l1param * l1 << "
+      // ";
+      f -= l1param * l1;
+      int nonzero = 0;
+      for (int j = 0; j < d; j++)
+        if (_vl[j] != 0) nonzero++;
+      //      cerr << " f = " << f << " l1 = " << l1 << " nonzero_features = "
+      // << nonzero << endl;
+    }
+    //    fprintf(stderr, "%4d  obj = %7.3f acc = %6.4f", iter+1, f,
+    // (double)ncorrect/ntotal);
+    //    fprintf(stderr, "%4d  obj = %f", iter+1, f);
+    fprintf(stderr, "%3d  obj(err) = %f (%6.4f)", iter + 1, f,
+            1 - (double)ncorrect / ntotal);
+
+    if (_nheldout > 0) {
+      double heldout_logl = heldout_likelihood();
+      //      fprintf(stderr, "  heldout_logl = %f  acc = %6.4f\n",
+      // heldout_logl, 1 - _heldout_error);
+      fprintf(stderr, "  heldout_logl(err) = %f (%6.4f)", heldout_logl,
+              _heldout_error);
+    }
+    fprintf(stderr, "\n");
+  }
+
+  return 0;
+}