SHOGUN  v1.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Tron.cpp
Go to the documentation of this file.
1 #include <math.h>
2 #include <stdio.h>
3 #include <string.h>
4 #include <stdarg.h>
5 
6 #include <shogun/lib/config.h>
7 #include <shogun/lib/Signal.h>
8 #include <shogun/lib/Time.h>
9 
10 #ifdef HAVE_LAPACK
13 
14 using namespace shogun;
15 
16 CTron::CTron(const function *f, float64_t e, int32_t it)
17 : CSGObject()
18 {
19  this->fun_obj=const_cast<function *>(f);
20  this->eps=e;
21  this->max_iter=it;
22 }
23 
25 {
26 }
27 
28 void CTron::tron(float64_t *w, float64_t max_train_time)
29 {
30  // Parameters for updating the iterates.
31  float64_t eta0 = 1e-4, eta1 = 0.25, eta2 = 0.75;
32 
33  // Parameters for updating the trust region size delta.
34  float64_t sigma1 = 0.25, sigma2 = 0.5, sigma3 = 4.;
35 
36  int32_t i, cg_iter;
37  float64_t delta, snorm, one=1.0;
38  float64_t alpha, f, fnew, prered, actred, gs;
39 
40  /* calling external lib */
41  int n = (int) fun_obj->get_nr_variable();
42  int search = 1, iter = 1, inc = 1;
43  double *s = SG_MALLOC(double, n);
44  double *r = SG_MALLOC(double, n);
45  double *w_new = SG_MALLOC(double, n);
46  double *g = SG_MALLOC(double, n);
47 
48  for (i=0; i<n; i++)
49  w[i] = 0;
50 
51  f = fun_obj->fun(w);
52  fun_obj->grad(w, g);
53  delta = cblas_dnrm2(n, g, inc);
54  float64_t gnorm1 = delta;
55  float64_t gnorm = gnorm1;
56 
57  if (gnorm <= eps*gnorm1)
58  search = 0;
59 
60  iter = 1;
61 
62  CSignal::clear_cancel();
63  CTime start_time;
64 
65  while (iter <= max_iter && search && (!CSignal::cancel_computations()))
66  {
67  if (max_train_time > 0 && start_time.cur_time_diff() > max_train_time)
68  break;
69 
70  cg_iter = trcg(delta, g, s, r);
71 
72  memcpy(w_new, w, sizeof(float64_t)*n);
73  cblas_daxpy(n, one, s, inc, w_new, inc);
74 
75  gs = cblas_ddot(n, g, inc, s, inc);
76  prered = -0.5*(gs-cblas_ddot(n, s, inc, r, inc));
77  fnew = fun_obj->fun(w_new);
78 
79  // Compute the actual reduction.
80  actred = f - fnew;
81 
82  // On the first iteration, adjust the initial step bound.
83  snorm = cblas_dnrm2(n, s, inc);
84  if (iter == 1)
85  delta = CMath::min(delta, snorm);
86 
87  // Compute prediction alpha*snorm of the step.
88  if (fnew - f - gs <= 0)
89  alpha = sigma3;
90  else
91  alpha = CMath::max(sigma1, -0.5*(gs/(fnew - f - gs)));
92 
93  // Update the trust region bound according to the ratio of actual to predicted reduction.
94  if (actred < eta0*prered)
95  delta = CMath::min(CMath::max(alpha, sigma1)*snorm, sigma2*delta);
96  else if (actred < eta1*prered)
97  delta = CMath::max(sigma1*delta, CMath::min(alpha*snorm, sigma2*delta));
98  else if (actred < eta2*prered)
99  delta = CMath::max(sigma1*delta, CMath::min(alpha*snorm, sigma3*delta));
100  else
101  delta = CMath::max(delta, CMath::min(alpha*snorm, sigma3*delta));
102 
103  SG_INFO("iter %2d act %5.3e pre %5.3e delta %5.3e f %5.3e |g| %5.3e CG %3d\n", iter, actred, prered, delta, f, gnorm, cg_iter);
104 
105  if (actred > eta0*prered)
106  {
107  iter++;
108  memcpy(w, w_new, sizeof(float64_t)*n);
109  f = fnew;
110  fun_obj->grad(w, g);
111 
112  gnorm = cblas_dnrm2(n, g, inc);
113  if (gnorm < eps*gnorm1)
114  break;
115  SG_SABS_PROGRESS(gnorm, -CMath::log10(gnorm), -CMath::log10(1), -CMath::log10(eps*gnorm1), 6);
116  }
117  if (f < -1.0e+32)
118  {
119  SG_WARNING("f < -1.0e+32\n");
120  break;
121  }
122  if (CMath::abs(actred) <= 0 && CMath::abs(prered) <= 0)
123  {
124  SG_WARNING("actred and prered <= 0\n");
125  break;
126  }
127  if (CMath::abs(actred) <= 1.0e-12*CMath::abs(f) &&
128  CMath::abs(prered) <= 1.0e-12*CMath::abs(f))
129  {
130  SG_WARNING("actred and prered too small\n");
131  break;
132  }
133  }
134 
135  SG_DONE();
136 
137  SG_FREE(g);
138  SG_FREE(r);
139  SG_FREE(w_new);
140  SG_FREE(s);
141 }
142 
143 int32_t CTron::trcg(float64_t delta, double* g, double* s, double* r)
144 {
145  /* calling external lib */
146  int i, cg_iter;
147  int n = (int) fun_obj->get_nr_variable();
148  int inc = 1;
149  double one = 1;
150  double *Hd = SG_MALLOC(double, n);
151  double *d = SG_MALLOC(double, n);
152  double rTr, rnewTrnew, alpha, beta, cgtol;
153 
154  for (i=0; i<n; i++)
155  {
156  s[i] = 0;
157  r[i] = -g[i];
158  d[i] = r[i];
159  }
160  cgtol = 0.1* cblas_dnrm2(n, g, inc);
161 
162  cg_iter = 0;
163  rTr = cblas_ddot(n, r, inc, r, inc);
164  while (1)
165  {
166  if (cblas_dnrm2(n, r, inc) <= cgtol)
167  break;
168  cg_iter++;
169  fun_obj->Hv(d, Hd);
170 
171  alpha = rTr/cblas_ddot(n, d, inc, Hd, inc);
172  cblas_daxpy(n, alpha, d, inc, s, inc);
173  if (cblas_dnrm2(n, s, inc) > delta)
174  {
175  SG_INFO("cg reaches trust region boundary\n");
176  alpha = -alpha;
177  cblas_daxpy(n, alpha, d, inc, s, inc);
178 
179  double std = cblas_ddot(n, s, inc, d, inc);
180  double sts = cblas_ddot(n, s, inc, s, inc);
181  double dtd = cblas_ddot(n, d, inc, d, inc);
182  double dsq = delta*delta;
183  double rad = sqrt(std*std + dtd*(dsq-sts));
184  if (std >= 0)
185  alpha = (dsq - sts)/(std + rad);
186  else
187  alpha = (rad - std)/dtd;
188  cblas_daxpy(n, alpha, d, inc, s, inc);
189  alpha = -alpha;
190  cblas_daxpy(n, alpha, Hd, inc, r, inc);
191  break;
192  }
193  alpha = -alpha;
194  cblas_daxpy(n, alpha, Hd, inc, r, inc);
195  rnewTrnew = cblas_ddot(n, r, inc, r, inc);
196  beta = rnewTrnew/rTr;
197  cblas_dscal(n, beta, d, inc);
198  cblas_daxpy(n, one, r, inc, d, inc);
199  rTr = rnewTrnew;
200  }
201 
202  SG_FREE(d);
203  SG_FREE(Hd);
204 
205  return(cg_iter);
206 }
207 
208 float64_t CTron::norm_inf(int32_t n, float64_t *x)
209 {
210  float64_t dmax = CMath::abs(x[0]);
211  for (int32_t i=1; i<n; i++)
212  if (CMath::abs(x[i]) >= dmax)
213  dmax = CMath::abs(x[i]);
214  return(dmax);
215 }
216 #endif //HAVE_LAPACK

SHOGUN Machine Learning Toolbox - Documentation