SHOGUN  v2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HessianLocallyLinearEmbedding.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) 2011 Sergey Lisitsyn
8  * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society
9  */
10 
12 #ifdef HAVE_LAPACK
14 #include <shogun/lib/common.h>
16 #include <shogun/io/SGIO.h>
18 #include <shogun/lib/Signal.h>
19 
20 #ifdef HAVE_PTHREAD
21 #include <pthread.h>
22 #endif
23 
24 using namespace shogun;
25 
26 #ifndef DOXYGEN_SHOULD_SKIP_THIS
27 struct HESSIANESTIMATION_THREAD_PARAM
28 {
30  int32_t idx_start;
32  int32_t idx_step;
34  int32_t idx_stop;
36  int32_t m_k;
38  int32_t target_dim;
40  const int32_t* neighborhood_matrix;
42  const float64_t* feature_matrix;
44  float64_t* local_feature_matrix;
46  float64_t* Yi_matrix;
48  float64_t* mean_vector;
50  float64_t* s_values_vector;
52  float64_t* tau;
54  int32_t tau_len;
56  float64_t* w_sum_vector;
58  float64_t* q_matrix;
60  float64_t* W_matrix;
62  int32_t N;
64  int32_t dim;
66  int32_t actual_k;
67 #ifdef HAVE_PTHREAD
68 
69  PTHREAD_LOCK_T* W_matrix_lock;
70 #endif
71 };
72 #endif
73 
76 {
77 }
78 
80 {
81 }
82 
84 {
85  return "HessianLocallyLinearEmbedding";
86 };
87 
89  SGMatrix<int32_t> neighborhood_matrix)
90 {
91  int32_t N = simple_features->get_num_vectors();
92  int32_t dim = simple_features->get_num_features();
93  int32_t dp = m_target_dim*(m_target_dim+1)/2;
94  if (m_k<(1+m_target_dim+dp))
95  SG_ERROR("K parameter should have value greater than 1+target dimensionality+dp.\n");
96 #ifdef HAVE_PTHREAD
97  int32_t t;
98  int32_t num_threads = parallel->get_num_threads();
99  ASSERT(num_threads>0);
100  // allocate threads and params
101  pthread_t* threads = SG_MALLOC(pthread_t, num_threads);
102  HESSIANESTIMATION_THREAD_PARAM* parameters = SG_MALLOC(HESSIANESTIMATION_THREAD_PARAM, num_threads);
103 
104 #else
105  int32_t num_threads = 1;
106 #endif
107 
108  // init matrices to be used
109  float64_t* local_feature_matrix = SG_MALLOC(float64_t, m_k*dim*num_threads);
110  float64_t* s_values_vector = SG_MALLOC(float64_t, dim*num_threads);
111  int32_t tau_len = CMath::min((1+m_target_dim+dp), m_k);
112  float64_t* tau = SG_MALLOC(float64_t, tau_len*num_threads);
113  float64_t* mean_vector = SG_MALLOC(float64_t, dim*num_threads);
114  float64_t* q_matrix = SG_MALLOC(float64_t, m_k*m_k*num_threads);
115  float64_t* w_sum_vector = SG_MALLOC(float64_t, dp*num_threads);
116  float64_t* Yi_matrix = SG_MALLOC(float64_t, m_k*(1+m_target_dim+dp)*num_threads);
117  // get feature matrix
118  SGMatrix<float64_t> feature_matrix = simple_features->get_feature_matrix();
119 
120 #ifdef HAVE_PTHREAD
121  PTHREAD_LOCK_T W_matrix_lock;
122  pthread_attr_t attr;
123  PTHREAD_LOCK_INIT(&W_matrix_lock);
124  pthread_attr_init(&attr);
125  pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
126 
127  for (t=0; t<num_threads; t++)
128  {
129  parameters[t].idx_start = t;
130  parameters[t].idx_step = num_threads;
131  parameters[t].idx_stop = N;
132  parameters[t].m_k = m_k;
133  parameters[t].target_dim = m_target_dim;
134  parameters[t].neighborhood_matrix = neighborhood_matrix.matrix;
135  parameters[t].feature_matrix = feature_matrix.matrix;
136  parameters[t].local_feature_matrix = local_feature_matrix + (m_k*dim)*t;
137  parameters[t].Yi_matrix = Yi_matrix + (m_k*(1+m_target_dim+dp))*t;
138  parameters[t].mean_vector = mean_vector + dim*t;
139  parameters[t].s_values_vector = s_values_vector + dim*t;
140  parameters[t].tau = tau+tau_len*t;
141  parameters[t].tau_len = tau_len;
142  parameters[t].w_sum_vector = w_sum_vector + dp*t;
143  parameters[t].q_matrix = q_matrix + (m_k*m_k)*t;
144  parameters[t].W_matrix = W_matrix;
145  parameters[t].W_matrix_lock = &W_matrix_lock;
146  parameters[t].N = N;
147  parameters[t].dim = dim;
148  parameters[t].actual_k = neighborhood_matrix.num_rows;
149  pthread_create(&threads[t], &attr, run_hessianestimation_thread, (void*)&parameters[t]);
150  }
151  for (t=0; t<num_threads; t++)
152  pthread_join(threads[t], NULL);
153  PTHREAD_LOCK_DESTROY(&W_matrix_lock);
154  SG_FREE(parameters);
155  SG_FREE(threads);
156 #else
157  HESSIANESTIMATION_THREAD_PARAM single_thread_param;
158  single_thread_param.idx_start = 0;
159  single_thread_param.idx_step = num_threads;
160  single_thread_param.idx_stop = N;
161  single_thread_param.m_k = m_k;
162  single_thread_param.neighborhood_matrix = neighborhood_matrix.matrix;
163  single_thread_param.feature_matrix = feature_matrix.matrix;
164  single_thread_param.local_feature_matrix = local_feature_matrix;
165  single_thread_param.Yi_matrix = Yi_matrix;
166  single_thread_param.mean_vector = mean_vector;
167  single_thread_param.s_values_vector = s_values_vector;
168  single_thread_param.tau = tau;
169  single_thread_param.tau_len = tau_len;
170  single_thread_param.w_sum_vector = w_sum_vector;
171  single_thread_param.q_matrix = q_matrix;
172  single_thread_param.W_matrix = W_matrix;
173  single_thread_param.N = N;
174  single_thread_param.dim = dim;
175  single_thread_param.actual_k = neighborhood_matrix.num_rows;
176  run_hessianestimation_thread((void*)&single_thread_param);
177 #endif
178 
179  // clean
180  SG_FREE(Yi_matrix);
181  SG_FREE(s_values_vector);
182  SG_FREE(mean_vector);
183  SG_FREE(tau);
184  SG_FREE(w_sum_vector);
185  SG_FREE(local_feature_matrix);
186  SG_FREE(q_matrix);
187 
188  return SGMatrix<float64_t>(W_matrix,N,N);
189 }
190 
192 {
193  HESSIANESTIMATION_THREAD_PARAM* parameters = (HESSIANESTIMATION_THREAD_PARAM*)p;
194  int32_t idx_start = parameters->idx_start;
195  int32_t idx_step = parameters->idx_step;
196  int32_t idx_stop = parameters->idx_stop;
197  int32_t m_k = parameters->m_k;
198  int32_t target_dim = parameters->target_dim;
199  const int32_t* neighborhood_matrix = parameters->neighborhood_matrix;
200  const float64_t* feature_matrix = parameters->feature_matrix;
201  float64_t* local_feature_matrix = parameters->local_feature_matrix;
202  float64_t* Yi_matrix = parameters->Yi_matrix;
203  float64_t* mean_vector = parameters->mean_vector;
204  float64_t* s_values_vector = parameters->s_values_vector;
205  float64_t* tau = parameters->tau;
206  int32_t tau_len = parameters->tau_len;
207  float64_t* w_sum_vector = parameters->w_sum_vector;
208  float64_t* q_matrix = parameters->q_matrix;
209  float64_t* W_matrix = parameters->W_matrix;
210 #ifdef HAVE_PTHREAD
211  PTHREAD_LOCK_T* W_matrix_lock = parameters->W_matrix_lock;
212 #endif
213 
214  int i,j,k,l;
215  int32_t dp = target_dim*(target_dim+1)/2;
216  int32_t actual_k = parameters->actual_k;
217  int32_t N = parameters->N;
218  int32_t dim = parameters->dim;
219 
220  for (i=idx_start; i<idx_stop; i+=idx_step)
221  {
222  // Yi(:,0) = 1
223  for (j=0; j<m_k; j++)
224  Yi_matrix[j] = 1.0;
225 
226  // fill mean vector with zeros
227  memset(mean_vector,0,sizeof(float64_t)*dim);
228 
229  // compute local feature matrix containing neighbors of i-th vector
230  for (j=0; j<m_k; j++)
231  {
232  for (k=0; k<dim; k++)
233  local_feature_matrix[j*dim+k] = feature_matrix[neighborhood_matrix[i*actual_k+j]*dim+k];
234 
235  cblas_daxpy(dim,1.0,local_feature_matrix+j*dim,1,mean_vector,1);
236  }
237 
238  // compute mean
239  cblas_dscal(dim,1.0/m_k,mean_vector,1);
240 
241  // center feature vectors by mean
242  for (j=0; j<m_k; j++)
243  cblas_daxpy(dim,-1.0,mean_vector,1,local_feature_matrix+j*dim,1);
244 
245  int32_t info = 0;
246  // find right eigenvectors of local_feature_matrix
247  wrap_dgesvd('N','O', dim,m_k,local_feature_matrix,dim,
248  s_values_vector,
249  NULL,1, NULL,1, &info);
250  ASSERT(info==0);
251 
252  // Yi(0:m_k,1:1+target_dim) = Vh(0:m_k, 0:target_dim)
253  for (j=0; j<target_dim; j++)
254  {
255  for (k=0; k<m_k; k++)
256  Yi_matrix[(j+1)*m_k+k] = local_feature_matrix[k*dim+j];
257  }
258 
259  int32_t ct = 0;
260 
261  // construct 2nd order hessian approx
262  for (j=0; j<target_dim; j++)
263  {
264  for (k=0; k<target_dim-j; k++)
265  {
266  for (l=0; l<m_k; l++)
267  {
268  Yi_matrix[(ct+k+1+target_dim)*m_k+l] = Yi_matrix[(j+1)*m_k+l]*Yi_matrix[(j+k+1)*m_k+l];
269  }
270  }
271  ct += ct + target_dim - j;
272  }
273 
274  // perform QR factorization
275  wrap_dgeqrf(m_k,(1+target_dim+dp),Yi_matrix,m_k,tau,&info);
276  ASSERT(info==0);
277  wrap_dorgqr(m_k,(1+target_dim+dp),tau_len,Yi_matrix,m_k,tau,&info);
278  ASSERT(info==0);
279 
280  float64_t* Pii = (Yi_matrix+m_k*(1+target_dim));
281 
282  for (j=0; j<dp; j++)
283  {
284  w_sum_vector[j] = 0.0;
285  for (k=0; k<m_k; k++)
286  {
287  w_sum_vector[j] += Pii[j*m_k+k];
288  }
289  if (w_sum_vector[j]<0.001)
290  w_sum_vector[j] = 1.0;
291  for (k=0; k<m_k; k++)
292  Pii[j*m_k+k] /= w_sum_vector[j];
293  }
294 
295  cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,
296  m_k,m_k,dp,
297  1.0,Pii,m_k,
298  Pii,m_k,
299  0.0,q_matrix,m_k);
300 #ifdef HAVE_PTHREAD
301  PTHREAD_LOCK(W_matrix_lock);
302 #endif
303  for (j=0; j<m_k; j++)
304  {
305  for (k=0; k<m_k; k++)
306  W_matrix[N*neighborhood_matrix[i*actual_k+k]+neighborhood_matrix[i*actual_k+j]] += q_matrix[j*m_k+k];
307  }
308 #ifdef HAVE_PTHREAD
309  PTHREAD_UNLOCK(W_matrix_lock);
310 #endif
311  }
312  return NULL;
313 }
314 #endif /* HAVE_LAPACK */

SHOGUN Machine Learning Toolbox - Documentation