Tapkee
laplacian_eigenmaps.hpp
Go to the documentation of this file.
00001 /* This software is distributed under BSD 3-clause license (see LICENSE file).
00002  *
00003  * Copyright (c) 2012-2013 Sergey Lisitsyn
00004  */
00005 
00006 #ifndef TAPKEE_LaplacianEigenmaps_H_
00007 #define TAPKEE_LaplacianEigenmaps_H_
00008 
00009 /* Tapkee includes */
00010 #include <tapkee/defines.hpp>
00011 #include <tapkee/utils/time.hpp>
00012 /* End of Tapkee includes */
00013 
00014 namespace tapkee
00015 {
00016 namespace tapkee_internal
00017 {
00018 
00038 template<class RandomAccessIterator, class DistanceCallback>
00039 Laplacian compute_laplacian(RandomAccessIterator begin, 
00040             RandomAccessIterator end,const Neighbors& neighbors, 
00041             DistanceCallback callback, ScalarType width)
00042 {
00043     SparseTriplets sparse_triplets;
00044 
00045     timed_context context("Laplacian computation");
00046     const IndexType k = neighbors[0].size();
00047     sparse_triplets.reserve((k+1)*(end-begin));
00048 
00049     DenseVector D = DenseVector::Zero(end-begin);
00050     for (RandomAccessIterator iter=begin; iter!=end; ++iter)
00051     {
00052         const LocalNeighbors& current_neighbors = neighbors[iter-begin];
00053 
00054         for (IndexType i=0; i<k; ++i)
00055         {
00056             ScalarType distance = callback.distance(*iter,begin[current_neighbors[i]]);
00057             ScalarType heat = exp(-distance*distance/width);
00058             D(iter-begin) += heat;
00059             D(current_neighbors[i]) += heat;
00060             sparse_triplets.push_back(SparseTriplet(current_neighbors[i],(iter-begin),-heat));
00061             sparse_triplets.push_back(SparseTriplet((iter-begin),current_neighbors[i],-heat));
00062         }
00063     }
00064     for (IndexType i=0; i<static_cast<IndexType>(end-begin); ++i)
00065         sparse_triplets.push_back(SparseTriplet(i,i,D(i)));
00066 
00067 #ifdef EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
00068     Eigen::DynamicSparseMatrix<ScalarType> dynamic_weight_matrix(end-begin,end-begin);
00069     dynamic_weight_matrix.reserve(sparse_triplets.size());
00070     for (SparseTriplets::const_iterator it=sparse_triplets.begin(); it!=sparse_triplets.end(); ++it)
00071         dynamic_weight_matrix.coeffRef(it->col(),it->row()) += it->value();
00072     SparseWeightMatrix weight_matrix(dynamic_weight_matrix);
00073 #else
00074     SparseWeightMatrix weight_matrix(end-begin,end-begin);
00075     weight_matrix.setFromTriplets(sparse_triplets.begin(),sparse_triplets.end());
00076 #endif
00077 
00078     return Laplacian(weight_matrix,DenseDiagonalMatrix(D));
00079 }
00080 
00081 template<class RandomAccessIterator, class FeatureVectorCallback>
00082 DenseSymmetricMatrixPair construct_locality_preserving_eigenproblem(SparseWeightMatrix& L,
00083         DenseDiagonalMatrix& D, RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback feature_vector_callback,
00084         IndexType dimension)
00085 {
00086     timed_context context("Constructing LPP eigenproblem");
00087 
00088     DenseSymmetricMatrix lhs = DenseSymmetricMatrix::Zero(dimension,dimension);
00089     DenseSymmetricMatrix rhs = DenseSymmetricMatrix::Zero(dimension,dimension);
00090 
00091     DenseVector rank_update_vector_i(dimension);
00092     DenseVector rank_update_vector_j(dimension);
00093     for (RandomAccessIterator iter=begin; iter!=end; ++iter)
00094     {
00095         feature_vector_callback.vector(*iter,rank_update_vector_i);
00096         rhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i,D.diagonal()(iter-begin));
00097     }
00098 
00099     for (int i=0; i<L.outerSize(); ++i)
00100     {
00101         for (SparseWeightMatrix::InnerIterator it(L,i); it; ++it)
00102         {
00103             feature_vector_callback.vector(begin[it.row()],rank_update_vector_i);
00104             feature_vector_callback.vector(begin[it.col()],rank_update_vector_j);
00105             lhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i, rank_update_vector_j, it.value());
00106         }
00107     }
00108 
00109     return DenseSymmetricMatrixPair(lhs,rhs);
00110 }
00111 
00112 }
00113 }
00114 
00115 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines