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 {
29  timed_context context("KLTSA weight matrix computation");
30  const IndexType k = neighbors[0].size();
31 
32  SparseTriplets sparse_triplets;
33  sparse_triplets.reserve((k*k+2*k+1)*(end-begin));
34 
35 #pragma omp parallel shared(begin,end,neighbors,callback,sparse_triplets) default(none)
36  {
37  IndexType index_iter;
38  DenseMatrix gram_matrix = DenseMatrix::Zero(k,k);
39  DenseVector rhs = DenseVector::Ones(k);
40  DenseMatrix G = DenseMatrix::Zero(k,target_dimension+1);
41  G.col(0).setConstant(1/sqrt(static_cast<ScalarType>(k)));
43  SparseTriplets local_triplets;
44  local_triplets.reserve(k*k+2*k+1);
45 
46 #pragma omp for nowait
47  for (index_iter=0; index_iter<static_cast<IndexType>(end-begin); index_iter++)
48  {
49  const LocalNeighbors& current_neighbors = neighbors[index_iter];
50 
51  for (IndexType i=0; i<k; ++i)
52  {
53  for (IndexType j=i; j<k; ++j)
54  {
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;
58  }
59  }
60 
61  centerMatrix(gram_matrix);
62 
63  //UNRESTRICT_ALLOC;
64  solver.compute(gram_matrix);
65  G.rightCols(target_dimension).noalias() = solver.eigenvectors().rightCols(target_dimension);
66  //RESTRICT_ALLOC;
67  gram_matrix.noalias() = G * G.transpose();
68 
69  SparseTriplet diagonal_triplet(index_iter,index_iter,shift);
70  local_triplets.push_back(diagonal_triplet);
71  for (IndexType i=0; i<k; ++i)
72  {
73  SparseTriplet neighborhood_diagonal_triplet(current_neighbors[i],current_neighbors[i],1.0);
74  local_triplets.push_back(neighborhood_diagonal_triplet);
75 
76  for (IndexType j=0; j<k; ++j)
77  {
78  SparseTriplet tangent_triplet(current_neighbors[i],current_neighbors[j],-gram_matrix(i,j));
79  local_triplets.push_back(tangent_triplet);
80  }
81  }
82 #pragma omp critical
83  {
84  copy(local_triplets.begin(),local_triplets.end(),std::back_inserter(sparse_triplets));
85  }
86 
87  local_triplets.clear();
88  }
89  }
90 
91  return sparse_matrix_from_triplets(sparse_triplets, end-begin, end-begin);
92 }
93 
94 template <class RandomAccessIterator, class PairwiseCallback>
95 SparseWeightMatrix linear_weight_matrix(const RandomAccessIterator& begin, const RandomAccessIterator& end,
96  const Neighbors& neighbors, PairwiseCallback callback,
97  const ScalarType shift, const ScalarType trace_shift)
98 {
99  timed_context context("KLLE weight computation");
100  const IndexType k = neighbors[0].size();
101 
102  SparseTriplets sparse_triplets;
103  sparse_triplets.reserve((k*k+2*k+1)*(end-begin));
104 
105 #pragma omp parallel shared(begin,end,neighbors,callback,sparse_triplets) default(none)
106  {
107  IndexType index_iter;
108  DenseMatrix gram_matrix = DenseMatrix::Zero(k,k);
109  DenseVector dots(k);
110  DenseVector rhs = DenseVector::Ones(k);
111  DenseVector weights;
112  SparseTriplets local_triplets;
113  local_triplets.reserve(k*k+2*k+1);
114 
115  //RESTRICT_ALLOC;
116 #pragma omp for nowait
117  for (index_iter=0; index_iter<static_cast<IndexType>(end-begin); index_iter++)
118  {
119  ScalarType kernel_value = callback.kernel(begin[index_iter],begin[index_iter]);
120  const LocalNeighbors& current_neighbors = neighbors[index_iter];
121 
122  for (IndexType i=0; i<k; ++i)
123  dots[i] = callback.kernel(begin[index_iter], begin[current_neighbors[i]]);
124 
125  for (IndexType i=0; i<k; ++i)
126  {
127  for (IndexType j=i; j<k; ++j)
128  gram_matrix(i,j) = kernel_value - dots(i) - dots(j) +
129  callback.kernel(begin[current_neighbors[i]],begin[current_neighbors[j]]);
130  }
131 
132  ScalarType trace = gram_matrix.trace();
133  gram_matrix.diagonal().array() += trace_shift*trace;
134  weights = gram_matrix.selfadjointView<Eigen::Upper>().ldlt().solve(rhs);
135  weights /= weights.sum();
136 
137  SparseTriplet diagonal_triplet(index_iter,index_iter,1.0+shift);
138  local_triplets.push_back(diagonal_triplet);
139  for (IndexType i=0; i<k; ++i)
140  {
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);
145  for (IndexType j=0; j<k; ++j)
146  {
147  SparseTriplet cross_triplet(current_neighbors[i],current_neighbors[j],weights(i)*weights(j));
148  local_triplets.push_back(cross_triplet);
149  }
150  }
151 
152 #pragma omp critical
153  {
154  copy(local_triplets.begin(),local_triplets.end(),std::back_inserter(sparse_triplets));
155  }
156 
157  local_triplets.clear();
158  }
159  //UNRESTRICT_ALLOC;
160  }
161 
162  return sparse_matrix_from_triplets(sparse_triplets, end-begin, end-begin);
163 }
164 
165 template <class RandomAccessIterator, class PairwiseCallback>
166 SparseWeightMatrix hessian_weight_matrix(RandomAccessIterator begin, RandomAccessIterator end,
167  const Neighbors& neighbors, PairwiseCallback callback,
168  const IndexType target_dimension)
169 {
170  timed_context context("Hessian weight matrix computation");
171  const IndexType k = neighbors[0].size();
172 
173  SparseTriplets sparse_triplets;
174  sparse_triplets.reserve(k*k*(end-begin));
175 
176  const IndexType dp = target_dimension*(target_dimension+1)/2;
177 
178 #pragma omp parallel shared(begin,end,neighbors,callback,sparse_triplets) default(none)
179  {
180  IndexType index_iter;
181  DenseMatrix gram_matrix = DenseMatrix::Zero(k,k);
182  DenseMatrix Yi(k,1+target_dimension+dp);
183 
184  SparseTriplets local_triplets;
185  local_triplets.reserve(k*k+2*k+1);
186 
187 #pragma omp for nowait
188  for (index_iter=0; index_iter<static_cast<IndexType>(end-begin); index_iter++)
189  {
190  const LocalNeighbors& current_neighbors = neighbors[index_iter];
191 
192  for (IndexType i=0; i<k; ++i)
193  {
194  for (IndexType j=i; j<k; ++j)
195  {
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;
199  }
200  }
201 
202  centerMatrix(gram_matrix);
203 
204  DenseSelfAdjointEigenSolver sae_solver;
205  sae_solver.compute(gram_matrix);
206 
207  Yi.col(0).setConstant(1.0);
208  Yi.block(0,1,k,target_dimension).noalias() = sae_solver.eigenvectors().rightCols(target_dimension);
209 
210  IndexType ct = 0;
211  for (IndexType j=0; j<target_dimension; ++j)
212  {
213  for (IndexType p=0; p<target_dimension-j; ++p)
214  {
215  Yi.col(ct+p+1+target_dimension).noalias() = Yi.col(j+1).cwiseProduct(Yi.col(j+p+1));
216  }
217  ct += ct + target_dimension - j;
218  }
219 
220  for (IndexType i=0; i<static_cast<IndexType>(Yi.cols()); i++)
221  {
222  for (IndexType j=0; j<i; j++)
223  {
224  ScalarType r = Yi.col(i).dot(Yi.col(j));
225  Yi.col(i) -= r*Yi.col(j);
226  }
227  ScalarType norm = Yi.col(i).norm();
228  Yi.col(i) *= (1.f / norm);
229  }
230  for (IndexType i=0; i<dp; i++)
231  {
232  ScalarType colsum = Yi.col(1+target_dimension+i).sum();
233  if (colsum > 1e-4)
234  Yi.col(1+target_dimension+i).array() /= colsum;
235  }
236 
237  // reuse gram matrix storage m'kay?
238  gram_matrix.noalias() = Yi.rightCols(dp)*(Yi.rightCols(dp).transpose());
239 
240  for (IndexType i=0; i<k; ++i)
241  {
242  for (IndexType j=0; j<k; ++j)
243  {
244  SparseTriplet hessian_triplet(current_neighbors[i],current_neighbors[j],gram_matrix(i,j));
245  local_triplets.push_back(hessian_triplet);
246  }
247  }
248 
249  #pragma omp critical
250  {
251  copy(local_triplets.begin(),local_triplets.end(),std::back_inserter(sparse_triplets));
252  }
253 
254  local_triplets.clear();
255  }
256  }
257 
258  return sparse_matrix_from_triplets(sparse_triplets, end-begin, end-begin);
259 }
260 
261 template<class RandomAccessIterator, class FeatureVectorCallback>
263  RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback feature_vector_callback,
264  IndexType dimension)
265 {
266  timed_context context("NPE eigenproblem construction");
267 
268  DenseSymmetricMatrix lhs = DenseSymmetricMatrix::Zero(dimension,dimension);
269  DenseSymmetricMatrix rhs = DenseSymmetricMatrix::Zero(dimension,dimension);
270 
271  DenseVector rank_update_vector_i(dimension);
272  DenseVector rank_update_vector_j(dimension);
273 
274  //RESTRICT_ALLOC;
275  for (RandomAccessIterator iter=begin; iter!=end; ++iter)
276  {
277  feature_vector_callback.vector(*iter,rank_update_vector_i);
278  rhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i);
279  }
280 
281  for (int i=0; i<W.outerSize(); ++i)
282  {
283  for (SparseWeightMatrix::InnerIterator it(W,i); it; ++it)
284  {
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());
288  }
289  }
290 
291  rhs += rhs.transpose().eval();
292  rhs /= 2;
293 
294  //UNRESTRICT_ALLOC;
295 
296  return DenseSymmetricMatrixPair(lhs,rhs);
297 }
298 
299 template<class RandomAccessIterator, class FeatureVectorCallback>
301  RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback feature_vector_callback,
302  IndexType dimension)
303 {
304  timed_context context("LLTSA eigenproblem construction");
305 
306  DenseSymmetricMatrix lhs = DenseSymmetricMatrix::Zero(dimension,dimension);
307  DenseSymmetricMatrix rhs = DenseSymmetricMatrix::Zero(dimension,dimension);
308 
309  DenseVector rank_update_vector_i(dimension);
310  DenseVector rank_update_vector_j(dimension);
311  DenseVector sum = DenseVector::Zero(dimension);
312 
313  //RESTRICT_ALLOC;
314  for (RandomAccessIterator iter=begin; iter!=end; ++iter)
315  {
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);
319  }
320  rhs.selfadjointView<Eigen::Upper>().rankUpdate(sum,-1./(end-begin));
321 
322  for (int i=0; i<W.outerSize(); ++i)
323  {
324  for (SparseWeightMatrix::InnerIterator it(W,i); it; ++it)
325  {
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());
329  }
330  }
331  lhs.selfadjointView<Eigen::Upper>().rankUpdate(sum,-1./(end-begin));
332 
333  rhs += rhs.transpose().eval();
334  rhs /= 2;
335 
336  //UNRESTRICT_ALLOC;
337 
338  return DenseSymmetricMatrixPair(lhs,rhs);
339 }
340 
341 } // End of namespace tapkee_internal
342 } // End of namespace tapkee
343 
344 #endif
void centerMatrix(DenseMatrix &matrix)
Definition: matrix.hpp:14
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)
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)
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)
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