Tapkee
eigendecomposition.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  * Randomized eigendecomposition code is inspired by the redsvd library
6  * code which is also distributed under BSD 3-clause license.
7  *
8  * Copyright (c) 2010-2013 Daisuke Okanohara
9  *
10  */
11 
12 #ifndef TAPKEE_EIGENDECOMPOSITION_H_
13 #define TAPKEE_EIGENDECOMPOSITION_H_
14 
15 /* Tapkee includes */
16 #ifdef TAPKEE_WITH_ARPACK
18 #endif
20 #include <tapkee/defines.hpp>
21 /* End of Tapkee includes */
22 
23 namespace tapkee
24 {
25 namespace tapkee_internal
26 {
27 
28 #ifdef TAPKEE_WITH_ARPACK
29 template <class MatrixType, class MatrixOperationType>
31 EigendecompositionResult eigendecomposition_impl_arpack(const MatrixType& wm, IndexType target_dimension, unsigned int skip)
32 {
33  timed_context context("ARPACK eigendecomposition");
34 
36  arpack(wm,target_dimension+skip,MatrixOperationType::ARPACK_CODE);
37 
38  if (arpack.info() == Eigen::Success)
39  {
40  std::string message = formatting::format("Took {} iterations.", arpack.getNbrIterations());
42  DenseMatrix selected_eigenvectors = arpack.eigenvectors().rightCols(target_dimension);
43  return EigendecompositionResult(selected_eigenvectors,arpack.eigenvalues().tail(target_dimension));
44  }
45  else
46  {
47  throw eigendecomposition_error("eigendecomposition failed");
48  }
49  return EigendecompositionResult();
50 }
51 #endif
52 
54 template <class MatrixType, class MatrixOperationType>
55 EigendecompositionResult eigendecomposition_impl_dense(const MatrixType& wm, IndexType target_dimension, unsigned int skip)
56 {
57  timed_context context("Eigen library dense eigendecomposition");
58 
59  DenseSymmetricMatrix dense_wm = wm;
60  DenseSelfAdjointEigenSolver solver(dense_wm.selfadjointView<Eigen::Upper>());
61 
62  if (solver.info() == Eigen::Success)
63  {
64  if (MatrixOperationType::largest)
65  {
66  assert(skip==0);
67  DenseMatrix selected_eigenvectors = solver.eigenvectors().rightCols(target_dimension);
68  return EigendecompositionResult(selected_eigenvectors,solver.eigenvalues().tail(target_dimension));
69  }
70  else
71  {
72  DenseMatrix selected_eigenvectors = solver.eigenvectors().leftCols(target_dimension+skip).rightCols(target_dimension);
73  return EigendecompositionResult(selected_eigenvectors,solver.eigenvalues().segment(skip,skip+target_dimension));
74  }
75  }
76  else
77  {
78  throw eigendecomposition_error("eigendecomposition failed");
79  }
80  return EigendecompositionResult();
81 }
82 
84 template <class MatrixType, class MatrixOperationType>
85 EigendecompositionResult eigendecomposition_impl_randomized(const MatrixType& wm, IndexType target_dimension, unsigned int skip)
86 {
87  timed_context context("Randomized eigendecomposition");
88 
89  DenseMatrix O(wm.rows(), target_dimension+skip);
90  for (IndexType i=0; i<O.rows(); ++i)
91  {
92  for (IndexType j=0; j<O.cols(); j++)
93  {
94  O(i,j) = tapkee::gaussian_random();
95  }
96  }
97  MatrixOperationType operation(wm);
98 
99  DenseMatrix Y = operation(O);
100  for (IndexType i=0; i<Y.cols(); i++)
101  {
102  for (IndexType j=0; j<i; j++)
103  {
104  ScalarType r = Y.col(i).dot(Y.col(j));
105  Y.col(i) -= r*Y.col(j);
106  }
107  ScalarType norm = Y.col(i).norm();
108  if (norm < 1e-4)
109  {
110  for (int k = i; k<Y.cols(); k++)
111  Y.col(k).setZero();
112  }
113  Y.col(i) *= (1.f / norm);
114  }
115 
116  DenseMatrix B1 = operation(Y);
117  DenseMatrix B = Y.householderQr().solve(B1);
118  DenseSelfAdjointEigenSolver eigenOfB(B);
119 
120  if (eigenOfB.info() == Eigen::Success)
121  {
122  if (MatrixOperationType::largest)
123  {
124  assert(skip==0);
125  DenseMatrix selected_eigenvectors = (Y*eigenOfB.eigenvectors()).rightCols(target_dimension);
126  return EigendecompositionResult(selected_eigenvectors,eigenOfB.eigenvalues());
127  }
128  else
129  {
130  DenseMatrix selected_eigenvectors = (Y*eigenOfB.eigenvectors()).leftCols(target_dimension+skip).rightCols(target_dimension);
131  return EigendecompositionResult(selected_eigenvectors,eigenOfB.eigenvalues());
132  }
133  }
134  else
135  {
136  throw eigendecomposition_error("eigendecomposition failed");
137  }
138  return EigendecompositionResult();
139 }
140 
141 template <typename MatrixType>
143 {
144 #ifdef TAPKEE_WITH_ARPACK
145  EigendecompositionResult arpack(const MatrixType& m, const ComputationStrategy& strategy,
146  const EigendecompositionStrategy& eigen_strategy,
147  IndexType target_dimension);
148 #endif
149  EigendecompositionResult dense(const MatrixType& m, const ComputationStrategy& strategy,
150  const EigendecompositionStrategy& eigen_strategy,
151  IndexType target_dimension);
152  EigendecompositionResult randomized(const MatrixType& m, const ComputationStrategy& strategy,
153  const EigendecompositionStrategy& eigen_strategy,
154  IndexType target_dimension);
155 };
156 
157 template <>
159 {
160 #ifdef TAPKEE_WITH_ARPACK
162  const EigendecompositionStrategy& eigen_strategy,
163  IndexType target_dimension)
164  {
165  if (strategy.is(HomogeneousCPUStrategy))
166  {
167  if (eigen_strategy.is(LargestEigenvalues))
168  return eigendecomposition_impl_arpack<DenseMatrix,DenseMatrixOperation>
169  (m,target_dimension,eigen_strategy.skip());
170  if (eigen_strategy.is(SquaredLargestEigenvalues))
171  return eigendecomposition_impl_arpack<DenseMatrix,DenseImplicitSquareMatrixOperation>
172  (m,target_dimension,eigen_strategy.skip());
173  if (eigen_strategy.is(SmallestEigenvalues))
174  return eigendecomposition_impl_arpack<DenseMatrix,DenseInverseMatrixOperation>
175  (m,target_dimension,eigen_strategy.skip());
176  unsupported();
177  }
178 #ifdef TAPKEE_WITH_VIENNACL
179  if (strategy.is(HeterogeneousOpenCLStrategy))
180  {
181  if (eigen_strategy.is(LargestEigenvalues))
182  return eigendecomposition_impl_arpack<DenseMatrix,GPUDenseMatrixOperation>
183  (m,target_dimension,eigen_strategy.skip());
184  if (eigen_strategy.is(SquaredLargestEigenvalues))
185  return eigendecomposition_impl_arpack<DenseMatrix,GPUDenseImplicitSquareMatrixOperation>
186  (m,target_dimension,eigen_strategy.skip());
187  unsupported();
188  }
189 #endif
190  unsupported();
191  return EigendecompositionResult();
192  }
193 #endif
195  const EigendecompositionStrategy& eigen_strategy,
196  IndexType target_dimension)
197  {
198  if(strategy.is(HomogeneousCPUStrategy))
199  {
200  if (eigen_strategy.is(LargestEigenvalues))
201  return eigendecomposition_impl_dense<DenseMatrix,DenseMatrixOperation>
202  (m,target_dimension,eigen_strategy.skip());
203  if (eigen_strategy.is(SquaredLargestEigenvalues))
204  return eigendecomposition_impl_dense<DenseMatrix,DenseMatrixOperation>
205  (m,target_dimension,eigen_strategy.skip());
206  if (eigen_strategy.is(SmallestEigenvalues))
207  return eigendecomposition_impl_dense<DenseMatrix,DenseInverseMatrixOperation>
208  (m,target_dimension,eigen_strategy.skip());
209  unsupported();
210  }
211  unsupported();
212  return EigendecompositionResult();
213  }
215  const EigendecompositionStrategy& eigen_strategy,
216  IndexType target_dimension)
217  {
218  if (strategy.is(HomogeneousCPUStrategy))
219  {
220  if (eigen_strategy.is(LargestEigenvalues))
221  return eigendecomposition_impl_randomized<DenseMatrix,DenseMatrixOperation>
222  (m,target_dimension,eigen_strategy.skip());
223  if (eigen_strategy.is(SquaredLargestEigenvalues))
224  return eigendecomposition_impl_randomized<DenseMatrix,DenseImplicitSquareMatrixOperation>
225  (m,target_dimension,eigen_strategy.skip());
226  if (eigen_strategy.is(SmallestEigenvalues))
227  return eigendecomposition_impl_randomized<DenseMatrix,DenseInverseMatrixOperation>
228  (m,target_dimension,eigen_strategy.skip());
229  unsupported();
230  }
231 #ifdef TAPKEE_WITH_VIENNACL
232  if (strategy.is(HeterogeneousOpenCLStrategy))
233  {
234  if (eigen_strategy.is(LargestEigenvalues))
235  return eigendecomposition_impl_randomized<DenseMatrix,GPUDenseMatrixOperation>
236  (m,target_dimension,eigen_strategy.skip());
237  if (eigen_strategy.is(SquaredLargestEigenvalues))
238  return eigendecomposition_impl_randomized<DenseMatrix,GPUDenseImplicitSquareMatrixOperation>
239  (m,target_dimension,eigen_strategy.skip());
240  unsupported();
241  }
242 #endif
243  unsupported();
244  return EigendecompositionResult();
245  }
246  inline void unsupported() const
247  {
248  throw unsupported_method_error("Unsupported method");
249  }
250 };
251 
252 template <>
254 {
255 #ifdef TAPKEE_WITH_ARPACK
257  const EigendecompositionStrategy& eigen_strategy,
258  IndexType target_dimension)
259  {
260  if (strategy.is(HomogeneousCPUStrategy))
261  {
262  if (eigen_strategy.is(SmallestEigenvalues))
263  return eigendecomposition_impl_arpack<SparseWeightMatrix,SparseInverseMatrixOperation>
264  (m,target_dimension,eigen_strategy.skip());
265  unsupported();
266  }
267  unsupported();
268  return EigendecompositionResult();
269  }
270 #endif
272  const EigendecompositionStrategy& eigen_strategy,
273  IndexType target_dimension)
274  {
275  if (strategy.is(HomogeneousCPUStrategy))
276  {
277  if (eigen_strategy.is(SmallestEigenvalues))
278  return eigendecomposition_impl_dense<SparseWeightMatrix,SparseInverseMatrixOperation>
279  (m,target_dimension,eigen_strategy.skip());
280  unsupported();
281  }
282  unsupported();
283  return EigendecompositionResult();
284  }
286  const EigendecompositionStrategy& eigen_strategy,
287  IndexType target_dimension)
288  {
289  if (strategy.is(HomogeneousCPUStrategy))
290  {
291  if (eigen_strategy.is(SmallestEigenvalues))
292  return eigendecomposition_impl_randomized<SparseWeightMatrix,SparseInverseMatrixOperation>
293  (m,target_dimension,eigen_strategy.skip());
294  unsupported();
295  }
296  unsupported();
297  return EigendecompositionResult();
298  }
299  inline void unsupported() const
300  {
301  throw unsupported_method_error("Unsupported method");
302  }
303 };
304 
333 template <class MatrixType>
335  const EigendecompositionStrategy& eigen_strategy,
336  const MatrixType& m, IndexType target_dimension)
337 {
338  LoggingSingleton::instance().message_info(formatting::format("Using the {} eigendecomposition method.",
339  get_eigen_method_name(method)));
340 #ifdef TAPKEE_WITH_ARPACK
341  if (method.is(Arpack))
342  return eigendecomposition_impl<MatrixType>().arpack(m,strategy,eigen_strategy,target_dimension);
343 #endif
344  if (method.is(Randomized))
345  return eigendecomposition_impl<MatrixType>().randomized(m,strategy,eigen_strategy,target_dimension);
346  if (method.is(Dense))
347  return eigendecomposition_impl<MatrixType>().dense(m,strategy,eigen_strategy,target_dimension);
348  return EigendecompositionResult();
349 }
350 
351 
352 } // End of namespace tapkee_internal
353 } // End of namespace tapkee
354 
355 #endif
EigendecompositionResult eigendecomposition_impl_randomized(const MatrixType &wm, IndexType target_dimension, unsigned int skip)
Randomized redsvd-like implementation of eigendecomposition-based embedding.
static const EigenMethod Dense("Dense")
Eigen library dense method (could be useful for debugging). Computes all eigenvectors thus can be ver...
static const EigendecompositionStrategy SmallestEigenvalues("Smallest eigenvalues", 1)
static const EigenMethod Randomized("Randomized")
Randomized method (implementation taken from the redsvd lib). Supports only standard but not generali...
EigendecompositionResult randomized(const SparseWeightMatrix &m, const ComputationStrategy &strategy, const EigendecompositionStrategy &eigen_strategy, IndexType target_dimension)
static const EigendecompositionStrategy SquaredLargestEigenvalues("Largest eigenvalues of squared matrix", 0)
const Matrix< Scalar, Dynamic, Dynamic > & eigenvectors() const
Returns the eigenvectors of given matrix.
Eigen::Matrix< tapkee::ScalarType, Eigen::Dynamic, Eigen::Dynamic > DenseMatrix
dense matrix type (non-overridable)
Definition: types.hpp:23
EigendecompositionResult eigendecomposition(const EigenMethod &method, const ComputationStrategy &strategy, const EigendecompositionStrategy &eigen_strategy, const MatrixType &m, IndexType target_dimension)
Multiple implementation handler method for various eigendecomposition methods.
EigendecompositionResult dense(const DenseMatrix &m, const ComputationStrategy &strategy, const EigendecompositionStrategy &eigen_strategy, IndexType target_dimension)
EigendecompositionResult randomized(const DenseMatrix &m, const ComputationStrategy &strategy, const EigendecompositionStrategy &eigen_strategy, IndexType target_dimension)
double ScalarType
default scalar value (can be overrided with TAPKEE_CUSTOM_INTERNAL_NUMTYPE define) ...
Definition: types.hpp:15
ComputationInfo info() const
Reports whether previous computation was successful.
EigendecompositionResult dense(const MatrixType &m, const ComputationStrategy &strategy, const EigendecompositionStrategy &eigen_strategy, IndexType target_dimension)
An exception type that is thrown when eigendecomposition is failed.
EigendecompositionResult randomized(const MatrixType &m, const ComputationStrategy &strategy, const EigendecompositionStrategy &eigen_strategy, IndexType target_dimension)
Eigen::SparseMatrix< tapkee::ScalarType > SparseWeightMatrix
sparse weight matrix type (non-overridable)
Definition: types.hpp:29
EigendecompositionResult arpack(const MatrixType &m, const ComputationStrategy &strategy, const EigendecompositionStrategy &eigen_strategy, IndexType target_dimension)
EigendecompositionResult arpack(const SparseWeightMatrix &m, const ComputationStrategy &strategy, const EigendecompositionStrategy &eigen_strategy, IndexType target_dimension)
static const EigenMethod Arpack("Arpack")
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.
std::string format(const std::string &fmt, const ValueWrapper a)
Definition: formatting.hpp:411
EigendecompositionResult dense(const SparseWeightMatrix &m, const ComputationStrategy &strategy, const EigendecompositionStrategy &eigen_strategy, IndexType target_dimension)
int IndexType
indexing type (non-overridable) set to int for compatibility with OpenMP 2.0
Definition: types.hpp:19
EigendecompositionResult eigendecomposition_impl_dense(const MatrixType &wm, IndexType target_dimension, unsigned int skip)
Eigen library dense implementation of eigendecomposition-based embedding.
static const ComputationStrategy HomogeneousCPUStrategy("CPU")
Eigen::SelfAdjointEigenSolver< tapkee::DenseMatrix > DenseSelfAdjointEigenSolver
selfadjoint solver (non-overridable)
Definition: types.hpp:33
static const EigendecompositionStrategy LargestEigenvalues("Largest eigenvalues", 0)
void message_info(const std::string &msg)
Definition: logging.hpp:115
std::string get_eigen_method_name(const EigenMethod &m)
Definition: naming.hpp:48
EigendecompositionResult arpack(const DenseMatrix &m, const ComputationStrategy &strategy, const EigendecompositionStrategy &eigen_strategy, IndexType target_dimension)
const Matrix< Scalar, Dynamic, 1 > & eigenvalues() const
Returns the eigenvalues of given matrix.
static LoggingSingleton & instance()
Definition: logging.hpp:102
EigendecompositionResult eigendecomposition_impl_arpack(const MatrixType &wm, IndexType target_dimension, unsigned int skip)
ARPACK implementation of eigendecomposition-based embedding.
ScalarType gaussian_random()
Definition: random.hpp:39
An exception type that is thrown when unsupported method is called.
TAPKEE_INTERNAL_PAIR< tapkee::DenseMatrix, tapkee::DenseVector > EigendecompositionResult
Definition: synonyms.hpp:41
tapkee::DenseMatrix DenseSymmetricMatrix
dense symmetric matrix (non-overridable, currently just dense matrix, can be improved later) ...
Definition: types.hpp:25
bool is(const M &m) const