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