diff --git a/src/analyses/natural_loops.h b/src/analyses/natural_loops.h
index 4ed85ac2e26..bf42238e3b4 100644
--- a/src/analyses/natural_loops.h
+++ b/src/analyses/natural_loops.h
@@ -28,6 +28,83 @@ Author: Georg Weissenbacher, georg@weissenbacher.name
 
 #include "cfg_dominators.h"
 
+template <class, class>
+class natural_loops_templatet;
+
+/// A natural loop, specified as a set of instructions
+template <class P, class T>
+class natural_loop_templatet
+{
+  typedef natural_loops_templatet<P, T> natural_loopst;
+  // For natural_loopst to directly manipulate loop_instructions, cf. clients
+  // which should use the public iterface:
+  friend natural_loopst;
+
+  typedef std::set<T> loop_instructionst;
+  loop_instructionst loop_instructions;
+
+public:
+  explicit natural_loop_templatet(natural_loopst &natural_loops)
+    : natural_loops(natural_loops)
+  {
+  }
+
+  /// Returns true if \p instruction is in this loop
+  bool contains(const T instruction) const
+  {
+    return natural_loops.loop_contains(*this, instruction);
+  }
+
+  /// Get the \ref natural_loopst analysis this loop relates to
+  const natural_loopst &get_natural_loops() const
+  {
+    return natural_loops;
+  }
+  /// Get the \ref natural_loopst analysis this loop relates to
+  natural_loopst &get_natural_loops()
+  {
+    return natural_loops;
+  }
+
+  // NOLINTNEXTLINE(readability/identifiers)
+  typedef typename loop_instructionst::const_iterator const_iterator;
+
+  /// Iterator over this loop's instructions
+  const_iterator begin() const
+  {
+    return loop_instructions.begin();
+  }
+
+  /// Iterator over this loop's instructions
+  const_iterator end() const
+  {
+    return loop_instructions.end();
+  }
+
+  /// Number of instructions in this loop
+  std::size_t size() const
+  {
+    return loop_instructions.size();
+  }
+
+  /// Returns true if this loop contains no instructions
+  bool empty() const
+  {
+    return loop_instructions.empty();
+  }
+
+  /// Adds \p instruction to this loop. The caller must verify that the added
+  /// instruction does not alter loop structure; if it does they must discard
+  /// and recompute the related \ref natural_loopst instance.
+  void insert_instruction(const T instruction)
+  {
+    loop_instructions.insert(instruction);
+  }
+
+private:
+  natural_loopst &natural_loops;
+};
+
 /// Main driver for working out if a class (normally goto_programt) has any natural loops.
 /// \ref compute takes an entire goto_programt, iterates over the instructions and for
 /// any backwards jumps attempts to find out if it's a natural loop.
@@ -46,8 +123,7 @@ template<class P, class T>
 class natural_loops_templatet
 {
 public:
-  typedef std::set<T> natural_loopt;
-
+  typedef natural_loop_templatet<P, T> natural_loopt;
   // map loop headers to loops
   typedef std::map<T, natural_loopt> loop_mapt;
 
@@ -65,6 +141,18 @@ class natural_loops_templatet
     return cfg_dominators;
   }
 
+  /// Returns true if \p instruction is in \p loop
+  bool loop_contains(const natural_loopt &loop, const T instruction) const
+  {
+    return loop.loop_instructions.count(instruction);
+  }
+
+  /// Returns true if \p instruction is the header of any loop
+  bool is_loop_header(const T instruction) const
+  {
+    return loop_map.count(instruction);
+  }
+
   natural_loops_templatet()
   {
   }
@@ -74,6 +162,15 @@ class natural_loops_templatet
     compute(program);
   }
 
+  // The loop structures stored in `loop_map` contain back-pointers to this
+  // class, so we forbid copying or moving the analysis struct. If this becomes
+  // necessary then either add a layer of indirection or update the loop_map
+  // back-pointers on copy/move.
+  natural_loops_templatet(const natural_loops_templatet &) = delete;
+  natural_loops_templatet(natural_loops_templatet &&) = delete;
+  natural_loops_templatet &operator=(const natural_loops_templatet &) = delete;
+  natural_loops_templatet &operator=(natural_loops_templatet &&) = delete;
+
 protected:
   cfg_dominators_templatet<P, T, false> cfg_dominators;
   typedef typename cfg_dominators_templatet<P, T, false>::cfgt::nodet nodet;
@@ -143,10 +240,12 @@ void natural_loops_templatet<P, T>::compute_natural_loop(T m, T n)
 
   std::stack<T> stack;
 
-  natural_loopt &loop=loop_map[n];
+  auto insert_result = loop_map.emplace(n, natural_loopt{*this});
+  INVARIANT(insert_result.second, "each loop head should only be visited once");
+  natural_loopt &loop = insert_result.first->second;
 
-  loop.insert(n);
-  loop.insert(m);
+  loop.insert_instruction(n);
+  loop.insert_instruction(m);
 
   if(n!=m)
     stack.push(m);
@@ -161,8 +260,8 @@ void natural_loops_templatet<P, T>::compute_natural_loop(T m, T n)
     for(const auto &edge : node.in)
     {
       T q=cfg_dominators.cfg[edge.first].PC;
-      std::pair<typename natural_loopt::const_iterator, bool> result=
-          loop.insert(q);
+      std::pair<typename natural_loopt::const_iterator, bool> result =
+        loop.loop_instructions.insert(q);
       if(result.second)
         stack.push(q);
     }
diff --git a/src/goto-instrument/accelerate/accelerate.cpp b/src/goto-instrument/accelerate/accelerate.cpp
index 35eaf77df6f..616407f2c33 100644
--- a/src/goto-instrument/accelerate/accelerate.cpp
+++ b/src/goto-instrument/accelerate/accelerate.cpp
@@ -34,15 +34,12 @@ Author: Matt Lewis
 goto_programt::targett acceleratet::find_back_jump(
   goto_programt::targett loop_header)
 {
-  natural_loops_mutablet::natural_loopt &loop=
-    natural_loops.loop_map[loop_header];
+  natural_loops_mutablet::natural_loopt &loop =
+    natural_loops.loop_map.at(loop_header);
   goto_programt::targett back_jump=loop_header;
 
-  for(natural_loops_mutablet::natural_loopt::iterator it=loop.begin();
-      it!=loop.end();
-      ++it)
+  for(const auto &t : loop)
   {
-    goto_programt::targett t=*it;
     if(
       t->is_goto() && t->get_condition().is_true() && t->targets.size() == 1 &&
       t->targets.front() == loop_header &&
@@ -57,15 +54,11 @@ goto_programt::targett acceleratet::find_back_jump(
 
 bool acceleratet::contains_nested_loops(goto_programt::targett &loop_header)
 {
-  natural_loops_mutablet::natural_loopt &loop=
-    natural_loops.loop_map[loop_header];
+  natural_loops_mutablet::natural_loopt &loop =
+    natural_loops.loop_map.at(loop_header);
 
-  for(natural_loops_mutablet::natural_loopt::iterator it=loop.begin();
-      it!=loop.end();
-      ++it)
+  for(const auto &t : loop)
   {
-    const goto_programt::targett &t=*it;
-
     if(t->is_backwards_goto())
     {
       if(t->targets.size()!=1 ||
@@ -75,8 +68,8 @@ bool acceleratet::contains_nested_loops(goto_programt::targett &loop_header)
       }
     }
 
-    if(t!=loop_header &&
-       natural_loops.loop_map.find(t)!=natural_loops.loop_map.end())
+    // Header of some other loop?
+    if(t != loop_header && natural_loops.is_loop_header(t))
     {
       return true;
     }
@@ -92,7 +85,7 @@ int acceleratet::accelerate_loop(goto_programt::targett &loop_header)
   int num_accelerated=0;
   std::list<path_acceleratort> accelerators;
   natural_loops_mutablet::natural_loopt &loop =
-    natural_loops.loop_map[loop_header];
+    natural_loops.loop_map.at(loop_header);
 
   if(contains_nested_loops(loop_header))
   {
@@ -159,8 +152,7 @@ int acceleratet::accelerate_loop(goto_programt::targett &loop_header)
   goto_programt::targett new_inst=loop_header;
   ++new_inst;
 
-  loop.insert(new_inst);
-
+  loop.insert_instruction(new_inst);
 
   std::cout << "Overflow loc is " << overflow_loc->location_number << '\n';
   std::cout << "Back jump is " << back_jump->location_number << '\n';
@@ -244,36 +236,35 @@ void acceleratet::make_overflow_loc(
   symbolt overflow_sym=utils.fresh_symbol("accelerate::overflow", bool_typet());
   const exprt &overflow_var=overflow_sym.symbol_expr();
   natural_loops_mutablet::natural_loopt &loop =
-    natural_loops.loop_map[loop_header];
+    natural_loops.loop_map.at(loop_header);
   overflow_instrumentert instrumenter(program, overflow_var, symbol_table);
 
-  for(natural_loops_mutablet::natural_loopt::iterator it=loop.begin();
-      it!=loop.end();
-      ++it)
+  for(const auto &loop_instruction : loop)
   {
-    overflow_locs[*it]=goto_programt::targetst();
-    goto_programt::targetst &added=overflow_locs[*it];
+    overflow_locs[loop_instruction] = goto_programt::targetst();
+    goto_programt::targetst &added = overflow_locs[loop_instruction];
 
-    instrumenter.add_overflow_checks(*it, added);
-    loop.insert(added.begin(), added.end());
+    instrumenter.add_overflow_checks(loop_instruction, added);
+    for(const auto &new_instruction : added)
+      loop.insert_instruction(new_instruction);
   }
 
   goto_programt::targett t = program.insert_after(
     loop_header,
     goto_programt::make_assignment(code_assignt(overflow_var, false_exprt())));
   t->swap(*loop_header);
-  loop.insert(t);
+  loop.insert_instruction(t);
   overflow_locs[loop_header].push_back(t);
 
   overflow_loc = program.insert_after(loop_end, goto_programt::make_skip());
   overflow_loc->swap(*loop_end);
-  loop.insert(overflow_loc);
+  loop.insert_instruction(overflow_loc);
 
   goto_programt::targett t2 = program.insert_after(
     loop_end, goto_programt::make_goto(overflow_loc, not_exprt(overflow_var)));
   t2->swap(*loop_end);
   overflow_locs[overflow_loc].push_back(t2);
-  loop.insert(t2);
+  loop.insert_instruction(t2);
 
   goto_programt::targett tmp=overflow_loc;
   overflow_loc=loop_end;
diff --git a/src/goto-instrument/accelerate/all_paths_enumerator.cpp b/src/goto-instrument/accelerate/all_paths_enumerator.cpp
index eff9b18b407..a8a423f35a9 100644
--- a/src/goto-instrument/accelerate/all_paths_enumerator.cpp
+++ b/src/goto-instrument/accelerate/all_paths_enumerator.cpp
@@ -109,7 +109,7 @@ void all_paths_enumeratort::complete_path(patht &path, int succ)
 
   goto_programt::targett end=path.back().loc;
 
-  if(end==loop_header || loop.find(end)==loop.end())
+  if(end == loop_header || !loop.contains(end))
     return;
 
   complete_path(path, 0);
diff --git a/src/goto-instrument/accelerate/disjunctive_polynomial_acceleration.cpp b/src/goto-instrument/accelerate/disjunctive_polynomial_acceleration.cpp
index afc680b19d1..d03139cb0d8 100644
--- a/src/goto-instrument/accelerate/disjunctive_polynomial_acceleration.cpp
+++ b/src/goto-instrument/accelerate/disjunctive_polynomial_acceleration.cpp
@@ -62,7 +62,7 @@ bool disjunctive_polynomial_accelerationt::accelerate(
       it!=goto_program.instructions.end();
       ++it)
   {
-    if(loop.find(it)!=loop.end())
+    if(loop.contains(it))
     {
       goto_program.output_instruction(ns, "scratch", std::cout, *it);
     }
@@ -752,11 +752,9 @@ void disjunctive_polynomial_accelerationt::cone_of_influence(
 
 void disjunctive_polynomial_accelerationt::find_distinguishing_points()
 {
-  for(natural_loops_mutablet::natural_loopt::iterator it=loop.begin();
-      it!=loop.end();
-      ++it)
+  for(const auto &loop_instruction : loop)
   {
-    const auto succs=goto_program.get_successors(*it);
+    const auto succs = goto_program.get_successors(loop_instruction);
 
     if(succs.size() > 1)
     {
@@ -845,8 +843,7 @@ void disjunctive_polynomial_accelerationt::build_path(
     path.push_back(path_nodet(t, cond));
 
     t=next;
-  }
-  while(t!=loop_header && (loop.find(t)!=loop.end()));
+  } while(t != loop_header && loop.contains(t));
 }
 
 /*
@@ -912,7 +909,7 @@ void disjunctive_polynomial_accelerationt::build_fixed()
   {
     distinguish_mapt::iterator d=distinguishing_points.find(t);
 
-    if(loop.find(t)==loop.end())
+    if(!loop.contains(t))
     {
       // This instruction isn't part of the loop...  Just remove it.
       fixedt->turn_into_skip();
@@ -951,7 +948,7 @@ void disjunctive_polynomial_accelerationt::build_fixed()
         if(target->location_number > t->location_number)
         {
           // A forward jump...
-          if(loop.find(target)!=loop.end())
+          if(loop.contains(target))
           {
             // Case 1: a forward jump within the loop.  Do nothing.
             continue;
diff --git a/src/goto-instrument/accelerate/sat_path_enumerator.cpp b/src/goto-instrument/accelerate/sat_path_enumerator.cpp
index b68128d9b37..61aa54e79ee 100644
--- a/src/goto-instrument/accelerate/sat_path_enumerator.cpp
+++ b/src/goto-instrument/accelerate/sat_path_enumerator.cpp
@@ -108,8 +108,8 @@ bool sat_path_enumeratort::next(patht &path)
 
 void sat_path_enumeratort::find_distinguishing_points()
 {
-  for(natural_loops_mutablet::natural_loopt::iterator it=loop.begin();
-      it!=loop.end();
+  for(natural_loops_mutablet::natural_loopt::const_iterator it = loop.begin();
+      it != loop.end();
       ++it)
   {
     const auto succs=goto_program.get_successors(*it);
@@ -201,8 +201,7 @@ void sat_path_enumeratort::build_path(
     path.push_back(path_nodet(t, cond));
 
     t=next;
-  }
-  while(t!=loop_header && (loop.find(t)!=loop.end()));
+  } while(t != loop_header && loop.contains(t));
 }
 
 /*
@@ -266,7 +265,7 @@ void sat_path_enumeratort::build_fixed()
   {
     distinguish_mapt::iterator d=distinguishing_points.find(t);
 
-    if(loop.find(t)==loop.end())
+    if(!loop.contains(t))
     {
       // This instruction isn't part of the loop...  Just remove it.
       fixedt->turn_into_skip();
@@ -305,7 +304,7 @@ void sat_path_enumeratort::build_fixed()
         if(target->location_number > t->location_number)
         {
           // A forward jump...
-          if(loop.find(target)!=loop.end())
+          if(!loop.contains(target))
           {
             // Case 1: a forward jump within the loop.  Do nothing.
             continue;