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,
33 sparse_triplets.reserve((k*k+2*k+1)*(end-begin));
35 #pragma omp parallel shared(begin,end,neighbors,callback,sparse_triplets) default(none)
40 DenseMatrix G = DenseMatrix::Zero(k,target_dimension+1);
41 G.col(0).setConstant(1/sqrt(static_cast<ScalarType>(k)));
44 local_triplets.reserve(k*k+2*k+1);
46 #pragma omp for nowait
47 for (index_iter=0; index_iter<static_cast<IndexType>(end-begin); index_iter++)
55 ScalarType kij = callback.kernel(begin[current_neighbors[i]],begin[current_neighbors[j]]);
56 gram_matrix(i,j) = kij;
57 gram_matrix(j,i) = kij;
64 solver.compute(gram_matrix);
65 G.rightCols(target_dimension).noalias() = solver.eigenvectors().rightCols(target_dimension);
67 gram_matrix.noalias() = G * G.transpose();
70 local_triplets.push_back(diagonal_triplet);
73 SparseTriplet neighborhood_diagonal_triplet(current_neighbors[i],current_neighbors[i],1.0);
74 local_triplets.push_back(neighborhood_diagonal_triplet);
78 SparseTriplet tangent_triplet(current_neighbors[i],current_neighbors[j],-gram_matrix(i,j));
79 local_triplets.push_back(tangent_triplet);
84 copy(local_triplets.begin(),local_triplets.end(),std::back_inserter(sparse_triplets));
87 local_triplets.clear();
94 template <
class RandomAccessIterator,
class PairwiseCallback>
96 const Neighbors& neighbors, PairwiseCallback callback,
103 sparse_triplets.reserve((k*k+2*k+1)*(end-begin));
105 #pragma omp parallel shared(begin,end,neighbors,callback,sparse_triplets) default(none)
113 local_triplets.reserve(k*k+2*k+1);
116 #pragma omp for nowait
117 for (index_iter=0; index_iter<static_cast<IndexType>(end-begin); index_iter++)
119 ScalarType kernel_value = callback.kernel(begin[index_iter],begin[index_iter]);
123 dots[i] = callback.kernel(begin[index_iter], begin[current_neighbors[i]]);
128 gram_matrix(i,j) = kernel_value - dots(i) - dots(j) +
129 callback.kernel(begin[current_neighbors[i]],begin[current_neighbors[j]]);
133 gram_matrix.diagonal().array() += trace_shift*trace;
134 weights = gram_matrix.selfadjointView<Eigen::Upper>().ldlt().solve(rhs);
135 weights /= weights.sum();
137 SparseTriplet diagonal_triplet(index_iter,index_iter,1.0+shift);
138 local_triplets.push_back(diagonal_triplet);
141 SparseTriplet row_side_triplet(current_neighbors[i],index_iter,-weights[i]);
142 SparseTriplet col_side_triplet(index_iter,current_neighbors[i],-weights[i]);
143 local_triplets.push_back(row_side_triplet);
144 local_triplets.push_back(col_side_triplet);
147 SparseTriplet cross_triplet(current_neighbors[i],current_neighbors[j],weights(i)*weights(j));
148 local_triplets.push_back(cross_triplet);
154 copy(local_triplets.begin(),local_triplets.end(),std::back_inserter(sparse_triplets));
157 local_triplets.clear();
165 template <
class RandomAccessIterator,
class PairwiseCallback>
167 const Neighbors& neighbors, PairwiseCallback callback,
174 sparse_triplets.reserve(k*k*(end-begin));
176 const IndexType dp = target_dimension*(target_dimension+1)/2;
178 #pragma omp parallel shared(begin,end,neighbors,callback,sparse_triplets) default(none)
185 local_triplets.reserve(k*k+2*k+1);
187 #pragma omp for nowait
188 for (index_iter=0; index_iter<static_cast<IndexType>(end-begin); index_iter++)
196 ScalarType kij = callback.kernel(begin[current_neighbors[i]],begin[current_neighbors[j]]);
197 gram_matrix(i,j) = kij;
198 gram_matrix(j,i) = kij;
205 sae_solver.compute(gram_matrix);
207 Yi.col(0).setConstant(1.0);
208 Yi.block(0,1,k,target_dimension).noalias() = sae_solver.eigenvectors().rightCols(target_dimension);
211 for (
IndexType j=0; j<target_dimension; ++j)
213 for (
IndexType p=0; p<target_dimension-j; ++p)
215 Yi.col(ct+p+1+target_dimension).noalias() = Yi.col(j+1).cwiseProduct(Yi.col(j+p+1));
217 ct += ct + target_dimension - j;
220 for (
IndexType i=0; i<static_cast<IndexType>(Yi.cols()); i++)
225 Yi.col(i) -= r*Yi.col(j);
228 Yi.col(i) *= (1.f / norm);
232 ScalarType colsum = Yi.col(1+target_dimension+i).sum();
234 Yi.col(1+target_dimension+i).array() /= colsum;
238 gram_matrix.noalias() = Yi.rightCols(dp)*(Yi.rightCols(dp).transpose());
244 SparseTriplet hessian_triplet(current_neighbors[i],current_neighbors[j],gram_matrix(i,j));
245 local_triplets.push_back(hessian_triplet);
251 copy(local_triplets.begin(),local_triplets.end(),std::back_inserter(sparse_triplets));
254 local_triplets.clear();
261 template<
class RandomAccessIterator,
class FeatureVectorCallback>
263 RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback feature_vector_callback,
275 for (RandomAccessIterator iter=begin; iter!=end; ++iter)
277 feature_vector_callback.vector(*iter,rank_update_vector_i);
278 rhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i);
281 for (
int i=0; i<W.outerSize(); ++i)
283 for (SparseWeightMatrix::InnerIterator it(W,i); it; ++it)
285 feature_vector_callback.vector(begin[it.row()],rank_update_vector_i);
286 feature_vector_callback.vector(begin[it.col()],rank_update_vector_j);
287 lhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i, rank_update_vector_j, it.value());
291 rhs += rhs.transpose().eval();
299 template<
class RandomAccessIterator,
class FeatureVectorCallback>
301 RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback feature_vector_callback,
314 for (RandomAccessIterator iter=begin; iter!=end; ++iter)
316 feature_vector_callback.vector(*iter,rank_update_vector_i);
317 sum += rank_update_vector_i;
318 rhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i);
320 rhs.selfadjointView<Eigen::Upper>().rankUpdate(sum,-1./(end-begin));
322 for (
int i=0; i<W.outerSize(); ++i)
324 for (SparseWeightMatrix::InnerIterator it(W,i); it; ++it)
326 feature_vector_callback.vector(begin[it.row()],rank_update_vector_i);
327 feature_vector_callback.vector(begin[it.col()],rank_update_vector_j);
328 lhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i, rank_update_vector_j, it.value());
331 lhs.selfadjointView<Eigen::Upper>().rankUpdate(sum,-1./(end-begin));
333 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)
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)
SparseWeightMatrix tangent_weight_matrix(RandomAccessIterator begin, RandomAccessIterator end, const Neighbors &neighbors, PairwiseCallback callback, const IndexType target_dimension, const ScalarType shift)
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) ...