Tapkee
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
locally_linear.hpp
Go to the documentation of this file.
1 /* This software is distributed under BSD 3-clause license (see LICENSE file).
2  *
3  * Copyright (c) 2012-2013 Sergey Lisitsyn
4  */
5 
6 #ifndef TAPKEE_LOCALLY_LINEAR_H_
7 #define TAPKEE_LOCALLY_LINEAR_H_
8 
9 /* Tapkee includes */
11 #include <tapkee/defines.hpp>
12 #include <tapkee/utils/matrix.hpp>
13 #include <tapkee/utils/time.hpp>
14 #include <tapkee/utils/sparse.hpp>
15 /* End of Tapkee includes */
16 
17 namespace tapkee
18 {
19 namespace tapkee_internal
20 {
21 
22 
23 
24 template <class RandomAccessIterator, class PairwiseCallback>
25 SparseWeightMatrix tangent_weight_matrix(RandomAccessIterator begin, RandomAccessIterator end,
26  const Neighbors& neighbors, PairwiseCallback callback,
27  const IndexType target_dimension, const ScalarType shift,
28  const bool partial_eigendecomposer=false)
29 {
30  timed_context context("KLTSA weight matrix computation");
31  const IndexType k = neighbors[0].size();
32 
33  SparseTriplets sparse_triplets;
34  sparse_triplets.reserve((k*k+2*k+1)*(end-begin));
35 
36 #pragma omp parallel shared(begin,end,neighbors,callback,sparse_triplets) default(none)
37  {
38  IndexType index_iter;
39  DenseMatrix gram_matrix = DenseMatrix::Zero(k,k);
40  DenseVector rhs = DenseVector::Ones(k);
41  DenseMatrix G = DenseMatrix::Zero(k,target_dimension+1);
42  G.col(0).setConstant(1/sqrt(static_cast<ScalarType>(k)));
44  SparseTriplets local_triplets;
45  local_triplets.reserve(k*k+2*k+1);
46 
47 #pragma omp for nowait
48  for (index_iter=0; index_iter<static_cast<IndexType>(end-begin); index_iter++)
49  {
50  const LocalNeighbors& current_neighbors = neighbors[index_iter];
51 
52  for (IndexType i=0; i<k; ++i)
53  {
54  for (IndexType j=i; j<k; ++j)
55  {
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;
59  }
60  }
61 
62  centerMatrix(gram_matrix);
63 
64  //UNRESTRICT_ALLOC;
65 #ifdef TAPKEE_WITH_ARPACK
66  if (partial_eigendecomposer)
67  {
68  G.rightCols(target_dimension).noalias() =
69  eigendecomposition<DenseMatrix,DenseMatrixOperation>(Arpack,gram_matrix,target_dimension,0).first;
70  }
71  else
72 #endif
73  {
74  solver.compute(gram_matrix);
75  G.rightCols(target_dimension).noalias() = solver.eigenvectors().rightCols(target_dimension);
76  }
77  //RESTRICT_ALLOC;
78  gram_matrix.noalias() = G * G.transpose();
79 
80  SparseTriplet diagonal_triplet(index_iter,index_iter,shift);
81  local_triplets.push_back(diagonal_triplet);
82  for (IndexType i=0; i<k; ++i)
83  {
84  SparseTriplet neighborhood_diagonal_triplet(current_neighbors[i],current_neighbors[i],1.0);
85  local_triplets.push_back(neighborhood_diagonal_triplet);
86 
87  for (IndexType j=0; j<k; ++j)
88  {
89  SparseTriplet tangent_triplet(current_neighbors[i],current_neighbors[j],-gram_matrix(i,j));
90  local_triplets.push_back(tangent_triplet);
91  }
92  }
93 #pragma omp critical
94  {
95  copy(local_triplets.begin(),local_triplets.end(),back_inserter(sparse_triplets));
96  }
97 
98  local_triplets.clear();
99  }
100  }
101 
102  return sparse_matrix_from_triplets(sparse_triplets, end-begin, end-begin);
103 }
104 
105 template <class RandomAccessIterator, class PairwiseCallback>
106 SparseWeightMatrix linear_weight_matrix(const RandomAccessIterator& begin, const RandomAccessIterator& end,
107  const Neighbors& neighbors, PairwiseCallback callback,
108  const ScalarType shift, const ScalarType trace_shift)
109 {
110  timed_context context("KLLE weight computation");
111  const IndexType k = neighbors[0].size();
112 
113  SparseTriplets sparse_triplets;
114  sparse_triplets.reserve((k*k+2*k+1)*(end-begin));
115 
116 #pragma omp parallel shared(begin,end,neighbors,callback,sparse_triplets) default(none)
117  {
118  IndexType index_iter;
119  DenseMatrix gram_matrix = DenseMatrix::Zero(k,k);
120  DenseVector dots(k);
121  DenseVector rhs = DenseVector::Ones(k);
122  DenseVector weights;
123  SparseTriplets local_triplets;
124  local_triplets.reserve(k*k+2*k+1);
125 
126  //RESTRICT_ALLOC;
127 #pragma omp for nowait
128  for (index_iter=0; index_iter<static_cast<IndexType>(end-begin); index_iter++)
129  {
130  ScalarType kernel_value = callback.kernel(begin[index_iter],begin[index_iter]);
131  const LocalNeighbors& current_neighbors = neighbors[index_iter];
132 
133  for (IndexType i=0; i<k; ++i)
134  dots[i] = callback.kernel(begin[index_iter], begin[current_neighbors[i]]);
135 
136  for (IndexType i=0; i<k; ++i)
137  {
138  for (IndexType j=i; j<k; ++j)
139  gram_matrix(i,j) = kernel_value - dots(i) - dots(j) +
140  callback.kernel(begin[current_neighbors[i]],begin[current_neighbors[j]]);
141  }
142 
143  ScalarType trace = gram_matrix.trace();
144  gram_matrix.diagonal().array() += trace_shift*trace;
145  weights = gram_matrix.selfadjointView<Eigen::Upper>().ldlt().solve(rhs);
146  weights /= weights.sum();
147 
148  SparseTriplet diagonal_triplet(index_iter,index_iter,1.0+shift);
149  local_triplets.push_back(diagonal_triplet);
150  for (IndexType i=0; i<k; ++i)
151  {
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);
156  for (IndexType j=0; j<k; ++j)
157  {
158  SparseTriplet cross_triplet(current_neighbors[i],current_neighbors[j],weights(i)*weights(j));
159  local_triplets.push_back(cross_triplet);
160  }
161  }
162 
163 #pragma omp critical
164  {
165  copy(local_triplets.begin(),local_triplets.end(),back_inserter(sparse_triplets));
166  }
167 
168  local_triplets.clear();
169  }
170  //UNRESTRICT_ALLOC;
171  }
172 
173  return sparse_matrix_from_triplets(sparse_triplets, end-begin, end-begin);
174 }
175 
176 template <class RandomAccessIterator, class PairwiseCallback>
177 SparseWeightMatrix hessian_weight_matrix(RandomAccessIterator begin, RandomAccessIterator end,
178  const Neighbors& neighbors, PairwiseCallback callback,
179  const IndexType target_dimension)
180 {
181  timed_context context("Hessian weight matrix computation");
182  const IndexType k = neighbors[0].size();
183 
184  SparseTriplets sparse_triplets;
185  sparse_triplets.reserve(k*k*(end-begin));
186 
187  const IndexType dp = target_dimension*(target_dimension+1)/2;
188 
189 #pragma omp parallel shared(begin,end,neighbors,callback,sparse_triplets) default(none)
190  {
191  IndexType index_iter;
192  DenseMatrix gram_matrix = DenseMatrix::Zero(k,k);
193  DenseMatrix Yi(k,1+target_dimension+dp);
194 
195  SparseTriplets local_triplets;
196  local_triplets.reserve(k*k+2*k+1);
197 
198 #pragma omp for nowait
199  for (index_iter=0; index_iter<static_cast<IndexType>(end-begin); index_iter++)
200  {
201  const LocalNeighbors& current_neighbors = neighbors[index_iter];
202 
203  for (IndexType i=0; i<k; ++i)
204  {
205  for (IndexType j=i; j<k; ++j)
206  {
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;
210  }
211  }
212 
213  centerMatrix(gram_matrix);
214 
215  DenseSelfAdjointEigenSolver sae_solver;
216  sae_solver.compute(gram_matrix);
217 
218  Yi.col(0).setConstant(1.0);
219  Yi.block(0,1,k,target_dimension).noalias() = sae_solver.eigenvectors().rightCols(target_dimension);
220 
221  IndexType ct = 0;
222  for (IndexType j=0; j<target_dimension; ++j)
223  {
224  for (IndexType p=0; p<target_dimension-j; ++p)
225  {
226  Yi.col(ct+p+1+target_dimension).noalias() = Yi.col(j+1).cwiseProduct(Yi.col(j+p+1));
227  }
228  ct += ct + target_dimension - j;
229  }
230 
231  for (IndexType i=0; i<static_cast<IndexType>(Yi.cols()); i++)
232  {
233  for (IndexType j=0; j<i; j++)
234  {
235  ScalarType r = Yi.col(i).dot(Yi.col(j));
236  Yi.col(i) -= r*Yi.col(j);
237  }
238  ScalarType norm = Yi.col(i).norm();
239  Yi.col(i) *= (1.f / norm);
240  }
241  for (IndexType i=0; i<dp; i++)
242  {
243  ScalarType colsum = Yi.col(1+target_dimension+i).sum();
244  if (colsum > 1e-4)
245  Yi.col(1+target_dimension+i).array() /= colsum;
246  }
247 
248  // reuse gram matrix storage m'kay?
249  gram_matrix.noalias() = Yi.rightCols(dp)*(Yi.rightCols(dp).transpose());
250 
251  for (IndexType i=0; i<k; ++i)
252  {
253  for (IndexType j=0; j<k; ++j)
254  {
255  SparseTriplet hessian_triplet(current_neighbors[i],current_neighbors[j],gram_matrix(i,j));
256  local_triplets.push_back(hessian_triplet);
257  }
258  }
259 
260  #pragma omp critical
261  {
262  copy(local_triplets.begin(),local_triplets.end(),back_inserter(sparse_triplets));
263  }
264 
265  local_triplets.clear();
266  }
267  }
268 
269  return sparse_matrix_from_triplets(sparse_triplets, end-begin, end-begin);
270 }
271 
272 template<class RandomAccessIterator, class FeatureVectorCallback>
274  RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback feature_vector_callback,
275  IndexType dimension)
276 {
277  timed_context context("NPE eigenproblem construction");
278 
279  DenseSymmetricMatrix lhs = DenseSymmetricMatrix::Zero(dimension,dimension);
280  DenseSymmetricMatrix rhs = DenseSymmetricMatrix::Zero(dimension,dimension);
281 
282  DenseVector rank_update_vector_i(dimension);
283  DenseVector rank_update_vector_j(dimension);
284 
285  //RESTRICT_ALLOC;
286  for (RandomAccessIterator iter=begin; iter!=end; ++iter)
287  {
288  feature_vector_callback.vector(*iter,rank_update_vector_i);
289  rhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i);
290  }
291 
292  for (int i=0; i<W.outerSize(); ++i)
293  {
294  for (SparseWeightMatrix::InnerIterator it(W,i); it; ++it)
295  {
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());
299  }
300  }
301 
302  rhs += rhs.transpose().eval();
303  rhs /= 2;
304 
305  //UNRESTRICT_ALLOC;
306 
307  return DenseSymmetricMatrixPair(lhs,rhs);
308 }
309 
310 template<class RandomAccessIterator, class FeatureVectorCallback>
312  RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback feature_vector_callback,
313  IndexType dimension)
314 {
315  timed_context context("LLTSA eigenproblem construction");
316 
317  DenseSymmetricMatrix lhs = DenseSymmetricMatrix::Zero(dimension,dimension);
318  DenseSymmetricMatrix rhs = DenseSymmetricMatrix::Zero(dimension,dimension);
319 
320  DenseVector rank_update_vector_i(dimension);
321  DenseVector rank_update_vector_j(dimension);
322  DenseVector sum = DenseVector::Zero(dimension);
323 
324  //RESTRICT_ALLOC;
325  for (RandomAccessIterator iter=begin; iter!=end; ++iter)
326  {
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);
330  }
331  rhs.selfadjointView<Eigen::Upper>().rankUpdate(sum,-1./(end-begin));
332 
333  for (int i=0; i<W.outerSize(); ++i)
334  {
335  for (SparseWeightMatrix::InnerIterator it(W,i); it; ++it)
336  {
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());
340  }
341  }
342  lhs.selfadjointView<Eigen::Upper>().rankUpdate(sum,-1./(end-begin));
343 
344  rhs += rhs.transpose().eval();
345  rhs /= 2;
346 
347  //UNRESTRICT_ALLOC;
348 
349  return DenseSymmetricMatrixPair(lhs,rhs);
350 }
351 
352 } // End of namespace tapkee_internal
353 } // End of namespace tapkee
354 
355 #endif
void centerMatrix(DenseMatrix &matrix)
Definition: matrix.hpp:9
Eigen::Matrix< tapkee::ScalarType, Eigen::Dynamic, Eigen::Dynamic > DenseMatrix
dense matrix type (non-overridable)
Definition: types.hpp:23
DenseSymmetricMatrixPair construct_lltsa_eigenproblem(SparseWeightMatrix W, RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback feature_vector_callback, IndexType dimension)
TAPKEE_INTERNAL_VECTOR< tapkee::tapkee_internal::SparseTriplet > SparseTriplets
Definition: synonyms.hpp:38
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) ...
Definition: types.hpp:15
TAPKEE_INTERNAL_PAIR< tapkee::DenseSymmetricMatrix, tapkee::DenseSymmetricMatrix > DenseSymmetricMatrixPair
Definition: synonyms.hpp:44
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)
Definition: types.hpp:29
Eigen::Matrix< tapkee::ScalarType, Eigen::Dynamic, 1 > DenseVector
dense vector type (non-overridable)
Definition: types.hpp:21
TAPKEE_INTERNAL_VECTOR< tapkee::IndexType > LocalNeighbors
Definition: synonyms.hpp:39
SparseMatrix sparse_matrix_from_triplets(const SparseTriplets &sparse_triplets, IndexType m, IndexType n)
Definition: sparse.hpp:18
int IndexType
indexing type (non-overridable) set to int for compatibility with OpenMP 2.0
Definition: types.hpp:19
Eigen::Triplet< tapkee::ScalarType > SparseTriplet
Definition: synonyms.hpp:35
Eigen::SelfAdjointEigenSolver< tapkee::DenseMatrix > DenseSelfAdjointEigenSolver
selfadjoint solver (non-overridable)
Definition: types.hpp:33
TAPKEE_INTERNAL_VECTOR< tapkee::tapkee_internal::LocalNeighbors > Neighbors
Definition: synonyms.hpp:40
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) ...
Definition: types.hpp:25
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)