diff --git a/cpp/mwmatching.h b/cpp/mwmatching.h
index 02995b5..cad8d19 100644
--- a/cpp/mwmatching.h
+++ b/cpp/mwmatching.h
@@ -10,6 +10,7 @@
 #include <cmath>
 #include <cstddef>
 #include <deque>
+#include <iterator>
 #include <limits>
 #include <list>
 #include <stack>
@@ -405,6 +406,9 @@ struct NonTrivialBlossom : public Blossom<WeightType>
         VertexPair edge;
     };
 
+    /** Iterator for this blossom in the list of non-trivial blossoms. */
+    typename std::list<NonTrivialBlossom>::iterator this_blossom_iterator;
+
     /**
      * List of sub-blossoms.
      * Ordered by their appearance in the alternating cycle.
@@ -585,10 +589,12 @@ public:
     std::vector<BlossomT> trivial_blossom;
 
     /**
-     * Non-trivial blossoms may be created and destroyed during
-     * the course of the algorithm.
+     * Linked list of non-trivial blossoms.
+     *
+     * These blossoms are created and destroyed during the course of the
+     * algorithm. Keep them in a linked list such that we can add and
+     * remove elements without invalidating pointers to other blossoms.
      */
-// NOTE - this MUST be a list, because we delete items from it while keeping pointers to other items
     std::list<NonTrivialBlossomT> nontrivial_blossom;
 
     /**
@@ -1255,8 +1261,12 @@ public:
         }
 
         // Create the new blossom object.
-        nontrivial_blossom.emplace_back(subblossoms, path.edges);
-        NonTrivialBlossomT* blossom = &nontrivial_blossom.back();
+        auto blossom_iterator = nontrivial_blossom.emplace(
+                                    nontrivial_blossom.end(),
+                                    subblossoms,
+                                    path.edges);
+        NonTrivialBlossomT* blossom = &(*blossom_iterator);
+        blossom->this_blossom_iterator = blossom_iterator;
 
         // Initialize modified blossom dual to set true dual to 0.
         blossom->dual_var = -2 * delta_sum;
@@ -1282,21 +1292,10 @@ public:
     /** Erase the specified non-trivial blossom. */
     void erase_nontrivial_blossom(NonTrivialBlossomT* blossom)
     {
-        auto blossom_it = std::find_if(
-            nontrivial_blossom.begin(),
-            nontrivial_blossom.end(),
-            [blossom](const NonTrivialBlossomT& b) {
-                return (&b == blossom);
-            });
-        assert(blossom_it != nontrivial_blossom.end());
-        nontrivial_blossom.erase(blossom_it);
+        nontrivial_blossom.erase(blossom->this_blossom_iterator);
     }
 
-    /**
-     * Expand the specified unlabeled blossom but do not yet delete it.
-     *
-     * This function takes time O(n).
-     */
+    /** Expand the specified unlabeled blossom but do not yet delete it. */
     void expand_unlabeled_blossom_core(NonTrivialBlossomT* blossom)
     {
         assert(blossom->parent == nullptr);
@@ -1329,22 +1328,14 @@ public:
         }
     }
 
-    /**
-     * Expand and delete the specified unlabeled blossom.
-     *
-     * This function takes time O(n).
-     */
+    /** Expand and delete the specified unlabeled blossom. */
     void expand_unlabeled_blossom(NonTrivialBlossomT* blossom)
     {
         expand_unlabeled_blossom_core(blossom);
         erase_nontrivial_blossom(blossom);
     }
 
-    /**
-     * Expand the specified T-blossom.
-     *
-     * This function takes time O(n).
-     */
+    /** Expand and delete the specified T-blossom. */
     void expand_t_blossom(NonTrivialBlossomT* blossom)
     {
         assert(blossom->parent == nullptr);