SHOGUN  6.1.3
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 #include <shogun/base/progress.h>
14 
15 using namespace shogun;
16 
17 double tron_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
18 {
19 #ifdef HAVE_LAPACK
20  return cblas_ddot(N,X,incX,Y,incY);
21 #else
22  double dot = 0.0;
23  for (int32_t i=0; i<N; i++)
24  dot += X[incX*i]*Y[incY*i];
25  return dot;
26 #endif
27 }
28 
29 double tron_dnrm2(const int N, const double *X, const int incX)
30 {
31 #ifdef HAVE_LAPACK
32  return cblas_dnrm2(N,X,incX);
33 #else
34  double dot = 0.0;
35  for (int32_t i=0; i<N; i++)
36  dot += X[incX*i]*X[incX*i];
37  return sqrt(dot);
38 #endif
39 }
40 
41 void tron_dscal(const int N, const double alpha, double *X, const int incX)
42 {
43 #ifdef HAVE_LAPACK
44  return cblas_dscal(N,alpha,X,incX);
45 #else
46  for (int32_t i=0; i<N; i++)
47  X[i]*= alpha;
48 #endif
49 }
50 
51 void tron_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
52 {
53 #ifdef HAVE_LAPACK
54  cblas_daxpy(N,alpha,X,incX,Y,incY);
55 #else
56  for (int32_t i=0; i<N; i++)
57  Y[i] += alpha*X[i];
58 #endif
59 }
60 
61 CTron::CTron(const function *f, float64_t e, int32_t it)
62 : CSGObject()
63 {
64  this->fun_obj=const_cast<function *>(f);
65  this->eps=e;
66  this->max_iter=it;
67 }
68 
70 {
71 }
72 
73 void CTron::tron(float64_t *w, float64_t max_train_time)
74 {
75  // Parameters for updating the iterates.
76  float64_t eta0 = 1e-4, eta1 = 0.25, eta2 = 0.75;
77 
78  // Parameters for updating the trust region size delta.
79  float64_t sigma1 = 0.25, sigma2 = 0.5, sigma3 = 4.;
80 
81  int32_t i, cg_iter;
82  float64_t delta, snorm, one=1.0;
83  float64_t alpha, f, fnew, prered, actred, gs;
84 
85  /* calling external lib */
86  int n = (int) fun_obj->get_nr_variable();
87  int search = 1, iter = 1, inc = 1;
88  double *s = SG_MALLOC(double, n);
89  double *r = SG_MALLOC(double, n);
90  double *w_new = SG_MALLOC(double, n);
91  double *g = SG_MALLOC(double, n);
92 
93  for (i=0; i<n; i++)
94  w[i] = 0;
95 
96  f = fun_obj->fun(w);
97  fun_obj->grad(w, g);
98  delta = tron_dnrm2(n, g, inc);
99  float64_t gnorm1 = delta;
100  float64_t gnorm = gnorm1;
101 
102  if (gnorm <= eps*gnorm1)
103  search = 0;
104 
105  iter = 1;
106 
107  CTime start_time;
108  auto pb = progress(range(10));
109 
110  // TODO: replace with new signal
111  // while (iter <= max_iter && search && (!CSignal::cancel_computations()))
112  while (iter <= max_iter && search)
113  {
114  if (max_train_time > 0 && start_time.cur_time_diff() > max_train_time)
115  break;
116 
117  cg_iter = trcg(delta, g, s, r);
118 
119  sg_memcpy(w_new, w, sizeof(float64_t)*n);
120  tron_daxpy(n, one, s, inc, w_new, inc);
121 
122  gs = tron_ddot(n, g, inc, s, inc);
123  prered = -0.5*(gs-tron_ddot(n, s, inc, r, inc));
124  fnew = fun_obj->fun(w_new);
125 
126  // Compute the actual reduction.
127  actred = f - fnew;
128 
129  // On the first iteration, adjust the initial step bound.
130  snorm = tron_dnrm2(n, s, inc);
131  if (iter == 1)
132  delta = CMath::min(delta, snorm);
133 
134  // Compute prediction alpha*snorm of the step.
135  if (fnew - f - gs <= 0)
136  alpha = sigma3;
137  else
138  alpha = CMath::max(sigma1, -0.5*(gs/(fnew - f - gs)));
139 
140  // Update the trust region bound according to the ratio of actual to predicted reduction.
141  if (actred < eta0*prered)
142  delta = CMath::min(CMath::max(alpha, sigma1)*snorm, sigma2*delta);
143  else if (actred < eta1*prered)
144  delta = CMath::max(sigma1*delta, CMath::min(alpha*snorm, sigma2*delta));
145  else if (actred < eta2*prered)
146  delta = CMath::max(sigma1*delta, CMath::min(alpha*snorm, sigma3*delta));
147  else
148  delta = CMath::max(delta, CMath::min(alpha*snorm, sigma3*delta));
149 
150  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)
151 
152  if (actred > eta0*prered)
153  {
154  iter++;
155  sg_memcpy(w, w_new, sizeof(float64_t)*n);
156  f = fnew;
157  fun_obj->grad(w, g);
158 
159  gnorm = tron_dnrm2(n, g, inc);
160  if (gnorm < eps*gnorm1)
161  break;
162  pb.print_absolute(
163  gnorm, -CMath::log10(gnorm), -CMath::log10(1),
164  -CMath::log10(eps * gnorm1));
165  }
166  if (f < -1.0e+32)
167  {
168  SG_WARNING("f < -1.0e+32\n")
169  break;
170  }
171  if (CMath::abs(actred) <= 0 && CMath::abs(prered) <= 0)
172  {
173  SG_WARNING("actred and prered <= 0\n")
174  break;
175  }
176  if (CMath::abs(actred) <= 1.0e-12*CMath::abs(f) &&
177  CMath::abs(prered) <= 1.0e-12*CMath::abs(f))
178  {
179  SG_WARNING("actred and prered too small\n")
180  break;
181  }
182  }
183 
184  pb.complete_absolute();
185 
186  SG_FREE(g);
187  SG_FREE(r);
188  SG_FREE(w_new);
189  SG_FREE(s);
190 }
191 
192 int32_t CTron::trcg(float64_t delta, double* g, double* s, double* r)
193 {
194  /* calling external lib */
195  int i, cg_iter;
196  int n = (int) fun_obj->get_nr_variable();
197  int inc = 1;
198  double one = 1;
199  double *Hd = SG_MALLOC(double, n);
200  double *d = SG_MALLOC(double, n);
201  double rTr, rnewTrnew, alpha, beta, cgtol;
202 
203  for (i=0; i<n; i++)
204  {
205  s[i] = 0;
206  r[i] = -g[i];
207  d[i] = r[i];
208  }
209  cgtol = 0.1* tron_dnrm2(n, g, inc);
210 
211  cg_iter = 0;
212  rTr = tron_ddot(n, r, inc, r, inc);
213  while (1)
214  {
215  if (tron_dnrm2(n, r, inc) <= cgtol)
216  break;
217  cg_iter++;
218  fun_obj->Hv(d, Hd);
219 
220  alpha = rTr/tron_ddot(n, d, inc, Hd, inc);
221  tron_daxpy(n, alpha, d, inc, s, inc);
222  if (tron_dnrm2(n, s, inc) > delta)
223  {
224  SG_INFO("cg reaches trust region boundary\n")
225  alpha = -alpha;
226  tron_daxpy(n, alpha, d, inc, s, inc);
227 
228  double std = tron_ddot(n, s, inc, d, inc);
229  double sts = tron_ddot(n, s, inc, s, inc);
230  double dtd = tron_ddot(n, d, inc, d, inc);
231  double dsq = delta*delta;
232  double rad = sqrt(std*std + dtd*(dsq-sts));
233  if (std >= 0)
234  alpha = (dsq - sts)/(std + rad);
235  else
236  alpha = (rad - std)/dtd;
237  tron_daxpy(n, alpha, d, inc, s, inc);
238  alpha = -alpha;
239  tron_daxpy(n, alpha, Hd, inc, r, inc);
240  break;
241  }
242  alpha = -alpha;
243  tron_daxpy(n, alpha, Hd, inc, r, inc);
244  rnewTrnew = tron_ddot(n, r, inc, r, inc);
245  beta = rnewTrnew/rTr;
246  tron_dscal(n, beta, d, inc);
247  tron_daxpy(n, one, r, inc, d, inc);
248  rTr = rnewTrnew;
249  }
250 
251  SG_FREE(d);
252  SG_FREE(Hd);
253 
254  return(cg_iter);
255 }
256 
257 float64_t CTron::norm_inf(int32_t n, float64_t *x)
258 {
259  float64_t dmax = CMath::abs(x[0]);
260  for (int32_t i=1; i<n; i++)
261  if (CMath::abs(x[i]) >= dmax)
262  dmax = CMath::abs(x[i]);
263  return(dmax);
264 }
Class Time that implements a stopwatch based on either cpu time or wall clock time.
Definition: Time.h:42
#define SG_INFO(...)
Definition: SGIO.h:117
virtual ~CTron()
Definition: tron.cpp:69
PRange< T > progress(Range< T > range, const SGIO &io, std::string prefix="PROGRESS: ", SG_PRG_MODE mode=UTF8, std::function< bool()> condition=[](){return true;})
Definition: progress.h:712
void tron(float64_t *w, float64_t max_train_time)
Definition: tron.cpp:73
Definition: basetag.h:132
double tron_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: tron.cpp:17
T dot(const SGVector< T > &a, const SGVector< T > &b)
Range< T > range(T rend)
Definition: range.h:136
float64_t cur_time_diff(bool verbose=false)
Definition: Time.cpp:68
Class SGObject is the base class of all shogun objects.
Definition: SGObject.h:124
void tron_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
Definition: tron.cpp:51
void tron_dscal(const int N, const double alpha, double *X, const int incX)
Definition: tron.cpp:41
double float64_t
Definition: common.h:60
CTron()
Definition: tron.h:58
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
#define SG_WARNING(...)
Definition: SGIO.h:127
double tron_dnrm2(const int N, const double *X, const int incX)
Definition: tron.cpp:29
T max(const Container< T > &a)

SHOGUN Machine Learning Toolbox - Documentation