SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
libp3bm.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * libppbm.h: Implementation of the Proximal Point BM solver for SO training
8  *
9  * Copyright (C) 2012 Michal Uricar, uricamic@cmp.felk.cvut.cz
10  *
11  * Implementation of the Proximal Point P-BMRM (p3bm)
12  *--------------------------------------------------------------------- */
13 
15 #include <shogun/lib/external/libqp.h>
16 #include <shogun/lib/Time.h>
17 
18 namespace shogun
19 {
20 static const uint32_t QPSolverMaxIter=0xFFFFFFFF;
21 static const float64_t epsilon=0.0;
22 
23 static float64_t *H, *H2;
24 
25 /*----------------------------------------------------------------------
26  Returns pointer at i-th column of Hessian matrix.
27  ----------------------------------------------------------------------*/
28 static const float64_t *get_col( uint32_t i)
29 {
30  return( &H2[ BufSize*i ] );
31 }
32 
34  CDualLibQPBMSOSVM *machine,
35  float64_t* W,
36  float64_t TolRel,
37  float64_t TolAbs,
38  float64_t _lambda,
39  uint32_t _BufSize,
40  bool cleanICP,
41  uint32_t cleanAfter,
42  float64_t K,
43  uint32_t Tmax,
44  uint32_t cp_models,
45  bool verbose)
46 {
47  BmrmStatistics p3bmrm;
48  libqp_state_T qp_exitflag={0, 0, 0, 0}, qp_exitflag_good={0, 0, 0, 0};
49  float64_t *b, *b2, *beta, *beta_good, *beta_start, *diag_H, *diag_H2;
50  float64_t R, *Rt, **subgrad_t, *A, QPSolverTolRel, *C=NULL;
51  float64_t *prevW, *wt, alpha, alpha_start, alpha_good=0.0, Fd_alpha0=0.0;
52  float64_t lastFp, wdist, gamma=0.0;
53  floatmax_t rsum, sq_norm_W, sq_norm_Wdiff, sq_norm_prevW, eps;
54  uint32_t *I, *I2, *I_start, *I_good;
55  uint8_t *S=NULL;
56  uint32_t qp_cnt=0;
57  bmrm_ll *CPList_head, *CPList_tail, *cp_ptr, *cp_ptr2, *cp_list=NULL;
58  float64_t *A_1=NULL;
59  bool *map=NULL, tuneAlpha=true, flag=true;
60  bool alphaChanged=false, isThereGoodSolution=false;
61  TMultipleCPinfo **info=NULL;
62  CStructuredModel* model=machine->get_model();
63  CSOSVMHelper* helper = NULL;
64  uint32_t nDim=model->get_dim();
65  uint32_t to=0, N=0, cp_i=0;
66 
67  CTime ttime;
68  float64_t tstart, tstop;
69 
70 
71  tstart=ttime.cur_time_diff(false);
72 
73  BufSize=_BufSize*cp_models;
74  QPSolverTolRel=1e-9;
75 
76  H=NULL;
77  b=NULL;
78  beta=NULL;
79  A=NULL;
80  subgrad_t=NULL;
81  diag_H=NULL;
82  I=NULL;
83  prevW=NULL;
84  wt=NULL;
85  diag_H2=NULL;
86  b2=NULL;
87  I2=NULL;
88  H2=NULL;
89  I_good=NULL;
90  I_start=NULL;
91  beta_start=NULL;
92  beta_good=NULL;
93 
94  alpha=0.0;
95 
97 
98  A= (float64_t*) LIBBMRM_CALLOC(nDim*BufSize, float64_t);
99 
100  b= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
101 
102  beta= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
103 
104  subgrad_t= (float64_t**) LIBBMRM_CALLOC(cp_models, float64_t*);
105 
106  Rt= (float64_t*) LIBBMRM_CALLOC(cp_models, float64_t);
107 
108  diag_H= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
109 
110  I= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
111 
112  cp_list= (bmrm_ll*) LIBBMRM_CALLOC(1, bmrm_ll);
113 
114  prevW= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
115 
116  wt= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
117 
118  C= (float64_t*) LIBBMRM_CALLOC(cp_models, float64_t);
119 
120  S= (uint8_t*) LIBBMRM_CALLOC(cp_models, uint8_t);
121 
122  info= (TMultipleCPinfo**) LIBBMRM_CALLOC(cp_models, TMultipleCPinfo*);
123 
124  CFeatures* features = model->get_features();
125  int32_t num_feats = features->get_num_vectors();
126  SG_UNREF(features);
127 
128  /* CP cleanup variables */
129  ICP_stats icp_stats;
130  icp_stats.maxCPs = BufSize;
131  icp_stats.ICPcounter= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
132  icp_stats.ICPs= (float64_t**) LIBBMRM_CALLOC(BufSize, float64_t*);
133  icp_stats.ACPs= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
134  icp_stats.H_buff= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t);
135 
136  if (H==NULL || A==NULL || b==NULL || beta==NULL || subgrad_t==NULL ||
137  diag_H==NULL || I==NULL || icp_stats.ICPcounter==NULL ||
138  icp_stats.ICPs==NULL || icp_stats.ACPs==NULL ||
139  cp_list==NULL || prevW==NULL || wt==NULL || Rt==NULL || C==NULL ||
140  S==NULL || info==NULL || icp_stats.H_buff==NULL)
141  {
142  p3bmrm.exitflag=-2;
143  goto cleanup;
144  }
145 
146  /* multiple cutting plane model init */
147 
148  to=0;
149  N= (uint32_t) round( (float64_t) ((float64_t)num_feats / (float64_t) cp_models));
150 
151  for (uint32_t p=0; p<cp_models; ++p)
152  {
153  S[p]=1;
154  C[p]=1.0;
156  subgrad_t[p]=(float64_t*)LIBBMRM_CALLOC(nDim, float64_t);
157 
158  if (subgrad_t[p]==NULL || info[p]==NULL)
159  {
160  p3bmrm.exitflag=-2;
161  goto cleanup;
162  }
163 
164  info[p]->m_from=to;
165  to=((p+1)*N > (uint32_t)num_feats) ? (uint32_t)num_feats : (p+1)*N;
166  info[p]->m_N=to-info[p]->m_from;
167  }
168 
169  map= (bool*) LIBBMRM_CALLOC(BufSize, bool);
170 
171  if (map==NULL)
172  {
173  p3bmrm.exitflag=-2;
174  goto cleanup;
175  }
176 
177  memset( (bool*) map, true, BufSize);
178 
179  /* Temporary buffers */
180  beta_start= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
181 
182  beta_good= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
183 
184  b2= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
185 
186  diag_H2= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
187 
188  H2= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t);
189 
190  I_start= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
191 
192  I_good= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
193 
194  I2= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
195 
196  if (beta_start==NULL || beta_good==NULL || b2==NULL || diag_H2==NULL ||
197  I_start==NULL || I_good==NULL || I2==NULL || H2==NULL)
198  {
199  p3bmrm.exitflag=-2;
200  goto cleanup;
201  }
202 
203  p3bmrm.hist_Fp.resize_vector(BufSize);
204  p3bmrm.hist_Fd.resize_vector(BufSize);
205  p3bmrm.hist_wdist.resize_vector(BufSize);
206 
207  /* Iinitial solution */
208  Rt[0] = machine->risk(subgrad_t[0], W, info[0]);
209 
210  p3bmrm.nCP=0;
211  p3bmrm.nIter=0;
212  p3bmrm.exitflag=0;
213 
214  b[0]=-Rt[0];
215 
216  /* Cutting plane auxiliary double linked list */
217  LIBBMRM_MEMCPY(A, subgrad_t[0], nDim*sizeof(float64_t));
218  map[0]=false;
219  cp_list->address=&A[0];
220  cp_list->idx=0;
221  cp_list->prev=NULL;
222  cp_list->next=NULL;
223  CPList_head=cp_list;
224  CPList_tail=cp_list;
225 
226  for (uint32_t p=1; p<cp_models; ++p)
227  {
228  Rt[p] = machine->risk(subgrad_t[p], W, info[p]);
229  b[p]=SGVector<float64_t>::dot(subgrad_t[p], W, nDim) - Rt[p];
230  add_cutting_plane(&CPList_tail, map, A, find_free_idx(map, BufSize), subgrad_t[p], nDim);
231  }
232 
233  /* Compute initial value of Fp, Fd, assuming that W is zero vector */
234  R=0.0;
235 
236  for (uint32_t p=0; p<cp_models; ++p)
237  R+=Rt[p];
238 
239  sq_norm_W=SGVector<float64_t>::dot(W, W, nDim);
240  sq_norm_Wdiff=0.0;
241 
242  for (uint32_t j=0; j<nDim; ++j)
243  {
244  sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
245  }
246 
247  wdist=CMath::sqrt(sq_norm_Wdiff);
248 
249  p3bmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
250  p3bmrm.Fd=-LIBBMRM_PLUS_INF;
251  lastFp=p3bmrm.Fp;
252 
253  /* if there is initial W, then set K to be 0.01 times its norm */
254  K = (sq_norm_W == 0.0) ? 0.4 : 0.01*CMath::sqrt(sq_norm_W);
255 
256  LIBBMRM_MEMCPY(prevW, W, nDim*sizeof(float64_t));
257 
258  tstop=ttime.cur_time_diff(false);
259 
260  /* Keep history of Fp, Fd, and wdist */
261  p3bmrm.hist_Fp[0]=p3bmrm.Fp;
262  p3bmrm.hist_Fd[0]=p3bmrm.Fd;
263  p3bmrm.hist_wdist[0]=wdist;
264 
265  /* Verbose output */
266  if (verbose)
267  SG_SDEBUG("%4d: tim=%.3lf, Fp=%lf, Fd=%lf, R=%lf, K=%lf, CPmodels=%d\n",
268  p3bmrm.nIter, tstop-tstart, p3bmrm.Fp, p3bmrm.Fd, R, K, cp_models);
269 
270  if (verbose)
271  helper = machine->get_helper();
272 
273  /* main loop */
274  while (p3bmrm.exitflag==0)
275  {
276  tstart=ttime.cur_time_diff(false);
277  p3bmrm.nIter++;
278 
279  /* Update H */
280  if (p3bmrm.nIter==1)
281  {
282  cp_ptr=CPList_head;
283 
284  for (cp_i=0; cp_i<cp_models; ++cp_i) /* for all cutting planes */
285  {
286  A_1=get_cutting_plane(cp_ptr);
287 
288  for (uint32_t p=0; p<cp_models; ++p)
289  {
290  rsum=SGVector<float64_t>::dot(A_1, subgrad_t[p], nDim);
291 
292  H[LIBBMRM_INDEX(p, cp_i, BufSize)]=rsum;
293  }
294 
295  cp_ptr=cp_ptr->next;
296  }
297  }
298  else
299  {
300  cp_ptr=CPList_head;
301 
302  for (cp_i=0; cp_i<p3bmrm.nCP+cp_models; ++cp_i) /* for all cutting planes */
303  {
304  A_1=get_cutting_plane(cp_ptr);
305 
306  for (uint32_t p=0; p<cp_models; ++p)
307  {
308  rsum=SGVector<float64_t>::dot(A_1, subgrad_t[p], nDim);
309 
310  H[LIBBMRM_INDEX(p3bmrm.nCP+p, cp_i, BufSize)]=rsum;
311  }
312 
313  cp_ptr=cp_ptr->next;
314  }
315 
316  for (uint32_t i=0; i<p3bmrm.nCP; ++i)
317  for (uint32_t j=0; j<cp_models; ++j)
318  H[LIBBMRM_INDEX(i, p3bmrm.nCP+j, BufSize)]=
319  H[LIBBMRM_INDEX(p3bmrm.nCP+j, i, BufSize)];
320  }
321 
322  for (uint32_t p=0; p<cp_models; ++p)
323  diag_H[p3bmrm.nCP+p]=H[LIBBMRM_INDEX(p3bmrm.nCP+p, p3bmrm.nCP+p, BufSize)];
324 
325  p3bmrm.nCP+=cp_models;
326 
327  /* tune alpha cycle */
328  /* ------------------------------------------------------------------------ */
329  flag=true;
330  isThereGoodSolution=false;
331 
332  for (uint32_t p=0; p<cp_models; ++p)
333  {
334  I[p3bmrm.nCP-cp_models+p]=p+1;
335  beta[p3bmrm.nCP-cp_models+p]=0.0;
336  }
337 
338  LIBBMRM_MEMCPY(beta_start, beta, p3bmrm.nCP*sizeof(float64_t));
339  LIBBMRM_MEMCPY(I_start, I, p3bmrm.nCP*sizeof(uint32_t));
340  qp_cnt=0;
341 
342  if (tuneAlpha)
343  {
344  alpha_start=alpha; alpha=0.0;
345  LIBBMRM_MEMCPY(I2, I_start, p3bmrm.nCP*sizeof(uint32_t));
346 
347  /* add alpha-dependent terms to H, diag_h and b */
348  cp_ptr=CPList_head;
349 
350  for (uint32_t i=0; i<p3bmrm.nCP; ++i)
351  {
352  A_1=get_cutting_plane(cp_ptr);
353  cp_ptr=cp_ptr->next;
354 
355  rsum = SGVector<float64_t>::dot(A_1, prevW, nDim);
356 
357  b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
358  diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
359 
360  for (uint32_t j=0; j<p3bmrm.nCP; ++j)
361  H2[LIBBMRM_INDEX(i, j, BufSize)]=
362  H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
363 
364  }
365 
366  /* solve QP with current alpha */
367  qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, C, I2, S, beta,
368  p3bmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
369  p3bmrm.qp_exitflag=qp_exitflag.exitflag;
370  qp_cnt++;
371  Fd_alpha0=-qp_exitflag.QP;
372 
373  /* obtain w_t and check if norm(w_{t+1} -w_t) <= K */
374  memset(wt, 0, sizeof(float64_t)*nDim);
375  SGVector<float64_t>::vec1_plus_scalar_times_vec2(wt, 2*alpha/(_lambda+2*alpha), prevW, nDim);
376  cp_ptr=CPList_head;
377  for (uint32_t j=0; j<p3bmrm.nCP; ++j)
378  {
379  A_1=get_cutting_plane(cp_ptr);
380  cp_ptr=cp_ptr->next;
381  SGVector<float64_t>::vec1_plus_scalar_times_vec2(wt, -beta[j]/(_lambda+2*alpha), A_1, nDim);
382  }
383 
384  sq_norm_Wdiff=0.0;
385 
386  for (uint32_t i=0; i<nDim; ++i)
387  sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
388 
389  if (CMath::sqrt(sq_norm_Wdiff) <= K)
390  {
391  flag=false;
392 
393  if (alpha!=alpha_start)
394  alphaChanged=true;
395  }
396  else
397  {
398  alpha=alpha_start;
399  }
400 
401  while(flag)
402  {
403  LIBBMRM_MEMCPY(I2, I_start, p3bmrm.nCP*sizeof(uint32_t));
404  LIBBMRM_MEMCPY(beta, beta_start, p3bmrm.nCP*sizeof(float64_t));
405 
406  /* add alpha-dependent terms to H, diag_h and b */
407  cp_ptr=CPList_head;
408 
409  for (uint32_t i=0; i<p3bmrm.nCP; ++i)
410  {
411  A_1=get_cutting_plane(cp_ptr);
412  cp_ptr=cp_ptr->next;
413 
414  rsum = SGVector<float64_t>::dot(A_1, prevW, nDim);
415 
416  b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
417  diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
418 
419  for (uint32_t j=0; j<p3bmrm.nCP; ++j)
420  H2[LIBBMRM_INDEX(i, j, BufSize)]=H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
421  }
422 
423  /* solve QP with current alpha */
424  qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, C, I2, S, beta,
425  p3bmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
426  p3bmrm.qp_exitflag=qp_exitflag.exitflag;
427  qp_cnt++;
428 
429  /* obtain w_t and check if norm(w_{t+1}-w_t) <= K */
430  memset(wt, 0, sizeof(float64_t)*nDim);
431  SGVector<float64_t>::vec1_plus_scalar_times_vec2(wt, 2*alpha/(_lambda+2*alpha), prevW, nDim);
432  cp_ptr=CPList_head;
433  for (uint32_t j=0; j<p3bmrm.nCP; ++j)
434  {
435  A_1=get_cutting_plane(cp_ptr);
436  cp_ptr=cp_ptr->next;
437  SGVector<float64_t>::vec1_plus_scalar_times_vec2(wt, -beta[j]/(_lambda+2*alpha), A_1, nDim);
438  }
439 
440  sq_norm_Wdiff=0.0;
441 
442  for (uint32_t i=0; i<nDim; ++i)
443  sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
444 
445  if (CMath::sqrt(sq_norm_Wdiff) > K)
446  {
447  /* if there is a record of some good solution (i.e. adjust alpha by division by 2) */
448 
449  if (isThereGoodSolution)
450  {
451  LIBBMRM_MEMCPY(beta, beta_good, p3bmrm.nCP*sizeof(float64_t));
452  LIBBMRM_MEMCPY(I2, I_good, p3bmrm.nCP*sizeof(uint32_t));
453  alpha=alpha_good;
454  qp_exitflag=qp_exitflag_good;
455  flag=false;
456  }
457  else
458  {
459  if (alpha == 0)
460  {
461  alpha=1.0;
462  alphaChanged=true;
463  }
464  else
465  {
466  alpha*=2;
467  alphaChanged=true;
468  }
469  }
470  }
471  else
472  {
473  if (alpha > 0)
474  {
475  /* keep good solution and try for alpha /= 2 if previous alpha was 1 */
476  LIBBMRM_MEMCPY(beta_good, beta, p3bmrm.nCP*sizeof(float64_t));
477  LIBBMRM_MEMCPY(I_good, I2, p3bmrm.nCP*sizeof(uint32_t));
478  alpha_good=alpha;
479  qp_exitflag_good=qp_exitflag;
480  isThereGoodSolution=true;
481 
482  if (alpha!=1.0)
483  {
484  alpha/=2.0;
485  alphaChanged=true;
486  }
487  else
488  {
489  alpha=0.0;
490  alphaChanged=true;
491  }
492  }
493  else
494  {
495  flag=false;
496  }
497  }
498  }
499  }
500  else
501  {
502  alphaChanged=false;
503  LIBBMRM_MEMCPY(I2, I_start, p3bmrm.nCP*sizeof(uint32_t));
504  LIBBMRM_MEMCPY(beta, beta_start, p3bmrm.nCP*sizeof(float64_t));
505 
506  /* add alpha-dependent terms to H, diag_h and b */
507  cp_ptr=CPList_head;
508 
509  for (uint32_t i=0; i<p3bmrm.nCP; ++i)
510  {
511  A_1=get_cutting_plane(cp_ptr);
512  cp_ptr=cp_ptr->next;
513 
514  rsum = SGVector<float64_t>::dot(A_1, prevW, nDim);
515 
516  b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
517  diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
518 
519  for (uint32_t j=0; j<p3bmrm.nCP; ++j)
520  H2[LIBBMRM_INDEX(i, j, BufSize)]=H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
521  }
522 
523  /* solve QP with current alpha */
524  qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, C, I2, S, beta,
525  p3bmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
526  p3bmrm.qp_exitflag=qp_exitflag.exitflag;
527  qp_cnt++;
528  }
529  /* ----------------------------------------------------------------------------------------------- */
530 
531  /* Update ICPcounter (add one to unused and reset used) + compute number of active CPs */
532  p3bmrm.nzA=0;
533 
534  for (uint32_t aaa=0; aaa<p3bmrm.nCP; ++aaa)
535  {
536  if (beta[aaa]>epsilon)
537  {
538  ++p3bmrm.nzA;
539  icp_stats.ICPcounter[aaa]=0;
540  }
541  else
542  {
543  icp_stats.ICPcounter[aaa]+=1;
544  }
545  }
546 
547  /* W update */
548  memset(W, 0, sizeof(float64_t)*nDim);
549  SGVector<float64_t>::vec1_plus_scalar_times_vec2(W, 2*alpha/(_lambda+2*alpha), prevW, nDim);
550  cp_ptr=CPList_head;
551  for (uint32_t j=0; j<p3bmrm.nCP; ++j)
552  {
553  A_1=get_cutting_plane(cp_ptr);
554  cp_ptr=cp_ptr->next;
555  SGVector<float64_t>::vec1_plus_scalar_times_vec2(W, -beta[j]/(_lambda+2*alpha), A_1, nDim);
556  }
557 
558  /* risk and subgradient computation */
559  R=0.0;
560 
561  for (uint32_t p=0; p<cp_models; ++p)
562  {
563  Rt[p] = machine->risk(subgrad_t[p], W, info[p]);
564  b[p3bmrm.nCP+p] = SGVector<float64_t>::dot(subgrad_t[p], W, nDim) - Rt[p];
565  add_cutting_plane(&CPList_tail, map, A, find_free_idx(map, BufSize), subgrad_t[p], nDim);
566  R+=Rt[p];
567  }
568 
569  sq_norm_W=SGVector<float64_t>::dot(W, W, nDim);
570  sq_norm_prevW=SGVector<float64_t>::dot(prevW, prevW, nDim);
571  sq_norm_Wdiff=0.0;
572 
573  for (uint32_t j=0; j<nDim; ++j)
574  {
575  sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
576  }
577 
578  /* compute Fp and Fd */
579  p3bmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
580  p3bmrm.Fd=-qp_exitflag.QP + ((alpha*_lambda)/(_lambda + 2*alpha))*sq_norm_prevW;
581 
582  /* gamma + tuneAlpha flag */
583  if (alphaChanged)
584  {
585  eps=1.0-(p3bmrm.Fd/p3bmrm.Fp);
586  gamma=(lastFp*(1-eps)-Fd_alpha0)/(Tmax*(1-eps));
587  }
588 
589  if ((lastFp-p3bmrm.Fp) <= gamma)
590  {
591  tuneAlpha=true;
592  }
593  else
594  {
595  tuneAlpha=false;
596  }
597 
598  /* Stopping conditions - set only with nonzero alpha */
599  if (alpha==0.0)
600  {
601  if (p3bmrm.Fp-p3bmrm.Fd<=TolRel*LIBBMRM_ABS(p3bmrm.Fp))
602  p3bmrm.exitflag=1;
603 
604  if (p3bmrm.Fp-p3bmrm.Fd<=TolAbs)
605  p3bmrm.exitflag=2;
606  }
607 
608  if (p3bmrm.nCP>=BufSize)
609  p3bmrm.exitflag=-1;
610 
611  tstop=ttime.cur_time_diff(false);
612 
613  /* compute wdist (= || W_{t+1} - W_{t} || ) */
614  sq_norm_Wdiff=0.0;
615 
616  for (uint32_t i=0; i<nDim; ++i)
617  {
618  sq_norm_Wdiff+=(W[i]-prevW[i])*(W[i]-prevW[i]);
619  }
620 
621  wdist=CMath::sqrt(sq_norm_Wdiff);
622 
623  /* Keep history of Fp, Fd and wdist */
624  p3bmrm.hist_Fp[p3bmrm.nIter]=p3bmrm.Fp;
625  p3bmrm.hist_Fd[p3bmrm.nIter]=p3bmrm.Fd;
626  p3bmrm.hist_wdist[p3bmrm.nIter]=wdist;
627 
628  /* Verbose output */
629  if (verbose)
630  SG_SDEBUG("%4d: tim=%.3lf, Fp=%lf, Fd=%lf, (Fp-Fd)=%lf, (Fp-Fd)/Fp=%lf, R=%lf, nCP=%d, nzA=%d, wdist=%lf, alpha=%lf, qp_cnt=%d, gamma=%lf, tuneAlpha=%d\n",
631  p3bmrm.nIter, tstop-tstart, p3bmrm.Fp, p3bmrm.Fd, p3bmrm.Fp-p3bmrm.Fd,
632  (p3bmrm.Fp-p3bmrm.Fd)/p3bmrm.Fp, R, p3bmrm.nCP, p3bmrm.nzA, wdist, alpha,
633  qp_cnt, gamma, tuneAlpha);
634 
635  /* Check size of Buffer */
636  if (p3bmrm.nCP>=BufSize)
637  {
638  p3bmrm.exitflag=-2;
639  SG_SERROR("Buffer exceeded.\n")
640  }
641 
642  /* keep w_t + Fp */
643  LIBBMRM_MEMCPY(prevW, W, nDim*sizeof(float64_t));
644  lastFp=p3bmrm.Fp;
645 
646  /* Inactive Cutting Planes (ICP) removal */
647  if (cleanICP)
648  {
649  clean_icp(&icp_stats, p3bmrm, &CPList_head,
650  &CPList_tail, H, diag_H, beta, map,
651  cleanAfter, b, I, cp_models);
652  }
653 
654  // next CP would exceed BufSize
655  if (p3bmrm.nCP+1 >= BufSize)
656  p3bmrm.exitflag=-1;
657 
658  /* Debug: compute objective and training error */
659  if (verbose)
660  {
661  SGVector<float64_t> w_debug(W, nDim, false);
662  float64_t primal = CSOSVMHelper::primal_objective(w_debug, model, _lambda);
663  float64_t train_error = CSOSVMHelper::average_loss(w_debug, model);
664  helper->add_debug_info(primal, p3bmrm.nIter, train_error);
665  }
666  } /* end of main loop */
667 
668  if (verbose)
669  {
670  helper->terminate();
671  SG_UNREF(helper);
672  }
673 
674  p3bmrm.hist_Fp.resize_vector(p3bmrm.nIter);
675  p3bmrm.hist_Fd.resize_vector(p3bmrm.nIter);
676  p3bmrm.hist_wdist.resize_vector(p3bmrm.nIter);
677 
678  cp_ptr=CPList_head;
679 
680  while(cp_ptr!=NULL)
681  {
682  cp_ptr2=cp_ptr;
683  cp_ptr=cp_ptr->next;
684  LIBBMRM_FREE(cp_ptr2);
685  cp_ptr2=NULL;
686  }
687 
688  cp_list=NULL;
689 
690 cleanup:
691 
692  LIBBMRM_FREE(H);
693  LIBBMRM_FREE(b);
694  LIBBMRM_FREE(beta);
695  LIBBMRM_FREE(A);
696  LIBBMRM_FREE(diag_H);
697  LIBBMRM_FREE(I);
698  LIBBMRM_FREE(icp_stats.ICPcounter);
699  LIBBMRM_FREE(icp_stats.ICPs);
700  LIBBMRM_FREE(icp_stats.ACPs);
701  LIBBMRM_FREE(icp_stats.H_buff);
702  LIBBMRM_FREE(map);
703  LIBBMRM_FREE(prevW);
704  LIBBMRM_FREE(wt);
705  LIBBMRM_FREE(beta_start);
706  LIBBMRM_FREE(beta_good);
707  LIBBMRM_FREE(I_start);
708  LIBBMRM_FREE(I_good);
709  LIBBMRM_FREE(I2);
710  LIBBMRM_FREE(b2);
711  LIBBMRM_FREE(diag_H2);
712  LIBBMRM_FREE(H2);
713  LIBBMRM_FREE(C);
714  LIBBMRM_FREE(S);
715  LIBBMRM_FREE(Rt);
716 
717  if (cp_list)
718  LIBBMRM_FREE(cp_list);
719 
720  for (uint32_t p=0; p<cp_models; ++p)
721  {
722  LIBBMRM_FREE(subgrad_t[p]);
723  LIBBMRM_FREE(info[p]);
724  }
725 
726  LIBBMRM_FREE(subgrad_t);
727  LIBBMRM_FREE(info);
728 
729  SG_UNREF(model);
730 
731  return(p3bmrm);
732 }
733 }

SHOGUN Machine Learning Toolbox - Documentation