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

SHOGUN Machine Learning Toolbox - Documentation