SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules
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 }
#define LIBBMRM_CALLOC(x, y)
Definition: libbmrm.h:24
Class Time that implements a stopwatch based on either cpu time or wall clock time.
Definition: Time.h:47
static float64_t * H
Definition: libbmrm.cpp:26
static const double * get_col(uint32_t j)
float64_t * H_buff
Definition: libbmrm.h:65
Class DualLibQPBMSOSVM that uses Bundle Methods for Regularized Risk Minimization algorithms for stru...
float64_t * get_cutting_plane(bmrm_ll *ptr)
Definition: libbmrm.h:120
uint32_t idx
Definition: libbmrm.h:46
static float64_t * H2
Definition: libp3bm.cpp:23
SGVector< float64_t > hist_wdist
bmrm_ll * next
Definition: libbmrm.h:42
virtual int32_t get_num_vectors() const =0
uint32_t * ACPs
Definition: libbmrm.h:62
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:128
static const float64_t epsilon
Definition: libbmrm.cpp:24
float64_t cur_time_diff(bool verbose=false)
Definition: Time.cpp:68
#define LIBBMRM_FREE(x)
Definition: libbmrm.h:26
#define LIBBMRM_ABS(A)
Definition: libbmrm.h:30
static const uint32_t QPSolverMaxIter
Definition: libbmrm.cpp:23
bmrm_ll * prev
Definition: libbmrm.h:40
double float64_t
Definition: common.h:50
long double floatmax_t
Definition: common.h:51
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:31
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:27
float64_t ** ICPs
Definition: libbmrm.h:59
static float64_t average_loss(SGVector< float64_t > w, CStructuredModel *model, bool is_ub=false)
Definition: SOSVMHelper.cpp:87
static void vec1_plus_scalar_times_vec2(T *vec1, const T scalar, const T *vec2, int32_t n)
x=x+alpha*y
Definition: SGVector.cpp:601
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:332
uint32_t maxCPs
Definition: libbmrm.h:53
Class CStructuredModel that represents the application specific model and contains most of the applic...
#define LIBBMRM_INDEX(ROW, COL, NUM_ROWS)
Definition: libbmrm.h:29
#define SG_UNREF(x)
Definition: SGObject.h:52
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
#define SG_SDEBUG(...)
Definition: SGIO.h:168
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:29
uint32_t * ICPcounter
Definition: libbmrm.h:56
The class Features is the base class of all feature objects.
Definition: Features.h:68
#define SG_SERROR(...)
Definition: SGIO.h:179
float64_t * address
Definition: libbmrm.h:44
BmrmStatistics svm_p3bm_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, uint32_t cp_models, bool verbose)
Definition: libp3bm.cpp:33
CStructuredModel * get_model() const
void resize_vector(int32_t n)
Definition: SGVector.cpp:328
static float32_t sqrt(float32_t x)
x^0.5
Definition: Math.h:401
#define LIBBMRM_PLUS_INF
Definition: libbmrm.h:23
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:92
uint32_t BufSize
Definition: libbmrm.cpp:27

SHOGUN Machine Learning Toolbox - Documentation