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