diff --git a/fastlib2/contrib/nvasil/non_convex_mvu/non_convex_mvu.h b/fastlib2/contrib/nvasil/non_convex_mvu/non_convex_mvu.h
index 8b08666329..b38b4d8276 100644
--- a/fastlib2/contrib/nvasil/non_convex_mvu/non_convex_mvu.h
+++ b/fastlib2/contrib/nvasil/non_convex_mvu/non_convex_mvu.h
@@ -31,7 +31,7 @@ class NonConvexMVU {
   void Init(std::string data_file, index_t knns);
   void Init(std::string data_file, index_t knns, index_t leaf_size);
   void ComputeLocalOptimum();
-  void ComputeLocalOptimumBFGS_();
+  void ComputeLocalOptimumBFGS();
   // eta < 1
   void set_eta(double eta);
   // gamma > 1
@@ -48,6 +48,10 @@ class NonConvexMVU {
    *  beta for armijo rule somewhere between 0.5 to 0.1
    */
   void set_armijo_beta(double armijo_beta);
+  /**
+   * Set the memory for the BFGS method
+   */
+  void set_mem_bfgs(index_t mem_bfgs);
   Matrix &coordinates();
 
  private:
@@ -79,7 +83,7 @@ class NonConvexMVU {
   // These parameters are used for limited BFGS
  
   //ro_k = 1/(y^T * s)
-  ArrayList<Matrix> ro_bfgs_; 
+  Vector ro_bfgs_; 
   // the memory of bfgs 
   index_t mem_bfgs_;
   // s_k = x_{k+1}-x_{k};
diff --git a/fastlib2/contrib/nvasil/non_convex_mvu/non_convex_mvu_impl.h b/fastlib2/contrib/nvasil/non_convex_mvu/non_convex_mvu_impl.h
index ebc5b80a8c..cd20d58f01 100644
--- a/fastlib2/contrib/nvasil/non_convex_mvu/non_convex_mvu_impl.h
+++ b/fastlib2/contrib/nvasil/non_convex_mvu/non_convex_mvu_impl.h
@@ -87,7 +87,7 @@ void NonConvexMVU::ComputeLocalOptimum() {
  
 }
 
-void NonConvexMVU::ComputeLocalOptimumBFGS_() {
+void NonConvexMVU::ComputeLocalOptimumBFGS() {
   double distance_constraint;
   double centering_constraint;
   double sum_of_dist_square = la::LengthEuclidean(distances_.size(), &distances_[0]);
@@ -107,7 +107,6 @@ void NonConvexMVU::ComputeLocalOptimumBFGS_() {
   for(index_t i=0; i<mem_bfgs_; i++) {
     s_bfgs_[i].Init(new_dimension_, num_of_points_);
     y_bfgs_[i].Init(new_dimension_, num_of_points_);
-    ro_bfgs_[i].Init(new_dimension_, 1);
   } 
   NOTIFY("Starting optimization ...\n");
   ComputeFeasibilityError_(&distance_constraint, &centering_constraint);
@@ -122,7 +121,8 @@ void NonConvexMVU::ComputeLocalOptimumBFGS_() {
     ComputeGradient_();
     la::SubOverwrite(coordinates_, previous_coordinates_, &s_bfgs_[i]);
     la::SubOverwrite(gradient_, previous_gradient_, &y_bfgs_[i]);
-    la::MulTransBOverwrite(s_bfgs_[i], y_bfgs_[i], &ro_bfgs_[i]);
+    ro_bfgs_[i] = la::Dot(num_of_points_ * new_dimension_, s_bfgs_[i].ptr(),
+                          y_bfgs_[i].ptr());
     previous_gradient_.CopyValues(gradient_);
     previous_coordinates_.CopyValues(coordinates_);
   } 
@@ -130,6 +130,8 @@ void NonConvexMVU::ComputeLocalOptimumBFGS_() {
   for(index_t it1=0; it1<max_iterations_; it1++) {  
     for(index_t it2=0; it2<max_iterations_; it2++) {
       ComputeBFGS_();
+      ComputeGradient_();
+      la::SubFrom(gradient_, &coordinates_);
       UpdateBFGS_();
       ComputeFeasibilityError_(&distance_constraint, &centering_constraint);
       NOTIFY("Iteration: %"LI"d : %"LI"d, feasibility error (dist)): %lg\n"
@@ -192,6 +194,10 @@ void NonConvexMVU::set_armijo_beta(double armijo_beta) {
   armijo_beta_ = armijo_beta;
 }
 
+void NonConvexMVU::set_mem_bfgs(index_t mem_bfgs) {
+  mem_bfgs_ = mem_bfgs;
+}
+
 void NonConvexMVU::InitOptimization_() {
   if (unlikely(new_dimension_<0)) {
     FATAL("You forgot to set the new dimension\n");
@@ -304,42 +310,50 @@ void NonConvexMVU::LocalSearch_(double *step) {
 }
 
 void NonConvexMVU::ComputeBFGS_() {
-  ArrayList<Matrix> alpha;
+  Vector alpha;
   alpha.Init(mem_bfgs_);
-  for(index_t i=0; i<mem_bfgs_; i++){
-    alpha[i].Init(new_dimension_, 1); 
-  }
   Matrix scaled_y;
-  scaled_y.Init(new_dimension_, 1);
+  scaled_y.Init(new_dimension_, num_of_points_);
   index_t num=0;
-  for(index_t i=index_bfgs_, num=0; num<mem_bfgs_; i=i%(mem_bfgs_+1), num++) {
-    la::MulTransBOverwrite(s_bfgs_[i], gradient_, &alpha[i]);
-    la::ScaleRows(ro_bfgs_[i], &alpha[i]);
+  for(index_t i=index_bfgs_, num=0; num<mem_bfgs_; i=(i+1)%mem_bfgs_, num++) {
+    alpha[i] = la::Dot(new_dimension_ * num_of_points_,
+                       s_bfgs_[i].ptr(), 
+                       gradient_.ptr());
+    alpha[i] *= ro_bfgs_[i];
     scaled_y.CopyValues(y_bfgs_[i]);
-    la::ScaleRows(alpha[i], &scaled_y);
+    la::Scale(alpha[i], &scaled_y);
     la::SubFrom(scaled_y, &gradient_);
   }
+  // We need to scale the gradient here
+  double norm_scale=1/ro_bfgs_[index_bfgs_]/la::Dot(num_of_points_ *
+                    new_dimension_, y_bfgs_[index_bfgs_].ptr(),
+                    y_bfgs_[index_bfgs_].ptr());
+  la::Scale(norm_scale, &gradient_);
+  
   Matrix scaled_s;
-  Matrix beta;
-  beta.Init(new_dimension_, 1);
+  double beta;
   scaled_s.Init(new_dimension_, num_of_points_);
   num=0;
-  for(index_t j=index_bfgs_, num=0; num<mem_bfgs_; 
-      num++, j=(j-1)%mem_bfgs_) {
-    la::MulTransBOverwrite(y_bfgs_[j], gradient_, &beta);
-    la::ScaleRows(ro_bfgs_[j], &beta);
-    la::SubFrom(alpha[j], &beta);
+  for(index_t j=(index_bfgs_-1+mem_bfgs_)%mem_bfgs_, num=0; num<mem_bfgs_; 
+      num++, j=(j-1+mem_bfgs_)%mem_bfgs_) {
+    beta = la::Dot(new_dimension_ * num_of_points_, 
+                   y_bfgs_[j].ptr(),
+                   gradient_.ptr());
+    beta *= ro_bfgs_[j];
     scaled_s.CopyValues(s_bfgs_[j]);
-    la::ScaleRows(beta, &scaled_s);
+    la::Scale(alpha[j]-beta, &scaled_s);
     la::AddTo(scaled_s, &gradient_);
   }
 }
 
 void NonConvexMVU::UpdateBFGS_() {
   // shift all values
-  index_bfgs_ = (index_bfgs_ - 1) % mem_bfgs_;
+  index_bfgs_ = (index_bfgs_ - 1 + mem_bfgs_ ) % mem_bfgs_;
   la::SubOverwrite(coordinates_, previous_coordinates_, &s_bfgs_[index_bfgs_]);
   la::SubOverwrite(gradient_, previous_gradient_, &y_bfgs_[index_bfgs_]);
+  ro_bfgs_[index_bfgs_] = 1.0/la::Dot(new_dimension_ * num_of_points_, 
+                                      s_bfgs_[index_bfgs_].ptr(),
+                                      y_bfgs_[index_bfgs_].ptr());
 }
 
 double NonConvexMVU::ComputeLagrangian_(Matrix &coord) {
diff --git a/fastlib2/contrib/nvasil/non_convex_mvu/non_convex_mvu_test.cc b/fastlib2/contrib/nvasil/non_convex_mvu/non_convex_mvu_test.cc
index dfe7c59a0f..44f7539682 100644
--- a/fastlib2/contrib/nvasil/non_convex_mvu/non_convex_mvu_test.cc
+++ b/fastlib2/contrib/nvasil/non_convex_mvu/non_convex_mvu_test.cc
@@ -34,6 +34,11 @@ class NonConvexMVUTest {
     engine_->Init("test_data_3_1000.csv", 5);
     engine_->coordinates_.Init(1, 1);
     engine_->gradient_.Init(1, 1);
+    engine_->previous_gradient_.Init(1, 1);
+    engine_->previous_coordinates_.Init(1, 1);
+    engine_->ro_bfgs_.Init(1);
+    engine_->s_bfgs_.Init();
+    engine_->y_bfgs_.Init();
     engine_->lagrange_mult_.Init(30);
     engine_->centering_lagrange_mult_.Init(20);
     NOTIFY("TestInit passed!!\n");
@@ -46,15 +51,27 @@ class NonConvexMVUTest {
     engine_->ComputeLocalOptimum();
     NOTIFY("TestComputeLocalOptimum() passed!!\n");
   }
-  
+  void TestComputeLocalOptimumBFGS() {
+    NOTIFY("Testing ComputeLocalOptimum() ...\n");
+    engine_->Init("test_data_3_1000.csv", 5);
+    engine_->set_new_dimension(3);
+    engine_->set_mem_bfgs(5);
+    engine_->ComputeLocalOptimumBFGS();
+    NOTIFY("TestComputeLocalOptimum() passed!!\n");
+  }
+ 
   void TestAll() {
     Init();
     TestInit();
     Destruct();
+   // Init();
+   // TestComputeLocalOptimum();
+   //  Destruct();
     Init();
-    TestComputeLocalOptimum();
+    TestComputeLocalOptimumBFGS();
     Destruct();
   }
+  
  private:
   NonConvexMVU *engine_;  
 };