6 #ifndef TAPKEE_LOCALLY_LINEAR_H_
7 #define TAPKEE_LOCALLY_LINEAR_H_
19 namespace tapkee_internal
24 template <
class RandomAccessIterator,
class PairwiseCallback>
26 const Neighbors& neighbors, PairwiseCallback callback,
28 const bool partial_eigendecomposer=
false)
34 sparse_triplets.reserve((k*k+2*k+1)*(end-begin));
36 #pragma omp parallel shared(begin,end,neighbors,callback,sparse_triplets) default(none)
41 DenseMatrix G = DenseMatrix::Zero(k,target_dimension+1);
42 G.col(0).setConstant(1/sqrt(static_cast<ScalarType>(k)));
45 local_triplets.reserve(k*k+2*k+1);
47 #pragma omp for nowait
48 for (index_iter=0; index_iter<static_cast<IndexType>(end-begin); index_iter++)
56 ScalarType kij = callback.kernel(begin[current_neighbors[i]],begin[current_neighbors[j]]);
57 gram_matrix(i,j) = kij;
58 gram_matrix(j,i) = kij;
65 #ifdef TAPKEE_WITH_ARPACK
66 if (partial_eigendecomposer)
68 G.rightCols(target_dimension).noalias() =
69 eigendecomposition<DenseMatrix,DenseMatrixOperation>(
Arpack,gram_matrix,target_dimension,0).first;
74 solver.compute(gram_matrix);
75 G.rightCols(target_dimension).noalias() = solver.eigenvectors().rightCols(target_dimension);
78 gram_matrix.noalias() = G * G.transpose();
81 local_triplets.push_back(diagonal_triplet);
84 SparseTriplet neighborhood_diagonal_triplet(current_neighbors[i],current_neighbors[i],1.0);
85 local_triplets.push_back(neighborhood_diagonal_triplet);
89 SparseTriplet tangent_triplet(current_neighbors[i],current_neighbors[j],-gram_matrix(i,j));
90 local_triplets.push_back(tangent_triplet);
95 copy(local_triplets.begin(),local_triplets.end(),back_inserter(sparse_triplets));
98 local_triplets.clear();
105 template <
class RandomAccessIterator,
class PairwiseCallback>
107 const Neighbors& neighbors, PairwiseCallback callback,
114 sparse_triplets.reserve((k*k+2*k+1)*(end-begin));
116 #pragma omp parallel shared(begin,end,neighbors,callback,sparse_triplets) default(none)
124 local_triplets.reserve(k*k+2*k+1);
127 #pragma omp for nowait
128 for (index_iter=0; index_iter<static_cast<IndexType>(end-begin); index_iter++)
130 ScalarType kernel_value = callback.kernel(begin[index_iter],begin[index_iter]);
134 dots[i] = callback.kernel(begin[index_iter], begin[current_neighbors[i]]);
139 gram_matrix(i,j) = kernel_value - dots(i) - dots(j) +
140 callback.kernel(begin[current_neighbors[i]],begin[current_neighbors[j]]);
144 gram_matrix.diagonal().array() += trace_shift*trace;
145 weights = gram_matrix.selfadjointView<Eigen::Upper>().ldlt().solve(rhs);
146 weights /= weights.sum();
148 SparseTriplet diagonal_triplet(index_iter,index_iter,1.0+shift);
149 local_triplets.push_back(diagonal_triplet);
152 SparseTriplet row_side_triplet(current_neighbors[i],index_iter,-weights[i]);
153 SparseTriplet col_side_triplet(index_iter,current_neighbors[i],-weights[i]);
154 local_triplets.push_back(row_side_triplet);
155 local_triplets.push_back(col_side_triplet);
158 SparseTriplet cross_triplet(current_neighbors[i],current_neighbors[j],weights(i)*weights(j));
159 local_triplets.push_back(cross_triplet);
165 copy(local_triplets.begin(),local_triplets.end(),back_inserter(sparse_triplets));
168 local_triplets.clear();
176 template <
class RandomAccessIterator,
class PairwiseCallback>
178 const Neighbors& neighbors, PairwiseCallback callback,
185 sparse_triplets.reserve(k*k*(end-begin));
187 const IndexType dp = target_dimension*(target_dimension+1)/2;
189 #pragma omp parallel shared(begin,end,neighbors,callback,sparse_triplets) default(none)
196 local_triplets.reserve(k*k+2*k+1);
198 #pragma omp for nowait
199 for (index_iter=0; index_iter<static_cast<IndexType>(end-begin); index_iter++)
207 ScalarType kij = callback.kernel(begin[current_neighbors[i]],begin[current_neighbors[j]]);
208 gram_matrix(i,j) = kij;
209 gram_matrix(j,i) = kij;
216 sae_solver.compute(gram_matrix);
218 Yi.col(0).setConstant(1.0);
219 Yi.block(0,1,k,target_dimension).noalias() = sae_solver.eigenvectors().rightCols(target_dimension);
222 for (
IndexType j=0; j<target_dimension; ++j)
224 for (
IndexType p=0; p<target_dimension-j; ++p)
226 Yi.col(ct+p+1+target_dimension).noalias() = Yi.col(j+1).cwiseProduct(Yi.col(j+p+1));
228 ct += ct + target_dimension - j;
231 for (
IndexType i=0; i<static_cast<IndexType>(Yi.cols()); i++)
236 Yi.col(i) -= r*Yi.col(j);
239 Yi.col(i) *= (1.f / norm);
243 ScalarType colsum = Yi.col(1+target_dimension+i).sum();
245 Yi.col(1+target_dimension+i).array() /= colsum;
249 gram_matrix.noalias() = Yi.rightCols(dp)*(Yi.rightCols(dp).transpose());
255 SparseTriplet hessian_triplet(current_neighbors[i],current_neighbors[j],gram_matrix(i,j));
256 local_triplets.push_back(hessian_triplet);
262 copy(local_triplets.begin(),local_triplets.end(),back_inserter(sparse_triplets));
265 local_triplets.clear();
272 template<
class RandomAccessIterator,
class FeatureVectorCallback>
274 RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback feature_vector_callback,
286 for (RandomAccessIterator iter=begin; iter!=end; ++iter)
288 feature_vector_callback.vector(*iter,rank_update_vector_i);
289 rhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i);
292 for (
int i=0; i<W.outerSize(); ++i)
294 for (SparseWeightMatrix::InnerIterator it(W,i); it; ++it)
296 feature_vector_callback.vector(begin[it.row()],rank_update_vector_i);
297 feature_vector_callback.vector(begin[it.col()],rank_update_vector_j);
298 lhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i, rank_update_vector_j, it.value());
302 rhs += rhs.transpose().eval();
310 template<
class RandomAccessIterator,
class FeatureVectorCallback>
312 RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback feature_vector_callback,
325 for (RandomAccessIterator iter=begin; iter!=end; ++iter)
327 feature_vector_callback.vector(*iter,rank_update_vector_i);
328 sum += rank_update_vector_i;
329 rhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i);
331 rhs.selfadjointView<Eigen::Upper>().rankUpdate(sum,-1./(end-begin));
333 for (
int i=0; i<W.outerSize(); ++i)
335 for (SparseWeightMatrix::InnerIterator it(W,i); it; ++it)
337 feature_vector_callback.vector(begin[it.row()],rank_update_vector_i);
338 feature_vector_callback.vector(begin[it.col()],rank_update_vector_j);
339 lhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i, rank_update_vector_j, it.value());
342 lhs.selfadjointView<Eigen::Upper>().rankUpdate(sum,-1./(end-begin));
344 rhs += rhs.transpose().eval();
void centerMatrix(DenseMatrix &matrix)
Eigen::Matrix< tapkee::ScalarType, Eigen::Dynamic, Eigen::Dynamic > DenseMatrix
dense matrix type (non-overridable)
DenseSymmetricMatrixPair construct_lltsa_eigenproblem(SparseWeightMatrix W, RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback feature_vector_callback, IndexType dimension)
TAPKEE_INTERNAL_VECTOR< tapkee::tapkee_internal::SparseTriplet > SparseTriplets
SparseWeightMatrix hessian_weight_matrix(RandomAccessIterator begin, RandomAccessIterator end, const Neighbors &neighbors, PairwiseCallback callback, const IndexType target_dimension)
ARPACK-based method (requires the ARPACK library binaries to be available around). Recommended to be used as a default method. Supports both generalized and standard eigenproblems.
double ScalarType
default scalar value (can be overrided with TAPKEE_CUSTOM_INTERNAL_NUMTYPE define) ...
TAPKEE_INTERNAL_PAIR< tapkee::DenseSymmetricMatrix, tapkee::DenseSymmetricMatrix > DenseSymmetricMatrixPair
DenseSymmetricMatrixPair construct_neighborhood_preserving_eigenproblem(SparseWeightMatrix W, RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback feature_vector_callback, IndexType dimension)
Eigen::SparseMatrix< tapkee::ScalarType > SparseWeightMatrix
sparse weight matrix type (non-overridable)
Eigen::Matrix< tapkee::ScalarType, Eigen::Dynamic, 1 > DenseVector
dense vector type (non-overridable)
TAPKEE_INTERNAL_VECTOR< tapkee::IndexType > LocalNeighbors
SparseMatrix sparse_matrix_from_triplets(const SparseTriplets &sparse_triplets, IndexType m, IndexType n)
int IndexType
indexing type (non-overridable) set to int for compatibility with OpenMP 2.0
Eigen::Triplet< tapkee::ScalarType > SparseTriplet
Eigen::SelfAdjointEigenSolver< tapkee::DenseMatrix > DenseSelfAdjointEigenSolver
selfadjoint solver (non-overridable)
TAPKEE_INTERNAL_VECTOR< tapkee::tapkee_internal::LocalNeighbors > Neighbors
SparseWeightMatrix linear_weight_matrix(const RandomAccessIterator &begin, const RandomAccessIterator &end, const Neighbors &neighbors, PairwiseCallback callback, const ScalarType shift, const ScalarType trace_shift)
tapkee::DenseMatrix DenseSymmetricMatrix
dense symmetric matrix (non-overridable, currently just dense matrix, can be improved later) ...
SparseWeightMatrix tangent_weight_matrix(RandomAccessIterator begin, RandomAccessIterator end, const Neighbors &neighbors, PairwiseCallback callback, const IndexType target_dimension, const ScalarType shift, const bool partial_eigendecomposer=false)