Still working on the dual manifold

2008-04-11 14:54:15 +00:00
parent fec6fd22c1
commit 9bc02c834e
5 changed files with 224 additions and 9 deletions
@@ -20,16 +20,47 @@
 #define DUAL_MANIFOLD_ENGINE_
 #include "dual_manifold_objective.h"
 #include "../l_bfgs/l_bfgs.h"
-
+#include "../l_bfgs/optimization_utils.h"
+/**
+ * DualManifoldEngine treats the problem of non-linear, non negative matrix factorization
+ * The classic problem of NMF is (\f& D \simeq WH\f&)
+ * For large scale problems D is sparse and it has a lot of zeros that represent
+ * missing or don't care data. This algorithm solves the following set of optimization
+ * problems:
+ * (\f$ \max trHH^T subject to A_i \bullet D = A_i \bullet (WH))
+ * (\f$ \max trWW^T subject to A_i \bullet D = A_i \bullet (WH))
+ * where (\f$\bullet \f$) is the matrix dot product
+ * (\f$ A_i \f$) is a selection matrix that selects the non-zero elements of (\f$D\f$)
+ * The matrix (\f& D \f&) can be real or non-negative.
+ * We can also restrict (\f$ W, H \f$) to be non-negative or sparse
+ * This depends on the definition of the OptimizedFunction
+ */
 template<typename OptimizedFunction>
 class DualManifoldEngine {
 public:
-  void Init(datanode *module);
+  /**
+   * pairs_to_consider are (row,column) indices from a given sparse D
+   * matrix, These are the elements that we care about in our factorization and 
+   * dot_prod_values are the values that we are trying to match
+   * 
+   */
+  void Init(datanode *module, 
+    // index pairs to consider from the matrix (row,column) pairs
+    ArrayList<std::pair<index_t, index_t> > &pairs_to_consider, 
+    // The values of the (row, column) values, also known as the dot products
+    ArrayList<double> &dot_prod_values);
  void Destruct();
+  void ComputeLocalOptimum(); 
+ 
+ private:
+  double feasibility_tolerance_;
+  double desired_feasibility_;
+  double norm_grad_tolerance_;
+  index_t iterations_;
+  index_t max_iterations_;
+  datanode *module;
    
  
- private:
-  datanode *module;
  LBfgs<OptimizedFunction> lbfgs1_;
  LBfgs<OptimizedFunction> lbfgs2_;
  OptimizedFunction optimized_function1_;
@@ -18,24 +18,39 @@

 template<typename OptimizedFunction>
 DualManifoldEngine<OptimizedFunction>::Init(datanode *module, 
-    ArrayList<std::pair<index_t, index_t> > &pairs_to_consider,
+    // index pairs to consider from the matrix (row,column) pairs
+    ArrayList<std::pair<index_t, index_t> > &pairs_to_consider, 
+    // The values of the (row, column) values, also known as the dot products
    ArrayList<double> &dot_prod_values) {

  module_=module;
-  lbfg1_.Init(&optimized_function1_, fx_param_node(module_, "lbfgs"));
-  lbfg2_.Init(&optimized_function2_, fx_param_node(module_, "lbfgs"));
+  l_bfg1_.Init(&optimized_function1_, fx_param_node(module_, "l_bfgs"));
+  l_bfg2_.Init(&optimized_function2_, fx_param_node(module_, "l_bfgs"));

  optimized_function1_.Init(fx_param_node(module_, "opt1"), 
-      lbfgs1_.coordinates(), 
+      l_bfgs1_.coordinates(), 
      pairs_to_consider,
      dot_prod_values);
  for(index_t i=0; i<pairs_to_consider.size(); i++) {
    std::swap(pairs_to_consider[i].first, pairs_to_consider[i].second);
  }
  optimized_function2_.Init(fx_param_node(module_, "opt2"), 
-      lbfgs1_.coordinates(), 
+      l_bfgs1_.coordinates(), 
      pairs_to_consider,
      dot_prod_values);
  
 }

+template<typename OptimizedFunction>
+void DualManifoldEngine<OptimizedFunction>::ComputeLocalOptimum() {
+  for(index_t i=0; i<max_iterations_; i++) {
+    l_bfgs1_.ComputeLocalOptimumBFGS();
+    l_bfgs2_.ComputeLocalOptimumBFGS();
+    double error;
+    optimoptimized_function1_.ComputeFeasibilityError(l_bfgs1_.coordinates(), &error);
+    if (error<desired_error_) {
+      break;
+    }
+  } 
+}
+
@@ -0,0 +1,51 @@
+/*
+ * =====================================================================================
+ * 
+ *       Filename:  mvu_dot_prod_objective.h
+ * 
+ *    Description:  
+ * 
+ *        Version:  1.0
+ *        Created:  04/09/2008 05:30:43 PM EDT
+ *       Revision:  none
+ *       Compiler:  gcc
+ * 
+ *         Author:  Nikolaos Vasiloglou (NV), nvasil@ieee.org
+ *        Company:  Georgia Tech Fastlab-ESP Lab
+ * 
+ * =====================================================================================
+ */
+
+#ifndef MVU_DOT_PROD_OBJECTIVE_H_
+#define MVU_DOT_PROD_OBJECTIVE_H_
+#include "fastlib/fastlib.h"
+#include "../l_bfgs/optimization_utils.h"
+class MVUDotProdObjective {
+ public:
+  void Init(datanode *module,
+      Matrix *coordinates, 
+      ArrayList<std::pair<index_t, index_t> > &pairs_to_consider, 
+      // The values of the (row, column) values, also known as the dot products
+      ArrayList<double> &dot_prod_values);
+  void ComputeGradient(Matrix &coordinates, Matrix *gradient);
+  void ComputeObjective(Matrix &coordinates, double *objective);
+  void ComputeFeasibilityError(Matrix &coordinates, double *error);
+  double ComputeLagrangian(Matrix &coordinates);
+  void UpdateLagrangeMult(Matrix &coordinates);
+  void Project(Matrix *coordinates);
+  void set_sigma(double sigma); 
+  bool IsDiverging(double objective); 
+
+ private:
+  datanode *module_;
+  Matrix *auxiliary_mat_;
+  ArrayList<std::pair<index_t, index_t> > pairs_to_consider_; 
+  ArrayList<double> dot_prod_values_;
+  Vector eq_lagrange_mult_;
+  double sigma_;
+  index_t num_of_constraints_;
+};
+
+
+#include "mvu_dot_prod_objective_impl.h"
+#endif // MVU_DOT_PROD_OBJECTIVE_H_
@@ -0,0 +1,115 @@
+/*
+ * =====================================================================================
+ * 
+ *       Filename:  mvu_dot_prod_objective_impl.h
+ * 
+ *    Description:  
+ * 
+ *        Version:  1.0
+ *        Created:  04/09/2008 06:24:01 PM EDT
+ *       Revision:  none
+ *       Compiler:  gcc
+ * 
+ *         Author:  Nikolaos Vasiloglou (NV), nvasil@ieee.org
+ *        Company:  Georgia Tech Fastlab-ESP Lab
+ * 
+ * =====================================================================================
+ */
+
+void MVUDotProdObjective::Init(datanode *module,
+    Matrix *coordinates, 
+    ArrayList<std::pair<index_t, index_t> > &pairs_to_consider, 
+    // The values of the (row, column) values, also known as the dot products
+    ArrayList<double> &dot_prod_values) {
+  
+  module_=module;
+  auxiliary_mat_=coordinates;
+  pairs_to_consider_.Copy(pairs_to_consider);
+  dot_prod_values_.Copy(dot_prod_values);
+  eq_lagrange_mult_.Init(dot_prod_values.size());
+  eq_lagrange_mult_.SetAll(0.0);
+  num_of_constraints_=dot_prod_values_.size();
+}
+
+void MVUDotProdObjective::ComputeGradient(Matrix &coordinates, Matrix *gradient) {
+  gradient->CopyValues(coordinates);
+  // we need to use -CRR^T because we want to maximize CRR^T
+  la::Scale(-1.0, gradient);
+  index_t dimension=auxiliary_mat_->n_rows();
+  Vector constant;
+  for (index_t i=0; i<num_of_constraints; i++) {
+    index_t ind1=pairs_to_consider_[i].first; 
+    index_t ind2=pairs_to_consider_[i].second;
+    double *p1=auxiliary_mat_->GetColumnPtr(ind1);
+    double *p2=coordinates.GetColumnPtr(ind2);
+    double dot_prod =la::Dot(dimension, p1, p2);
+    double diff=dot_prod-dot_prod_values_[i];
+    la::AddExpert(dimension,  
+        -eq_lagrange_mult_[i]+sigma_*diff,
+        p1,
+        gradient->GetColumnPtr(ind1));
+  } 
+}
+
+void MVUDotProdObjective::ComputeObjective(Matrix &coordinates, double *objective) {
+  *objective=0;
+  index_t dimension = coordinates.n_rows();
+  for(index_t i=0; i< coordinates.n_cols(); i++) {
+     *objective-=la::Dot(dimension, 
+                        coordinates.GetColumnPtr(i),
+                        coordinates.GetColumnPtr(i));
+  } 
+}
+
+void MVUDotProdObjective::ComputeFeasibilityError(Matrix &coordinates, double *error) {
+  *error=0;
+  for(index_t i=0; i<num_of_constraints_; i++) {
+    index_t ind1=pairs_to_consider_[i].first; 
+    index_t ind2=pairs_to_consider_[i].second;
+    double *p1=auxiliary_mat_->GetColumnPtr(ind1);
+    double *p2=coordinates.GetColumnPtr(ind2);
+    double dot_prod =la::Dot(dimension, p1, p2);
+    double diff=dot_prod-dot_prod_values_[i];
+    error +=diff*diff;
+  }
+}
+
+double MVUDotProdObjective::ComputeLagrangian(Matrix &coordinates) {
+  double lagrangian=0;
+  ComputeObjective(coordinates, &lagrangian);
+  for(index_t i=0; i<num_of_constriants_; i++) {
+    index_t ind1=pairs_to_consider_[i].first; 
+    index_t ind2=pairs_to_consider_[i].second;
+    double *p1=auxiliary_mat_->GetColumnPtr(ind1);
+    double *p2=coordinates.GetColumnPtr(ind2);
+    double dot_prod =la::Dot(dimension, p1, p2);
+    double diff=dot_prod-dot_prod_values_[i];
+    lagrangian+= -eq_lagrange_mult_[i]*diff + sigma_*diff*diff;
+  }
+}
+
+void MVUDotProdObjective::UpdateLagrangeMult(Matrix &coordinates) {
+  index_t dimension=coordinates.n_rows();
+  for(index_t i=0; i<num_of_constraints_; i++) {
+    index_t ind1=pairs_to_consider_[i].first; 
+    index_t ind2=pairs_to_consider_[i].second;
+    double *p1=auxiliary_mat_->GetColumnPtr(ind1);
+    double *p2=coordinates.GetColumnPtr(ind2);
+    double dot_prod =la::Dot(dimension, p1, p2);
+    double diff=dot_prod-dot_prod_values_[i];
+    eq_lagrange_mult_[i]-=sigma_*diff;
+  }
+}
+
+void MVUDotProdObjective::Project(Matrix *coordinates) {
+  OptUtils::RemoveMean(coordinates);
+}
+
+void MVUDotProdObjective::set_sigma(double sigma) {
+   sigma_=sigma;
+} 
+
+bool MVUDotProdObjective::IsDiverging(double objective) {
+  return false;
+}
+
@@ -33,6 +33,9 @@ class LBfgs {
  void GetResults(Matrix *result);
  void Reset(); 
  void set_coordinates(Matrix &coordinates);
+  void set_desired_feasibility(double desired_feasibility);
+  void set_feasibility_tolerance(double feasibility_tolerance);
+  void set_norm_grad_tolerance(double norm_grad_tolerance);
  Matrix *coordinates();
  
 private: