Skip to content

Commit

Permalink
Optimized LLEs using reformulation with alignment
Browse files Browse the repository at this point in the history
  • Loading branch information
@lisitsyn
lisitsyn committed Nov 6, 2011
1 parent a7bec52 commit 5d580e0
Show file tree
Hide file tree
Showing 4 changed files with 53 additions and 208 deletions.
101 changes: 16 additions & 85 deletions src/shogun/converter/KernelLocallyLinearEmbedding.cpp
View file
Expand Up @@ -68,24 +68,6 @@ struct K_NEIGHBORHOOD_THREAD_PARAM
/// matrix containing neighbors indexes
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 */

CKernelLocallyLinearEmbedding::CKernelLocallyLinearEmbedding() :
Expand Down Expand Up @@ -142,7 +124,7 @@ SGMatrix<float64_t> CKernelLocallyLinearEmbedding::construct_weight_matrix(SGMat
{
int32_t N = kernel_matrix.num_cols;
// loop variables
int32_t i,j,t;
int32_t t;
#ifdef HAVE_PTHREAD
int32_t num_threads = parallel->get_num_threads();
ASSERT(num_threads>0);
Expand Down Expand Up @@ -199,70 +181,7 @@ SGMatrix<float64_t> CKernelLocallyLinearEmbedding::construct_weight_matrix(SGMat
SG_FREE(id_vector);
SG_FREE(local_gram_matrix);

// diag(W) = 1
for (i=0; i<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(kernel_matrix.matrix,N,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[i*N+j] = M_matrix[j*N+i];
}
}
SG_FREE(nz_idxs);
SG_FREE(W_matrix);
return M_matrix;
return SGMatrix<float64_t>(W_matrix,N,N);
}

void* CKernelLocallyLinearEmbedding::run_linearreconstruction_thread(void* p)
Expand Down Expand Up @@ -313,11 +232,23 @@ void* CKernelLocallyLinearEmbedding::run_linearreconstruction_thread(void* p)
for (j=0; j<m_k; j++)
norming += id_vector[j];

cblas_dscal(m_k,-1.0/norming,id_vector,1);
cblas_dscal(m_k,1.0/norming,id_vector,1);

memset(local_gram_matrix,0,sizeof(float64_t)*m_k*m_k);
cblas_dger(CblasColMajor,m_k,m_k,1.0,id_vector,1,id_vector,1,local_gram_matrix,m_k);

// put weights into W matrix
W_matrix[N*i+i] += 1.0;
for (j=0; j<m_k; j++)
W_matrix[N*neighborhood_matrix[j*N+i]+i]=id_vector[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];
}
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];
}
}
return NULL;
}
Expand Down
12 changes: 10 additions & 2 deletions src/shogun/converter/KernelLocallyLinearEmbedding.h
View file
Expand Up @@ -27,9 +27,17 @@ class CKernel;
* data using kernel extension of Locally Linear Embedding algorithm as
* described in
*
* Kayo, O. (2006). Locally linear embedding algorithm. Extensions and applications. October.
* Retrieved from: http://herkules.oulu.fi/isbn9514280415/isbn9514280415.pdf
* Decoste, D. (2001).
* Visualizing Mercer Kernel Feature Spaces Via Kernelized Locally-Linear Embeddings.
* The 8th International Conference on Neural Information Processing ICONIP2001
*
* It is optimized with alignment formulation as described in
*
* Zhao, D. (2006).
* Formulating LLE using alignment technique.
* Pattern Recognition, 39(11), 2233-2235.
* Retrieved from http://linkinghub.elsevier.com/retrieve/pii/S0031320306002160
*
*/
class CKernelLocallyLinearEmbedding: public CLocallyLinearEmbedding
{
Expand Down
135 changes: 19 additions & 116 deletions src/shogun/converter/LocallyLinearEmbedding.cpp
View file
Expand Up @@ -78,24 +78,6 @@ 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() :
Expand Down Expand Up @@ -342,7 +324,7 @@ SGMatrix<float64_t> CLocallyLinearEmbedding::construct_weight_matrix(CSimpleFeat
{
int32_t N = simple_features->get_num_vectors();
int32_t dim = simple_features->get_num_features();
int32_t i,j,t;
int32_t t;
#ifdef HAVE_PTHREAD
int32_t num_threads = parallel->get_num_threads();
ASSERT(num_threads>0);
Expand Down Expand Up @@ -409,70 +391,7 @@ SGMatrix<float64_t> CLocallyLinearEmbedding::construct_weight_matrix(CSimpleFeat
SG_FREE(z_matrix);
SG_FREE(covariance_matrix);

// diag(W)=1
for (i=0; i<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);
#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[i*N+j] = M_matrix[j*N+i];
}
}
SG_FREE(nz_idxs);
SG_FREE(W_matrix);
return M_matrix;
return SGMatrix<float64_t>(W_matrix,N,N);
}

SGMatrix<float64_t> CLocallyLinearEmbedding::construct_embedding(CFeatures* features,SGMatrix<float64_t> matrix,int dimension)
Expand All @@ -488,9 +407,9 @@ SGMatrix<float64_t> CLocallyLinearEmbedding::construct_embedding(CFeatures* feat
float64_t* nullspace_features = NULL;
if (m_use_arpack)
{
#ifndef HAVE_ARPACK
#ifndef HAVE_ARPACK
SG_ERROR("ARPACK is not supported in this configuration.\n");
#endif
#endif
// using ARPACK (faster)
eigenvalues_vector = SG_MALLOC(float64_t, dimension+1);
#ifdef HAVE_ARPACK
Expand Down Expand Up @@ -545,7 +464,7 @@ void* CLocallyLinearEmbedding::run_linearreconstruction_thread(void* p)
float64_t* W_matrix = parameters->W_matrix;
float64_t m_reconstruction_shift = parameters->m_reconstruction_shift;

int32_t i,j;
int32_t i,j,k;
float64_t norming,trace;

for (i=idx_start; i<idx_stop; i+=idx_step)
Expand Down Expand Up @@ -589,11 +508,23 @@ void* CLocallyLinearEmbedding::run_linearreconstruction_thread(void* p)
for (j=0; j<m_k; j++)
norming += id_vector[j];

cblas_dscal(m_k,-1.0/norming,id_vector,1);
cblas_dscal(m_k,1.0/norming,id_vector,1);

memset(covariance_matrix,0,sizeof(float64_t)*m_k*m_k);
cblas_dger(CblasColMajor,m_k,m_k,1.0,id_vector,1,id_vector,1,covariance_matrix,m_k);

// put weights into W matrix
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];
}
for (j=0; j<m_k; j++)
W_matrix[N*neighborhood_matrix[j*N+i]+i]=id_vector[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];
}
}
return NULL;
}
Expand Down Expand Up @@ -657,34 +588,6 @@ SGMatrix<int32_t> CLocallyLinearEmbedding::get_neighborhood_matrix(SGMatrix<floa
return SGMatrix<int32_t>(neighborhood_matrix,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)
{
NEIGHBORHOOD_THREAD_PARAM* parameters = (NEIGHBORHOOD_THREAD_PARAM*)p;
Expand Down
13 changes: 8 additions & 5 deletions src/shogun/converter/LocallyLinearEmbedding.h
View file
Expand Up @@ -51,6 +51,14 @@ class CDistance;
* is being done using ternary search of minima of
* the mean reconstruction error. The reconstruction error is
* considered to have only one global minimum in this mode.
*
* It is optimized with alignment formulation as described in
*
* Zhao, D. (2006).
* Formulating LLE using alignment technique.
* Pattern Recognition, 39(11), 2233-2235.
* Retrieved from http://linkinghub.elsevier.com/retrieve/pii/S0031320306002160
*
*/
class CLocallyLinearEmbedding: public CEmbeddingConverter
{
Expand Down Expand Up @@ -217,11 +225,6 @@ class CLocallyLinearEmbedding: public CEmbeddingConverter
*/
static void* run_linearreconstruction_thread(void* p);

/** runs sparse matrix-matrix multiplication thread
* @param p thread params
*/
static void* run_sparsedot_thread(void* p);

};
}

Expand Down

0 comments on commit 5d580e0

Please sign in to comment.