diff --git a/src/feature_creation/feature_space/FeatureSpace.cpp b/src/feature_creation/feature_space/FeatureSpace.cpp
index 938abef653bd56874b85b21f7326206d439cdd6a..6a4900af52fad2fef2af0b82904fb1f3c09a30e4 100644
--- a/src/feature_creation/feature_space/FeatureSpace.cpp
+++ b/src/feature_creation/feature_space/FeatureSpace.cpp
@@ -923,7 +923,9 @@ void FeatureSpace::sis(const std::vector<double>& prop)
         // Collect the best features from all ranks
         std::vector<node_sc_pair> local_sel(_n_sis_select, std::make_tuple<node_ptr, double>(nullptr, std::numeric_limits<double>::max()));
         for(int ff = 0; ff < phi_sel.size(); ++ff)
-            local_sel.push_back(mpi_reduce_op::make_node_sc_pair(phi_sel[ff], scores_sel[ff]));
+        {
+            local_sel[ff] = mpi_reduce_op::make_node_sc_pair(phi_sel[ff], scores_sel[ff]);
+        }
 
         std::vector<node_sc_pair> selected(_n_sis_select);
         mpi::all_reduce(
diff --git a/src/mpi_interface/MPI_ops.cpp b/src/mpi_interface/MPI_ops.cpp
index d04260084c45f8dc4f780ee53ba4ca6e77b1a908..ff2f11f178010aa5d25ff8dda2352e160685890b 100644
--- a/src/mpi_interface/MPI_ops.cpp
+++ b/src/mpi_interface/MPI_ops.cpp
@@ -28,12 +28,14 @@ void mpi_reduce_op::set_op(std::string project_type, double cross_cor_max, int n
 
 std::vector<node_sc_pair> mpi_reduce_op::select_top_feats(std::vector<node_sc_pair> in_vec_1, std::vector<node_sc_pair> in_vec_2)
 {
-    std::vector<node_sc_pair> out_vec;
-    out_vec.reserve(N_SIS_SELECT);
-    in_vec_2.insert(in_vec_2.end(), in_vec_1.begin(), in_vec_1.end());
+    // Set up an output vector
+    std::vector<node_sc_pair> out_vec(N_SIS_SELECT);
 
+    // Merge input vectors and sort
+    in_vec_2.insert(in_vec_2.end(), in_vec_1.begin(), in_vec_1.end());
     std::sort(in_vec_2.begin(), in_vec_2.end(), my_sorter);
 
+    // Populate the output vector
     int ff = 0;
     int out_ind = 0;
     while((out_ind < N_SIS_SELECT) && (ff < in_vec_2.size()))
@@ -41,7 +43,7 @@ std::vector<node_sc_pair> mpi_reduce_op::select_top_feats(std::vector<node_sc_pa
         const node_ptr cur_node = std::get<0>(in_vec_2[ff]);
         if(cur_node && IS_VALID(cur_node->value_ptr(), cur_node->n_samp(), CROSS_COR_MAX, out_vec, std::get<1>(in_vec_2[ff])))
         {
-            out_vec.push_back(in_vec_2[ff]);
+            out_vec[out_ind] = in_vec_2[ff];
             ++out_ind;
         }
         ++ff;