Improved performance of locally linear embedding

shogun-toolbox · Sep 10, 2011 · 613d9dc · 613d9dc
1 parent f8944e0
commit 613d9dc
Show file tree

Hide file tree

Showing 2 changed files with 113 additions and 10 deletions.
diff --git a/src/shogun/preprocessor/LocallyLinearEmbedding.cpp b/src/shogun/preprocessor/LocallyLinearEmbedding.cpp
@@ -15,6 +15,8 @@
 #include <shogun/mathematics/lapack.h>
 #include <shogun/lib/FibonacciHeap.h>
 #include <shogun/lib/common.h>
+#include <shogun/lib/Time.h>
+#include <list>
 #include <shogun/mathematics/Math.h>
 #include <shogun/io/SGIO.h>
 #include <shogun/distance/EuclidianDistance.h>
@@ -74,6 +76,24 @@ struct NEIGHBORHOOD_THREAD_PARAM
 	/// matrix containing indexes of neighbors of ith object in ith column
 	int32_t* neighborhood_matrix;
 };
+
+struct SPARSEDOT_THREAD_PARAM
+{
+	/// starting index of loop
+	int32_t idx_start;
+	/// step of loop
+	int32_t idx_step;
+	/// end index of loop
+	int32_t idx_stop;
+	/// number of vectors
+	int32_t N;
+	/// weight matrix
+	const float64_t* W_matrix;
+	/// result matrix
+	float64_t* M_matrix;
+	/// non zero indexes dynamic array
+	DynArray<int32_t>** nz_idxs;
+};
 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
 
 CLocallyLinearEmbedding::CLocallyLinearEmbedding() :
@@ -202,26 +222,77 @@ SGMatrix<float64_t> CLocallyLinearEmbedding::apply_to_feature_matrix(CFeatures*
 #endif
 
 	// clean
-	SG_FREE(id_vector);
 	neighborhood_matrix.destroy_matrix();
+	SG_FREE(id_vector);
 	SG_FREE(z_matrix);
 	SG_FREE(covariance_matrix);
 
 	// W=I-W
+	for (i=0; i<N*N; i++)
+	{
+		W_matrix[i] *= -1.0;
+	}
 	for (i=0; i<N; i++)
 	{
-		for (j=0; j<N; j++)
-			W_matrix[j*N+i] = (i==j) ? 1.0-W_matrix[j*N+i] : -W_matrix[j*N+i];
+		W_matrix[i*N+i] = 1.0;
 	}
 
 	// compute M=(W-I)'*(W-I)
+	DynArray<int32_t>** nz_idxs = SG_MALLOC(DynArray<int32_t>*,N);
+	for (i=0; i<N; i++)
+	{
+		nz_idxs[i] = new DynArray<int32_t>(m_k,false);
+		for (j=0; j<N; j++)
+		{
+			if (W_matrix[i*N+j]!=0.0)
+				nz_idxs[i]->push_back(j);
+		}
+	}
 	SGMatrix<float64_t> M_matrix(N,N);
-	cblas_dgemm(CblasColMajor,CblasTrans, CblasNoTrans,
-	            N,N,N,
-	            1.0,W_matrix,N,
-	            W_matrix,N,
-	            0.0,M_matrix.matrix,N);
-
+#ifdef HAVE_PTHREAD
+	// allocate threads
+	threads = SG_MALLOC(pthread_t, num_threads);
+	SPARSEDOT_THREAD_PARAM* parameters_ = SG_MALLOC(SPARSEDOT_THREAD_PARAM, num_threads);
+	pthread_attr_t attr_;
+	pthread_attr_init(&attr_);
+	pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
+
+	for (t=0; t<num_threads; t++)
+	{
+		parameters_[t].idx_start = t;
+		parameters_[t].idx_step = num_threads;
+		parameters_[t].idx_stop = N;
+		parameters_[t].N = N;
+		parameters_[t].W_matrix = W_matrix;
+		parameters_[t].M_matrix = M_matrix.matrix;
+		parameters_[t].nz_idxs = nz_idxs;
+		pthread_create(&threads[t], &attr_, run_sparsedot_thread, (void*)&parameters_[t]);
+	}
+	for (t=0; t<num_threads; t++)
+		pthread_join(threads[t], NULL);
+	pthread_attr_destroy(&attr_);
+	SG_FREE(parameters_);
+	SG_FREE(threads);
+#else
+	SPARSEDOT_THREAD_PARAM single_thread_param;
+	single_thread_param.idx_start = 0;
+	single_thread_param.idx_step = 1;
+	single_thread_param.idx_stop = N;
+	single_thread_param.N = N;
+	single_thread_param.W_matrix = W_matrix;
+	single_thread_param.M_matrix = M_matrix.matrix;
+	single_thread_param.nz_idxs = nz_idxs;
+	run_sparsedot_thread((void*)single_thread_param);
+#endif
+	for (i=0; i<N; i++)
+	{
+		delete nz_idxs[i];
+		for (j=0; j<i; j++)
+		{
+			M_matrix.matrix[i*N+j] = M_matrix.matrix[j*N+i];
+		}
+	}
+	SG_FREE(nz_idxs);
 	SG_FREE(W_matrix);
 
 	simple_features->set_feature_matrix(find_null_space(M_matrix,m_target_dim,false));
@@ -260,7 +331,7 @@ SGMatrix<float64_t> CLocallyLinearEmbedding::find_null_space(SGMatrix<float64_t>
 		// using ARPACK (faster)
 		eigenvalues_vector = SG_MALLOC(float64_t, dimension+1);
 		#ifdef HAVE_ARPACK
-		arpack_dsaeupd_wrap(matrix.matrix,NULL,N,dimension+1,"LA",3,m_posdef,-1e-6, 0.0,
+		arpack_dsaeupd_wrap(matrix.matrix,NULL,N,dimension+1,"LA",3,m_posdef,0.0, 0.0,
 		                    eigenvalues_vector,matrix.matrix,eigenproblem_status);
 		#endif
 	}
@@ -441,6 +512,33 @@ SGMatrix<int32_t> CLocallyLinearEmbedding::get_neighborhood_matrix(CDistance* di
 	return SGMatrix<int32_t>(neighborhood_matrix,m_k,N);
 }
 
+void* CLocallyLinearEmbedding::run_sparsedot_thread(void* p)
+{
+	SPARSEDOT_THREAD_PARAM* parameters = (SPARSEDOT_THREAD_PARAM*)p;
+	int32_t idx_start = parameters->idx_start;
+	int32_t idx_step = parameters->idx_step;
+	int32_t idx_stop = parameters->idx_stop;
+	int32_t N = parameters->N;
+	const float64_t* W_matrix = parameters->W_matrix;
+	float64_t* M_matrix = parameters->M_matrix;
+	DynArray<int32_t>** nz_idxs = parameters->nz_idxs;
+
+	int i,j,k;
+	for (i=idx_start; i<idx_stop; i+=idx_step)
+	{
+		for (j=i; j<N; j++)
+		{
+			M_matrix[i*N+j] = 0.0;
+			for (k=0; k<nz_idxs[i]->get_num_elements(); k++)
+			{
+				M_matrix[i*N+j] += W_matrix[i*N+(*nz_idxs[i])[k]]*W_matrix[j*N+(*nz_idxs[i])[k]];
+			}
+
+		}
+	}
+
+	return NULL;
+}
 
 void* CLocallyLinearEmbedding::run_neighborhood_thread(void* p)
 {

diff --git a/src/shogun/preprocessor/LocallyLinearEmbedding.h b/src/shogun/preprocessor/LocallyLinearEmbedding.h
@@ -129,6 +129,11 @@ class CLocallyLinearEmbedding: public CDimensionReductionPreprocessor
 	 */
 	static void* run_linearreconstruction_thread(void* p);
 
+	/** runs sparse matrix-matrix multiplication thread
+	 * @param p thread params
+	 */
+	static void* run_sparsedot_thread(void* p);
+
 	/** find null space of given matrix 
 	 * @param matrix given matrix
 	 * @param dimension dimension of null space to be computed