SHOGUN  v1.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
LocallyLinearEmbedding.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 #include <shogun/lib/config.h>
13 #ifdef HAVE_LAPACK
18 #include <shogun/base/DynArray.h>
20 #include <shogun/io/SGIO.h>
21 #include <shogun/lib/Time.h>
23 #include <shogun/lib/Signal.h>
24 
25 #ifdef HAVE_PTHREAD
26 #include <pthread.h>
27 #endif
28 
29 using namespace shogun;
30 
31 #ifndef DOXYGEN_SHOULD_SKIP_THIS
32 struct LINRECONSTRUCTION_THREAD_PARAM
33 {
35  int32_t idx_start;
37  int32_t idx_stop;
39  int32_t idx_step;
41  int32_t m_k;
43  int32_t dim;
45  int32_t N;
47  const int32_t* neighborhood_matrix;
49  const float64_t* feature_matrix;
51  float64_t* z_matrix;
53  float64_t* covariance_matrix;
55  float64_t* id_vector;
57  float64_t* W_matrix;
59  float64_t m_reconstruction_shift;
60 };
61 
62 struct NEIGHBORHOOD_THREAD_PARAM
63 {
65  int32_t idx_start;
67  int32_t idx_step;
69  int32_t idx_stop;
71  int32_t N;
73  int32_t m_k;
75  CFibonacciHeap* heap;
77  const float64_t* distance_matrix;
79  int32_t* neighborhood_matrix;
80 };
81 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
82 
85 {
86  m_k = 10;
87  m_max_k = 60;
88  m_auto_k = false;
89  m_nullspace_shift = -1e-9;
91 #ifdef HAVE_ARPACK
92  m_use_arpack = true;
93 #else
94  m_use_arpack = false;
95 #endif
96  init();
97 }
98 
100 {
101  m_parameters->add(&m_auto_k, "auto_k",
102  "whether k should be determined automatically in range");
103  m_parameters->add(&m_k, "k", "number of neighbors");
104  m_parameters->add(&m_max_k, "max_k",
105  "maximum number of neighbors used to compute optimal one");
106  m_parameters->add(&m_nullspace_shift, "nullspace_shift",
107  "nullspace finding regularization shift");
108  m_parameters->add(&m_reconstruction_shift, "reconstruction_shift",
109  "shift used to regularize reconstruction step");
110  m_parameters->add(&m_use_arpack, "use_arpack",
111  "whether arpack is being used or not");
112 }
113 
114 
116 {
117 }
118 
120 {
121  ASSERT(k>0);
122  m_k = k;
123 }
124 
126 {
127  return m_k;
128 }
129 
131 {
132  ASSERT(max_k>=m_k);
133  m_max_k = max_k;
134 }
135 
137 {
138  return m_max_k;
139 }
140 
142 {
143  m_auto_k = auto_k;
144 }
145 
147 {
148  return m_auto_k;
149 }
150 
152 {
153  m_nullspace_shift = nullspace_shift;
154 }
155 
157 {
158  return m_nullspace_shift;
159 }
160 
162 {
163  m_reconstruction_shift = reconstruction_shift;
164 }
165 
167 {
168  return m_reconstruction_shift;
169 }
170 
172 {
173  m_use_arpack = use_arpack;
174 }
175 
177 {
178  return m_use_arpack;
179 }
180 
182 {
183  return "LocallyLinearEmbedding";
184 }
185 
187 {
188  ASSERT(features);
189  // check features
190  if (!(features->get_feature_class()==C_SIMPLE &&
191  features->get_feature_type()==F_DREAL))
192  {
193  SG_ERROR("Given features are not of SimpleRealFeatures type.\n");
194  }
195  // shorthand for simplefeatures
196  CSimpleFeatures<float64_t>* simple_features = (CSimpleFeatures<float64_t>*) features;
197  SG_REF(features);
198 
199  // get and check number of vectors
200  int32_t N = simple_features->get_num_vectors();
201  if (m_k>=N)
202  SG_ERROR("Number of neighbors (%d) should be less than number of objects (%d).\n",
203  m_k, N);
204 
205  // compute distance matrix
206  SG_DEBUG("Computing distance matrix\n");
208  m_distance->init(simple_features,simple_features);
211  SG_DEBUG("Calculating neighborhood matrix\n");
212  SGMatrix<int32_t> neighborhood_matrix;
213 
214  if (m_auto_k)
215  {
216  neighborhood_matrix = get_neighborhood_matrix(distance_matrix,m_max_k);
217  m_k = estimate_k(simple_features,neighborhood_matrix);
218  SG_DEBUG("Estimated k with value of %d\n",m_k);
219  }
220  else
221  neighborhood_matrix = get_neighborhood_matrix(distance_matrix,m_k);
222 
223  // init W (weight) matrix
224  float64_t* W_matrix = distance_matrix.matrix;
225  memset(W_matrix,0,sizeof(float64_t)*N*N);
226 
227  // construct weight matrix
228  SG_DEBUG("Constructing weight matrix\n");
229  SGMatrix<float64_t> weight_matrix = construct_weight_matrix(simple_features,W_matrix,neighborhood_matrix);
230  neighborhood_matrix.destroy_matrix();
231 
232  // find null space of weight matrix
233  SG_DEBUG("Finding nullspace\n");
234  SGMatrix<float64_t> new_feature_matrix = construct_embedding(weight_matrix,m_target_dim);
235  weight_matrix.destroy_matrix();
236 
237  SG_UNREF(features);
238  return (CFeatures*)(new CSimpleFeatures<float64_t>(new_feature_matrix));
239 }
240 
242 {
243  int32_t right = m_max_k;
244  int32_t left = m_k;
245  int32_t left_third;
246  int32_t right_third;
247  ASSERT(right>=left);
248  if (right==left) return left;
249  int32_t dim;
250  int32_t N;
251  float64_t* feature_matrix= simple_features->get_feature_matrix(dim,N);
252  float64_t* z_matrix = SG_MALLOC(float64_t,right*dim);
253  float64_t* covariance_matrix = SG_MALLOC(float64_t,right*right);
254  float64_t* resid_vector = SG_MALLOC(float64_t, right);
255  float64_t* id_vector = SG_MALLOC(float64_t, right);
256  while (right-left>2)
257  {
258  left_third = (left*2+right)/3;
259  right_third = (right*2+left)/3;
260  float64_t left_val = compute_reconstruction_error(left_third,dim,N,feature_matrix,z_matrix,
261  covariance_matrix,resid_vector,
262  id_vector,neighborhood_matrix);
263  float64_t right_val = compute_reconstruction_error(right_third,dim,N,feature_matrix,z_matrix,
264  covariance_matrix,resid_vector,
265  id_vector,neighborhood_matrix);
266  if (left_val<right_val)
267  right = right_third;
268  else
269  left = left_third;
270  }
271  SG_FREE(z_matrix);
272  SG_FREE(covariance_matrix);
273  SG_FREE(resid_vector);
274  SG_FREE(id_vector);
275  return right;
276 }
277 
279  float64_t* z_matrix, float64_t* covariance_matrix,
280  float64_t* resid_vector, float64_t* id_vector,
281  SGMatrix<int32_t> neighborhood_matrix)
282 {
283  // todo parse params
284  int32_t i,j;
285  float64_t total_residual_norm = 0.0;
286  for (i=CMath::random(0,20); i<N; i+=N/20)
287  {
288  for (j=0; j<k; j++)
289  {
290  cblas_dcopy(dim,feature_matrix+neighborhood_matrix[j*N+i]*dim,1,z_matrix+j*dim,1);
291  cblas_daxpy(dim,-1.0,feature_matrix+i*dim,1,z_matrix+j*dim,1);
292  }
293  cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans,
294  k,k,dim,
295  1.0,z_matrix,dim,
296  z_matrix,dim,
297  0.0,covariance_matrix,k);
298  for (j=0; j<k; j++)
299  {
300  resid_vector[j] = 1.0;
301  id_vector[j] = 1.0;
302  }
303  if (k>dim)
304  {
305  float64_t trace = 0.0;
306  for (j=0; j<k; j++)
307  trace += covariance_matrix[j*k+j];
308  for (j=0; j<m_k; j++)
309  covariance_matrix[j*k+j] += m_reconstruction_shift*trace;
310  }
311  clapack_dposv(CblasColMajor,CblasLower,k,1,covariance_matrix,k,id_vector,k);
312  float64_t norming=0.0;
313  for (j=0; j<k; j++)
314  norming += id_vector[j];
315  cblas_dscal(k,-1.0/norming,id_vector,1);
316  cblas_dsymv(CblasColMajor,CblasLower,k,-1.0,covariance_matrix,k,id_vector,1,1.0,resid_vector,1);
317  total_residual_norm += cblas_dnrm2(k,resid_vector,1);
318  }
319  return total_residual_norm/k;
320 }
321 
323  float64_t* W_matrix, SGMatrix<int32_t> neighborhood_matrix)
324 {
325  int32_t N = simple_features->get_num_vectors();
326  int32_t dim = simple_features->get_num_features();
327  int32_t t;
328 #ifdef HAVE_PTHREAD
329  int32_t num_threads = parallel->get_num_threads();
330  ASSERT(num_threads>0);
331  // allocate threads
332  pthread_t* threads = SG_MALLOC(pthread_t, num_threads);
333  LINRECONSTRUCTION_THREAD_PARAM* parameters = SG_MALLOC(LINRECONSTRUCTION_THREAD_PARAM, num_threads);
334  pthread_attr_t attr;
335  pthread_attr_init(&attr);
336  pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
337 #else
338  int32_t num_threads = 1;
339 #endif
340  // init storages to be used
341  float64_t* z_matrix = SG_MALLOC(float64_t, m_k*dim*num_threads);
342  float64_t* covariance_matrix = SG_MALLOC(float64_t, m_k*m_k*num_threads);
343  float64_t* id_vector = SG_MALLOC(float64_t, m_k*num_threads);
344 
345  // get feature matrix
346  SGMatrix<float64_t> feature_matrix = simple_features->get_feature_matrix();
347 
348 #ifdef HAVE_PTHREAD
349  for (t=0; t<num_threads; t++)
350  {
351  parameters[t].idx_start = t;
352  parameters[t].idx_step = num_threads;
353  parameters[t].idx_stop = N;
354  parameters[t].m_k = m_k;
355  parameters[t].dim = dim;
356  parameters[t].N = N;
357  parameters[t].neighborhood_matrix = neighborhood_matrix.matrix;
358  parameters[t].z_matrix = z_matrix+(m_k*dim)*t;
359  parameters[t].feature_matrix = feature_matrix.matrix;
360  parameters[t].covariance_matrix = covariance_matrix+(m_k*m_k)*t;
361  parameters[t].id_vector = id_vector+m_k*t;
362  parameters[t].W_matrix = W_matrix;
363  parameters[t].m_reconstruction_shift = m_reconstruction_shift;
364  pthread_create(&threads[t], &attr, run_linearreconstruction_thread, (void*)&parameters[t]);
365  }
366  for (t=0; t<num_threads; t++)
367  pthread_join(threads[t], NULL);
368  pthread_attr_destroy(&attr);
369  SG_FREE(parameters);
370  SG_FREE(threads);
371 #else
372  LINRECONSTRUCTION_THREAD_PARAM single_thread_param;
373  single_thread_param.idx_start = 0;
374  single_thread_param.idx_step = 1;
375  single_thread_param.idx_stop = N;
376  single_thread_param.m_k = m_k;
377  single_thread_param.dim = dim;
378  single_thread_param.N = N;
379  single_thread_param.neighborhood_matrix = neighborhood_matrix.matrix;
380  single_thread_param.z_matrix = z_matrix;
381  single_thread_param.feature_matrix = feature_matrix.matrix;
382  single_thread_param.covariance_matrix = covariance_matrix;
383  single_thread_param.id_vector = id_vector;
384  single_thread_param.W_matrix = W_matrix;
385  single_thread_param.m_reconstruction_shift = m_reconstruction_shift;
386  run_linearreconstruction_thread((void*)single_thread_param);
387 #endif
388 
389  // clean
390  SG_FREE(id_vector);
391  SG_FREE(z_matrix);
392  SG_FREE(covariance_matrix);
393 
394  return SGMatrix<float64_t>(W_matrix,N,N);
395 }
396 
398 {
399  int i,j;
400  ASSERT(matrix.num_cols==matrix.num_rows);
401  int N = matrix.num_cols;
402  // get eigenvectors with ARPACK or LAPACK
403  int eigenproblem_status = 0;
404 
405  float64_t* eigenvalues_vector = NULL;
406  float64_t* eigenvectors = NULL;
407  float64_t* nullspace_features = NULL;
408  if (m_use_arpack)
409  {
410 #ifndef HAVE_ARPACK
411  SG_ERROR("ARPACK is not supported in this configuration.\n");
412 #endif
413  // using ARPACK (faster)
414  eigenvalues_vector = SG_MALLOC(float64_t, dimension+1);
415 #ifdef HAVE_ARPACK
416  arpack_dsxupd(matrix.matrix,NULL,false,N,dimension+1,"LA",true,3,true,m_nullspace_shift,0.0,
417  eigenvalues_vector,matrix.matrix,eigenproblem_status);
418 #endif
419  if (eigenproblem_status)
420  SG_ERROR("ARPACK failed with code: %d", eigenproblem_status);
421  nullspace_features = SG_MALLOC(float64_t, N*dimension);
422  for (i=0; i<dimension; i++)
423  {
424  for (j=0; j<N; j++)
425  nullspace_features[j*dimension+i] = matrix[j*(dimension+1)+i+1];
426  }
427  SG_FREE(eigenvalues_vector);
428  }
429  else
430  {
431  // using LAPACK (slower)
432  eigenvalues_vector = SG_MALLOC(float64_t, N);
433  eigenvectors = SG_MALLOC(float64_t,(dimension+1)*N);
434  wrap_dsyevr('V','U',N,matrix.matrix,N,2,dimension+2,eigenvalues_vector,eigenvectors,&eigenproblem_status);
435  if (eigenproblem_status)
436  SG_ERROR("LAPACK failed with code: %d", eigenproblem_status);
437  nullspace_features = SG_MALLOC(float64_t, N*dimension);
438  // LAPACKed eigenvectors
439  for (i=0; i<dimension; i++)
440  {
441  for (j=0; j<N; j++)
442  nullspace_features[j*dimension+i] = eigenvectors[i*N+j];
443  }
444  SG_FREE(eigenvectors);
445  SG_FREE(eigenvalues_vector);
446  }
447  return SGMatrix<float64_t>(nullspace_features,dimension,N);
448 }
449 
451 {
452  LINRECONSTRUCTION_THREAD_PARAM* parameters = (LINRECONSTRUCTION_THREAD_PARAM*)p;
453  int32_t idx_start = parameters->idx_start;
454  int32_t idx_step = parameters->idx_step;
455  int32_t idx_stop = parameters->idx_stop;
456  int32_t m_k = parameters->m_k;
457  int32_t dim = parameters->dim;
458  int32_t N = parameters->N;
459  const int32_t* neighborhood_matrix = parameters->neighborhood_matrix;
460  float64_t* z_matrix = parameters->z_matrix;
461  const float64_t* feature_matrix = parameters->feature_matrix;
462  float64_t* covariance_matrix = parameters->covariance_matrix;
463  float64_t* id_vector = parameters->id_vector;
464  float64_t* W_matrix = parameters->W_matrix;
465  float64_t m_reconstruction_shift = parameters->m_reconstruction_shift;
466 
467  int32_t i,j,k;
468  float64_t norming,trace;
469 
470  for (i=idx_start; i<idx_stop; i+=idx_step)
471  {
472  // compute local feature matrix containing neighbors of i-th vector
473  // center features by subtracting i-th feature column
474  for (j=0; j<m_k; j++)
475  {
476  cblas_dcopy(dim,feature_matrix+neighborhood_matrix[j*N+i]*dim,1,z_matrix+j*dim,1);
477  cblas_daxpy(dim,-1.0,feature_matrix+i*dim,1,z_matrix+j*dim,1);
478  }
479 
480  // compute local covariance matrix
481  cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans,
482  m_k,m_k,dim,
483  1.0,z_matrix,dim,
484  z_matrix,dim,
485  0.0,covariance_matrix,m_k);
486 
487  for (j=0; j<m_k; j++)
488  id_vector[j] = 1.0;
489 
490  // regularize in case of ill-posed system
491  if (m_k>dim)
492  {
493  // compute tr(C)
494  trace = 0.0;
495  for (j=0; j<m_k; j++)
496  trace += covariance_matrix[j*m_k+j];
497 
498  for (j=0; j<m_k; j++)
499  covariance_matrix[j*m_k+j] += m_reconstruction_shift*trace;
500  }
501 
502  // solve system of linear equations: covariance_matrix * X = 1
503  // covariance_matrix is a pos-def matrix
504  clapack_dposv(CblasColMajor,CblasLower,m_k,1,covariance_matrix,m_k,id_vector,m_k);
505 
506  // normalize weights
507  norming=0.0;
508  for (j=0; j<m_k; j++)
509  norming += id_vector[j];
510 
511  cblas_dscal(m_k,1.0/norming,id_vector,1);
512 
513  memset(covariance_matrix,0,sizeof(float64_t)*m_k*m_k);
514  cblas_dger(CblasColMajor,m_k,m_k,1.0,id_vector,1,id_vector,1,covariance_matrix,m_k);
515 
516  // put weights into W matrix
517  W_matrix[N*i+i] += 1.0;
518  for (j=0; j<m_k; j++)
519  {
520  W_matrix[N*i+neighborhood_matrix[j*N+i]] -= id_vector[j];
521  W_matrix[N*neighborhood_matrix[j*N+i]+i] -= id_vector[j];
522  }
523  for (j=0; j<m_k; j++)
524  {
525  for (k=0; k<m_k; k++)
526  W_matrix[N*neighborhood_matrix[j*N+i]+neighborhood_matrix[k*N+i]]+=covariance_matrix[j*m_k+k];
527  }
528  }
529  return NULL;
530 }
531 
533 {
534  int32_t t;
535  int32_t N = distance_matrix.num_rows;
536  // init matrix and heap to be used
537  int32_t* neighborhood_matrix = SG_MALLOC(int32_t, N*k);
538 #ifdef HAVE_PTHREAD
539  int32_t num_threads = parallel->get_num_threads();
540  ASSERT(num_threads>0);
541  NEIGHBORHOOD_THREAD_PARAM* parameters = SG_MALLOC(NEIGHBORHOOD_THREAD_PARAM, num_threads);
542  pthread_t* threads = SG_MALLOC(pthread_t, num_threads);
543  pthread_attr_t attr;
544  pthread_attr_init(&attr);
545  pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
546 #else
547  int32_t num_threads = 1;
548 #endif
549  CFibonacciHeap** heaps = SG_MALLOC(CFibonacciHeap*, num_threads);
550  for (t=0; t<num_threads; t++)
551  heaps[t] = new CFibonacciHeap(N);
552 
553 #ifdef HAVE_PTHREAD
554  for (t=0; t<num_threads; t++)
555  {
556  parameters[t].idx_start = t;
557  parameters[t].idx_step = num_threads;
558  parameters[t].idx_stop = N;
559  parameters[t].m_k = k;
560  parameters[t].N = N;
561  parameters[t].heap = heaps[t];
562  parameters[t].neighborhood_matrix = neighborhood_matrix;
563  parameters[t].distance_matrix = distance_matrix.matrix;
564  pthread_create(&threads[t], &attr, run_neighborhood_thread, (void*)&parameters[t]);
565  }
566  for (t=0; t<num_threads; t++)
567  pthread_join(threads[t], NULL);
568  pthread_attr_destroy(&attr);
569  SG_FREE(threads);
570  SG_FREE(parameters);
571 #else
572  NEIGHBORHOOD_THREAD_PARAM single_thread_param;
573  single_thread_param.idx_start = 0;
574  single_thread_param.idx_step = 1;
575  single_thread_param.idx_stop = N;
576  single_thread_param.m_k = k;
577  single_thread_param.N = N;
578  single_thread_param.heap = heaps[0]
579  single_thread_param.neighborhood_matrix = neighborhood_matrix;
580  single_thread_param.distance_matrix = distance_matrix.matrix;
581  run_neighborhood_thread((void*)&single_thread_param);
582 #endif
583 
584  for (t=0; t<num_threads; t++)
585  delete heaps[t];
586  SG_FREE(heaps);
587 
588  return SGMatrix<int32_t>(neighborhood_matrix,k,N);
589 }
590 
592 {
593  NEIGHBORHOOD_THREAD_PARAM* parameters = (NEIGHBORHOOD_THREAD_PARAM*)p;
594  int32_t idx_start = parameters->idx_start;
595  int32_t idx_step = parameters->idx_step;
596  int32_t idx_stop = parameters->idx_stop;
597  int32_t N = parameters->N;
598  int32_t m_k = parameters->m_k;
599  CFibonacciHeap* heap = parameters->heap;
600  const float64_t* distance_matrix = parameters->distance_matrix;
601  int32_t* neighborhood_matrix = parameters->neighborhood_matrix;
602 
603  int32_t i,j;
604  float64_t tmp;
605  for (i=idx_start; i<idx_stop; i+=idx_step)
606  {
607  for (j=0; j<N; j++)
608  {
609  heap->insert(j,distance_matrix[i*N+j]);
610  }
611 
612  heap->extract_min(tmp);
613 
614  for (j=0; j<m_k; j++)
615  neighborhood_matrix[j*N+i] = heap->extract_min(tmp);
616 
617  heap->clear();
618  }
619 
620  return NULL;
621 }
622 #endif /* HAVE_LAPACK */

SHOGUN Machine Learning Toolbox - Documentation