Tapkee
locally_linear.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) 2012-2013 Sergey Lisitsyn
00004  */
00005 
00006 #ifndef TAPKEE_LOCALLY_LINEAR_H_
00007 #define TAPKEE_LOCALLY_LINEAR_H_
00008 
00009 /* Tapkee includes */
00010 #include <tapkee/routines/eigendecomposition.hpp>
00011 #include <tapkee/defines.hpp>
00012 #include <tapkee/utils/matrix.hpp>
00013 #include <tapkee/utils/time.hpp>
00014 #include <tapkee/utils/sparse.hpp>
00015 /* End of Tapkee includes */
00016 
00017 namespace tapkee
00018 {
00019 namespace tapkee_internal
00020 {
00021 
00022 
00023 
00024 template <class RandomAccessIterator, class PairwiseCallback>
00025 SparseWeightMatrix tangent_weight_matrix(RandomAccessIterator begin, RandomAccessIterator end, 
00026                                          const Neighbors& neighbors, PairwiseCallback callback, 
00027                                          const IndexType target_dimension, const ScalarType shift,
00028                                          const bool partial_eigendecomposer=false)
00029 {
00030     timed_context context("KLTSA weight matrix computation");
00031     const IndexType k = neighbors[0].size();
00032 
00033     SparseTriplets sparse_triplets;
00034     sparse_triplets.reserve((k*k+2*k+1)*(end-begin));
00035 
00036 #pragma omp parallel shared(begin,end,neighbors,callback,sparse_triplets) default(none)
00037     {
00038         IndexType index_iter;
00039         DenseMatrix gram_matrix = DenseMatrix::Zero(k,k);
00040         DenseVector rhs = DenseVector::Ones(k);
00041         DenseMatrix G = DenseMatrix::Zero(k,target_dimension+1);
00042         G.col(0).setConstant(1/sqrt(static_cast<ScalarType>(k)));
00043         DenseSelfAdjointEigenSolver solver;
00044         SparseTriplets local_triplets;
00045         local_triplets.reserve(k*k+2*k+1);
00046 
00047 #pragma omp for nowait
00048         for (index_iter=0; index_iter<static_cast<IndexType>(end-begin); index_iter++)
00049         {
00050             const LocalNeighbors& current_neighbors = neighbors[index_iter];
00051         
00052             for (IndexType i=0; i<k; ++i)
00053             {
00054                 for (IndexType j=i; j<k; ++j)
00055                 {
00056                     ScalarType kij = callback.kernel(begin[current_neighbors[i]],begin[current_neighbors[j]]);
00057                     gram_matrix(i,j) = kij;
00058                     gram_matrix(j,i) = kij;
00059                 }
00060             }
00061             
00062             centerMatrix(gram_matrix);
00063 
00064             //UNRESTRICT_ALLOC;
00065 #ifdef TAPKEE_WITH_ARPACK
00066             if (partial_eigendecomposer)
00067             {
00068                 G.rightCols(target_dimension).noalias() =
00069                     eigendecomposition<DenseMatrix,DenseMatrixOperation>(Arpack,gram_matrix,target_dimension,0).first;
00070             }
00071             else
00072 #endif
00073             {
00074                 solver.compute(gram_matrix);
00075                 G.rightCols(target_dimension).noalias() = solver.eigenvectors().rightCols(target_dimension);
00076             }
00077             //RESTRICT_ALLOC;
00078             gram_matrix.noalias() = G * G.transpose();
00079             
00080             SparseTriplet diagonal_triplet(index_iter,index_iter,shift);
00081             local_triplets.push_back(diagonal_triplet);
00082             for (IndexType i=0; i<k; ++i)
00083             {
00084                 SparseTriplet neighborhood_diagonal_triplet(current_neighbors[i],current_neighbors[i],1.0);
00085                 local_triplets.push_back(neighborhood_diagonal_triplet);
00086 
00087                 for (IndexType j=0; j<k; ++j) 
00088                 {
00089                     SparseTriplet tangent_triplet(current_neighbors[i],current_neighbors[j],-gram_matrix(i,j));
00090                     local_triplets.push_back(tangent_triplet);
00091                 }
00092             }
00093 #pragma omp critical
00094             {
00095                 copy(local_triplets.begin(),local_triplets.end(),back_inserter(sparse_triplets));
00096             }
00097             
00098             local_triplets.clear();
00099         }
00100     }
00101 
00102     return sparse_matrix_from_triplets(sparse_triplets, end-begin, end-begin);
00103 }
00104 
00105 template <class RandomAccessIterator, class PairwiseCallback>
00106 SparseWeightMatrix linear_weight_matrix(const RandomAccessIterator& begin, const RandomAccessIterator& end, 
00107                                         const Neighbors& neighbors, PairwiseCallback callback,
00108                                         const ScalarType shift, const ScalarType trace_shift)
00109 {
00110     timed_context context("KLLE weight computation");
00111     const IndexType k = neighbors[0].size();
00112 
00113     SparseTriplets sparse_triplets;
00114     sparse_triplets.reserve((k*k+2*k+1)*(end-begin));
00115 
00116 #pragma omp parallel shared(begin,end,neighbors,callback,sparse_triplets) default(none)
00117     {
00118         IndexType index_iter;
00119         DenseMatrix gram_matrix = DenseMatrix::Zero(k,k);
00120         DenseVector dots(k);
00121         DenseVector rhs = DenseVector::Ones(k);
00122         DenseVector weights;
00123         SparseTriplets local_triplets;
00124         local_triplets.reserve(k*k+2*k+1);
00125         
00126         //RESTRICT_ALLOC;
00127 #pragma omp for nowait
00128         for (index_iter=0; index_iter<static_cast<IndexType>(end-begin); index_iter++)
00129         {
00130             ScalarType kernel_value = callback.kernel(begin[index_iter],begin[index_iter]);
00131             const LocalNeighbors& current_neighbors = neighbors[index_iter];
00132             
00133             for (IndexType i=0; i<k; ++i)
00134                 dots[i] = callback.kernel(begin[index_iter], begin[current_neighbors[i]]);
00135 
00136             for (IndexType i=0; i<k; ++i)
00137             {
00138                 for (IndexType j=i; j<k; ++j)
00139                     gram_matrix(i,j) = kernel_value - dots(i) - dots(j) + 
00140                                        callback.kernel(begin[current_neighbors[i]],begin[current_neighbors[j]]);
00141             }
00142             
00143             ScalarType trace = gram_matrix.trace();
00144             gram_matrix.diagonal().array() += trace_shift*trace;
00145             weights = gram_matrix.selfadjointView<Eigen::Upper>().ldlt().solve(rhs);
00146             weights /= weights.sum();
00147 
00148             SparseTriplet diagonal_triplet(index_iter,index_iter,1.0+shift);
00149             local_triplets.push_back(diagonal_triplet);
00150             for (IndexType i=0; i<k; ++i)
00151             {
00152                 SparseTriplet row_side_triplet(current_neighbors[i],index_iter,-weights[i]);
00153                 SparseTriplet col_side_triplet(index_iter,current_neighbors[i],-weights[i]);
00154                 local_triplets.push_back(row_side_triplet);
00155                 local_triplets.push_back(col_side_triplet);
00156                 for (IndexType j=0; j<k; ++j)
00157                 {
00158                     SparseTriplet cross_triplet(current_neighbors[i],current_neighbors[j],weights(i)*weights(j));
00159                     local_triplets.push_back(cross_triplet);
00160                 }
00161             }
00162 
00163 #pragma omp critical
00164             {
00165                 copy(local_triplets.begin(),local_triplets.end(),back_inserter(sparse_triplets));
00166             }
00167             
00168             local_triplets.clear();
00169         }
00170         //UNRESTRICT_ALLOC;
00171     }
00172 
00173     return sparse_matrix_from_triplets(sparse_triplets, end-begin, end-begin);
00174 }
00175 
00176 template <class RandomAccessIterator, class PairwiseCallback>
00177 SparseWeightMatrix hessian_weight_matrix(RandomAccessIterator begin, RandomAccessIterator end, 
00178                                          const Neighbors& neighbors, PairwiseCallback callback, 
00179                                          const IndexType target_dimension)
00180 {
00181     timed_context context("Hessian weight matrix computation");
00182     const IndexType k = neighbors[0].size();
00183 
00184     SparseTriplets sparse_triplets;
00185     sparse_triplets.reserve(k*k*(end-begin));
00186 
00187     const IndexType dp = target_dimension*(target_dimension+1)/2;
00188 
00189 #pragma omp parallel shared(begin,end,neighbors,callback,sparse_triplets) default(none)
00190     {
00191         IndexType index_iter;
00192         DenseMatrix gram_matrix = DenseMatrix::Zero(k,k);
00193         DenseMatrix Yi(k,1+target_dimension+dp);
00194 
00195         SparseTriplets local_triplets;
00196         local_triplets.reserve(k*k+2*k+1);
00197 
00198 #pragma omp for nowait
00199         for (index_iter=0; index_iter<static_cast<IndexType>(end-begin); index_iter++)
00200         {
00201             const LocalNeighbors& current_neighbors = neighbors[index_iter];
00202         
00203             for (IndexType i=0; i<k; ++i)
00204             {
00205                 for (IndexType j=i; j<k; ++j)
00206                 {
00207                     ScalarType kij = callback.kernel(begin[current_neighbors[i]],begin[current_neighbors[j]]);
00208                     gram_matrix(i,j) = kij;
00209                     gram_matrix(j,i) = kij;
00210                 }
00211             }
00212             
00213             centerMatrix(gram_matrix);
00214             
00215             DenseSelfAdjointEigenSolver sae_solver;
00216             sae_solver.compute(gram_matrix);
00217 
00218             Yi.col(0).setConstant(1.0);
00219             Yi.block(0,1,k,target_dimension).noalias() = sae_solver.eigenvectors().rightCols(target_dimension);
00220 
00221             IndexType ct = 0;
00222             for (IndexType j=0; j<target_dimension; ++j)
00223             {
00224                 for (IndexType p=0; p<target_dimension-j; ++p)
00225                 {
00226                     Yi.col(ct+p+1+target_dimension).noalias() = Yi.col(j+1).cwiseProduct(Yi.col(j+p+1));
00227                 }
00228                 ct += ct + target_dimension - j;
00229             }
00230             
00231             for (IndexType i=0; i<static_cast<IndexType>(Yi.cols()); i++)
00232             {
00233                 for (IndexType j=0; j<i; j++)
00234                 {
00235                     ScalarType r = Yi.col(i).dot(Yi.col(j));
00236                     Yi.col(i) -= r*Yi.col(j);
00237                 }
00238                 ScalarType norm = Yi.col(i).norm();
00239                 Yi.col(i) *= (1.f / norm);
00240             }
00241             for (IndexType i=0; i<dp; i++)
00242             {
00243                 ScalarType colsum = Yi.col(1+target_dimension+i).sum();
00244                 if (colsum > 1e-4)
00245                     Yi.col(1+target_dimension+i).array() /= colsum;
00246             }
00247 
00248             // reuse gram matrix storage m'kay?
00249             gram_matrix.noalias() = Yi.rightCols(dp)*(Yi.rightCols(dp).transpose());
00250 
00251             for (IndexType i=0; i<k; ++i)
00252             {
00253                 for (IndexType j=0; j<k; ++j)
00254                 {
00255                     SparseTriplet hessian_triplet(current_neighbors[i],current_neighbors[j],gram_matrix(i,j));
00256                     local_triplets.push_back(hessian_triplet);
00257                 }
00258             }
00259 
00260             #pragma omp critical
00261             {
00262                 copy(local_triplets.begin(),local_triplets.end(),back_inserter(sparse_triplets));
00263             }
00264 
00265             local_triplets.clear();
00266         }
00267     }
00268 
00269     return sparse_matrix_from_triplets(sparse_triplets, end-begin, end-begin);
00270 }
00271 
00272 template<class RandomAccessIterator, class FeatureVectorCallback>
00273 DenseSymmetricMatrixPair construct_neighborhood_preserving_eigenproblem(SparseWeightMatrix W,
00274         RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback feature_vector_callback,
00275         IndexType dimension)
00276 {
00277     timed_context context("NPE eigenproblem construction");
00278     
00279     DenseSymmetricMatrix lhs = DenseSymmetricMatrix::Zero(dimension,dimension);
00280     DenseSymmetricMatrix rhs = DenseSymmetricMatrix::Zero(dimension,dimension);
00281 
00282     DenseVector rank_update_vector_i(dimension);
00283     DenseVector rank_update_vector_j(dimension);
00284 
00285     //RESTRICT_ALLOC;
00286     for (RandomAccessIterator iter=begin; iter!=end; ++iter)
00287     {
00288         feature_vector_callback.vector(*iter,rank_update_vector_i);
00289         rhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i);
00290     }
00291 
00292     for (int i=0; i<W.outerSize(); ++i)
00293     {
00294         for (SparseWeightMatrix::InnerIterator it(W,i); it; ++it)
00295         {
00296             feature_vector_callback.vector(begin[it.row()],rank_update_vector_i);
00297             feature_vector_callback.vector(begin[it.col()],rank_update_vector_j);
00298             lhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i, rank_update_vector_j, it.value());
00299         }
00300     }
00301     
00302     rhs += rhs.transpose().eval();
00303     rhs /= 2;
00304 
00305     //UNRESTRICT_ALLOC;
00306 
00307     return DenseSymmetricMatrixPair(lhs,rhs);
00308 }
00309 
00310 template<class RandomAccessIterator, class FeatureVectorCallback>
00311 DenseSymmetricMatrixPair construct_lltsa_eigenproblem(SparseWeightMatrix W,
00312         RandomAccessIterator begin, RandomAccessIterator end, FeatureVectorCallback feature_vector_callback,
00313         IndexType dimension)
00314 {
00315     timed_context context("LLTSA eigenproblem construction");
00316 
00317     DenseSymmetricMatrix lhs = DenseSymmetricMatrix::Zero(dimension,dimension);
00318     DenseSymmetricMatrix rhs = DenseSymmetricMatrix::Zero(dimension,dimension);
00319 
00320     DenseVector rank_update_vector_i(dimension);
00321     DenseVector rank_update_vector_j(dimension);
00322     DenseVector sum = DenseVector::Zero(dimension);
00323     
00324     //RESTRICT_ALLOC;
00325     for (RandomAccessIterator iter=begin; iter!=end; ++iter)
00326     {
00327         feature_vector_callback.vector(*iter,rank_update_vector_i);
00328         sum += rank_update_vector_i;
00329         rhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i);
00330     }
00331     rhs.selfadjointView<Eigen::Upper>().rankUpdate(sum,-1./(end-begin));
00332 
00333     for (int i=0; i<W.outerSize(); ++i)
00334     {
00335         for (SparseWeightMatrix::InnerIterator it(W,i); it; ++it)
00336         {
00337             feature_vector_callback.vector(begin[it.row()],rank_update_vector_i);
00338             feature_vector_callback.vector(begin[it.col()],rank_update_vector_j);
00339             lhs.selfadjointView<Eigen::Upper>().rankUpdate(rank_update_vector_i, rank_update_vector_j, it.value());
00340         }
00341     }
00342     lhs.selfadjointView<Eigen::Upper>().rankUpdate(sum,-1./(end-begin));
00343 
00344     rhs += rhs.transpose().eval();
00345     rhs /= 2;
00346 
00347     //UNRESTRICT_ALLOC;
00348 
00349     return DenseSymmetricMatrixPair(lhs,rhs);
00350 }
00351 
00352 } // End of namespace tapkee_internal
00353 } // End of namespace tapkee
00354 
00355 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines