Tapkee
fa.hpp
Go to the documentation of this file.
00001 /* This software is distributed under BSD 3-clause license (see LICENSE file).
00002  *
00003  * Copyright (c) 2013 Sergey Lisitsyn, Fernando J. Iglesias Garcia
00004  */
00005 
00006 #ifndef TAPKEE_FA_H_
00007 #define TAPKEE_FA_H_
00008 
00009 namespace tapkee
00010 {
00011 namespace tapkee_internal
00012 {
00013 
00014 template <class RandomAccessIterator, class FeatureVectorCallback>
00015 DenseMatrix project(RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback callback,
00016         IndexType dimension, const IndexType max_iter, const ScalarType epsilon,
00017         const IndexType target_dimension, const DenseVector& mean_vector)
00018 {
00019     timed_context context("Data projection");
00020 
00021     // The number of data points
00022     const IndexType n = end-begin;
00023 
00024     // Dense representation of the data points
00025 
00026     DenseVector current_vector(dimension);
00027 
00028     DenseMatrix X = DenseMatrix::Zero(dimension,n);
00029 
00030     for (RandomAccessIterator iter=begin; iter!=end; ++iter)
00031     {
00032         callback.vector(*iter,current_vector);
00033         X.col(iter-begin) = current_vector - mean_vector;
00034     }
00035 
00036     // Initialize FA model
00037 
00038     // Initial variances
00039     DenseMatrix sig = DenseMatrix::Identity(dimension,dimension);
00040     // Initial linear mapping
00041     DenseMatrix A = DenseMatrix::Random(dimension, target_dimension).cwiseAbs();
00042 
00043     // Main loop
00044     IndexType iter = 0;
00045     DenseMatrix invC,M,SC;
00046     ScalarType ll = 0, newll = 0;
00047     while (iter < max_iter)
00048     {
00049         ++iter;
00050 
00051         // Perform E-step
00052 
00053         // Compute the inverse of the covariance matrix
00054         invC = (A*A.transpose() + sig).inverse();
00055         M = A.transpose()*invC*X;
00056         SC = n*(DenseMatrix::Identity(target_dimension,target_dimension) - A.transpose()*invC*A) + M*M.transpose();
00057 
00058         // Perform M-step
00059         A = (X*M.transpose())*SC.inverse();
00060         sig = DenseMatrix(((X*X.transpose() - A*M*X.transpose()).diagonal()/n).asDiagonal()).array() + epsilon;
00061 
00062         // Compute log-likelihood of FA model
00063         newll = 0.5*(log(invC.determinant()) - (invC*X).cwiseProduct(X).sum()/n);
00064 
00065         // Check for convergence
00066         if ((iter > 1) && (fabs(newll - ll) < epsilon))
00067             break;
00068 
00069         ll = newll;
00070     }
00071 
00072     return X.transpose()*A;
00073 }
00074 
00075 } /* namespace tapkee */
00076 } /* namespace tapkee_internal */
00077 
00078 #endif /* TAPKEE_FA_H_ */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines