SHOGUN
v1.1.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2011 Sergey Lisitsyn 00008 * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society 00009 */ 00010 00011 #include <shogun/converter/HessianLocallyLinearEmbedding.h> 00012 #ifdef HAVE_LAPACK 00013 #include <shogun/mathematics/lapack.h> 00014 #include <shogun/lib/common.h> 00015 #include <shogun/mathematics/Math.h> 00016 #include <shogun/io/SGIO.h> 00017 #include <shogun/distance/Distance.h> 00018 #include <shogun/lib/Signal.h> 00019 00020 #ifdef HAVE_PTHREAD 00021 #include <pthread.h> 00022 #endif 00023 00024 using namespace shogun; 00025 00026 #ifndef DOXYGEN_SHOULD_SKIP_THIS 00027 struct HESSIANESTIMATION_THREAD_PARAM 00028 { 00030 int32_t idx_start; 00032 int32_t idx_step; 00034 int32_t idx_stop; 00036 int32_t m_k; 00038 int32_t dim; 00040 int32_t N; 00042 int32_t dp; 00044 int32_t target_dim; 00046 const int32_t* neighborhood_matrix; 00048 const float64_t* feature_matrix; 00050 float64_t* local_feature_matrix; 00052 float64_t* Yi_matrix; 00054 float64_t* mean_vector; 00056 float64_t* s_values_vector; 00058 float64_t* tau; 00060 int32_t tau_len; 00062 float64_t* w_sum_vector; 00064 float64_t* q_matrix; 00066 float64_t* W_matrix; 00067 #ifdef HAVE_PTHREAD 00068 00069 PTHREAD_LOCK_T* W_matrix_lock; 00070 #endif 00071 }; 00072 #endif 00073 00074 CHessianLocallyLinearEmbedding::CHessianLocallyLinearEmbedding() : 00075 CLocallyLinearEmbedding() 00076 { 00077 } 00078 00079 CHessianLocallyLinearEmbedding::~CHessianLocallyLinearEmbedding() 00080 { 00081 } 00082 00083 const char* CHessianLocallyLinearEmbedding::get_name() const 00084 { 00085 return "HessianLocallyLinearEmbedding"; 00086 }; 00087 00088 SGMatrix<float64_t> CHessianLocallyLinearEmbedding::construct_weight_matrix(CSimpleFeatures<float64_t>* simple_features,float64_t* W_matrix, 00089 SGMatrix<int32_t> neighborhood_matrix) 00090 { 00091 int32_t N = simple_features->get_num_vectors(); 00092 int32_t dim = simple_features->get_num_features(); 00093 int32_t dp = m_target_dim*(m_target_dim+1)/2; 00094 if (m_k<(1+m_target_dim+dp)) 00095 SG_ERROR("K parameter should have value greater than 1+target dimensionality+dp.\n"); 00096 int32_t t; 00097 #ifdef HAVE_PTHREAD 00098 int32_t num_threads = parallel->get_num_threads(); 00099 ASSERT(num_threads>0); 00100 // allocate threads and params 00101 pthread_t* threads = SG_MALLOC(pthread_t, num_threads); 00102 HESSIANESTIMATION_THREAD_PARAM* parameters = SG_MALLOC(HESSIANESTIMATION_THREAD_PARAM, num_threads); 00103 00104 #else 00105 int32_t num_threads = 1; 00106 #endif 00107 00108 // init matrices to be used 00109 float64_t* local_feature_matrix = SG_MALLOC(float64_t, m_k*dim*num_threads); 00110 float64_t* s_values_vector = SG_MALLOC(float64_t, dim*num_threads); 00111 int32_t tau_len = CMath::min((1+m_target_dim+dp), m_k); 00112 float64_t* tau = SG_MALLOC(float64_t, tau_len*num_threads); 00113 float64_t* mean_vector = SG_MALLOC(float64_t, dim*num_threads); 00114 float64_t* q_matrix = SG_MALLOC(float64_t, m_k*m_k*num_threads); 00115 float64_t* w_sum_vector = SG_MALLOC(float64_t, dp*num_threads); 00116 float64_t* Yi_matrix = SG_MALLOC(float64_t, m_k*(1+m_target_dim+dp)*num_threads); 00117 // get feature matrix 00118 SGMatrix<float64_t> feature_matrix = simple_features->get_feature_matrix(); 00119 00120 #ifdef HAVE_PTHREAD 00121 PTHREAD_LOCK_T W_matrix_lock; 00122 pthread_attr_t attr; 00123 PTHREAD_LOCK_INIT(&W_matrix_lock); 00124 pthread_attr_init(&attr); 00125 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE); 00126 00127 for (t=0; t<num_threads; t++) 00128 { 00129 parameters[t].idx_start = t; 00130 parameters[t].idx_step = num_threads; 00131 parameters[t].idx_stop = N; 00132 parameters[t].m_k = m_k; 00133 parameters[t].dim = dim; 00134 parameters[t].target_dim = m_target_dim; 00135 parameters[t].N = N; 00136 parameters[t].dp = dp; 00137 parameters[t].neighborhood_matrix = neighborhood_matrix.matrix; 00138 parameters[t].feature_matrix = feature_matrix.matrix; 00139 parameters[t].local_feature_matrix = local_feature_matrix + (m_k*dim)*t; 00140 parameters[t].Yi_matrix = Yi_matrix + (m_k*(1+m_target_dim+dp))*t; 00141 parameters[t].mean_vector = mean_vector + dim*t; 00142 parameters[t].s_values_vector = s_values_vector + dim*t; 00143 parameters[t].tau = tau+tau_len*t; 00144 parameters[t].tau_len = tau_len; 00145 parameters[t].w_sum_vector = w_sum_vector + dp*t; 00146 parameters[t].q_matrix = q_matrix + (m_k*m_k)*t; 00147 parameters[t].W_matrix = W_matrix; 00148 parameters[t].W_matrix_lock = &W_matrix_lock; 00149 pthread_create(&threads[t], &attr, run_hessianestimation_thread, (void*)¶meters[t]); 00150 } 00151 for (t=0; t<num_threads; t++) 00152 pthread_join(threads[t], NULL); 00153 PTHREAD_LOCK_DESTROY(&W_matrix_lock); 00154 SG_FREE(parameters); 00155 SG_FREE(threads); 00156 #else 00157 HESSIANESTIMATION_THREAD_PARAM single_thread_param; 00158 single_thread_param.idx_start = t; 00159 single_thread_param.idx_step = num_threads; 00160 single_thread_param.idx_stop = N; 00161 single_thread_param.m_k = m_k; 00162 single_thread_param.dim = dim; 00163 single_thread_param.target_dim = m_target_dim; 00164 single_thread_param.N = N; 00165 single_thread_param.dp = dp; 00166 single_thread_param.neighborhood_matrix = neighborhood_matrix.matrix; 00167 single_thread_param.feature_matrix = feature_matrix.matrix; 00168 single_thread_param.local_feature_matrix = local_feature_matrix; 00169 single_thread_param.Yi_matrix = Yi_matrix; 00170 single_thread_param.mean_vector = mean_vector; 00171 single_thread_param.s_values_vector = s_values_vector; 00172 single_thread_param.tau = tau; 00173 single_thread_param.tau_len = tau_len; 00174 single_thread_param.w_sum_vector = w_sum_vector; 00175 single_thread_param.q_matrix = q_matrix; 00176 single_thread_param.W_matrix = W_matrix; 00177 run_hessianestimation_thread((void*)&single_thread_param); 00178 #endif 00179 00180 // clean 00181 SG_FREE(Yi_matrix); 00182 SG_FREE(s_values_vector); 00183 SG_FREE(mean_vector); 00184 SG_FREE(tau); 00185 SG_FREE(w_sum_vector); 00186 SG_FREE(local_feature_matrix); 00187 SG_FREE(q_matrix); 00188 00189 return SGMatrix<float64_t>(W_matrix,N,N); 00190 } 00191 00192 void* CHessianLocallyLinearEmbedding::run_hessianestimation_thread(void* p) 00193 { 00194 HESSIANESTIMATION_THREAD_PARAM* parameters = (HESSIANESTIMATION_THREAD_PARAM*)p; 00195 int32_t idx_start = parameters->idx_start; 00196 int32_t idx_step = parameters->idx_step; 00197 int32_t idx_stop = parameters->idx_stop; 00198 int32_t m_k = parameters->m_k; 00199 int32_t dim = parameters->dim; 00200 int32_t N = parameters->N; 00201 int32_t dp = parameters->dp; 00202 int32_t target_dim = parameters->target_dim; 00203 const int32_t* neighborhood_matrix = parameters->neighborhood_matrix; 00204 const float64_t* feature_matrix = parameters->feature_matrix; 00205 float64_t* local_feature_matrix = parameters->local_feature_matrix; 00206 float64_t* Yi_matrix = parameters->Yi_matrix; 00207 float64_t* mean_vector = parameters->mean_vector; 00208 float64_t* s_values_vector = parameters->s_values_vector; 00209 float64_t* tau = parameters->tau; 00210 int32_t tau_len = parameters->tau_len; 00211 float64_t* w_sum_vector = parameters->w_sum_vector; 00212 float64_t* q_matrix = parameters->q_matrix; 00213 float64_t* W_matrix = parameters->W_matrix; 00214 #ifdef HAVE_PTHREAD 00215 PTHREAD_LOCK_T* W_matrix_lock = parameters->W_matrix_lock; 00216 #endif 00217 00218 int i,j,k,l; 00219 00220 for (i=idx_start; i<idx_stop; i+=idx_step) 00221 { 00222 // Yi(:,0) = 1 00223 for (j=0; j<m_k; j++) 00224 Yi_matrix[j] = 1.0; 00225 00226 // fill mean vector with zeros 00227 memset(mean_vector,0,sizeof(float64_t)*dim); 00228 00229 // compute local feature matrix containing neighbors of i-th vector 00230 for (j=0; j<m_k; j++) 00231 { 00232 for (k=0; k<dim; k++) 00233 local_feature_matrix[j*dim+k] = feature_matrix[neighborhood_matrix[j*N+i]*dim+k]; 00234 00235 cblas_daxpy(dim,1.0,local_feature_matrix+j*dim,1,mean_vector,1); 00236 } 00237 00238 // compute mean 00239 cblas_dscal(dim,1.0/m_k,mean_vector,1); 00240 00241 // center feature vectors by mean 00242 for (j=0; j<m_k; j++) 00243 cblas_daxpy(dim,-1.0,mean_vector,1,local_feature_matrix+j*dim,1); 00244 00245 int32_t info = 0; 00246 // find right eigenvectors of local_feature_matrix 00247 wrap_dgesvd('N','O', dim,m_k,local_feature_matrix,dim, 00248 s_values_vector, 00249 NULL,1, NULL,1, &info); 00250 ASSERT(info==0); 00251 00252 // Yi(0:m_k,1:1+target_dim) = Vh(0:m_k, 0:target_dim) 00253 for (j=0; j<target_dim; j++) 00254 { 00255 for (k=0; k<m_k; k++) 00256 Yi_matrix[(j+1)*m_k+k] = local_feature_matrix[k*dim+j]; 00257 } 00258 00259 int32_t ct = 0; 00260 00261 // construct 2nd order hessian approx 00262 for (j=0; j<target_dim; j++) 00263 { 00264 for (k=0; k<target_dim-j; k++) 00265 { 00266 for (l=0; l<m_k; l++) 00267 { 00268 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]; 00269 } 00270 } 00271 ct += ct + target_dim - j; 00272 } 00273 00274 // perform QR factorization 00275 wrap_dgeqrf(m_k,(1+target_dim+dp),Yi_matrix,m_k,tau,&info); 00276 ASSERT(info==0); 00277 wrap_dorgqr(m_k,(1+target_dim+dp),tau_len,Yi_matrix,m_k,tau,&info); 00278 ASSERT(info==0); 00279 00280 float64_t* Pii = (Yi_matrix+m_k*(1+target_dim)); 00281 00282 for (j=0; j<dp; j++) 00283 { 00284 w_sum_vector[j] = 0.0; 00285 for (k=0; k<m_k; k++) 00286 { 00287 w_sum_vector[j] += Pii[j*m_k+k]; 00288 } 00289 if (w_sum_vector[j]<0.001) 00290 w_sum_vector[j] = 1.0; 00291 for (k=0; k<m_k; k++) 00292 Pii[j*m_k+k] /= w_sum_vector[j]; 00293 } 00294 00295 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans, 00296 m_k,m_k,dp, 00297 1.0,Pii,m_k, 00298 Pii,m_k, 00299 0.0,q_matrix,m_k); 00300 #ifdef HAVE_PTHREAD 00301 PTHREAD_LOCK(W_matrix_lock); 00302 #endif 00303 for (j=0; j<m_k; j++) 00304 { 00305 for (k=0; k<m_k; k++) 00306 W_matrix[N*neighborhood_matrix[k*N+i]+neighborhood_matrix[j*N+i]] += q_matrix[j*m_k+k]; 00307 } 00308 #ifdef HAVE_PTHREAD 00309 PTHREAD_UNLOCK(W_matrix_lock); 00310 #endif 00311 } 00312 return NULL; 00313 } 00314 #endif /* HAVE_LAPACK */