00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <shogun/converter/LinearLocalTangentSpaceAlignment.h>
00012 #ifdef HAVE_LAPACK
00013 #include <shogun/mathematics/arpack.h>
00014 #include <shogun/mathematics/lapack.h>
00015 #include <shogun/mathematics/Math.h>
00016 #include <shogun/io/SGIO.h>
00017 #include <shogun/lib/Time.h>
00018 #include <shogun/distance/Distance.h>
00019 #include <shogun/lib/Signal.h>
00020
00021 using namespace shogun;
00022
00023 CLinearLocalTangentSpaceAlignment::CLinearLocalTangentSpaceAlignment() :
00024 CLocalTangentSpaceAlignment()
00025 {
00026 }
00027
00028 CLinearLocalTangentSpaceAlignment::~CLinearLocalTangentSpaceAlignment()
00029 {
00030 }
00031
00032 const char* CLinearLocalTangentSpaceAlignment::get_name() const
00033 {
00034 return "LinearLocalTangentSpaceAlignment";
00035 }
00036
00037 SGMatrix<float64_t> CLinearLocalTangentSpaceAlignment::construct_embedding(CFeatures* features, SGMatrix<float64_t> matrix, int dimension)
00038 {
00039 CDenseFeatures<float64_t>* simple_features = (CDenseFeatures<float64_t>*)features;
00040 ASSERT(simple_features);
00041 int i,j;
00042
00043 SGMatrix<float64_t> feature_matrix = simple_features->get_feature_matrix().clone();
00044 int N=feature_matrix.num_cols;
00045 int dim=feature_matrix.num_rows;
00046 ASSERT(dimension<=dim);
00047 float64_t* XTM = SG_MALLOC(float64_t, dim*N);
00048 float64_t* lhs_M = SG_MALLOC(float64_t, dim*dim);
00049 float64_t* rhs_M = SG_MALLOC(float64_t, dim*dim);
00050 SGMatrix<float64_t>::center_matrix(matrix.matrix,N,N);
00051
00052 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,dim,N,N,1.0,feature_matrix.matrix,dim,matrix.matrix,N,0.0,XTM,dim);
00053 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,dim,dim,N,1.0,XTM,dim,feature_matrix.matrix,dim,0.0,lhs_M,dim);
00054
00055 float64_t* mean_vector = SG_CALLOC(float64_t, dim);
00056 for (i=0; i<N; i++)
00057 cblas_daxpy(dim,1.0,feature_matrix.matrix+i*dim,1,mean_vector,1);
00058
00059 cblas_dscal(dim,1.0/N,mean_vector,1);
00060
00061 for (i=0; i<N; i++)
00062 cblas_daxpy(dim,-1.0,mean_vector,1,feature_matrix.matrix+i*dim,1);
00063
00064 SG_FREE(mean_vector);
00065
00066 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,dim,dim,N,1.0,feature_matrix.matrix,dim,feature_matrix.matrix,dim,0.0,rhs_M,dim);
00067
00068 for (i=0; i<dim; i++) rhs_M[i*dim+i] += 1e-6;
00069
00070 float64_t* evals = SG_MALLOC(float64_t, dim);
00071 float64_t* evectors = SG_MALLOC(float64_t, dimension*dim);
00072 int32_t info = 0;
00073 #ifdef HAVE_ARPACK
00074 arpack_dsxupd(lhs_M,rhs_M,false,dim,dimension,"LA",false,3,true,false,m_nullspace_shift,0.0,
00075 evals,evectors,info);
00076 #else
00077 wrap_dsygvx(1,'V','U',dim,lhs_M,dim,rhs_M,dim,dim-dimension+1,dim,evals,evectors,&info);
00078 #endif
00079 SG_FREE(lhs_M);
00080 SG_FREE(rhs_M);
00081 SG_FREE(evals);
00082 if (info!=0) SG_ERROR("Failed to solve eigenproblem (%d)\n",info);
00083
00084 for (i=0; i<dimension/2; i++)
00085 {
00086 cblas_dswap(dim,evectors+i*dim,1,evectors+(dimension-i-1)*dim,1);
00087 }
00088
00089 cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans,N,dimension,dim,1.0,feature_matrix.matrix,dim,evectors,dim,0.0,XTM,N);
00090 SG_FREE(evectors);
00091
00092 float64_t* new_features = SG_MALLOC(float64_t, dimension*N);
00093 for (i=0; i<dimension; i++)
00094 {
00095 for (j=0; j<N; j++)
00096 new_features[j*dimension+i] = XTM[i*N+j];
00097 }
00098 SG_FREE(XTM);
00099 return SGMatrix<float64_t>(new_features,dimension,N);
00100 }
00101
00102 #endif