SHOGUN  v3.0.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
libppbm.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  * libp3bm.h: Implementation of the Proximal Point P-BMRM 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
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 static uint32_t BufSize;
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  bool verbose)
46 {
47  BmrmStatistics ppbmrm;
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, *subgrad, *A, QPSolverTolRel, C=1.0;
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=1;
56  CStructuredModel* model=machine->get_model();
57  uint32_t nDim=model->get_dim();
58  CSOSVMHelper* helper = NULL;
59  uint32_t qp_cnt=0;
60  bmrm_ll *CPList_head, *CPList_tail, *cp_ptr, *cp_ptr2, *cp_list=NULL;
61  float64_t *A_1=NULL, *A_2=NULL;
62  bool *map=NULL, tuneAlpha=true, flag=true, alphaChanged=false, isThereGoodSolution=false;
63 
64  CTime ttime;
65  float64_t tstart, tstop;
66 
67 
68  tstart=ttime.cur_time_diff(false);
69 
70  BufSize=_BufSize;
71  QPSolverTolRel=1e-9;
72 
73  H=NULL;
74  b=NULL;
75  beta=NULL;
76  A=NULL;
77  subgrad=NULL;
78  diag_H=NULL;
79  I=NULL;
80  prevW=NULL;
81  wt=NULL;
82  diag_H2=NULL;
83  b2=NULL;
84  I2=NULL;
85  H2=NULL;
86  I_good=NULL;
87  I_start=NULL;
88  beta_start=NULL;
89  beta_good=NULL;
90 
91  alpha=0.0;
92 
94 
95  A= (float64_t*) LIBBMRM_CALLOC(nDim*BufSize, float64_t);
96 
97  b= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
98 
99  beta= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
100 
101  subgrad= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
102 
103  diag_H= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
104 
105  I= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
106 
107  /* structure for maintaining inactive CPs info */
108  ICP_stats icp_stats;
109  icp_stats.maxCPs = BufSize;
110  icp_stats.ICPcounter= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
111  icp_stats.ICPs= (float64_t**) LIBBMRM_CALLOC(BufSize, float64_t*);
112  icp_stats.ACPs= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
113 
114  cp_list= (bmrm_ll*) LIBBMRM_CALLOC(1, bmrm_ll);
115 
116  prevW= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
117 
118  wt= (float64_t*) LIBBMRM_CALLOC(nDim, float64_t);
119 
120  if (H==NULL || A==NULL || b==NULL || beta==NULL || subgrad==NULL ||
121  diag_H==NULL || I==NULL || icp_stats.ICPcounter==NULL ||
122  icp_stats.ICPs==NULL || icp_stats.ACPs==NULL ||
123  cp_list==NULL || prevW==NULL || wt==NULL)
124  {
125  ppbmrm.exitflag=-2;
126  goto cleanup;
127  }
128 
129  map= (bool*) LIBBMRM_CALLOC(BufSize, bool);
130 
131  if (map==NULL)
132  {
133  ppbmrm.exitflag=-2;
134  goto cleanup;
135  }
136 
137  memset( (bool*) map, true, BufSize);
138 
139  /* Temporary buffers for ICP removal */
140  icp_stats.H_buff= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t);
141 
142  if (icp_stats.H_buff==NULL)
143  {
144  ppbmrm.exitflag=-2;
145  goto cleanup;
146  }
147 
148  /* Temporary buffers */
149  beta_start= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
150 
151  beta_good= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
152 
153  b2= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
154 
155  diag_H2= (float64_t*) LIBBMRM_CALLOC(BufSize, float64_t);
156 
157  H2= (float64_t*) LIBBMRM_CALLOC(BufSize*BufSize, float64_t);
158 
159  I_start= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
160 
161  I_good= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
162 
163  I2= (uint32_t*) LIBBMRM_CALLOC(BufSize, uint32_t);
164 
165  if (beta_start==NULL || beta_good==NULL || b2==NULL || diag_H2==NULL ||
166  I_start==NULL || I_good==NULL || I2==NULL || H2==NULL)
167  {
168  ppbmrm.exitflag=-2;
169  goto cleanup;
170  }
171 
172  ppbmrm.hist_Fp.resize_vector(BufSize);
173  ppbmrm.hist_Fd.resize_vector(BufSize);
174  ppbmrm.hist_wdist.resize_vector(BufSize);
175 
176  /* Iinitial solution */
177  R = machine->risk(subgrad, W);
178 
179  ppbmrm.nCP=0;
180  ppbmrm.nIter=0;
181  ppbmrm.exitflag=0;
182 
183  b[0]=-R;
184 
185  /* Cutting plane auxiliary double linked list */
186  LIBBMRM_MEMCPY(A, subgrad, nDim*sizeof(float64_t));
187  map[0]=false;
188  cp_list->address=&A[0];
189  cp_list->idx=0;
190  cp_list->prev=NULL;
191  cp_list->next=NULL;
192  CPList_head=cp_list;
193  CPList_tail=cp_list;
194 
195  /* Compute initial value of Fp, Fd, assuming that W is zero vector */
196  sq_norm_Wdiff=0.0;
197 
198  b[0] = SGVector<float64_t>::dot(subgrad, W, nDim);
199  sq_norm_W = SGVector<float64_t>::dot(W, W, nDim);
200  for (uint32_t j=0; j<nDim; ++j)
201  {
202  sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
203  }
204 
205  ppbmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
206  ppbmrm.Fd=-LIBBMRM_PLUS_INF;
207  lastFp=ppbmrm.Fp;
208  wdist=CMath::sqrt(sq_norm_Wdiff);
209 
210  K = (sq_norm_W == 0.0) ? 0.4 : 0.01*CMath::sqrt(sq_norm_W);
211 
212  LIBBMRM_MEMCPY(prevW, W, nDim*sizeof(float64_t));
213 
214  tstop=ttime.cur_time_diff(false);
215 
216  /* Keep history of Fp, Fd, wdist */
217  ppbmrm.hist_Fp[0]=ppbmrm.Fp;
218  ppbmrm.hist_Fd[0]=ppbmrm.Fd;
219  ppbmrm.hist_wdist[0]=wdist;
220 
221  /* Verbose output */
222 
223  if (verbose)
224  SG_SDEBUG("%4d: tim=%.3lf, Fp=%lf, Fd=%lf, R=%lf, K=%lf\n",
225  ppbmrm.nIter, tstop-tstart, ppbmrm.Fp, ppbmrm.Fd, R, K);
226 
227  if (verbose)
228  helper = machine->get_helper();
229 
230  /* main loop */
231 
232  while (ppbmrm.exitflag==0)
233  {
234  tstart=ttime.cur_time_diff(false);
235  ppbmrm.nIter++;
236 
237  /* Update H */
238 
239  if (ppbmrm.nCP>0)
240  {
241  A_2=get_cutting_plane(CPList_tail);
242  cp_ptr=CPList_head;
243 
244  for (uint32_t i=0; i<ppbmrm.nCP; ++i)
245  {
246  A_1=get_cutting_plane(cp_ptr);
247  cp_ptr=cp_ptr->next;
248  rsum=SGVector<float64_t>::dot(A_1, A_2, nDim);
249 
250  H[LIBBMRM_INDEX(ppbmrm.nCP, i, BufSize)]
251  = H[LIBBMRM_INDEX(i, ppbmrm.nCP, BufSize)]
252  = rsum;
253  }
254  }
255 
256  A_2=get_cutting_plane(CPList_tail);
257  rsum = SGVector<float64_t>::dot(A_2, A_2, nDim);
258 
259  H[LIBBMRM_INDEX(ppbmrm.nCP, ppbmrm.nCP, BufSize)]=rsum;
260 
261  diag_H[ppbmrm.nCP]=H[LIBBMRM_INDEX(ppbmrm.nCP, ppbmrm.nCP, BufSize)];
262  I[ppbmrm.nCP]=1;
263 
264  ppbmrm.nCP++;
265  beta[ppbmrm.nCP]=0.0; // [beta; 0]
266 
267  /* tune alpha cycle */
268  /* ---------------------------------------------------------------------- */
269 
270  flag=true;
271  isThereGoodSolution=false;
272  LIBBMRM_MEMCPY(beta_start, beta, ppbmrm.nCP*sizeof(float64_t));
273  LIBBMRM_MEMCPY(I_start, I, ppbmrm.nCP*sizeof(uint32_t));
274  qp_cnt=0;
275  alpha_good=alpha;
276 
277  if (tuneAlpha)
278  {
279  alpha_start=alpha; alpha=0.0;
280  beta[ppbmrm.nCP]=0.0;
281  LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t));
282  I2[ppbmrm.nCP]=1;
283 
284  /* add alpha-dependent terms to H, diag_h and b */
285  cp_ptr=CPList_head;
286 
287  for (uint32_t i=0; i<ppbmrm.nCP; ++i)
288  {
289  A_1=get_cutting_plane(cp_ptr);
290  cp_ptr=cp_ptr->next;
291 
292  rsum = SGVector<float64_t>::dot(A_1, prevW, nDim);
293 
294  b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
295  diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
296 
297  for (uint32_t j=0; j<ppbmrm.nCP; ++j)
298  H2[LIBBMRM_INDEX(i, j, BufSize)]=
299  H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
300  }
301 
302  /* solve QP with current alpha */
303  qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta,
304  ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
305  ppbmrm.qp_exitflag=qp_exitflag.exitflag;
306  qp_cnt++;
307  Fd_alpha0=-qp_exitflag.QP;
308 
309  /* obtain w_t and check if norm(w_{t+1} -w_t) <= K */
310  for (uint32_t i=0; i<nDim; ++i)
311  {
312  rsum=0.0;
313  cp_ptr=CPList_head;
314 
315  for (uint32_t j=0; j<ppbmrm.nCP; ++j)
316  {
317  A_1=get_cutting_plane(cp_ptr);
318  cp_ptr=cp_ptr->next;
319  rsum+=A_1[i]*beta[j];
320  }
321 
322  wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha);
323  }
324 
325  sq_norm_Wdiff=0.0;
326 
327  for (uint32_t i=0; i<nDim; ++i)
328  sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
329 
330  if (CMath::sqrt(sq_norm_Wdiff) <= K)
331  {
332  flag=false;
333 
334  if (alpha!=alpha_start)
335  alphaChanged=true;
336  }
337  else
338  {
339  alpha=alpha_start;
340  }
341 
342  while(flag)
343  {
344  LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t));
345  LIBBMRM_MEMCPY(beta, beta_start, ppbmrm.nCP*sizeof(float64_t));
346  I2[ppbmrm.nCP]=1;
347  beta[ppbmrm.nCP]=0.0;
348 
349  /* add alpha-dependent terms to H, diag_h and b */
350  cp_ptr=CPList_head;
351 
352  for (uint32_t i=0; i<ppbmrm.nCP; ++i)
353  {
354  A_1=get_cutting_plane(cp_ptr);
355  cp_ptr=cp_ptr->next;
356 
357  rsum = SGVector<float64_t>::dot(A_1, prevW, nDim);
358 
359  b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
360  diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
361 
362  for (uint32_t j=0; j<ppbmrm.nCP; ++j)
363  H2[LIBBMRM_INDEX(i, j, BufSize)]=H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
364  }
365 
366  /* solve QP with current alpha */
367  qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta,
368  ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
369  ppbmrm.qp_exitflag=qp_exitflag.exitflag;
370  qp_cnt++;
371 
372  /* obtain w_t and check if norm(w_{t+1}-w_t) <= K */
373  for (uint32_t i=0; i<nDim; ++i)
374  {
375  rsum=0.0;
376  cp_ptr=CPList_head;
377 
378  for (uint32_t j=0; j<ppbmrm.nCP; ++j)
379  {
380  A_1=get_cutting_plane(cp_ptr);
381  cp_ptr=cp_ptr->next;
382  rsum+=A_1[i]*beta[j];
383  }
384 
385  wt[i]=(2*alpha*prevW[i] - rsum)/(_lambda+2*alpha);
386  }
387 
388  sq_norm_Wdiff=0.0;
389  for (uint32_t i=0; i<nDim; ++i)
390  sq_norm_Wdiff+=(wt[i]-prevW[i])*(wt[i]-prevW[i]);
391 
392  if (CMath::sqrt(sq_norm_Wdiff) > K)
393  {
394  /* if there is a record of some good solution
395  * (i.e. adjust alpha by division by 2) */
396 
397  if (isThereGoodSolution)
398  {
399  LIBBMRM_MEMCPY(beta, beta_good, ppbmrm.nCP*sizeof(float64_t));
400  LIBBMRM_MEMCPY(I2, I_good, ppbmrm.nCP*sizeof(uint32_t));
401  alpha=alpha_good;
402  qp_exitflag=qp_exitflag_good;
403  flag=false;
404  }
405  else
406  {
407  if (alpha == 0)
408  {
409  alpha=1.0;
410  alphaChanged=true;
411  }
412  else
413  {
414  alpha*=2;
415  alphaChanged=true;
416  }
417  }
418  }
419  else
420  {
421  if (alpha > 0)
422  {
423  /* keep good solution and try for alpha /= 2 if previous alpha was 1 */
424  LIBBMRM_MEMCPY(beta_good, beta, ppbmrm.nCP*sizeof(float64_t));
425  LIBBMRM_MEMCPY(I_good, I2, ppbmrm.nCP*sizeof(uint32_t));
426  alpha_good=alpha;
427  qp_exitflag_good=qp_exitflag;
428  isThereGoodSolution=true;
429 
430  if (alpha!=1.0)
431  {
432  alpha/=2.0;
433  alphaChanged=true;
434  }
435  else
436  {
437  alpha=0.0;
438  alphaChanged=true;
439  }
440  }
441  else
442  {
443  flag=false;
444  }
445  }
446  }
447  }
448  else
449  {
450  alphaChanged=false;
451  LIBBMRM_MEMCPY(I2, I_start, ppbmrm.nCP*sizeof(uint32_t));
452  LIBBMRM_MEMCPY(beta, beta_start, ppbmrm.nCP*sizeof(float64_t));
453 
454  /* add alpha-dependent terms to H, diag_h and b */
455  cp_ptr=CPList_head;
456 
457  for (uint32_t i=0; i<ppbmrm.nCP; ++i)
458  {
459  A_1=get_cutting_plane(cp_ptr);
460  cp_ptr=cp_ptr->next;
461 
462  rsum = SGVector<float64_t>::dot(A_1, prevW, nDim);
463 
464  b2[i]=b[i]-((2*alpha)/(_lambda+2*alpha))*rsum;
465  diag_H2[i]=diag_H[i]/(_lambda+2*alpha);
466 
467  for (uint32_t j=0; j<ppbmrm.nCP; ++j)
468  H2[LIBBMRM_INDEX(i, j, BufSize)]=
469  H[LIBBMRM_INDEX(i, j, BufSize)]/(_lambda+2*alpha);
470  }
471  /* solve QP with current alpha */
472  qp_exitflag=libqp_splx_solver(&get_col, diag_H2, b2, &C, I2, &S, beta,
473  ppbmrm.nCP, QPSolverMaxIter, 0.0, QPSolverTolRel, -LIBBMRM_PLUS_INF, 0);
474  ppbmrm.qp_exitflag=qp_exitflag.exitflag;
475  qp_cnt++;
476  }
477 
478  /* ----------------------------------------------------------------------------------------------- */
479 
480  /* Update ICPcounter (add one to unused and reset used) + compute number of active CPs */
481  ppbmrm.nzA=0;
482 
483  for (uint32_t aaa=0; aaa<ppbmrm.nCP; ++aaa)
484  {
485  if (beta[aaa]>epsilon)
486  {
487  ++ppbmrm.nzA;
488  icp_stats.ICPcounter[aaa]=0;
489  }
490  else
491  {
492  icp_stats.ICPcounter[aaa]+=1;
493  }
494  }
495 
496  /* W update */
497  for (uint32_t i=0; i<nDim; ++i)
498  {
499  rsum=0.0;
500  cp_ptr=CPList_head;
501 
502  for (uint32_t j=0; j<ppbmrm.nCP; ++j)
503  {
504  A_1=get_cutting_plane(cp_ptr);
505  cp_ptr=cp_ptr->next;
506  rsum+=A_1[i]*beta[j];
507  }
508 
509  W[i]=(2*alpha*prevW[i]-rsum)/(_lambda+2*alpha);
510  }
511 
512  /* risk and subgradient computation */
513  R = machine->risk(subgrad, W);
514  add_cutting_plane(&CPList_tail, map, A,
515  find_free_idx(map, BufSize), subgrad, nDim);
516 
517  sq_norm_W=SGVector<float64_t>::dot(W, W, nDim);
518  sq_norm_prevW=SGVector<float64_t>::dot(prevW, prevW, nDim);
519  b[ppbmrm.nCP]=SGVector<float64_t>::dot(subgrad, W, nDim) - R;
520 
521  sq_norm_Wdiff=0.0;
522  for (uint32_t j=0; j<nDim; ++j)
523  {
524  sq_norm_Wdiff+=(W[j]-prevW[j])*(W[j]-prevW[j]);
525  }
526 
527  /* compute Fp and Fd */
528  ppbmrm.Fp=R+0.5*_lambda*sq_norm_W + alpha*sq_norm_Wdiff;
529  ppbmrm.Fd=-qp_exitflag.QP+((alpha*_lambda)/(_lambda + 2*alpha))*sq_norm_prevW;
530 
531  /* gamma + tuneAlpha flag */
532  if (alphaChanged)
533  {
534  eps=1.0-(ppbmrm.Fd/ppbmrm.Fp);
535  gamma=(lastFp*(1-eps)-Fd_alpha0)/(Tmax*(1-eps));
536  }
537 
538  if ((lastFp-ppbmrm.Fp) <= gamma)
539  {
540  tuneAlpha=true;
541  }
542  else
543  {
544  tuneAlpha=false;
545  }
546 
547  /* Stopping conditions - set only with nonzero alpha */
548  if (alpha==0.0)
549  {
550  if (ppbmrm.Fp-ppbmrm.Fd<=TolRel*LIBBMRM_ABS(ppbmrm.Fp))
551  ppbmrm.exitflag=1;
552 
553  if (ppbmrm.Fp-ppbmrm.Fd<=TolAbs)
554  ppbmrm.exitflag=2;
555  }
556 
557  if (ppbmrm.nCP>=BufSize)
558  ppbmrm.exitflag=-1;
559 
560  tstop=ttime.cur_time_diff(false);
561 
562  /* compute wdist (= || W_{t+1} - W_{t} || ) */
563  sq_norm_Wdiff=0.0;
564 
565  for (uint32_t i=0; i<nDim; ++i)
566  {
567  sq_norm_Wdiff+=(W[i]-prevW[i])*(W[i]-prevW[i]);
568  }
569 
570  wdist=CMath::sqrt(sq_norm_Wdiff);
571 
572  /* Keep history of Fp, Fd, wdist */
573  ppbmrm.hist_Fp[ppbmrm.nIter]=ppbmrm.Fp;
574  ppbmrm.hist_Fd[ppbmrm.nIter]=ppbmrm.Fd;
575  ppbmrm.hist_wdist[ppbmrm.nIter]=wdist;
576 
577  /* Verbose output */
578  if (verbose)
579  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",
580  ppbmrm.nIter, tstop-tstart, ppbmrm.Fp, ppbmrm.Fd, ppbmrm.Fp-ppbmrm.Fd,
581  (ppbmrm.Fp-ppbmrm.Fd)/ppbmrm.Fp, R, ppbmrm.nCP, ppbmrm.nzA, wdist, alpha,
582  qp_cnt, gamma, tuneAlpha);
583 
584  /* Check size of Buffer */
585  if (ppbmrm.nCP>=BufSize)
586  {
587  ppbmrm.exitflag=-2;
588  SG_SERROR("Buffer exceeded.\n")
589  }
590 
591  /* keep w_t + Fp */
592  LIBBMRM_MEMCPY(prevW, W, nDim*sizeof(float64_t));
593  lastFp=ppbmrm.Fp;
594 
595  /* Inactive Cutting Planes (ICP) removal */
596  if (cleanICP)
597  {
598  clean_icp(&icp_stats, ppbmrm, &CPList_head, &CPList_tail, H, diag_H, beta, map, cleanAfter, b, I);
599  }
600 
601  /* Debug: compute objective and training error */
602  if (verbose)
603  {
604  SGVector<float64_t> w_debug(W, nDim, false);
605  float64_t primal = CSOSVMHelper::primal_objective(w_debug, model, _lambda);
606  float64_t train_error = CSOSVMHelper::average_loss(w_debug, model);
607  helper->add_debug_info(primal, ppbmrm.nIter, train_error);
608  }
609  } /* end of main loop */
610 
611  if (verbose)
612  {
613  helper->terminate();
614  SG_UNREF(helper);
615  }
616 
617  ppbmrm.hist_Fp.resize_vector(ppbmrm.nIter);
618  ppbmrm.hist_Fd.resize_vector(ppbmrm.nIter);
619  ppbmrm.hist_wdist.resize_vector(ppbmrm.nIter);
620 
621  cp_ptr=CPList_head;
622 
623  while(cp_ptr!=NULL)
624  {
625  cp_ptr2=cp_ptr;
626  cp_ptr=cp_ptr->next;
627  LIBBMRM_FREE(cp_ptr2);
628  cp_ptr2=NULL;
629  }
630 
631  cp_list=NULL;
632 
633 cleanup:
634 
635  LIBBMRM_FREE(H);
636  LIBBMRM_FREE(b);
637  LIBBMRM_FREE(beta);
638  LIBBMRM_FREE(A);
639  LIBBMRM_FREE(subgrad);
640  LIBBMRM_FREE(diag_H);
641  LIBBMRM_FREE(I);
642  LIBBMRM_FREE(icp_stats.ICPcounter);
643  LIBBMRM_FREE(icp_stats.ICPs);
644  LIBBMRM_FREE(icp_stats.ACPs);
645  LIBBMRM_FREE(icp_stats.H_buff);
646  LIBBMRM_FREE(map);
647  LIBBMRM_FREE(prevW);
648  LIBBMRM_FREE(wt);
649  LIBBMRM_FREE(beta_start);
650  LIBBMRM_FREE(beta_good);
651  LIBBMRM_FREE(I_start);
652  LIBBMRM_FREE(I_good);
653  LIBBMRM_FREE(I2);
654  LIBBMRM_FREE(b2);
655  LIBBMRM_FREE(diag_H2);
656  LIBBMRM_FREE(H2);
657 
658  if (cp_list)
659  LIBBMRM_FREE(cp_list);
660 
661  SG_UNREF(model);
662 
663  return(ppbmrm);
664 }
665 }
#define LIBBMRM_CALLOC(x, y)
Definition: libbmrm.h:23
Class Time that implements a stopwatch based on either cpu time or wall clock time.
Definition: Time.h:46
static float64_t dot(const bool *v1, const bool *v2, int32_t n)
compute dot product between v1 and v2 (blas optimized)
Definition: SGVector.h:344
static float64_t * H
Definition: libbmrm.cpp:24
static const double * get_col(uint32_t j)
float64_t * H_buff
Definition: libbmrm.h:61
Class DualLibQPBMSOSVM that uses Bundle Methods for Regularized Risk Minimization algorithms for stru...
float64_t * get_cutting_plane(bmrm_ll *ptr)
Definition: libbmrm.h:116
uint32_t idx
Definition: libbmrm.h:42
static float64_t * H2
Definition: libp3bm.cpp:23
SGVector< float64_t > hist_wdist
bmrm_ll * next
Definition: libbmrm.h:38
uint32_t * ACPs
Definition: libbmrm.h:58
SGVector< float64_t > hist_Fd
static float64_t primal_objective(SGVector< float64_t > w, CStructuredModel *model, float64_t lbda)
Definition: SOSVMHelper.cpp:56
virtual int32_t get_dim() const =0
uint32_t find_free_idx(bool *map, uint32_t size)
Definition: libbmrm.h:124
static float64_t average_loss(SGVector< float64_t > w, CStructuredModel *model)
Definition: SOSVMHelper.cpp:85
static const float64_t epsilon
Definition: libbmrm.cpp:22
float64_t cur_time_diff(bool verbose=false)
Definition: Time.cpp:67
#define LIBBMRM_FREE(x)
Definition: libbmrm.h:25
#define LIBBMRM_ABS(A)
Definition: libbmrm.h:29
static const uint32_t QPSolverMaxIter
Definition: libbmrm.cpp:21
BmrmStatistics svm_ppbm_solver(CDualLibQPBMSOSVM *machine, float64_t *W, float64_t TolRel, float64_t TolAbs, float64_t _lambda, uint32_t _BufSize, bool cleanICP, uint32_t cleanAfter, float64_t K, uint32_t Tmax, bool verbose)
Definition: libppbm.cpp:34
bmrm_ll * prev
Definition: libbmrm.h:36
double float64_t
Definition: common.h:48
long double floatmax_t
Definition: common.h:49
class CSOSVMHelper contains helper functions to compute primal objectives, dual objectives, average training losses, duality gaps etc. These values will be recorded to check convergence. This class is inspired by the matlab implementation of the block coordinate Frank-Wolfe SOSVM solver [1].
Definition: SOSVMHelper.h:29
virtual float64_t risk(float64_t *subgrad, float64_t *W, TMultipleCPinfo *info=0, EStructRiskType rtype=N_SLACK_MARGIN_RESCALING)
SGVector< float64_t > hist_Fp
#define LIBBMRM_MEMCPY(x, y, z)
Definition: libbmrm.h:26
float64_t ** ICPs
Definition: libbmrm.h:55
uint32_t maxCPs
Definition: libbmrm.h:49
Class CStructuredModel that represents the application specific model and contains most of the applic...
#define LIBBMRM_INDEX(ROW, COL, NUM_ROWS)
Definition: libbmrm.h:28
#define SG_UNREF(x)
Definition: SGObject.h:54
#define SG_SDEBUG(...)
Definition: SGIO.h:170
virtual void add_debug_info(float64_t primal, float64_t eff_pass, float64_t train_error, float64_t dual=-1, float64_t dgap=-1)
void add_cutting_plane(bmrm_ll **tail, bool *map, float64_t *A, uint32_t free_idx, float64_t *cp_data, uint32_t dim)
Definition: libbmrm.cpp:27
uint32_t * ICPcounter
Definition: libbmrm.h:52
#define SG_SERROR(...)
Definition: SGIO.h:181
float64_t * address
Definition: libbmrm.h:40
CStructuredModel * get_model() const
void resize_vector(int32_t n)
Definition: SGVector.cpp:307
static float32_t sqrt(float32_t x)
x^0.5
Definition: Math.h:252
#define LIBBMRM_PLUS_INF
Definition: libbmrm.h:22
void clean_icp(ICP_stats *icp_stats, BmrmStatistics &bmrm, bmrm_ll **head, bmrm_ll **tail, float64_t *&Hmat, float64_t *&diag_H, float64_t *&beta, bool *&map, uint32_t cleanAfter, float64_t *&b, uint32_t *&I, uint32_t cp_models)
Definition: libbmrm.cpp:89
static uint32_t BufSize
Definition: libbmrm.cpp:25

SHOGUN Machine Learning Toolbox - Documentation