diff --git a/fastlib2/contrib/nvasil/convex_nmf/main.cc b/fastlib2/contrib/nvasil/convex_nmf/main.cc
index b1eff0220d..d3bbf2a7dd 100644
--- a/fastlib2/contrib/nvasil/convex_nmf/main.cc
+++ b/fastlib2/contrib/nvasil/convex_nmf/main.cc
@@ -5,14 +5,14 @@ int main(int argc, char *argv[]) {
   fx_module *fx_root;
 	fx_root=fx_init(argc, argv, NULL);
 	fx_module *nmf_module=fx_submodule(fx_root, "/engine");
-	NmfEngine engine;
+	NmfEngine<BigSdpNmfObjectiveMinVar> engine;
 	fx_set_param_str(nmf_module, "data_file", 
-      "/net/hg200/nvasil/dataset/orl_faces/orl_test_faces_100.csv");
-	fx_set_param_int(nmf_module, "sdp_rank",1 );
+      "/net/hg200/nvasil/dataset/orl_faces/orl_test_faces_474.csv");
+	fx_set_param_int(nmf_module, "sdp_rank", 4);
   fx_set_param_int(nmf_module, "new_dim", 20);
-	fx_set_param_double(nmf_module, "l_bfgs/sigma", 0.010);
+	fx_set_param_double(nmf_module, "l_bfgs/sigma", 0.1);
 	fx_set_param_double(nmf_module, "l_bfgs/norm_grad_tolerance", 10);
-	fx_set_param_double(nmf_module, "l_bfgs/desired_feasibility", 1);
+	fx_set_param_double(nmf_module, "l_bfgs/desired_feasibility", 5);
   fx_set_param_double(nmf_module, "l_bfgs/feasibility_tolerance", 0.01);
   fx_set_param_int(nmf_module, "l_bfgs/mem_bfgs", 5);
   engine.Init(nmf_module);
diff --git a/fastlib2/contrib/nvasil/convex_nmf/nmf_engine.h b/fastlib2/contrib/nvasil/convex_nmf/nmf_engine.h
index 24ce696f2d..59cdd24272 100644
--- a/fastlib2/contrib/nvasil/convex_nmf/nmf_engine.h
+++ b/fastlib2/contrib/nvasil/convex_nmf/nmf_engine.h
@@ -4,6 +4,7 @@
 #include "../l_bfgs/l_bfgs.h"
 #include "nmf_objectives.h"
 
+template<typename NmfObjective>
 class NmfEngine {
  public:
 	 void  Init(fx_module *module) {
@@ -59,44 +60,55 @@ class NmfEngine {
         FATAL("Svd failed soething is wrong...\n");
       }
       bool negative_flag=false;
+      index_t positives=0;
+      index_t negatives=0;
+      index_t zeros=0;
       for (index_t i=0; i<vt_mat.n_cols(); i++) {
         if (vt_mat.get(0, i)<0) {
           negative_flag=true;
+          negatives++;
+        }  
+        if (vt_mat.get(0,i)>0) {
+          positives++;
         }
-        if (negative_flag==true && vt_mat.get(0, i)>0) {
-          NONFATAL("Method failed, first eigenvector has positive and negative elements");
-          break;          
-        }
-      } 
-      if (negative_flag==true) {
-        la::Scale(-1, &vt_mat);
-      }
-      opt_function_.Project(&vt_mat);
-      for(index_t i=0; i<num_of_rows_; i++) {
-			  for(index_t j=0; j<new_dim_; j++) {
-				  w_mat_.set(i, j, s[0]*vt_mat.get(0, i*new_dim_+j));
-				}
-			}
-			index_t offset_h=num_of_rows_*new_dim_;
-			for(index_t i=0; i<num_of_columns_; i++) {
-			  for(index_t j=0; j<new_dim_; j++) {
-				   h_mat_.set(j, i , s[0]*vt_mat.get(0, offset_h+i*new_dim_+j ));
-				}
-			}
-      // now compute reconstruction error
-      Matrix v_rec;
-      la::MulInit(w_mat_, h_mat_, &v_rec);
-      double error=0;
-      double v_sum=0;
-      for(index_t i=0; i<values_.size(); i++) {
-        index_t r=rows_[i];
-        index_t c=columns_[i];
-        error+=fabs(v_rec.get(r, c)-values_[i]);
-        v_sum+=values_[i];
-      }
-      data::Save("result.csv", result);
-      NOTIFY("Reconstruction error: %lg%%\n", error*100/v_sum);
+         if (vt_mat.get(0, i)==0) {
+          zeros++;
+        }     
+     } 
+      
+     data::Save("result.csv", result);
+     data::Save("vt_mat.csv", vt_mat);
+     NOTIFY("We found %i positives, %i negatives and %i zeros", 
+         positives, negatives, zeros);
+     if (negatives>=positives) {
+       la::Scale(-1, &vt_mat);
+     }
+     OptUtils::NonNegativeProjection(&vt_mat);
+     for(index_t i=0; i<num_of_rows_; i++) {
+		   for(index_t j=0; j<new_dim_; j++) {
+			   w_mat_.set(i, j, s[0]*vt_mat.get(0, i*new_dim_+j));
+			 }
+		 }
+		 index_t offset_h=num_of_rows_*new_dim_;
+		 for(index_t i=0; i<num_of_columns_; i++) {
+		   for(index_t j=0; j<new_dim_; j++) {
+			   h_mat_.set(j, i , s[0]*vt_mat.get(0, offset_h+i*new_dim_+j));
+			 }
+		 }
+     // now compute reconstruction error
+     Matrix v_rec;
+     la::MulInit(w_mat_, h_mat_, &v_rec);
+     double error=0;
+     double v_sum=0;
+     for(index_t i=0; i<values_.size(); i++) {
+       index_t r=rows_[i];
+       index_t c=columns_[i];
+       error+=fabs(v_rec.get(r, c)-values_[i]);
+       v_sum+=values_[i];
+     }
+     NOTIFY("Reconstruction error: %lg%%\n", error*100/v_sum);
 	 }
+   
 	 void GetW(Matrix *w_mat) {
 	   w_mat->Copy(w_mat_);
 	 }
@@ -106,8 +118,8 @@ class NmfEngine {
 	 
  private:
 	fx_module *module_;
-  LBfgs<BigSdpNmfObjective> engine_;
-	BigSdpNmfObjective opt_function_;
+  LBfgs<NmfObjective> engine_;
+	NmfObjective opt_function_;
 	ArrayList<index_t> rows_;
 	ArrayList<index_t> columns_;
 	ArrayList<double> values_;
diff --git a/fastlib2/contrib/nvasil/convex_nmf/nmf_objectives.h b/fastlib2/contrib/nvasil/convex_nmf/nmf_objectives.h
index 789968e8fc..4873cd24aa 100644
--- a/fastlib2/contrib/nvasil/convex_nmf/nmf_objectives.h
+++ b/fastlib2/contrib/nvasil/convex_nmf/nmf_objectives.h
@@ -3,7 +3,7 @@
 #include "fastlib/fastlib.h"
 #include "../l_bfgs/optimization_utils.h"
 
-class BigSdpNmfObjective {
+class BigSdpNmfObjectiveMaxVar {
  public:
 	void Init(fx_module *module, 
 			ArrayList<index_t> &rows,
@@ -33,7 +33,71 @@ class BigSdpNmfObjective {
 	ArrayList<index_t> columns_;
   ArrayList<double>  values_;
   Vector eq_lagrange_mult_; 
+};
 
+class BigSdpNmfObjectiveMinVar {
+ public:
+	void Init(fx_module *module, 
+			ArrayList<index_t> &rows,
+			ArrayList<index_t> &columns,
+      ArrayList<double>  &values);
+  void Destruct();
+  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); 
+  void GiveInitMatrix(Matrix *init_data);
+	bool IsDiverging(double objective); 
+ 
+ private:
+  fx_module *module_;
+  double sigma_;
+	index_t num_of_columns_;
+	index_t num_of_rows_;
+	index_t new_dim_;
+	index_t rank_;
+	index_t offset_h_;
+	index_t offset_h_mat_; // we need this to know wher h matrix starts 
+	ArrayList<index_t> rows_;
+	ArrayList<index_t> columns_;
+  ArrayList<double>  values_;
+  Vector eq_lagrange_mult_; 
+};
+
+
+class BigSdpNmfObjectiveMinVarIneq {
+ public:
+	void Init(fx_module *module, 
+			ArrayList<index_t> &rows,
+			ArrayList<index_t> &columns,
+      ArrayList<double>  &values);
+  void Destruct();
+  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); 
+  void GiveInitMatrix(Matrix *init_data);
+	bool IsDiverging(double objective); 
+ 
+ private:
+  fx_module *module_;
+  double sigma_;
+	index_t num_of_columns_;
+	index_t num_of_rows_;
+	index_t new_dim_;
+	index_t rank_;
+	index_t offset_h_;
+	index_t offset_h_mat_; // we need this to know wher h matrix starts 
+	ArrayList<index_t> rows_;
+	ArrayList<index_t> columns_;
+  ArrayList<double>  values_;
+  Vector ineqeq_lagrange_mult_; 
 };
 
 #include "nmf_objectives_impl.h"
diff --git a/fastlib2/contrib/nvasil/convex_nmf/nmf_objectives_impl.h b/fastlib2/contrib/nvasil/convex_nmf/nmf_objectives_impl.h
index 8462b0900a..a311de34db 100644
--- a/fastlib2/contrib/nvasil/convex_nmf/nmf_objectives_impl.h
+++ b/fastlib2/contrib/nvasil/convex_nmf/nmf_objectives_impl.h
@@ -1,33 +1,33 @@
 
-void BigSdpNmfObjective::Init(fx_module *module, 
+void BigSdpNmfObjectiveMaxVar::Init(fx_module *module, 
 			ArrayList<index_t> &rows,
 			ArrayList<index_t> &columns,
       ArrayList<double>  &values) {
   module_=module;
-	rows_.InitCopy(rows);
+	new_dim_=fx_param_int(module_, "new_dim", 5); 
+  rank_=fx_param_int(module_, "rank", 3);
+  rows_.InitCopy(rows);
 	columns_.InitCopy(columns);
 	values_.InitCopy(values);
 	num_of_rows_=*std::max_element(rows_.begin(), rows_.end())+1;
 	num_of_columns_=*std::max_element(columns_.begin(), columns_.end())+1;
 	eq_lagrange_mult_.Init(values_.size());
-	for(index_t i=0; i<eq_lagrange_mult_.length(); i++) {
-    eq_lagrange_mult_[i]=math::Random(-0.0, 1.0);
-  }
-  rank_=fx_param_int(module_, "rank", 3);
-	new_dim_=fx_param_int(module_, "new_dim", 5);
+  eq_lagrange_mult_.SetAll(0);
 	offset_h_ = num_of_rows_*new_dim_;
 }
 
-void BigSdpNmfObjective::Destruct() {
+void BigSdpNmfObjectiveMaxVar::Destruct() {
   rows_.Renew();
 	columns_.Renew();
   values_.Renew();
 	eq_lagrange_mult_.Destruct();
 }
 
-void BigSdpNmfObjective::ComputeGradient(Matrix &coordinates, Matrix *gradient) {
+void BigSdpNmfObjectiveMaxVar::ComputeGradient(Matrix &coordinates, Matrix *gradient) {
+  // gradient of the variance
   gradient->CopyValues(coordinates);
-	la::Scale(2.0, gradient);
+	la::Scale(-2.0, gradient);
+  // Gradient on the equalities and inequalities  
 	for(index_t i=0; i<values_.size(); i++) {
     index_t w =rows_[i]*new_dim_;
     index_t h =offset_h_+columns_[i]*new_dim_;
@@ -38,19 +38,18 @@ void BigSdpNmfObjective::ComputeGradient(Matrix &coordinates, Matrix *gradient)
 				 coordinates.GetColumnPtr(w), gradient->GetColumnPtr(h));	
 		la::AddExpert(rank_*new_dim_, -eq_lagrange_mult_[i]+2*sigma_*diff, 
 				coordinates.GetColumnPtr(h), gradient->GetColumnPtr(w));
-	}
-	
+  }	
 }
 
-void BigSdpNmfObjective::ComputeObjective(Matrix &coordinates, double *objective) {
+void BigSdpNmfObjectiveMaxVar::ComputeObjective(Matrix &coordinates, double *objective) {
   *objective=0;	
-	for(index_t i=0; i<coordinates.n_cols(); i++) {
-	   *objective+=la::Dot(coordinates.n_rows(), 
+  for(index_t i=0; i<coordinates.n_cols(); i++) {
+	   *objective-=la::Dot(coordinates.n_rows(), 
 				 coordinates.GetColumnPtr(i), coordinates.GetColumnPtr(i));
 	}
 }
 
-void BigSdpNmfObjective::ComputeFeasibilityError(Matrix &coordinates, double *error) {
+void BigSdpNmfObjectiveMaxVar::ComputeFeasibilityError(Matrix &coordinates, double *error) {
   *error=0;
 	for(index_t i=0; i<values_.size(); i++) {
     index_t w =rows_[i]*new_dim_;
@@ -59,14 +58,14 @@ void BigSdpNmfObjective::ComputeFeasibilityError(Matrix &coordinates, double *er
 				                coordinates.GetColumnPtr(w),
 												coordinates.GetColumnPtr(h))-values_[i];
 		*error+=diff*diff;
-   		
-	}
+ 	}
 }
 
-double BigSdpNmfObjective::ComputeLagrangian(Matrix &coordinates) {
+double BigSdpNmfObjectiveMaxVar::ComputeLagrangian(Matrix &coordinates) {
   double lagrangian=0;
   ComputeObjective(coordinates, &lagrangian);
   for(index_t i=0; i<values_.size(); i++) {
+    // equalities
     index_t w =rows_[i]*new_dim_;
     index_t h =offset_h_+columns_[i]*new_dim_;
     double diff=la::Dot(rank_*new_dim_,
@@ -77,70 +76,27 @@ double BigSdpNmfObjective::ComputeLagrangian(Matrix &coordinates) {
 	return lagrangian;
 }
 
-void BigSdpNmfObjective::UpdateLagrangeMult(Matrix &coordinates) {
+void BigSdpNmfObjectiveMaxVar::UpdateLagrangeMult(Matrix &coordinates) {
   for(index_t i=0; i<values_.size(); i++) {
     index_t w =rows_[i]*new_dim_;
     index_t h =offset_h_+columns_[i]*new_dim_;
+    // equalities
     double diff=la::Dot(rank_*new_dim_,
 				                coordinates.GetColumnPtr(w),
 												coordinates.GetColumnPtr(h))-values_[i];
 		eq_lagrange_mult_[i]-=sigma_*diff;
-	}
+ 	}
 }
 
-void BigSdpNmfObjective::Project(Matrix *coordinates) {
- for(index_t i=0; i<coordinates->n_cols(); i++) {
-    for(index_t j=0; j< coordinates->n_rows(); j++) {
-	    if (coordinates->get(j, i)<0.0) {
-		    coordinates->set(j, i, 0.0);
-		  }
-    }
-	}
- /* Matrix u_mat, vt_mat;
-  Vector s;
-  success_t success=la::SVDInit(*coordinates, &s, &u_mat, &vt_mat);
-  if (success==SUCCESS_FAIL) {
-    FATAL("Svd failed...\n");
-  }
-  index_t positives=0;
-  index_t negatives=0;
-  for(index_t i=0; i<vt_mat.n_cols(); i++) {
-    if (vt_mat.get(0, i)<0) {
-     negatives++;
-    } else {
-      positives++;
-    }
-  }
-  if (positives>=negatives) {
-    for(index_t i=0; i<vt_mat.n_cols(); i++) {
-      if (vt_mat.get(0, i)<0) {
-        vt_mat.set(0, i, vt_mat.get(0, i));
-      }
-    }
-  } else {
-    for(index_t i=0; i<vt_mat.n_cols(); i++) {
-      if (vt_mat.get(0, i)>0) {
-        vt_mat.set(0, i, vt_mat.get(0, i));
-      }
-    }
-  }
-
-  Matrix s_mat;
-  s_mat.Init(s.length(), s.length());
-  s_mat.SetDiagonal(s);
-  Matrix tmp;
-  coordinates->Destruct();
-  la::MulInit(u_mat, s_mat, &tmp);
-  la::MulInit(tmp, vt_mat, coordinates);  
-*/
-  
+void BigSdpNmfObjectiveMaxVar::Project(Matrix *coordinates) {
+  OptUtils::NonNegativeProjection(coordinates);
 }
 
-void BigSdpNmfObjective::set_sigma(double sigma) {
+void BigSdpNmfObjectiveMaxVar::set_sigma(double sigma) {
   sigma_=sigma;
 }
 
-void BigSdpNmfObjective::GiveInitMatrix(Matrix *init_data) {
+void BigSdpNmfObjectiveMaxVar::GiveInitMatrix(Matrix *init_data) {
   init_data->Init(rank_, (num_of_rows_+num_of_columns_)*new_dim_);
 	for(index_t i=0; i<init_data->n_rows(); i++) {
 	  for(index_t j=0; j<init_data->n_cols(); j++) {
@@ -150,7 +106,240 @@ void BigSdpNmfObjective::GiveInitMatrix(Matrix *init_data) {
 }
 
 
-bool BigSdpNmfObjective::IsDiverging(double objective) {
+bool BigSdpNmfObjectiveMaxVar::IsDiverging(double objective) {
   return false;
 } 
 
+///////////////////////////////////////////////////////////////////////////
+
+void BigSdpNmfObjectiveMinVar::Init(fx_module *module, 
+			ArrayList<index_t> &rows,
+			ArrayList<index_t> &columns,
+      ArrayList<double>  &values) {
+  module_=module;
+	new_dim_=fx_param_int(module_, "new_dim", 5); 
+  rank_=fx_param_int(module_, "rank", 3);
+  rows_.InitCopy(rows);
+	columns_.InitCopy(columns);
+	values_.InitCopy(values);
+	num_of_rows_=*std::max_element(rows_.begin(), rows_.end())+1;
+	num_of_columns_=*std::max_element(columns_.begin(), columns_.end())+1;
+	eq_lagrange_mult_.Init(values_.size());
+  eq_lagrange_mult_.SetAll(0);
+	offset_h_ = num_of_rows_*new_dim_;
+}
+
+void BigSdpNmfObjectiveMinVar::Destruct() {
+  rows_.Renew();
+	columns_.Renew();
+  values_.Renew();
+	eq_lagrange_mult_.Destruct();
+}
+
+void BigSdpNmfObjectiveMinVar::ComputeGradient(Matrix &coordinates, 
+    Matrix *gradient) {
+  // gradient of the variance
+  gradient->CopyValues(coordinates);
+	la::Scale(2.0, gradient);
+  // Gradient on the equalities and inequalities  
+	for(index_t i=0; i<values_.size(); i++) {
+    index_t w =rows_[i]*new_dim_;
+    index_t h =offset_h_+columns_[i]*new_dim_;
+    double diff=la::Dot(rank_*new_dim_,
+				                coordinates.GetColumnPtr(w),
+												coordinates.GetColumnPtr(h))-values_[i];
+	  la::AddExpert(rank_*new_dim_, -eq_lagrange_mult_[i]+2*sigma_*diff,
+				 coordinates.GetColumnPtr(w), gradient->GetColumnPtr(h));	
+		la::AddExpert(rank_*new_dim_, -eq_lagrange_mult_[i]+2*sigma_*diff, 
+				coordinates.GetColumnPtr(h), gradient->GetColumnPtr(w));
+  }	
+}
+
+void BigSdpNmfObjectiveMinVar::ComputeObjective(Matrix &coordinates, 
+    double *objective) {
+  *objective=0;	
+  for(index_t i=0; i<coordinates.n_cols(); i++) {
+	   *objective+=la::Dot(coordinates.n_rows(), 
+				 coordinates.GetColumnPtr(i), coordinates.GetColumnPtr(i));
+	}
+}
+
+void BigSdpNmfObjectiveMinVar::ComputeFeasibilityError(Matrix &coordinates, 
+    double *error) {
+  *error=0;
+	for(index_t i=0; i<values_.size(); i++) {
+    index_t w =rows_[i]*new_dim_;
+    index_t h =offset_h_+columns_[i]*new_dim_;
+    double diff=la::Dot(rank_*new_dim_,
+				                coordinates.GetColumnPtr(w),
+												coordinates.GetColumnPtr(h))-values_[i];
+		*error+=diff*diff;
+ 	}
+}
+
+double BigSdpNmfObjectiveMinVar::ComputeLagrangian(Matrix &coordinates) {
+  double lagrangian=0;
+  ComputeObjective(coordinates, &lagrangian);
+  for(index_t i=0; i<values_.size(); i++) {
+    // equalities
+    index_t w =rows_[i]*new_dim_;
+    index_t h =offset_h_+columns_[i]*new_dim_;
+    double diff=la::Dot(rank_*new_dim_,
+				                coordinates.GetColumnPtr(w),
+												coordinates.GetColumnPtr(h))-values_[i];
+		lagrangian+=(sigma_*diff-eq_lagrange_mult_[i])*diff;
+	}
+	return lagrangian;
+}
+
+void BigSdpNmfObjectiveMinVar::UpdateLagrangeMult(Matrix &coordinates) {
+  for(index_t i=0; i<values_.size(); i++) {
+    index_t w =rows_[i]*new_dim_;
+    index_t h =offset_h_+columns_[i]*new_dim_;
+    // equalities
+    double diff=la::Dot(rank_*new_dim_,
+				                coordinates.GetColumnPtr(w),
+												coordinates.GetColumnPtr(h))-values_[i];
+		eq_lagrange_mult_[i]-=sigma_*diff;
+ 	}
+}
+
+void BigSdpNmfObjectiveMinVar::Project(Matrix *coordinates) {
+  OptUtils::NonNegativeProjection(coordinates);
+}
+
+void BigSdpNmfObjectiveMinVar::set_sigma(double sigma) {
+  sigma_=sigma;
+}
+
+void BigSdpNmfObjectiveMinVar::GiveInitMatrix(Matrix *init_data) {
+  init_data->Init(rank_, (num_of_rows_+num_of_columns_)*new_dim_);
+	for(index_t i=0; i<init_data->n_rows(); i++) {
+	  for(index_t j=0; j<init_data->n_cols(); j++) {
+		  init_data->set(i, j, math::Random());
+		}
+	}
+}
+
+
+bool BigSdpNmfObjectiveMinVar::IsDiverging(double objective) {
+  return false;
+} 
+
+///////////////////////////////////////////////////////////////////////////////
+void BigSdpNmfObjectiveMinVarIneq::Init(fx_module *module, 
+			ArrayList<index_t> &rows,
+			ArrayList<index_t> &columns,
+      ArrayList<double>  &values) {
+  module_=module;
+	new_dim_=fx_param_int(module_, "new_dim", 5); 
+  rank_=fx_param_int(module_, "rank", 3);
+  rows_.InitCopy(rows);
+	columns_.InitCopy(columns);
+	values_.InitCopy(values);
+	num_of_rows_=*std::max_element(rows_.begin(), rows_.end())+1;
+	num_of_columns_=*std::max_element(columns_.begin(), columns_.end())+1;
+	eq_lagrange_mult_.Init(values_.size());
+  eq_lagrange_mult_.SetAll(0);
+	offset_h_ = num_of_rows_*new_dim_;
+}
+
+void BigSdpNmfObjectiveMinVar::Destruct() {
+  rows_.Renew();
+	columns_.Renew();
+  values_.Renew();
+	eq_lagrange_mult_.Destruct();
+}
+
+void BigSdpNmfObjectiveMinVarIneq::ComputeGradient(Matrix &coordinates, 
+    Matrix *gradient) {
+  // gradient of the variance
+  gradient->CopyValues(coordinates);
+	la::Scale(2.0, gradient);
+  // Gradient on the equalities and inequalities  
+	for(index_t i=0; i<values_.size(); i++) {
+    index_t w =rows_[i]*new_dim_;
+    index_t h =offset_h_+columns_[i]*new_dim_;
+    double diff=la::Dot(rank_*new_dim_,
+				                coordinates.GetColumnPtr(w),
+												coordinates.GetColumnPtr(h))-values_[i];
+	  la::AddExpert(rank_*new_dim_, -eq_lagrange_mult_[i]+2*sigma_*diff,
+				 coordinates.GetColumnPtr(w), gradient->GetColumnPtr(h));	
+		la::AddExpert(rank_*new_dim_, -eq_lagrange_mult_[i]+2*sigma_*diff, 
+				coordinates.GetColumnPtr(h), gradient->GetColumnPtr(w));
+  }	
+}
+
+void BigSdpNmfObjectiveMinVarIneq::ComputeObjective(Matrix &coordinates, 
+    double *objective) {
+  *objective=0;	
+  for(index_t i=0; i<coordinates.n_cols(); i++) {
+	   *objective+=la::Dot(coordinates.n_rows(), 
+				 coordinates.GetColumnPtr(i), coordinates.GetColumnPtr(i));
+	}
+}
+
+void BigSdpNmfObjectiveMinVarIneq::ComputeFeasibilityError(Matrix &coordinates, 
+    double *error) {
+  *error=0;
+	for(index_t i=0; i<values_.size(); i++) {
+    index_t w =rows_[i]*new_dim_;
+    index_t h =offset_h_+columns_[i]*new_dim_;
+    double diff=la::Dot(rank_*new_dim_,
+				                coordinates.GetColumnPtr(w),
+												coordinates.GetColumnPtr(h))-values_[i];
+		*error+=diff*diff;
+ 	}
+}
+
+
+double BigSdpNmfObjectiveMinVarIneq::ComputeLagrangian(Matrix &coordinates) {
+  double lagrangian=0;
+  ComputeObjective(coordinates, &lagrangian);
+  for(index_t i=0; i<values_.size(); i++) {
+    // equalities
+    index_t w =rows_[i]*new_dim_;
+    index_t h =offset_h_+columns_[i]*new_dim_;
+    double diff=la::Dot(rank_*new_dim_,
+				                coordinates.GetColumnPtr(w),
+												coordinates.GetColumnPtr(h))-values_[i];
+		lagrangian+=(sigma_*diff-eq_lagrange_mult_[i])*diff;
+	}
+	return lagrangian;
+}
+
+void BigSdpNmfObjectiveMinVarIneq::UpdateLagrangeMult(Matrix &coordinates) {
+  for(index_t i=0; i<values_.size(); i++) {
+    index_t w =rows_[i]*new_dim_;
+    index_t h =offset_h_+columns_[i]*new_dim_;
+    // equalities
+    double diff=la::Dot(rank_*new_dim_,
+				                coordinates.GetColumnPtr(w),
+												coordinates.GetColumnPtr(h))-values_[i];
+		eq_lagrange_mult_[i]-=sigma_*diff;
+ 	}
+}
+
+void BigSdpNmfObjectiveMinVarIneq::Project(Matrix *coordinates) {
+  OptUtils::NonNegativeProjection(coordinates);
+}
+
+void BigSdpNmfObjectiveMinVarIneq::set_sigma(double sigma) {
+  sigma_=sigma;
+}
+
+void BigSdpNmfObjectiveMinVarIneq::GiveInitMatrix(Matrix *init_data) {
+  init_data->Init(rank_, (num_of_rows_+num_of_columns_)*new_dim_);
+	for(index_t i=0; i<init_data->n_rows(); i++) {
+	  for(index_t j=0; j<init_data->n_cols(); j++) {
+		  init_data->set(i, j, math::Random());
+		}
+	}
+}
+
+
+bool BigSdpNmfObjectiveMinVarIneq::IsDiverging(double objective) {
+  return false;
+} 
+
+