Tapkee
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
fa.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) 2013 Sergey Lisitsyn, Fernando J. Iglesias Garcia
4  */
5 
6 #ifndef TAPKEE_FA_H_
7 #define TAPKEE_FA_H_
8 
9 namespace tapkee
10 {
11 namespace tapkee_internal
12 {
13 
14 template <class RandomAccessIterator, class FeatureVectorCallback>
15 DenseMatrix project(RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback callback,
16  IndexType dimension, const IndexType max_iter, const ScalarType epsilon,
17  const IndexType target_dimension, const DenseVector& mean_vector)
18 {
19  timed_context context("Data projection");
20 
21  // The number of data points
22  const IndexType n = end-begin;
23 
24  // Dense representation of the data points
25 
26  DenseVector current_vector(dimension);
27 
28  DenseMatrix X = DenseMatrix::Zero(dimension,n);
29 
30  for (RandomAccessIterator iter=begin; iter!=end; ++iter)
31  {
32  callback.vector(*iter,current_vector);
33  X.col(iter-begin) = current_vector - mean_vector;
34  }
35 
36  // Initialize FA model
37 
38  // Initial variances
39  DenseMatrix sig = DenseMatrix::Identity(dimension,dimension);
40  // Initial linear mapping
41  DenseMatrix A = DenseMatrix::Random(dimension, target_dimension).cwiseAbs();
42 
43  // Main loop
44  IndexType iter = 0;
45  DenseMatrix invC,M,SC;
46  ScalarType ll = 0, newll = 0;
47  while (iter < max_iter)
48  {
49  ++iter;
50 
51  // Perform E-step
52 
53  // Compute the inverse of the covariance matrix
54  invC = (A*A.transpose() + sig).inverse();
55  M = A.transpose()*invC*X;
56  SC = n*(DenseMatrix::Identity(target_dimension,target_dimension) - A.transpose()*invC*A) + M*M.transpose();
57 
58  // Perform M-step
59  A = (X*M.transpose())*SC.inverse();
60  sig = DenseMatrix(((X*X.transpose() - A*M*X.transpose()).diagonal()/n).asDiagonal()).array() + epsilon;
61 
62  // Compute log-likelihood of FA model
63  newll = 0.5*(log(invC.determinant()) - (invC*X).cwiseProduct(X).sum()/n);
64 
65  // Check for convergence
66  if ((iter > 1) && (fabs(newll - ll) < epsilon))
67  break;
68 
69  ll = newll;
70  }
71 
72  return X.transpose()*A;
73 }
74 
75 } /* namespace tapkee */
76 } /* namespace tapkee_internal */
77 
78 #endif /* TAPKEE_FA_H_ */
Eigen::Matrix< tapkee::ScalarType, Eigen::Dynamic, Eigen::Dynamic > DenseMatrix
dense matrix type (non-overridable)
Definition: types.hpp:23
double ScalarType
default scalar value (can be overrided with TAPKEE_CUSTOM_INTERNAL_NUMTYPE define) ...
Definition: types.hpp:15
Eigen::Matrix< tapkee::ScalarType, Eigen::Dynamic, 1 > DenseVector
dense vector type (non-overridable)
Definition: types.hpp:21
int IndexType
indexing type (non-overridable) set to int for compatibility with OpenMP 2.0
Definition: types.hpp:19
DenseMatrix project(RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback callback, IndexType dimension, const IndexType max_iter, const ScalarType epsilon, const IndexType target_dimension, const DenseVector &mean_vector)
Definition: fa.hpp:15