diff --git a/fastlib2/contrib/nvasil/convex_nmf/main.cc b/fastlib2/contrib/nvasil/convex_nmf/main.cc
index 2d21313e8b..d108df94b0 100644
--- a/fastlib2/contrib/nvasil/convex_nmf/main.cc
+++ b/fastlib2/contrib/nvasil/convex_nmf/main.cc
@@ -24,15 +24,16 @@ int main(int argc, char *argv[]) {
 	fx_root=fx_init(argc, argv, NULL);
 	fx_module *nmf_module=fx_submodule(fx_root, "/engine");
   SdpNmfEngine<SmallSdpNmf> engine;
-  fx_set_param_str(nmf_module, "data_file", "1.csv");
+  fx_set_param_str(nmf_module, "data_file", "5.csv");
   //    "/net/hg200/nvasil/dataset/orl_faces/orl_test_faces_100.csv");
-  fx_set_param_int(nmf_module, "new_dim", 1);
+  fx_set_param_int(nmf_module, "new_dim", 2);
 	fx_set_param_double(nmf_module, "l_bfgs/sigma", 1);
-  fx_set_param_double(nmf_module, "l_bfgs/gamma", 1.7);
-  fx_set_param_double(nmf_module, "l_bfgs/wolfe_sigma2", 0.9999);
+  fx_set_param_double(nmf_module, "l_bfgs/gamma", 2);
+  fx_set_param_double(nmf_module, "l_bfgs/wolfe_sigma2", 0.99);
+  fx_set_param_double(nmf_module, "l_bfgs/wolfe_sigma1", 1e-4);
   fx_set_param_double(nmf_module, "l_bfgs/use_default_termination", false);
-  fx_set_param_int(nmf_module, "l_bfgs/mem_bfgs", 10); 
-  fx_set_param_double(nmf_module, "optfun/desired_duality_gap", 1e-3); 
+  fx_set_param_int(nmf_module, "l_bfgs/mem_bfgs", 4); 
+  fx_set_param_double(nmf_module, "optfun/desired_duality_gap", 1e-4); 
 	fx_set_param_double(nmf_module, "optfun/gradient_tolerance", 1e-2);
   engine.Init(nmf_module);
   engine.ComputeNmf();
@@ -40,6 +41,6 @@ int main(int argc, char *argv[]) {
 	Matrix h_mat;
   engine.GetW(&w_mat);
 	engine.GetH(&h_mat);
-  data::Save("W.csv", w_mat);
-	data::Save("H.csv", h_mat);
+  data::Save("w.csv", w_mat);
+	data::Save("h.csv", h_mat);
 }
diff --git a/fastlib2/contrib/nvasil/convex_nmf/sdp_nmf_engine.h b/fastlib2/contrib/nvasil/convex_nmf/sdp_nmf_engine.h
index 20be290f6d..93303d9f0f 100644
--- a/fastlib2/contrib/nvasil/convex_nmf/sdp_nmf_engine.h
+++ b/fastlib2/contrib/nvasil/convex_nmf/sdp_nmf_engine.h
@@ -25,6 +25,14 @@
 template<typename SdpNmfObjective>
 class SdpNmfEngine {
  public:
+   SdpNmfEngine() {
+     engine_=NULL;
+   }
+   ~SdpNmfEngine() {
+     if (engine_!=NULL) {
+       delete engine_;
+     }
+   }
 	 void  Init(fx_module *module) {
  		 module_=module;
 	   std::string data_file=fx_param_str_req(module_, "data_file");
@@ -40,21 +48,36 @@ class SdpNmfEngine {
 		 fx_module *opt_function_module=fx_submodule(module_, "optfun");
 		 fx_set_param_int(opt_function_module, "new_dim", new_dim_);
 		 opt_function_.Init(opt_function_module, rows_, columns_, values_);
-		 fx_module *l_bfgs_module=fx_submodule(module_, "l_bfgs");
-     Matrix init_data;
-		 opt_function_.GiveInitMatrix(&init_data);
-	   fx_set_param_int(l_bfgs_module, "num_of_points", init_data.n_cols());
-		 fx_set_param_int(l_bfgs_module, "new_dimension", init_data.n_rows());
-     engine_.Init(&opt_function_, l_bfgs_module);
-	 	 engine_.set_coordinates(init_data);
+		 l_bfgs_module_=fx_submodule(module_, "l_bfgs");
 	 }
 	 void Destruct() {
-	   
+     if (engine_!=NULL) {
+	     delete engine_;
+       engine_ = NULL;
+     }
+     SdpNmfObjective  opt_function_;
+	   rows_.Renew();
+	   columns_.Renew();
+	   values_.Renew();
+	   w_mat_.Destruct();
+	   h_mat_.Destruct();
+
 	 };
 	 void ComputeNmf() {
-     engine_.ComputeLocalOptimumBFGS();
      Matrix result;
-     engine_.GetResults(&result);
+		 opt_function_.GiveInitMatrix(&result);
+	   fx_set_param_int(l_bfgs_module_, "num_of_points", result.n_cols());
+		 fx_set_param_int(l_bfgs_module_, "new_dimension", result.n_rows());
+
+     for(double sigma=1.0; sigma<=1e4; sigma*=4) {
+       engine_=new LBfgs<SdpNmfObjective>();
+       engine_->Init(&opt_function_, l_bfgs_module_);
+       engine_->set_coordinates(result);
+       result.Destruct();
+       engine_->set_sigma(sigma);
+       engine_->ComputeLocalOptimumBFGS();
+       engine_->GetResults(&result);
+     }   
 		 w_mat_.Init(new_dim_, num_of_rows_);
 		 h_mat_.Init(new_dim_, num_of_columns_);
      w_mat_.CopyColumnFromMat(0, 0, num_of_rows_, result); 
@@ -72,8 +95,6 @@ class SdpNmfEngine {
        error+=fabs(v_rec.get(r, c)-values_[i]);
        v_sum+=values_[i];
      }
-     data::Save("w.csv", w_mat_);
-     data::Save("h.csv", h_mat_);
      NOTIFY("Reconstruction error: %lg%%\n", error*100/v_sum);
 	 }
    
@@ -86,7 +107,8 @@ class SdpNmfEngine {
 	 
  private:
 	fx_module *module_;
-  LBfgs<SdpNmfObjective> engine_;
+  fx_module *l_bfgs_module_;
+  LBfgs<SdpNmfObjective> *engine_;
 	SdpNmfObjective  opt_function_;
 	ArrayList<index_t> rows_;
 	ArrayList<index_t> columns_;
diff --git a/fastlib2/contrib/nvasil/convex_nmf/sdp_objectives_impl.h b/fastlib2/contrib/nvasil/convex_nmf/sdp_objectives_impl.h
index b990657012..5dc404b39b 100644
--- a/fastlib2/contrib/nvasil/convex_nmf/sdp_objectives_impl.h
+++ b/fastlib2/contrib/nvasil/convex_nmf/sdp_objectives_impl.h
@@ -41,32 +41,32 @@ void SmallSdpNmf::Init(fx_module *module,
   // W
   for(index_t i=0; i<num_of_rows_; i++) {
     for(index_t j=0; j<new_dim_; j++) {
-      objective_factor_.set(j, i, -1);
+      objective_factor_.set(j, i, -1.0);
     }
   }
   // H
   for(index_t i=offset_h_; i<offset_h_+num_of_columns_; i++) {
     for(index_t j=0; j<new_dim_; j++) {
-      objective_factor_.set(j, i, -1);
+      objective_factor_.set(j, i, -1.0);
     }
   }
   // tw
   for(index_t i=offset_tw_; i<offset_tw_+num_of_rows_; i++) {
     for(index_t j=0; j<new_dim_; j++) {
-      objective_factor_.set(j, i, 1);
+      objective_factor_.set(j, i, 1.0);
     }
   }
   // th
   for(index_t i=offset_th_; i<offset_th_+num_of_rows_; i++) {
     for(index_t j=0; j<new_dim_; j++) {
-      objective_factor_.set(j, i, 1);
+      objective_factor_.set(j, i, 1.0);
     }
   }
  
   // v
   for(index_t i=offset_v_; i<offset_v_+values_.size(); i++) {
     for(index_t j=0; j<new_dim_; j++) {
-      objective_factor_.set(j, i, 1);
+      objective_factor_.set(j, i, 1.0);
     }
   }
  
@@ -84,6 +84,7 @@ void SmallSdpNmf::ComputeGradient(Matrix &coordinates, Matrix *gradient) {
       diff+=coordinates.get(j, v_i);
     }
     diff-=values_[i];
+    //NOTIFY("diff:%lg", diff);
     for(index_t j=0; j<new_dim_; j++) {
       gradient->set(j, v_i, gradient->get(j, v_i)-1.0/diff);
     }     
@@ -110,16 +111,18 @@ void SmallSdpNmf::ComputeGradient(Matrix &coordinates, Matrix *gradient) {
       double t1_minus_ww=(t1-w*w);
       double t2_minus_hh=(t2-h*h);
       double wh_minus_v=(w*h-v);
-      double determinant=t1_minus_ww*t2_minus_hh-math::Pow<2,1>(wh_minus_v);
-      NOTIFY("determinant:%lg", determinant);
-      DEBUG_ERR_MSG_IF(determinant==0.0, "Determinant equal to zero");
-      if (determinant==0.0) {
+      double determinant=t1_minus_ww*t2_minus_hh
+        -math::Pow<2,1>(wh_minus_v);
+      
+      //NOTIFY("determinant:%lg", determinant);
+      if (determinant<=0.0) {
         gradient->set(j, w_i, -DBL_MAX);
         gradient->set(j, h_i, -DBL_MAX);
         gradient->set(j, t1_i, -DBL_MAX);
         gradient->set(j, t2_i, -DBL_MAX);
         gradient->set(j, v_i, -DBL_MAX);  
-        NONFATAL("Determinant is zero");    
+        NONFATAL("Determinant is less or equal to zero");
+        return;    
       } else {     
         double dw=(-2*w*(t2_minus_hh)-2*h*(wh_minus_v))/determinant;
         double dh=(-2*h*(t1_minus_ww)-2*w*(wh_minus_v))/determinant;
@@ -129,7 +132,7 @@ void SmallSdpNmf::ComputeGradient(Matrix &coordinates, Matrix *gradient) {
         gradient->set(j, w_i, gradient->get(j, w_i)-dw);
         gradient->set(j, h_i, gradient->get(j, h_i)-dh);
         gradient->set(j, t1_i, gradient->get(j, t1_i)-dt1);
-        gradient->set(j, t2_i, gradient->get(j, t1_i)-dt2);
+        gradient->set(j, t2_i, gradient->get(j, t2_i)-dt2);
         gradient->set(j, v_i, gradient->get(j, v_i)-dv);  
       }
     }  
@@ -144,6 +147,10 @@ void SmallSdpNmf::ComputeObjective(Matrix &coordinates, double *objective) {
 void SmallSdpNmf::ComputeFeasibilityError(Matrix &coordinates, double *error) {
   // return duality gap instead
   *error=number_of_cones_/sigma_;
+  double lagrangian=ComputeLagrangian(coordinates);
+  double objective;
+  ComputeObjective(coordinates, &objective);
+  NOTIFY("sum of log barriers:%lg ", lagrangian-sigma_*objective);
 }
 
 double SmallSdpNmf::ComputeLagrangian(Matrix &coordinates) {
@@ -152,7 +159,9 @@ double SmallSdpNmf::ComputeLagrangian(Matrix &coordinates) {
   ComputeObjective(coordinates, &lagrangian);
   lagrangian*=sigma_;
   // from the LP cones
+  double temp_prod=1.0;
   for(index_t i=0; i<values_.size(); i++) {
+    
     index_t v_i=offset_v_+i;
     double diff=0;
     for(index_t j=0; j<new_dim_; j++) {
@@ -162,13 +171,12 @@ double SmallSdpNmf::ComputeLagrangian(Matrix &coordinates) {
     if unlikely(diff<0) {
       return DBL_MAX;
     }
-//    DEBUG_ERR_MSG_IF(diff<=0, "LP cone is invalid, you are "
-//       " out of the feasible region, constraint %i, diff %lg",
-//       i, diff);
-    lagrangian-=log(diff+1e-100);   
-  } 
+    temp_prod*=diff;
+  }
+  lagrangian-=log(temp_prod);   
   // from the SDP cones
   // determinant=(t1-w*w)*(t2-h*h)-math::Pow<2,1>(w*h-v);
+  temp_prod=1.0;
  for(index_t i=0; i<values_.size(); i++) {
     index_t w_i=rows_[i];
     index_t h_i=offset_h_+columns_[i];
@@ -183,16 +191,25 @@ double SmallSdpNmf::ComputeLagrangian(Matrix &coordinates) {
       double v=coordinates.get(j, v_i);
       double t1_minus_ww=(t1-w*w);
       double t2_minus_hh=(t2-h*h);
+      // we need this, because during Wolfe step
+      // we might accidently cross the barrier 
+      // and never realize it since the determinant
+      // will stay positive
+      if (t1_minus_ww<0 || t2_minus_hh<0) {
+        return DBL_MAX;
+      }
       double wh_minus_v=(w*h-v);
-      double determinant=t1_minus_ww*t2_minus_hh-math::Pow<2,1>(wh_minus_v);
-      if (unlikely(determinant<0)) {
+      double determinant=t1_minus_ww*t2_minus_hh
+        -math::Pow<2,1>(wh_minus_v);
+      if (unlikely(determinant<=0)) {
         return DBL_MAX;
       }
 //      DEBUG_ERR_MSG_IF(determinant<=0, "SDP cone is invalid, you are "
 //       " out of the feasible region");
-      lagrangian-=log(determinant+1e-100);
+      temp_prod*=determinant;
     }  
   }
+  lagrangian-=log(temp_prod);
   return lagrangian;
 }
 
@@ -201,7 +218,7 @@ void SmallSdpNmf::UpdateLagrangeMult(Matrix &coordinates) {
 }
 
 void SmallSdpNmf::Project(Matrix *coordinates) {
-  OptUtils::NonNegativeProjection(coordinates);
+  //OptUtils::NonNegativeProjection(coordinates);
 }
 
 void SmallSdpNmf::set_sigma(double sigma) {
@@ -213,7 +230,7 @@ void SmallSdpNmf::GiveInitMatrix(Matrix *init_data) {
       2*(num_of_rows_+num_of_columns_)+values_.size());
   for(index_t i=0; i<num_of_rows_+num_of_columns_; i++) {
     for(index_t j=0; j<new_dim_; j++) {
-      init_data->set(j, i, math::Random(0.0, 1.0));
+      init_data->set(j, i, 0*math::Random(0.0, 1.0));
     }
   }
   for(index_t i=0; i<values_.size(); i++) {
@@ -244,7 +261,7 @@ bool SmallSdpNmf::IsDiverging(double objective) {
 
 bool SmallSdpNmf::IsOptimizationOver(Matrix &coordinates, 
     Matrix &gradient, double step) {
-  double norm_gradient = la::Dot(gradient.n_elements(), 
+/*  double norm_gradient = la::Dot(gradient.n_elements(), 
       gradient.ptr(), gradient.ptr());
   //  one of our barriers is zero
   if (norm_gradient>=DBL_MAX) {
@@ -255,6 +272,8 @@ bool SmallSdpNmf::IsOptimizationOver(Matrix &coordinates,
   } else {
     return false;
   }
+  */
+  return true;
 }
 
 bool SmallSdpNmf::IsIntermediateStepOver(Matrix &coordinates, 
diff --git a/fastlib2/contrib/nvasil/l_bfgs/l_bfgs_impl.h b/fastlib2/contrib/nvasil/l_bfgs/l_bfgs_impl.h
index 536459b660..56e73b3674 100644
--- a/fastlib2/contrib/nvasil/l_bfgs/l_bfgs_impl.h
+++ b/fastlib2/contrib/nvasil/l_bfgs/l_bfgs_impl.h
@@ -40,7 +40,7 @@ void LBfgs<OptimizedFunction>::Init(OptimizedFunction *optimized_function,
   }
   DEBUG_ASSERT(wolfe_sigma1_>0);
   DEBUG_ASSERT(wolfe_sigma1_<1);
-  DEBUG_ASSERT(wolfe_sigma2_>0);
+  DEBUG_ASSERT(wolfe_sigma2_>wolfe_sigma1_);
   DEBUG_ASSERT(wolfe_sigma2_<1);
   if (unlikely(wolfe_sigma1_>=wolfe_sigma2_)) {
     FATAL("Wolfe sigma1 %lg should be less than sigma2 %lg", 
@@ -48,7 +48,7 @@ void LBfgs<OptimizedFunction>::Init(OptimizedFunction *optimized_function,
   }
   DEBUG_ASSERT(wolfe_sigma1_>0);
   DEBUG_ASSERT(wolfe_sigma1_<1);
-  DEBUG_ASSERT(wolfe_sigma2_>0);
+  DEBUG_ASSERT(wolfe_sigma2_>wolfe_sigma1_);
   DEBUG_ASSERT(wolfe_sigma2_<1);
   norm_grad_tolerance_ = fx_param_double(module_, "norm_grad_tolerance", 0.1); 
   wolfe_beta_   = fx_param_double(module_, "wolfe_beta", 0.8);
@@ -108,7 +108,7 @@ void LBfgs<OptimizedFunction>::ComputeLocalOptimumBFGS() {
     success_t success=ComputeBFGS_(&step_, gradient_, i);
     if (success==SUCCESS_FAIL) {
       NOTIFY("LBFGS failed to find a direction, continuing with gradient descent\n");
-      ComputeWolfeStep_(&step_, gradient_);
+      //ComputeWolfeStep_(&step_, gradient_);
     }
     optimized_function_->ComputeGradient(coordinates_, &gradient_);
     UpdateBFGS_();
@@ -118,7 +118,7 @@ void LBfgs<OptimizedFunction>::ComputeLocalOptimumBFGS() {
     if (silent_==false) {
         ReportProgressFile_();
     }
-    if (use_default_termination_== true) {
+/*    if (use_default_termination_== true) {
        optimized_function_->ComputeFeasibilityError(coordinates_,
           &feasibility_error); 
       if (feasibility_error < desired_feasibility_) {
@@ -132,6 +132,7 @@ void LBfgs<OptimizedFunction>::ComputeLocalOptimumBFGS() {
         return;
       }
     }
+*/
   }
  
   NOTIFY("Now starting optimizing with BFGS...\n");
@@ -145,8 +146,11 @@ void LBfgs<OptimizedFunction>::ComputeLocalOptimumBFGS() {
           gradient_.ptr(), gradient_.ptr());
       num_of_iterations_++;
       if (success_bfgs==SUCCESS_FAIL){
-        NOTIFY("LBFGS failed to find a direction, continuing with gradient descent\n");
-        ComputeWolfeStep_(&step_, gradient_);
+         NOTIFY("LBFGS failed to find a direction, continuing with gradient descent\n");
+        // if  (ComputeWolfeStep_(&step_, gradient_)==SUCCESS_FAIL) {
+        //   NONFATAL("Gradient descent failed too");
+        // }
+        break;
       }
       if (silent_==false) {
         ReportProgressFile_();
@@ -357,10 +361,10 @@ success_t LBfgs<OptimizedFunction>::ComputeBFGS_(double *step, Matrix &grad, ind
   double y_y = la::Dot(num_of_points_* 
                    new_dimension_, y_bfgs_[index_bfgs_].ptr(),
                    y_bfgs_[index_bfgs_].ptr());
-  if (unlikely(y_y<1e-10)){
+  if (unlikely(y_y<1e-40)){
     NONFATAL("Gradient differences close to singular...norm=%lg\n", y_y);
   } 
-  double norm_scale=s_y/(y_y+1e-10);
+  double norm_scale=s_y/(y_y+1e-40);
   la::Scale(norm_scale, &temp_direction);
   Matrix scaled_s;
   double beta;