granicus.if.org Git - liblinear/blob - newton.cpp

   1 #include <math.h>
   2 #include <stdio.h>
   3 #include <string.h>
   4 #include <stdarg.h>
   5 #include "newton.h"
   6
   7 #ifndef min
   8 template <class T> static inline T min(T x,T y) { return (x<y)?x:y; }
   9 #endif
  10
  11 #ifndef max
  12 template <class T> static inline T max(T x,T y) { return (x>y)?x:y; }
  13 #endif
  14
  15 #ifdef __cplusplus
  16 extern "C" {
  17 #endif
  18
  19 extern double dnrm2_(int *, double *, int *);
  20 extern double ddot_(int *, double *, int *, double *, int *);
  21 extern int daxpy_(int *, double *, double *, int *, double *, int *);
  22 extern int dscal_(int *, double *, double *, int *);
  23
  24 #ifdef __cplusplus
  25 }
  26 #endif
  27
  28 static void default_print(const char *buf)
  29 {
  30         fputs(buf,stdout);
  31         fflush(stdout);
  32 }
  33
  34 // On entry *f must be the function value of w
  35 // On exit w is updated and *f is the new function value
  36 double function::linesearch_and_update(double *w, double *s, double *f, double *g, double alpha)
  37 {
  38         double gTs = 0;
  39         double eta = 0.01;
  40         int n = get_nr_variable();
  41         int max_num_linesearch = 20;
  42         double *w_new = new double[n];
  43         double fold = *f;
  44
  45         for (int i=0;i<n;i++)
  46                 gTs += s[i] * g[i];
  47
  48         int num_linesearch = 0;
  49         for(num_linesearch=0; num_linesearch < max_num_linesearch; num_linesearch++)
  50         {
  51                 for (int i=0;i<n;i++)
  52                         w_new[i] = w[i] + alpha*s[i];
  53                 *f = fun(w_new);
  54                 if (*f - fold <= eta * alpha * gTs)
  55                         break;
  56                 else
  57                         alpha *= 0.5;
  58         }
  59
  60         if (num_linesearch >= max_num_linesearch)
  61         {
  62                 *f = fold;
  63                 return 0;
  64         }
  65         else
  66                 memcpy(w, w_new, sizeof(double)*n);
  67
  68         delete [] w_new;
  69         return alpha;
  70 }
  71
  72 void NEWTON::info(const char *fmt,...)
  73 {
  74         char buf[BUFSIZ];
  75         va_list ap;
  76         va_start(ap,fmt);
  77         vsprintf(buf,fmt,ap);
  78         va_end(ap);
  79         (*newton_print_string)(buf);
  80 }
  81
  82 NEWTON::NEWTON(const function *fun_obj, double eps, double eps_cg, int max_iter)
  83 {
  84         this->fun_obj=const_cast<function *>(fun_obj);
  85         this->eps=eps;
  86         this->eps_cg=eps_cg;
  87         this->max_iter=max_iter;
  88         newton_print_string = default_print;
  89 }
  90
  91 NEWTON::~NEWTON()
  92 {
  93 }
  94
  95 void NEWTON::newton(double *w)
  96 {
  97         int n = fun_obj->get_nr_variable();
  98         int i, cg_iter;
  99         double step_size;
 100         double f, fold, actred;
 101         double init_step_size = 1;
 102         int search = 1, iter = 1, inc = 1;
 103         double *s = new double[n];
 104         double *r = new double[n];
 105         double *g = new double[n];
 106
 107         const double alpha_pcg = 0.01;
 108         double *M = new double[n];
 109
 110         // calculate gradient norm at w=0 for stopping condition.
 111         double *w0 = new double[n];
 112         for (i=0; i<n; i++)
 113                 w0[i] = 0;
 114         fun_obj->fun(w0);
 115         fun_obj->grad(w0, g);
 116         double gnorm0 = dnrm2_(&n, g, &inc);
 117         delete [] w0;
 118
 119         f = fun_obj->fun(w);
 120         fun_obj->grad(w, g);
 121         double gnorm = dnrm2_(&n, g, &inc);
 122         info("init f %5.3e |g| %5.3e\n", f, gnorm);
 123
 124         if (gnorm <= eps*gnorm0)
 125                 search = 0;
 126
 127         while (iter <= max_iter && search)
 128         {
 129                 fun_obj->get_diag_preconditioner(M);
 130                 for(i=0; i<n; i++)
 131                         M[i] = (1-alpha_pcg) + alpha_pcg*M[i];
 132                 cg_iter = pcg(g, M, s, r);
 133
 134                 fold = f;
 135                 step_size = fun_obj->linesearch_and_update(w, s, &f, g, init_step_size);
 136
 137                 if (step_size == 0)
 138                 {
 139                         info("WARNING: line search fails\n");
 140                         break;
 141                 }
 142
 143                 fun_obj->grad(w, g);
 144                 gnorm = dnrm2_(&n, g, &inc);
 145
 146                 info("iter %2d f %5.3e |g| %5.3e CG %3d step_size %4.2e \n", iter, f, gnorm, cg_iter, step_size);
 147
 148                 if (gnorm <= eps*gnorm0)
 149                         break;
 150                 if (f < -1.0e+32)
 151                 {
 152                         info("WARNING: f < -1.0e+32\n");
 153                         break;
 154                 }
 155                 actred = fold - f;
 156                 if (fabs(actred) <= 1.0e-12*fabs(f))
 157                 {
 158                         info("WARNING: actred too small\n");
 159                         break;
 160                 }
 161
 162                 iter++;
 163         }
 164
 165         if(iter >= max_iter)
 166                 info("\nWARNING: reaching max number of Newton iterations\n");
 167
 168         delete[] g;
 169         delete[] r;
 170         delete[] s;
 171         delete[] M;
 172 }
 173
 174 int NEWTON::pcg(double *g, double *M, double *s, double *r)
 175 {
 176         int i, inc = 1;
 177         int n = fun_obj->get_nr_variable();
 178         double one = 1;
 179         double *d = new double[n];
 180         double *Hd = new double[n];
 181         double zTr, znewTrnew, alpha, beta, cgtol, dHd;
 182         double *z = new double[n];
 183         double Q = 0, newQ, Qdiff;
 184
 185         for (i=0; i<n; i++)
 186         {
 187                 s[i] = 0;
 188                 r[i] = -g[i];
 189                 z[i] = r[i] / M[i];
 190                 d[i] = z[i];
 191         }
 192
 193         zTr = ddot_(&n, z, &inc, r, &inc);
 194         double gMinv_norm = sqrt(zTr);
 195         cgtol = min(eps_cg, sqrt(gMinv_norm));
 196         int cg_iter = 0;
 197         int max_cg_iter = max(n, 5);
 198
 199         while (cg_iter < max_cg_iter)
 200         {
 201                 cg_iter++;
 202
 203                 fun_obj->Hv(d, Hd);
 204                 dHd = ddot_(&n, d, &inc, Hd, &inc);
 205                 // avoid 0/0 in getting alpha
 206                 if (dHd <= 1.0e-16)
 207                         break;
 208
 209                 alpha = zTr/dHd;
 210                 daxpy_(&n, &alpha, d, &inc, s, &inc);
 211                 alpha = -alpha;
 212                 daxpy_(&n, &alpha, Hd, &inc, r, &inc);
 213
 214                 // Using quadratic approximation as CG stopping criterion
 215                 newQ = -0.5*(ddot_(&n, s, &inc, r, &inc) - ddot_(&n, s, &inc, g, &inc));
 216                 Qdiff = newQ - Q;
 217                 if (newQ <= 0 && Qdiff <= 0)
 218                 {
 219                         if (cg_iter * Qdiff >= cgtol * newQ)
 220                                 break;
 221                 }
 222                 else
 223                 {
 224                         info("WARNING: quadratic approximation > 0 or increasing in CG\n");
 225                         break;
 226                 }
 227                 Q = newQ;
 228
 229                 for (i=0; i<n; i++)
 230                         z[i] = r[i] / M[i];
 231                 znewTrnew = ddot_(&n, z, &inc, r, &inc);
 232                 beta = znewTrnew/zTr;
 233                 dscal_(&n, &beta, d, &inc);
 234                 daxpy_(&n, &one, z, &inc, d, &inc);
 235                 zTr = znewTrnew;
 236         }
 237
 238         if (cg_iter == max_cg_iter)
 239                 info("WARNING: reaching maximal number of CG steps\n");
 240
 241         delete[] d;
 242         delete[] Hd;
 243         delete[] z;
 244
 245         return cg_iter;
 246 }
 247
 248 void NEWTON::set_print_string(void (*print_string) (const char *buf))
 249 {
 250         newton_print_string = print_string;
 251 }