Fixes and improvements of EDRT

shogun-toolbox · Jan 30, 2012 · 689b5ed · 689b5ed
1 parent bb5c66e
commit 689b5ed
Show file tree

Hide file tree

Showing 6 changed files with 65 additions and 31 deletions.
diff --git a/src/shogun/converter/DiffusionMaps.cpp b/src/shogun/converter/DiffusionMaps.cpp
@@ -122,32 +122,34 @@ CSimpleFeatures<float64_t>* CDiffusionMaps::embed_kernel(CKernel* kernel)
 	float64_t* s_values = SG_MALLOC(float64_t, N);
 
 	int32_t info = 0;
-	SGMatrix<float64_t> new_feature_matrix;
+	SGMatrix<float64_t> new_feature_matrix = SGMatrix<float64_t>(m_target_dim,N);
 
 	if (use_arpack)
 	{
 #ifdef HAVE_ARPACK
 		arpack_dsxupd(kernel_matrix.matrix,NULL,false,N,m_target_dim,"LA",false,1,false,true,0.0,0.0,s_values,kernel_matrix.matrix,info);
 #endif /* HAVE_ARPACK */
-		SG_REALLOC(float64_t,kernel_matrix.matrix,N*m_target_dim);
-		new_feature_matrix = SGMatrix<float64_t>(kernel_matrix.matrix,m_target_dim,N);
+		for (i=0; i<m_target_dim; i++)
+		{
+			for (j=0; j<N; j++)
+				new_feature_matrix[j*m_target_dim+i] = kernel_matrix[j*m_target_dim+i];
+		}
 	}
 	else 
 	{
 		SG_WARNING("LAPACK does not provide efficient routines to construct embedding (this may take time). Consider installing ARPACK.");
 		wrap_dgesvd('O','N',N,N,kernel_matrix.matrix,N,s_values,NULL,1,NULL,1,&info);
-		new_feature_matrix = SGMatrix<float64_t>(m_target_dim,N);
 		for (i=0; i<m_target_dim; i++)
 		{
 			for (j=0; j<N; j++)
 				new_feature_matrix[j*m_target_dim+i] = 
 				    kernel_matrix[(m_target_dim-i-1)*N+j];
 		}
-		kernel_matrix.destroy_matrix();
 	}
 	if (info)
 		SG_ERROR("Eigenproblem solving  failed with %d code", info);
 
+	kernel_matrix.destroy_matrix();
 	SG_FREE(s_values);
 
 	return new CSimpleFeatures<float64_t>(new_feature_matrix);

diff --git a/src/shogun/converter/HessianLocallyLinearEmbedding.cpp b/src/shogun/converter/HessianLocallyLinearEmbedding.cpp
@@ -230,7 +230,7 @@ void* CHessianLocallyLinearEmbedding::run_hessianestimation_thread(void* p)
 		for (j=0; j<m_k; j++)
 		{
 			for (k=0; k<dim; k++)
-				local_feature_matrix[j*dim+k] = feature_matrix[neighborhood_matrix[j*N+i]*dim+k];
+				local_feature_matrix[j*dim+k] = feature_matrix[neighborhood_matrix[i*m_k+j]*dim+k];
 
 			cblas_daxpy(dim,1.0,local_feature_matrix+j*dim,1,mean_vector,1);
 		}
@@ -303,7 +303,7 @@ void* CHessianLocallyLinearEmbedding::run_hessianestimation_thread(void* p)
 		for (j=0; j<m_k; j++)
 		{
 			for (k=0; k<m_k; k++)
-				W_matrix[N*neighborhood_matrix[k*N+i]+neighborhood_matrix[j*N+i]] += q_matrix[j*m_k+k];
+				W_matrix[N*neighborhood_matrix[i*m_k+k]+neighborhood_matrix[i*m_k+j]] += q_matrix[j*m_k+k];
 		}
 #ifdef HAVE_PTHREAD
 		PTHREAD_UNLOCK(W_matrix_lock);

diff --git a/src/shogun/converter/KernelLocalTangentSpaceAlignment.cpp b/src/shogun/converter/KernelLocalTangentSpaceAlignment.cpp
@@ -131,12 +131,11 @@ SGMatrix<float64_t> CKernelLocalTangentSpaceAlignment::construct_weight_matrix(S
 	SG_FREE(local_gram_matrix);
 	SG_FREE(ev_vector);
 	SG_FREE(G_matrix);
-	kernel_matrix.destroy_matrix();
 
 	for (int32_t i=0; i<N; i++)
 	{
 		for (int32_t j=0; j<m_k; j++)
-			W_matrix[N*neighborhood_matrix[j*N+i]+neighborhood_matrix[j*N+i]] += 1.0;
+			W_matrix[N*neighborhood_matrix[i*m_k+j]+neighborhood_matrix[i*m_k+j]] += 1.0;
 	}
 
 	return SGMatrix<float64_t>(W_matrix,N,N);
@@ -171,7 +170,7 @@ void* CKernelLocalTangentSpaceAlignment::run_kltsa_thread(void* p)
 		for (j=0; j<m_k; j++)
 		{
 			for (k=0; k<m_k; k++)
-				local_gram_matrix[j*m_k+k] = kernel_matrix[neighborhood_matrix[j*N+i]*N+neighborhood_matrix[k*N+i]];
+				local_gram_matrix[j*m_k+k] = kernel_matrix[neighborhood_matrix[i*m_k+j]*N+neighborhood_matrix[i*m_k+k]];
 		}
 
 		CMath::center_matrix(local_gram_matrix,m_k,m_k);
@@ -192,7 +191,7 @@ void* CKernelLocalTangentSpaceAlignment::run_kltsa_thread(void* p)
 		for (j=0; j<m_k; j++)
 		{
 			for (k=0; k<m_k; k++)
-				W_matrix[N*neighborhood_matrix[k*N+i]+neighborhood_matrix[j*N+i]] -= local_gram_matrix[j*m_k+k];
+				W_matrix[N*neighborhood_matrix[i*m_k+k]+neighborhood_matrix[i*m_k+j]] -= local_gram_matrix[j*m_k+k];
 		}
 #ifdef HAVE_PTHREAD
 		PTHREAD_UNLOCK(W_matrix_lock);

diff --git a/src/shogun/converter/KernelLocallyLinearEmbedding.cpp b/src/shogun/converter/KernelLocallyLinearEmbedding.cpp
@@ -16,6 +16,7 @@
 #include <shogun/base/DynArray.h>
 #include <shogun/mathematics/Math.h>
 #include <shogun/io/SGIO.h>
+#include <shogun/lib/Time.h>
 #include <shogun/lib/Signal.h>
 
 #ifdef HAVE_PTHREAD
@@ -67,6 +68,8 @@ struct K_NEIGHBORHOOD_THREAD_PARAM
 	CFibonacciHeap* heap;
 	/// kernel matrix
 	const float64_t* kernel_matrix;
+	/// kernel diag
+	const float64_t* kernel_diag;
 	/// matrix containing neighbors indexes
 	int32_t* neighborhood_matrix;
 };
@@ -114,15 +117,29 @@ CFeatures* CKernelLocallyLinearEmbedding::apply(CFeatures* features)
 
 CSimpleFeatures<float64_t>* CKernelLocallyLinearEmbedding::embed_kernel(CKernel* kernel)
 {
+	CTime* time = new CTime();
+
+	time->start();
 	SGMatrix<float64_t> kernel_matrix = kernel->get_kernel_matrix();
+	SG_DEBUG("Kernel matrix computation took %fs\n",time->cur_time_diff());
+
+	time->start();
 	SGMatrix<int32_t> neighborhood_matrix = get_neighborhood_matrix(kernel_matrix,m_k);
-
+	SG_DEBUG("Neighbors finding took %fs\n",time->cur_time_diff());
+
+	time->start();
 	SGMatrix<float64_t> M_matrix = construct_weight_matrix(kernel_matrix,neighborhood_matrix);
+	SG_DEBUG("Weights computation took %fs\n",time->cur_time_diff());
+	kernel_matrix.destroy_matrix();
 	neighborhood_matrix.destroy_matrix();
 
+	time->start();
 	SGMatrix<float64_t> nullspace = construct_embedding(M_matrix,m_target_dim);
+	SG_DEBUG("Embedding construction took %fs\n",time->cur_time_diff());
 	M_matrix.destroy_matrix();
 
+	delete time;
+
 	return new CSimpleFeatures<float64_t>(nullspace);
 }
 
@@ -217,9 +234,9 @@ void* CKernelLocallyLinearEmbedding::run_linearreconstruction_thread(void* p)
 			for (k=0; k<m_k; k++)
 				local_gram_matrix[j*m_k+k] = 
 					kernel_matrix[i*N+i] -
-					kernel_matrix[i*N+neighborhood_matrix[j*N+i]] -
-					kernel_matrix[i*N+neighborhood_matrix[k*N+i]] +
-					kernel_matrix[neighborhood_matrix[j*N+i]*N+neighborhood_matrix[k*N+i]];
+					kernel_matrix[i*N+neighborhood_matrix[i*m_k+j]] -
+					kernel_matrix[i*N+neighborhood_matrix[i*m_k+k]] +
+					kernel_matrix[neighborhood_matrix[i*m_k+j]*N+neighborhood_matrix[i*m_k+k]];
 		}
 
 		for (j=0; j<m_k; j++)
@@ -250,13 +267,13 @@ void* CKernelLocallyLinearEmbedding::run_linearreconstruction_thread(void* p)
 		W_matrix[N*i+i] += 1.0;
 		for (j=0; j<m_k; j++)
 		{
-			W_matrix[N*i+neighborhood_matrix[j*N+i]] -= id_vector[j];
-			W_matrix[N*neighborhood_matrix[j*N+i]+i] -= id_vector[j];
+			W_matrix[N*i+neighborhood_matrix[i*m_k+j]] -= id_vector[j];
+			W_matrix[N*neighborhood_matrix[i*m_k+j]+i] -= id_vector[j];
 		}
 		for (j=0; j<m_k; j++)
 		{
 			for (k=0; k<m_k; k++)
-				W_matrix[N*neighborhood_matrix[j*N+i]+neighborhood_matrix[k*N+i]]+=local_gram_matrix[j*m_k+k];
+				W_matrix[N*neighborhood_matrix[i*m_k+j]+neighborhood_matrix[i*m_k+k]]+=local_gram_matrix[j*m_k+k];
 		}
 	}
 	return NULL;
@@ -268,6 +285,11 @@ SGMatrix<int32_t> CKernelLocallyLinearEmbedding::get_neighborhood_matrix(SGMatri
 	int32_t N = kernel_matrix.num_cols;
 	// init matrix and heap to be used
 	int32_t* neighborhood_matrix = SG_MALLOC(int32_t, N*k);
+	float64_t* kernel_diag = SG_MALLOC(float64_t, N);
+	// dirty
+	for (t=0; t<N; t++)
+		kernel_diag[t] = kernel_matrix.matrix[t*N+t];
+
 #ifdef HAVE_PTHREAD
 	int32_t num_threads = parallel->get_num_threads();
 	ASSERT(num_threads>0);
@@ -294,6 +316,7 @@ SGMatrix<int32_t> CKernelLocallyLinearEmbedding::get_neighborhood_matrix(SGMatri
 		parameters[t].heap = heaps[t];
 		parameters[t].neighborhood_matrix = neighborhood_matrix;
 		parameters[t].kernel_matrix = kernel_matrix.matrix;
+		parameters[t].kernel_diag = kernel_diag;
 		pthread_create(&threads[t], &attr, run_neighborhood_thread, (void*)&parameters[t]);
 	}
 	for (t=0; t<num_threads; t++)
@@ -311,9 +334,11 @@ SGMatrix<int32_t> CKernelLocallyLinearEmbedding::get_neighborhood_matrix(SGMatri
 	single_thread_param.heap = heaps[0]
 	single_thread_param.neighborhood_matrix = neighborhood_matrix;
 	single_thread_param.kernel_matrix = kernel_matrix.matrix;
+	single_thread_param.kernel_diag = kernel_diag;
 	run_neighborhood_thread((void*)&single_thread_param);
 #endif
 
+	SG_FREE(kernel_diag);
 	for (t=0; t<num_threads; t++)
 		delete heaps[t];
 	SG_FREE(heaps);
@@ -331,21 +356,23 @@ void* CKernelLocallyLinearEmbedding::run_neighborhood_thread(void* p)
 	int32_t m_k = parameters->m_k;
 	CFibonacciHeap* heap = parameters->heap;
 	const float64_t* kernel_matrix = parameters->kernel_matrix;
+	const float64_t* kernel_diag = parameters->kernel_diag;
 	int32_t* neighborhood_matrix = parameters->neighborhood_matrix;
 
 	int32_t i,j;
-	float64_t tmp;
+	float64_t tmp,k_diag_i,k_diag_j;
 	for (i=idx_start; i<idx_stop; i+=idx_step)
 	{
+		k_diag_i = kernel_diag[i];
 		for (j=0; j<N; j++)
 		{
-			heap->insert(j,kernel_matrix[i*N+i]-2*kernel_matrix[i*N+j]+kernel_matrix[j*N+j]);
+			heap->insert(j,k_diag_i+kernel_diag[j]-2.0*kernel_matrix[i*N+j]);
 		}
 
 		heap->extract_min(tmp);
 
 		for (j=0; j<m_k; j++)
-			neighborhood_matrix[j*N+i] = heap->extract_min(tmp);
+			neighborhood_matrix[i*m_k+j] = heap->extract_min(tmp);
 
 		heap->clear();
 	}

diff --git a/src/shogun/converter/LocalTangentSpaceAlignment.cpp b/src/shogun/converter/LocalTangentSpaceAlignment.cpp
@@ -166,7 +166,7 @@ SGMatrix<float64_t> CLocalTangentSpaceAlignment::construct_weight_matrix(CSimple
 	for (int32_t i=0; i<N; i++)
 	{
 		for (int32_t j=0; j<m_k; j++)
-			W_matrix[N*neighborhood_matrix[j*N+i]+neighborhood_matrix[j*N+i]] += 1.0;
+			W_matrix[N*neighborhood_matrix[i*m_k+j]+neighborhood_matrix[i*m_k+j]] += 1.0;
 	}
 
 	return SGMatrix<float64_t>(W_matrix,N,N);
@@ -208,7 +208,7 @@ void* CLocalTangentSpaceAlignment::run_ltsa_thread(void* p)
 		for (j=0; j<m_k; j++)
 		{
 			for (k=0; k<dim; k++)
-				local_feature_matrix[j*dim+k] = feature_matrix[neighborhood_matrix[j*N+i]*dim+k];
+				local_feature_matrix[j*dim+k] = feature_matrix[neighborhood_matrix[i*m_k+j]*dim+k];
 
 			cblas_daxpy(dim,1.0,local_feature_matrix+j*dim,1,mean_vector,1);
 		}
@@ -246,7 +246,7 @@ void* CLocalTangentSpaceAlignment::run_ltsa_thread(void* p)
 		for (j=0; j<m_k; j++)
 		{
 			for (k=0; k<m_k; k++)
-				W_matrix[N*neighborhood_matrix[k*N+i]+neighborhood_matrix[j*N+i]] -= q_matrix[j*m_k+k];
+				W_matrix[N*neighborhood_matrix[i*m_k+k]+neighborhood_matrix[i*m_k+j]] -= q_matrix[j*m_k+k];
 		}
 #ifdef HAVE_PTHREAD
 		PTHREAD_UNLOCK(W_matrix_lock);

diff --git a/src/shogun/converter/LocallyLinearEmbedding.cpp b/src/shogun/converter/LocallyLinearEmbedding.cpp
@@ -205,12 +205,16 @@ CFeatures* CLocallyLinearEmbedding::apply(CFeatures* features)
 	// compute distance matrix
 	SG_DEBUG("Computing distance matrix\n");
 	ASSERT(m_distance);
+	CTime* time = new CTime();
+	time->start();
 	m_distance->init(simple_features,simple_features);
 	SGMatrix<float64_t> distance_matrix = m_distance->get_distance_matrix();
 	m_distance->remove_lhs_and_rhs();
+	SG_DEBUG("Distance matrix computation took %fs\n",time->cur_time_diff());
 	SG_DEBUG("Calculating neighborhood matrix\n");
 	SGMatrix<int32_t> neighborhood_matrix;
 
+	time->start();
 	if (m_auto_k) 
 	{
 		neighborhood_matrix = get_neighborhood_matrix(distance_matrix,m_max_k);
@@ -219,14 +223,16 @@ CFeatures* CLocallyLinearEmbedding::apply(CFeatures* features)
 	} 
 	else
 		neighborhood_matrix = get_neighborhood_matrix(distance_matrix,m_k);
+
+	SG_DEBUG("Neighbors finding took %fs\n",time->cur_time_diff());
 
 	// init W (weight) matrix
 	float64_t* W_matrix = distance_matrix.matrix;
 	memset(W_matrix,0,sizeof(float64_t)*N*N);
 
 	// construct weight matrix
 	SG_DEBUG("Constructing weight matrix\n");
-	CTime* time = new CTime();
+	time->start();
 	SGMatrix<float64_t> weight_matrix = construct_weight_matrix(simple_features,W_matrix,neighborhood_matrix);
 	SG_DEBUG("Weight matrix construction took %.5fs\n", time->cur_time_diff());
 	neighborhood_matrix.destroy_matrix();
@@ -292,7 +298,7 @@ float64_t CLocallyLinearEmbedding::compute_reconstruction_error(int32_t k, int d
 	{
 		for (j=0; j<k; j++)
 		{
-			cblas_dcopy(dim,feature_matrix+neighborhood_matrix[j*N+i]*dim,1,z_matrix+j*dim,1);   
+			cblas_dcopy(dim,feature_matrix+neighborhood_matrix[i*m_k+j]*dim,1,z_matrix+j*dim,1);   
 			cblas_daxpy(dim,-1.0,feature_matrix+i*dim,1,z_matrix+j*dim,1);
 		}
 		cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans,
@@ -479,7 +485,7 @@ void* CLocallyLinearEmbedding::run_linearreconstruction_thread(void* p)
 		// center features by subtracting i-th feature column
 		for (j=0; j<m_k; j++)
 		{
-			cblas_dcopy(dim,feature_matrix+neighborhood_matrix[j*N+i]*dim,1,z_matrix+j*dim,1);
+			cblas_dcopy(dim,feature_matrix+neighborhood_matrix[i*m_k+j]*dim,1,z_matrix+j*dim,1);
 			cblas_daxpy(dim,-1.0,feature_matrix+i*dim,1,z_matrix+j*dim,1);
 		}
 
@@ -523,13 +529,13 @@ void* CLocallyLinearEmbedding::run_linearreconstruction_thread(void* p)
 		W_matrix[N*i+i] += 1.0;
 		for (j=0; j<m_k; j++)
 		{
-			W_matrix[N*i+neighborhood_matrix[j*N+i]] -= id_vector[j];
-			W_matrix[N*neighborhood_matrix[j*N+i]+i] -= id_vector[j];
+			W_matrix[N*i+neighborhood_matrix[i*m_k+j]] -= id_vector[j];
+			W_matrix[N*neighborhood_matrix[i*m_k+j]+i] -= id_vector[j];
 		}
 		for (j=0; j<m_k; j++)
 		{
 			for (k=0; k<m_k; k++)
-				W_matrix[N*neighborhood_matrix[j*N+i]+neighborhood_matrix[k*N+i]]+=covariance_matrix[j*m_k+k];
+				W_matrix[N*neighborhood_matrix[i*m_k+j]+neighborhood_matrix[i*m_k+k]]+=covariance_matrix[j*m_k+k];
 		}
 	}
 	return NULL;
@@ -618,7 +624,7 @@ void* CLocallyLinearEmbedding::run_neighborhood_thread(void* p)
 		heap->extract_min(tmp);
 
 		for (j=0; j<m_k; j++)
-			neighborhood_matrix[j*N+i] = heap->extract_min(tmp);
+			neighborhood_matrix[i*m_k+j] = heap->extract_min(tmp);
 
 		heap->clear();
 	}