SHOGUN  4.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
malsar_clustered.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  * Written (W) 2012 Sergey Lisitsyn
8  * Copyright (C) 2012 Jiayu Zhou and Jieping Ye
9  */
10 
12 #ifdef HAVE_EIGEN3
13 #ifndef HAVE_CXX11
16 #include <iostream>
17 #include <shogun/lib/external/libqp.h>
18 
19 using namespace Eigen;
20 using namespace std;
21 
22 namespace shogun
23 {
24 
25 static double* H_diag_matrix;
26 static int H_diag_matrix_ld;
27 
28 static const double* get_col(uint32_t j)
29 {
30  return H_diag_matrix + j*H_diag_matrix_ld;
31 }
32 
33 malsar_result_t malsar_clustered(
34  CDotFeatures* features,
35  double* y,
36  double rho1,
37  double rho2,
38  const malsar_options& options)
39 {
40  int task;
41  int n_feats = features->get_dim_feature_space();
42  SG_SDEBUG("n feats = %d\n", n_feats)
43  int n_vecs = features->get_num_vectors();
44  SG_SDEBUG("n vecs = %d\n", n_vecs)
45  int n_tasks = options.n_tasks;
46  SG_SDEBUG("n tasks = %d\n", n_tasks)
47 
48  H_diag_matrix = SG_CALLOC(double, n_tasks*n_tasks);
49  for (int i=0; i<n_tasks; i++)
50  H_diag_matrix[i*n_tasks+i] = 2.0;
51  H_diag_matrix_ld = n_tasks;
52 
53  int iter = 0;
54 
55  // initialize weight vector and bias for each task
56  MatrixXd Ws = MatrixXd::Zero(n_feats, n_tasks);
57  VectorXd Cs = VectorXd::Zero(n_tasks);
58  MatrixXd Ms = MatrixXd::Identity(n_tasks, n_tasks)*options.n_clusters/n_tasks;
59  MatrixXd IM = Ms;
60  MatrixXd IMsqinv = Ms;
61  MatrixXd invEtaMWt = Ms;
62 
63  MatrixXd Wz=Ws, Wzp=Ws, Wz_old=Ws, delta_Wzp=Ws, gWs=Ws;
64  VectorXd Cz=Cs, Czp=Cs, Cz_old=Cs, delta_Czp=Cs, gCs=Cs;
65  MatrixXd Mz=Ms, Mzp=Ms, Mz_old=Ms, delta_Mzp=Ms, gMs=Ms;
66  MatrixXd Mzp_Pz;
67 
68  double eta = rho2/rho1;
69  double c = rho1*eta*(1+eta);
70 
71  double t=1, t_old=0;
72  double gamma=1, gamma_inc=2;
73  double obj=0.0, obj_old=0.0;
74 
75  double* diag_H = SG_MALLOC(double, n_tasks);
76  double* f = SG_MALLOC(double, n_tasks);
77  double* a = SG_MALLOC(double, n_tasks);
78  double* lb = SG_MALLOC(double, n_tasks);
79  double* ub = SG_MALLOC(double, n_tasks);
80  double* x = SG_CALLOC(double, n_tasks);
81 
82  //internal::set_is_malloc_allowed(false);
83  bool done = false;
84  while (!done && iter <= options.max_iter)
85  {
86  double alpha = double(t_old - 1)/t;
87  SG_SDEBUG("alpha=%f\n",alpha)
88 
89  // compute search point
90  Ws = (1+alpha)*Wz - alpha*Wz_old;
91  Cs = (1+alpha)*Cz - alpha*Cz_old;
92  Ms = (1+alpha)*Mz - alpha*Mz_old;
93 
94  // zero gradient
95  gWs.setZero();
96  gCs.setZero();
97  //internal::set_is_malloc_allowed(true);
98  SG_SDEBUG("Computing gradient\n")
99  IM = (eta*MatrixXd::Identity(n_tasks,n_tasks)+Ms);
100  IMsqinv = (IM*IM).inverse();
101  invEtaMWt = IM.inverse()*Ws.transpose();
102  gMs.noalias() = -c*(Ws.transpose()*Ws)*IMsqinv;
103  gWs.noalias() += 2*c*invEtaMWt.transpose();
104  //internal::set_is_malloc_allowed(false);
105 
106  // compute gradient and objective at search point
107  double Fs = 0;
108  for (task=0; task<n_tasks; task++)
109  {
110  SGVector<index_t> task_idx = options.tasks_indices[task];
111  int n_vecs_task = task_idx.vlen;
112  for (int i=0; i<n_vecs_task; i++)
113  {
114  double aa = -y[task_idx[i]]*(features->dense_dot(task_idx[i], Ws.col(task).data(), n_feats)+Cs[task]);
115  double bb = CMath::max(aa,0.0);
116 
117  // avoid underflow when computing exponential loss
118  Fs += (CMath::log(CMath::exp(-bb) + CMath::exp(aa-bb)) + bb)/n_vecs_task;
119  double b = -y[task_idx[i]]*(1 - 1/(1+CMath::exp(aa)))/n_vecs_task;
120  gCs[task] += b;
121  features->add_to_dense_vec(b, task_idx[i], gWs.col(task).data(), n_feats);
122  }
123  }
124  SG_SDEBUG("Computed gradient\n")
125 
126  // add regularizer
127  Fs += c*(Ws*invEtaMWt).trace();
128  SG_SDEBUG("Fs = %f \n", Fs)
129 
130  double Fzp = 0.0;
131 
132  int inner_iter = 0;
133  // line search, Armijo-Goldstein scheme
134  while (inner_iter <= 1)
135  {
136  Wzp = Ws - gWs/gamma;
137  Czp = Cs - gCs/gamma;
138  // compute singular projection of Ms - gMs/gamma with k
139  //internal::set_is_malloc_allowed(true);
140  EigenSolver<MatrixXd> eigensolver;
141  eigensolver.compute(Ms-gMs/gamma, true);
142  if (eigensolver.info()!=Eigen::Success)
143  SG_SERROR("Eigendecomposition failed")
144 
145  // solve problem
146  // min sum_i (s_i - s*_i)^2 s.t. sum_i s_i = k, 0<=s_i<=1
147  for (int i=0; i<n_tasks; i++)
148  {
149  diag_H[i] = 2.0;
150  // TODO fails with C++11
151  //std::complex<MatrixXd::Scalar> eigenvalue = eigensolver.eigenvalues()[i];
152  //cout << "eigenvalue " << eigenvalue << "=" << std::real(eigenvalue) << "+i" << std::imag(eigenvalue) << endl;
153  f[i] = -2*eigensolver.eigenvalues()[i].real();
154  if (f[i]!=f[i])
155  SG_SERROR("NaN %d eigenvalue", i)
156 
157  SG_SDEBUG("%dth eigenvalue %f\n",i,eigensolver.eigenvalues()[i].real())
158  a[i] = 1.0;
159  lb[i] = 0.0;
160  ub[i] = 1.0;
161  x[i] = double(options.n_clusters)/n_tasks;
162  }
163  double b = options.n_clusters;//eigensolver.eigenvalues().sum().real();
164  SG_SDEBUG("b = %f\n", b)
165  SG_SDEBUG("Calling libqp\n")
166  libqp_state_T problem_state = libqp_gsmo_solver(&get_col,diag_H,f,a,b,lb,ub,x,n_tasks,1000,1e-6,NULL);
167  SG_SDEBUG("Libqp objective = %f\n",problem_state.QP)
168  SG_SDEBUG("Exit code = %d\n",problem_state.exitflag)
169  SG_SDEBUG("%d iteration passed\n",problem_state.nIter)
170  SG_SDEBUG("Solution is \n [ ")
171  for (int i=0; i<n_tasks; i++)
172  SG_SDEBUG("%f ", x[i])
173  SG_SDEBUG("]\n")
174  Map<VectorXd> Mzp_DiagSigz(x,n_tasks);
175  Mzp_Pz = eigensolver.eigenvectors().real();
176  Mzp = Mzp_Pz*Mzp_DiagSigz.asDiagonal()*Mzp_Pz.transpose();
177  //internal::set_is_malloc_allowed(false);
178  // walk in direction of antigradient
179  for (int i=0; i<n_tasks; i++)
180  Mzp_DiagSigz[i] += eta;
181  //internal::set_is_malloc_allowed(true);
182  invEtaMWt = (Mzp_Pz*
183  (Mzp_DiagSigz.cwiseInverse().asDiagonal())*
184  Mzp_Pz.transpose())*
185  Wzp.transpose();
186  //internal::set_is_malloc_allowed(false);
187  // compute objective at line search point
188  Fzp = 0.0;
189  for (task=0; task<n_tasks; task++)
190  {
191  SGVector<index_t> task_idx = options.tasks_indices[task];
192  int n_vecs_task = task_idx.vlen;
193  for (int i=0; i<n_vecs_task; i++)
194  {
195  double aa = -y[task_idx[i]]*(features->dense_dot(task_idx[i], Wzp.col(task).data(), n_feats)+Cs[task]);
196  double bb = CMath::max(aa,0.0);
197 
198  Fzp += (CMath::log(CMath::exp(-bb) + CMath::exp(aa-bb)) + bb)/n_vecs_task;
199  }
200  }
201  Fzp += c*(Wzp*invEtaMWt).trace();
202 
203  // compute delta between line search point and search point
204  delta_Wzp = Wzp - Ws;
205  delta_Czp = Czp - Cs;
206  delta_Mzp = Mzp - Ms;
207 
208  // norms of delta
209  double nrm_delta_Wzp = delta_Wzp.squaredNorm();
210  double nrm_delta_Czp = delta_Czp.squaredNorm();
211  double nrm_delta_Mzp = delta_Mzp.squaredNorm();
212 
213  double r_sum = (nrm_delta_Wzp + nrm_delta_Czp + nrm_delta_Mzp)/3;
214 
215  double Fzp_gamma = 0.0;
216  if (n_feats > n_tasks)
217  {
218  Fzp_gamma = Fs + (delta_Wzp.transpose()*gWs).trace() +
219  (delta_Czp.transpose()*gCs).trace() +
220  (delta_Mzp.transpose()*gMs).trace() +
221  (gamma/2)*nrm_delta_Wzp +
222  (gamma/2)*nrm_delta_Czp +
223  (gamma/2)*nrm_delta_Mzp;
224  }
225  else
226  {
227  Fzp_gamma = Fs + (gWs.transpose()*delta_Wzp).trace() +
228  (gCs.transpose()*delta_Czp).trace() +
229  (gMs.transpose()*delta_Mzp).trace() +
230  (gamma/2)*nrm_delta_Wzp +
231  (gamma/2)*nrm_delta_Czp +
232  (gamma/2)*nrm_delta_Mzp;
233  }
234 
235  // break if delta is getting too small
236  if (r_sum <= 1e-20)
237  {
238  done = true;
239  break;
240  }
241 
242  // break if objective at line searc point is smaller than Fzp_gamma
243  if (Fzp <= Fzp_gamma)
244  break;
245  else
246  gamma *= gamma_inc;
247 
248  inner_iter++;
249  }
250 
251  Wz_old = Wz;
252  Cz_old = Cz;
253  Mz_old = Mz;
254  Wz = Wzp;
255  Cz = Czp;
256  Mz = Mzp;
257 
258  // compute objective value
259  obj_old = obj;
260  obj = Fzp;
261 
262  // check if process should be terminated
263  switch (options.termination)
264  {
265  case 0:
266  if (iter>=2)
267  {
268  if ( CMath::abs(obj-obj_old) <= options.tolerance )
269  done = true;
270  }
271  break;
272  case 1:
273  if (iter>=2)
274  {
275  if ( CMath::abs(obj-obj_old) <= options.tolerance*CMath::abs(obj_old))
276  done = true;
277  }
278  break;
279  case 2:
280  if (CMath::abs(obj) <= options.tolerance)
281  done = true;
282  break;
283  case 3:
284  if (iter>=options.max_iter)
285  done = true;
286  break;
287  }
288 
289  iter++;
290  t_old = t;
291  t = 0.5 * (1 + CMath::sqrt(1.0 + 4*t*t));
292  }
293  //internal::set_is_malloc_allowed(true);
294  SG_SDEBUG("%d iteration passed, objective = %f\n",iter,obj)
295 
296  SG_FREE(H_diag_matrix);
297  SG_FREE(diag_H);
298  SG_FREE(f);
299  SG_FREE(a);
300  SG_FREE(lb);
301  SG_FREE(ub);
302  SG_FREE(x);
303 
304  SGMatrix<float64_t> tasks_w(n_feats, n_tasks);
305  for (int i=0; i<n_feats; i++)
306  {
307  for (task=0; task<n_tasks; task++)
308  tasks_w(i,task) = Wzp(i,task);
309  }
310  //tasks_w.display_matrix();
311  SGVector<float64_t> tasks_c(n_tasks);
312  for (int i=0; i<n_tasks; i++) tasks_c[i] = Czp[i];
313  return malsar_result_t(tasks_w, tasks_c);
314 };
315 };
316 #endif
317 #endif
static const double * get_col(uint32_t j)
virtual float64_t dense_dot(int32_t vec_idx1, const float64_t *vec2, int32_t vec2_len)=0
Definition: SGMatrix.h:20
virtual int32_t get_num_vectors() const =0
virtual void add_to_dense_vec(float64_t alpha, int32_t vec_idx1, float64_t *vec2, int32_t vec2_len, bool abs_val=false)=0
static double * H_diag_matrix
Features that support dot products among other operations.
Definition: DotFeatures.h:44
virtual int32_t get_dim_feature_space() const =0
static int H_diag_matrix_ld
index_t vlen
Definition: SGVector.h:494
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
#define SG_SDEBUG(...)
Definition: SGIO.h:168
malsar_result_t malsar_clustered(CDotFeatures *features, double *y, double rho1, double rho2, const malsar_options &options)
#define SG_SERROR(...)
Definition: SGIO.h:179
Matrix::Scalar max(Matrix m)
Definition: Redux.h:66
T * data() const
Definition: SGVector.h:118

SHOGUN Machine Learning Toolbox - Documentation