diff --git a/dwave/optimization/_model.pyi b/dwave/optimization/_model.pyi
index 53adc213..2afa1d23 100644
--- a/dwave/optimization/_model.pyi
+++ b/dwave/optimization/_model.pyi
@@ -71,6 +71,7 @@ class _Graph:
     def num_inputs(self) -> int: ...
     def num_nodes(self) -> int: ...
     def num_symbols(self) -> int: ...
+    def remove_redundant_symbols(self, time_limit_s: None | float = None) -> int: ...
     def remove_unused_symbols(self) -> int: ...
     def state_size(self) -> int: ...
     def unlock(self): ...
diff --git a/dwave/optimization/_model.pyx b/dwave/optimization/_model.pyx
index 2eb45de0..43f20e5a 100644
--- a/dwave/optimization/_model.pyx
+++ b/dwave/optimization/_model.pyx
@@ -928,6 +928,46 @@ cdef class _Graph:
         """
         return self.num_nodes()
 
+    def remove_redundant_symbols(self, time_limit_s=None):
+        """Remove redundant symbols from the model.
+
+        Symbols are redundant if they are the same type, they share the
+        same predecessors, and they encode the same operation.
+
+        Args:
+            time_limit_s (float): The maximum amount of time to spend trying
+                to remove redundant symbols. If not provided, the method will
+                run until there are no redundant symbols to remove.
+
+        Returns:
+            int: Number of symbols removed.
+
+        Examples:
+
+            >>> from dwave.optimization import Model
+            ...
+            >>> model = Model()
+            ...
+            >>> x = model.integer(5)
+            ...
+            >>> -(x + 1).sum()  # doctest: +ELLIPSIS
+            <dwave.optimization.symbols.unaryop.Negative at ...>
+            >>> -(x + 1).sum()  # doctest: +ELLIPSIS
+            <dwave.optimization.symbols.unaryop.Negative at ...>
+            >>> model.num_symbols()
+            8
+            >>> model.remove_redundant_symbols()
+            3
+            >>> model.num_symbols()
+            5
+        """
+        if self.is_locked():
+            raise ValueError("cannot remove symbols from a locked model")
+        if time_limit_s is None:
+            return self._graph.remove_redundant_nodes(False)
+        else:
+            return self._graph.remove_redundant_nodes(False, time_limit_s)
+
     def remove_unused_symbols(self):
         """Remove unused symbols from the model.
 
@@ -981,7 +1021,7 @@ cdef class _Graph:
         """
         if self.is_locked():
             raise ValueError("cannot remove symbols from a locked model")
-        return self._graph.remove_unused_nodes()
+        return self._graph.remove_unused_nodes(False)
 
     def state_size(self):
         r"""Return an estimate of the size, in bytes, of a model's state.
diff --git a/dwave/optimization/include/dwave-optimization/graph.hpp b/dwave/optimization/include/dwave-optimization/graph.hpp
index 22fbc3d6..e52b56df 100644
--- a/dwave/optimization/include/dwave-optimization/graph.hpp
+++ b/dwave/optimization/include/dwave-optimization/graph.hpp
@@ -156,6 +156,17 @@ class Graph {
     /// Reset the state of the given node and all successors recursively.
     static void recursive_reset(State& state, const Node* ptr);
 
+    /// Remove redundant nodes from the graph.
+    ///
+    /// Redundant nodes are ones that are Node::equal_to() another node in the
+    /// graph.
+    ///
+    /// Returns the number of nodes removed from the graph.
+    ssize_t remove_redundant_nodes(
+        bool ignore_listeners = false,
+        double time_limit_s = std::numeric_limits<double>::infinity()
+    );
+
     /// Remove unused nodes from the graph.
     ///
     /// This method will reset the topological sort if there is one.
@@ -263,6 +274,12 @@ class Node {
     /// derived from its predecessors. Defaults to `true`, except for decisions.
     virtual bool deterministic_state() const { return true; }
 
+    /// Test whether two nodes are equal. Each node class defines equality for
+    /// itself but nodes *must* share the same set of
+    /// predecessors (permutations are sometimes allowed) and they *must* be the
+    /// same type. Also that they can have false negatives in some cases.
+    virtual bool equal_to(const Node& rhs) const = 0;
+
     /// Return a shared pointer to a bool value. When the node is destructed
     /// the bool will be set to True
     std::shared_ptr<const bool> expired_ptr() const { return expired_ptr_; }
@@ -305,6 +322,11 @@ class Node {
     /// Return the successors of the node.
     const std::vector<SuccessorView>& successors() const { return successors_; }
 
+    /// Take the successors from `from` and make them successors of the node.
+    /// Note that this bypasses most checks! Use with extreme caution as it
+    /// make create undefined models.
+    void take_successors(Node& from);
+
     /// The current topological index. Will be negative if unsorted.
     ssize_t topological_index() const { return topological_index_; }
 
@@ -322,6 +344,7 @@ class Node {
     friend void Graph::reset_topological_sort();
     template <class NodeType, class... Args>
     friend NodeType* Graph::emplace_node(Args&&...);
+    friend ssize_t Graph::remove_redundant_nodes(bool, double);
     friend ssize_t Graph::remove_unused_nodes(bool);
 
  protected:
@@ -380,6 +403,10 @@ class Node {
     // Remove a successor. *Does not* remove itself from it's successor's predecessors.
     ssize_t remove_successor_(const Node* ptr);
 
+    // Replace the predecessor at `index` with the node specified by `node_ptr`.
+    // Does not update `node_ptr`.
+    virtual void replace_predecessor_(ssize_t index, Node* node_ptr);
+
  private:
     ssize_t topological_index_ = -1;  // negative is unset
 
@@ -431,6 +458,10 @@ class DecisionNode : public Decision, public virtual Node {
     /// Decision nodes by definition do not have a deterministic state.
     bool deterministic_state() const final { return false; }
 
+    /// Decision can only ever be equal to themselves because they are
+    /// independent variables.
+    bool equal_to(const Node& rhs) const final { return static_cast<const Node*>(this) == &rhs; }
+
     /// Decisions don't have predecessors so no one should be calling update().
     /// Always throws a logic_error.
     [[noreturn]] void update(State& state, int index) const override;
@@ -440,4 +471,44 @@ class DecisionNode : public Decision, public virtual Node {
     bool removable_() const override { return false; }
 };
 
+/// Provide an implementation of equal_to() that defers to a type-specific
+/// implementation.
+template<std::derived_from<Node> Base, typename Derived = void>
+struct EqualityMixin : Base {
+    /// @copydoc Node::equal_to()
+    bool equal_to(const Node& rhs) const final {
+        // Every node is equal with itself
+        if (&rhs == static_cast<const Node*>(this)) return true;
+
+        // A node cannot be equal if they are not of the same type.
+        const auto* rhs_ptr = dynamic_cast<const Derived*>(&rhs);
+        if (rhs_ptr == nullptr) return false;
+
+        // lhs and rhs are the same type, so let's go to our type-specific
+        // method.
+        return equal_to(*rhs_ptr);
+    }
+
+    // Nodes must implement `equal_to()` for their own type.
+    virtual bool equal_to(const Derived& rhs) const = 0;
+};
+
+/// Overload for EqualityMixin that does not require a class-specific equal_to()
+/// method overload. This simply checks that they have the same type and same
+/// predecessors.
+template<std::derived_from<Node> Base>
+struct EqualityMixin<Base, void> : Base {
+    /// @copydoc Node::equal_to()
+    bool equal_to(const Node& rhs) const final {
+        // Every node is equal with itself
+        if (&rhs == static_cast<const Node*>(this)) return true;
+
+        // Nodes that are not the same type cannot be equal
+        if (typeid(*this) != typeid(rhs)) return false;
+
+        // Nodes that are of the same type must have the same predecessors
+        return std::ranges::equal(this->predecessors(), rhs.predecessors());
+    }
+};
+
 }  // namespace dwave::optimization
diff --git a/dwave/optimization/include/dwave-optimization/nodes/binaryop.hpp b/dwave/optimization/include/dwave-optimization/nodes/binaryop.hpp
index fb0a62e4..3e778561 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/binaryop.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/binaryop.hpp
@@ -29,7 +29,7 @@
 namespace dwave::optimization {
 
 template <class BinaryOp>
-class BinaryOpNode : public ArrayOutputMixin<ArrayNode> {
+class BinaryOpNode : public ArrayOutputMixin<EqualityMixin<ArrayNode, BinaryOpNode<BinaryOp>>> {
  public:
     // We need at least two nodes, and they must be the same shape
     BinaryOpNode(ArrayNode* a_ptr, ArrayNode* b_ptr);
@@ -37,6 +37,11 @@ class BinaryOpNode : public ArrayOutputMixin<ArrayNode> {
     double const* buff(const State& state) const override;
     std::span<const Update> diff(const State& state) const override;
 
+    /// Two BinaryOpNodes are equal if they have the same operation and the
+    /// same predecessors. If the BinaryOp is commutative, then permutations
+    /// of predecessors are allowed.
+    bool equal_to(const BinaryOpNode& rhs) const override;
+
     /// @copydoc Array::integral()
     bool integral() const override;
 
@@ -46,10 +51,10 @@ class BinaryOpNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Array::max()
     double max() const override;
 
-    using ArrayOutputMixin::shape;
+    using ArrayOutputMixin<EqualityMixin<ArrayNode, BinaryOpNode<BinaryOp>>>::shape;
     std::span<const ssize_t> shape(const State& state) const override;
 
-    using ArrayOutputMixin::size;
+    using ArrayOutputMixin<EqualityMixin<ArrayNode, BinaryOpNode<BinaryOp>>>::size;
     ssize_t size(const State& state) const override;
 
     ssize_t size_diff(const State& state) const override;
@@ -63,20 +68,22 @@ class BinaryOpNode : public ArrayOutputMixin<ArrayNode> {
 
     // The predecessors of the operation, as Array*.
     std::span<Array* const> operands() {
-        assert(predecessors().size() == operands_.size());
+        assert(this->predecessors().size() == operands_.size());
         return operands_;
     }
     std::span<const Array* const> operands() const {
-        assert(predecessors().size() == operands_.size());
+        assert(this->predecessors().size() == operands_.size());
         return operands_;
     }
 
  private:
+    void replace_predecessor_(ssize_t previous_index, Node* node_ptr) override;
+
     BinaryOp op;
 
     // There are redundant, because we could dynamic_cast each time from
     // predecessors(), but this is more performant
-    std::array<Array* const, 2> operands_;
+    std::array<Array*, 2> operands_;
 
     const ValuesInfo values_info_;
     const SizeInfo sizeinfo_;
diff --git a/dwave/optimization/include/dwave-optimization/nodes/collections.hpp b/dwave/optimization/include/dwave-optimization/nodes/collections.hpp
index 48b20119..8dadd26e 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/collections.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/collections.hpp
@@ -138,6 +138,9 @@ class DisjointBitSetNode : public ArrayOutputMixin<ArrayNode> {
 
     std::span<const Update> diff(const State& state) const override;
 
+    // A DisjointBitSetNode can only ever be equal to itself.
+    bool equal_to(const Node& rhs) const override;
+
     /// @copydoc Array::integral()
     bool integral() const override;
 
@@ -228,6 +231,9 @@ class DisjointListNode : public ArrayOutputMixin<ArrayNode> {
 
     std::span<const Update> diff(const State& state) const override;
 
+    // A DisjointListNode can only ever be equal to itself.
+    bool equal_to(const Node& rhs) const override;
+
     /// @copydoc Array::integral()
     bool integral() const override;
 
diff --git a/dwave/optimization/include/dwave-optimization/nodes/constants.hpp b/dwave/optimization/include/dwave-optimization/nodes/constants.hpp
index 2fc0f194..57efe7a9 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/constants.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/constants.hpp
@@ -26,7 +26,7 @@
 namespace dwave::optimization {
 
 /// A contiguous block of numbers.
-class ConstantNode : public ArrayOutputMixin<ArrayNode> {
+class ConstantNode : public ArrayOutputMixin<EqualityMixin<ArrayNode, ConstantNode>> {
  public:
     struct DataSource {
         DataSource() = default;
@@ -100,6 +100,9 @@ class ConstantNode : public ArrayOutputMixin<ArrayNode> {
 
     // Overloads required by the Node ABC *************************************
 
+    /// @copydoc Node::equal_to()
+    bool equal_to(const ConstantNode& rhs) const override;
+
     // This node never needs to update its successors
     void propagate(State&) const noexcept override {}
 
diff --git a/dwave/optimization/include/dwave-optimization/nodes/creation.hpp b/dwave/optimization/include/dwave-optimization/nodes/creation.hpp
index 876c1a10..5e77e6ff 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/creation.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/creation.hpp
@@ -24,7 +24,7 @@
 namespace dwave::optimization {
 
 /// A node encoding evenly spaced values within a given interval.
-class ARangeNode : public ArrayOutputMixin<ArrayNode> {
+class ARangeNode : public ArrayOutputMixin<EqualityMixin<ArrayNode, ARangeNode>> {
  public:
     using array_or_int = std::variant<const Array*, ssize_t>;
 
@@ -83,6 +83,9 @@ class ARangeNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Array::diff()
     std::span<const Update> diff(const State& state) const override;
 
+    /// ARangeNode is equal to another ARangeNode encoding the same range.
+    bool equal_to(const ARangeNode& rhs) const override;
+
     /// @copydoc Node::initialize_state()
     void initialize_state(State& state) const override;
 
@@ -126,10 +129,13 @@ class ARangeNode : public ArrayOutputMixin<ArrayNode> {
     /// The value or array defining the spacing between values.
     array_or_int step() const { return step_; }
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
-    const array_or_int start_;
-    const array_or_int stop_;
-    const array_or_int step_;
+    array_or_int start_;
+    array_or_int stop_;
+    array_or_int step_;
 
     const std::pair<double, double> values_minmax_;
     const SizeInfo sizeinfo_;
diff --git a/dwave/optimization/include/dwave-optimization/nodes/flow.hpp b/dwave/optimization/include/dwave-optimization/nodes/flow.hpp
index 0487636a..a8e01cc9 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/flow.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/flow.hpp
@@ -26,7 +26,7 @@ namespace dwave::optimization {
 ///
 /// `condition` and `arr` must be the same size. This always outputs a
 /// 1d array.
-class ExtractNode : public ArrayOutputMixin<ArrayNode> {
+class ExtractNode : public ArrayOutputMixin<EqualityMixin<ArrayNode>> {
  public:
     ExtractNode(ArrayNode* condition_ptr, ArrayNode* arr_ptr);
 
@@ -73,6 +73,9 @@ class ExtractNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Array::sizeinfo()
     SizeInfo sizeinfo() const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     // these are redundant, but convenient
     const Array* condition_ptr_;
@@ -87,13 +90,14 @@ class ExtractNode : public ArrayOutputMixin<ArrayNode> {
 /// `condition` must be either a scalar array or the same shape as `x` and `y`.
 /// `x` and `y` must have the same shape, including dynamic.
 /// dynamically sized `condition`s are not allowed.
-class WhereNode : public ArrayOutputMixin<ArrayNode> {
+class WhereNode : public ArrayOutputMixin<EqualityMixin<ArrayNode>> {
  public:
     WhereNode(ArrayNode* condition_ptr, ArrayNode* x_ptr, ArrayNode* y_ptr);
 
     double const* buff(const State& state) const override;
     void commit(State& state) const override;
     std::span<const Update> diff(const State& state) const override;
+
     void initialize_state(State& state) const override;
 
     /// @copydoc Array::integral()
@@ -116,6 +120,9 @@ class WhereNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Array::sizeinfo()
     SizeInfo sizeinfo() const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     // these are redundant, but convenient
     const Array* condition_ptr_;
diff --git a/dwave/optimization/include/dwave-optimization/nodes/indexing.hpp b/dwave/optimization/include/dwave-optimization/nodes/indexing.hpp
index c3968d0d..e61cce7e 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/indexing.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/indexing.hpp
@@ -41,7 +41,7 @@ namespace dwave::optimization {
 //      https://numpy.org/doc/stable/user/basics.indexing.html#combining-advanced-and-basic-indexing
 //      A combination of int/slice/array indices. The rules are... more complicated.
 
-class AdvancedIndexingNode : public ArrayNode {
+class AdvancedIndexingNode : public EqualityMixin<ArrayNode, AdvancedIndexingNode> {
  public:
     // Indices are some combination of arrays and slices
     using array_or_slice = std::variant<ArrayNode*, Slice>;
@@ -80,6 +80,7 @@ class AdvancedIndexingNode : public ArrayNode {
 
     // Node overloads
     void commit(State& state) const override;
+    bool equal_to(const AdvancedIndexingNode& rhs) const override;
     void initialize_state(State& state) const override;
     void propagate(State& state) const override;
     void revert(State& state) const override;
@@ -87,6 +88,9 @@ class AdvancedIndexingNode : public ArrayNode {
     // AdvancedIndexingHelper methods
     std::span<const array_or_slice> indices() const { return indices_; }
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     struct IndexParser_;
     AdvancedIndexingNode(ArrayNode* array_ptr, IndexParser_&& parser);
@@ -156,7 +160,7 @@ class AdvancedIndexingNode : public ArrayNode {
 
     const ssize_t size_;
 
-    const std::vector<array_or_slice> indices_;
+    std::vector<array_or_slice> indices_;
     const ssize_t indexing_arrays_ndim_;
 
     // "Grouped-indexers mode" refers to NumPy combined indexing where all array indexers
@@ -179,7 +183,7 @@ class AdvancedIndexingNode : public ArrayNode {
 // See https://numpy.org/doc/stable/user/basics.indexing.html#basic-indexing
 // BasicIndex nodes always have exactly one predecessor representing the array.
 // to be indexed.
-class BasicIndexingNode : public ArrayNode {
+class BasicIndexingNode : public EqualityMixin<ArrayNode, BasicIndexingNode> {
  public:
     // Indices are some combination of slices and integers.
     using slice_or_int = std::variant<Slice, ssize_t>;
@@ -229,6 +233,8 @@ class BasicIndexingNode : public ArrayNode {
 
     // don't need to overload update()
 
+    bool equal_to(const BasicIndexingNode& rhs) const override;
+
     void propagate(State& state) const override;
 
     void initialize_state(State& state) const override;
@@ -242,6 +248,9 @@ class BasicIndexingNode : public ArrayNode {
     // Infer the indices used to create the node.
     std::vector<slice_or_int> infer_indices() const;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     // Private constructor using an intermediate object
     struct IndexParser_;
@@ -293,7 +302,7 @@ class BasicIndexingNode : public ArrayNode {
     const SizeInfo sizeinfo_;
 };
 
-class PermutationNode : public ArrayOutputMixin<ArrayNode> {
+class PermutationNode : public ArrayOutputMixin<EqualityMixin<ArrayNode, PermutationNode>> {
  public:
     // We use this style rather than a template to support Cython later
     PermutationNode(ArrayNode* array_ptr, ArrayNode* order_ptr);
@@ -301,6 +310,9 @@ class PermutationNode : public ArrayOutputMixin<ArrayNode> {
     double const* buff(const State& state) const override;
     std::span<const Update> diff(const State& state) const override;
 
+    /// @copydoc Node::equal_to()
+    bool equal_to(const PermutationNode& rhs) const override;
+
     /// @copydoc Array::integral()
     bool integral() const override;
 
@@ -316,6 +328,9 @@ class PermutationNode : public ArrayOutputMixin<ArrayNode> {
 
     void propagate(State& state) const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     // we could dynamically cast each time, but it's easier to just keep separate
     // pointers to the "array" part of the two predecessors
diff --git a/dwave/optimization/include/dwave-optimization/nodes/inputs.hpp b/dwave/optimization/include/dwave-optimization/nodes/inputs.hpp
index 6b13f35d..578799f6 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/inputs.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/inputs.hpp
@@ -76,6 +76,9 @@ class InputNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Array::diff()
     std::span<const Update> diff(const State& state) const noexcept override;
 
+    /// InputNodes are never equal to other nodes
+    bool equal_to(const Node& rhs) const override { return static_cast<const Node*>(this) == &rhs; }
+
     /// @copydoc Array::integral()
     bool integral() const override { return values_info_.integral; };
 
@@ -102,6 +105,9 @@ class InputNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Node::revert()
     void revert(State& state) const noexcept override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     void check_values(std::span<const double> new_values) const;
 
diff --git a/dwave/optimization/include/dwave-optimization/nodes/interpolation.hpp b/dwave/optimization/include/dwave-optimization/nodes/interpolation.hpp
index 928ce9e5..d77073be 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/interpolation.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/interpolation.hpp
@@ -24,7 +24,7 @@
 #include "dwave-optimization/state.hpp"
 
 namespace dwave::optimization {
-class BSplineNode : public ArrayOutputMixin<ArrayNode> {
+class BSplineNode : public ArrayOutputMixin<EqualityMixin<ArrayNode, BSplineNode>> {
  public:
     explicit BSplineNode(
         ArrayNode* array_ptr,
@@ -36,6 +36,10 @@ class BSplineNode : public ArrayOutputMixin<ArrayNode> {
     double const* buff(const State& state) const override;
     void commit(State& state) const override;
     std::span<const Update> diff(const State&) const override;
+
+    /// @copydoc Node::equal_to()
+    bool equal_to(const BSplineNode& rhs) const override;
+
     void initialize_state(State& state) const override;
 
     /// @copydoc Array::integral()
@@ -58,6 +62,9 @@ class BSplineNode : public ArrayOutputMixin<ArrayNode> {
     using Array::size;
     ssize_t size(const State& state) const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     const Array* array_ptr_;
 
diff --git a/dwave/optimization/include/dwave-optimization/nodes/lambda.hpp b/dwave/optimization/include/dwave-optimization/nodes/lambda.hpp
index 4c85618e..bb1b4dc5 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/lambda.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/lambda.hpp
@@ -36,7 +36,7 @@ namespace dwave::optimization {
 // `std::accumulate()`, the special input should be the first `InputNode`
 // on the given `Graph`, with the remaining inputs used for the values of
 // the operands.
-class AccumulateZipNode : public ArrayOutputMixin<ArrayNode> {
+class AccumulateZipNode : public ArrayOutputMixin<EqualityMixin<ArrayNode, AccumulateZipNode>> {
  public:
     // Initial value can either be a double or another node
     using array_or_double = std::variant<ArrayNode*, double>;
@@ -97,10 +97,15 @@ class AccumulateZipNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Array::diff()
     std::span<const Update> diff(const State& state) const override;
 
+    /// @copydoc Node::equal_to()
+    bool equal_to(const AccumulateZipNode& rhs) const override;
+
     /// Access the underlying shared_ptr holding the Graph.
     /// Modifying the Graph leads to undefined behavior.
     std::shared_ptr<Graph>& expression_ptr() { return expression_ptr_; }
 
+    const array_or_double& initial() const { return initial_; }
+
     /// @copydoc Node::initialize_state()
     void initialize_state(State& state) const override;
 
@@ -125,7 +130,8 @@ class AccumulateZipNode : public ArrayOutputMixin<ArrayNode> {
     ssize_t size_diff(const State& state) const override;
     SizeInfo sizeinfo() const override;
 
-    const array_or_double initial;
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
 
  private:
     double evaluate_expression(State& register_) const;
@@ -136,7 +142,8 @@ class AccumulateZipNode : public ArrayOutputMixin<ArrayNode> {
     const InputNode* const accumulate_input() const;
 
     std::shared_ptr<Graph> expression_ptr_;
-    const std::vector<ArrayNode*> operands_;
+    array_or_double initial_;
+    std::vector<ArrayNode*> operands_;
 
     const SizeInfo sizeinfo_;
 };
diff --git a/dwave/optimization/include/dwave-optimization/nodes/linear_algebra.hpp b/dwave/optimization/include/dwave-optimization/nodes/linear_algebra.hpp
index 1f6a0c38..ec6260f6 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/linear_algebra.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/linear_algebra.hpp
@@ -14,15 +14,15 @@
 
 #pragma once
 
+#include <array>
 #include <string>
-#include <vector>
 
 #include "dwave-optimization/array.hpp"
 #include "dwave-optimization/graph.hpp"
 
 namespace dwave::optimization {
 
-class MatrixMultiplyNode : public ArrayOutputMixin<ArrayNode> {
+class MatrixMultiplyNode : public ArrayOutputMixin<EqualityMixin<ArrayNode>> {
  public:
     MatrixMultiplyNode(ArrayNode* x_ptr, ArrayNode* y_ptr);
 
@@ -47,6 +47,9 @@ class MatrixMultiplyNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Array::min()
     double min() const override;
 
+    /// The predecessors, as ArrayNode*
+    std::span<const ArrayNode* const, 2> operands() const { return operands_; }
+
     /// @copydoc Node::propagate()
     void propagate(State& state) const override;
 
@@ -71,12 +74,14 @@ class MatrixMultiplyNode : public ArrayOutputMixin<ArrayNode> {
     /// Either `"blas"` or `"fallback"`.
     static std::string implementation;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     void matmul_(State& state, std::span<double> out, std::span<const ssize_t> out_shape) const;
     void update_shape_(State& state) const;
 
-    const ArrayNode* x_ptr_;
-    const ArrayNode* y_ptr_;
+    std::array<const ArrayNode*, 2> operands_;
 
     const SizeInfo sizeinfo_;
     const ValuesInfo values_info_;
diff --git a/dwave/optimization/include/dwave-optimization/nodes/lp.hpp b/dwave/optimization/include/dwave-optimization/nodes/lp.hpp
index 9143140b..8c67313f 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/lp.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/lp.hpp
@@ -25,7 +25,7 @@ namespace dwave::optimization {
 class LinearProgramNodeBase;
 
 /// A logical node that propagates whether or not its predecessor LinearProgram is feasible.
-class LinearProgramFeasibleNode : public ScalarOutputMixin<ArrayNode, true> {
+class LinearProgramFeasibleNode : public ScalarOutputMixin<EqualityMixin<ArrayNode>, true> {
  public:
     explicit LinearProgramFeasibleNode(LinearProgramNodeBase* lp_ptr);
 
@@ -44,6 +44,9 @@ class LinearProgramFeasibleNode : public ScalarOutputMixin<ArrayNode, true> {
     /// @copydoc Node::propagate()
     void propagate(State& state) const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     const LinearProgramNodeBase* lp_ptr_;
 };
@@ -97,7 +100,7 @@ class LinearProgramNodeBase : public Node {
 /// Following https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html
 /// linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, bounds=(0, None),
 ///         callback=None, options=None, x0=None, integrality=None)
-class LinearProgramNode : public LinearProgramNodeBase {
+class LinearProgramNode : public EqualityMixin<LinearProgramNodeBase, LinearProgramNode> {
  public:
     /// Construct a LinearProgramNode
     ///
@@ -119,6 +122,9 @@ class LinearProgramNode : public LinearProgramNodeBase {
     /// The LP node's state is potentially degenerate, and therefore not deterministic.
     bool deterministic_state() const override;
 
+    /// @copydoc Node::equal_to()
+    bool equal_to(const LinearProgramNode& rhs) const override;
+
     /// @copydoc LinearProgramNodeBase::feasible()
     bool feasible(const State& state) const override;
 
@@ -149,6 +155,9 @@ class LinearProgramNode : public LinearProgramNodeBase {
     /// @copydoc LinearProgramNodeBase::variables_shape()
     std::span<const ssize_t> variables_shape() const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     /// Read out each of the predecessor arrays and copy the data to `lp`, where
     /// it can be easily passed in to linprog (which expects (contiguous) vectors).
@@ -157,16 +166,19 @@ class LinearProgramNode : public LinearProgramNodeBase {
 
     // The coefficients of the linear objective function to be minimized.
     // minimize c @ x
+    // Should never be nullptr
     const ArrayNode* c_ptr_;
 
     // The inequality constraint matrix and vector.
     // b_lb <= A @ x <= b_ub
+    // Either all are nullptr or A and at least one of b_lb/b_ub exist
     const ArrayNode* b_lb_ptr_;
     const ArrayNode* A_ptr_;
     const ArrayNode* b_ub_ptr_;
 
     // The equality constaint matrix and vector.
     // A_eq @ x == b_eq
+    // Either both are nullptr or neither
     const ArrayNode* A_eq_ptr_;
     const ArrayNode* b_eq_ptr_;
 
@@ -180,7 +192,7 @@ class LinearProgramNode : public LinearProgramNodeBase {
 
 /// A scalar node that propagates the objective value of the solution found by the
 /// LinearProgramNode. Note that the output is undefined if the solution is not feasible.
-class LinearProgramObjectiveValueNode : public ScalarOutputMixin<ArrayNode, true> {
+class LinearProgramObjectiveValueNode : public ScalarOutputMixin<EqualityMixin<ArrayNode>, true> {
  public:
     explicit LinearProgramObjectiveValueNode(LinearProgramNodeBase* lp_ptr);
 
@@ -199,13 +211,16 @@ class LinearProgramObjectiveValueNode : public ScalarOutputMixin<ArrayNode, true
     /// @copydoc Node::propagate()
     void propagate(State& state) const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     const LinearProgramNodeBase* lp_ptr_;
 };
 
 /// An array node that propagates the solution found by the LinearProgramNode. Note that the
 /// solution may not be feasible or optimial.
-class LinearProgramSolutionNode : public ArrayOutputMixin<ArrayNode> {
+class LinearProgramSolutionNode : public ArrayOutputMixin<EqualityMixin<ArrayNode>> {
  public:
     explicit LinearProgramSolutionNode(LinearProgramNodeBase* lp_ptr);
 
@@ -240,6 +255,9 @@ class LinearProgramSolutionNode : public ArrayOutputMixin<ArrayNode> {
 
     using ArrayOutputMixin::size;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     const LinearProgramNodeBase* lp_ptr_;
 };
diff --git a/dwave/optimization/include/dwave-optimization/nodes/manipulation.hpp b/dwave/optimization/include/dwave-optimization/nodes/manipulation.hpp
index 19ee2094..21cfc549 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/manipulation.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/manipulation.hpp
@@ -26,7 +26,7 @@
 
 namespace dwave::optimization {
 
-class BroadcastToNode : public ArrayNode {
+class BroadcastToNode : public EqualityMixin<ArrayNode, BroadcastToNode> {
  public:
     BroadcastToNode(ArrayNode* array_ptr, std::initializer_list<ssize_t> shape);
     BroadcastToNode(ArrayNode* array_ptr, std::span<const ssize_t> shape);
@@ -43,6 +43,9 @@ class BroadcastToNode : public ArrayNode {
     /// @copydoc Array::diff()
     std::span<const Update> diff(const State& state) const override;
 
+    /// @copydoc Node::equal_to()
+    bool equal_to(const BroadcastToNode& rhs) const override;
+
     /// @copydoc Node::initialize_state()
     void initialize_state(State& state) const override;
 
@@ -78,6 +81,9 @@ class BroadcastToNode : public ArrayNode {
     /// @copydoc Array::strides()
     std::span<const ssize_t> strides() const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     /// Translate a linear index of the predecessor into a linear index of the
     /// BroadcastToNode.
@@ -94,7 +100,7 @@ class BroadcastToNode : public ArrayNode {
     const ValuesInfo values_info_;
 };
 
-class ConcatenateNode : public ArrayOutputMixin<ArrayNode> {
+class ConcatenateNode : public ArrayOutputMixin<EqualityMixin<ArrayNode, ConcatenateNode>> {
  public:
     explicit ConcatenateNode(std::span<ArrayNode*> array_ptrs, ssize_t axis);
     explicit ConcatenateNode(std::ranges::contiguous_range auto&& array_ptrs, ssize_t axis) :
@@ -103,6 +109,10 @@ class ConcatenateNode : public ArrayOutputMixin<ArrayNode> {
     double const* buff(const State& statfe) const override;
     void commit(State& state) const override;
     std::span<const Update> diff(const State& state) const override;
+
+    /// @copydoc Node::equal_to()
+    bool equal_to(const ConcatenateNode& rhs) const override;
+
     void initialize_state(State& state) const override;
 
     /// @copydoc Array::integral()
@@ -119,6 +129,9 @@ class ConcatenateNode : public ArrayOutputMixin<ArrayNode> {
 
     ssize_t axis() const { return axis_; }
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     ssize_t axis_;
     std::vector<ArrayNode*> array_ptrs_;
@@ -128,7 +141,7 @@ class ConcatenateNode : public ArrayOutputMixin<ArrayNode> {
 };
 
 /// An array node that is a contiguous copy of its predecessor.
-class CopyNode : public ArrayOutputMixin<ArrayNode> {
+class CopyNode : public ArrayOutputMixin<EqualityMixin<ArrayNode>> {
  public:
     explicit CopyNode(ArrayNode* array_ptr);
 
@@ -172,6 +185,9 @@ class CopyNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Array::size_diff()
     ssize_t size_diff(const State& state) const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     const Array* array_ptr_;
 
@@ -190,7 +206,7 @@ class CopyNode : public ArrayOutputMixin<ArrayNode> {
 /// @endcode
 ///
 /// In the case of duplicate indices, the most recent update will be propagated.
-class PutNode : public ArrayOutputMixin<ArrayNode> {
+class PutNode : public ArrayOutputMixin<EqualityMixin<ArrayNode>> {
  public:
     /// Constructor for PutNode.
     ///
@@ -230,6 +246,9 @@ class PutNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Node::revert()
     void revert(State& state) const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     const Array* array_ptr_;
     const Array* indices_ptr_;
@@ -239,7 +258,7 @@ class PutNode : public ArrayOutputMixin<ArrayNode> {
 };
 
 /// Propagates the values of its predecessor, interpreted into a different shape.
-class ReshapeNode : public ArrayOutputMixin<ArrayNode> {
+class ReshapeNode : public ArrayOutputMixin<EqualityMixin<ArrayNode, ReshapeNode>> {
  public:
     /// Constructor for ReshapeNode.
     ///
@@ -264,6 +283,9 @@ class ReshapeNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Array::diff()
     std::span<const Update> diff(const State& state) const override;
 
+    /// @copydoc Node::equal_to()
+    bool equal_to(const ReshapeNode& rhs) const override;
+
     /// @copydoc Node::initialize_state()
     void initialize_state(State& state) const override;
 
@@ -298,6 +320,9 @@ class ReshapeNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Array::size_diff()
     ssize_t size_diff(const State& state) const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     // we could dynamically cast each time, but it's easier to just keep separate
     // pointer to the "array" part of the predecessor
@@ -309,7 +334,7 @@ class ReshapeNode : public ArrayOutputMixin<ArrayNode> {
 
 /// Reshape a node to a specific non-dynamic shape. Use fill_value for any missing
 /// values.
-class ResizeNode : public ArrayOutputMixin<ArrayNode> {
+class ResizeNode : public ArrayOutputMixin<EqualityMixin<ArrayNode, ResizeNode>> {
  public:
     /// Constructor for ResizeNode.
     ///
@@ -336,6 +361,8 @@ class ResizeNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Array::diff()
     std::span<const Update> diff(const State& state) const override;
 
+    bool equal_to(const ResizeNode& rhs) const override;
+
     /// The fill value.
     double fill_value() const { return fill_value_; }
 
@@ -357,6 +384,9 @@ class ResizeNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Node::revert()
     void revert(State& state) const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     const Array* array_ptr_;
 
@@ -365,7 +395,7 @@ class ResizeNode : public ArrayOutputMixin<ArrayNode> {
     const ValuesInfo values_info_;
 };
 
-class RollNode : public ArrayOutputMixin<ArrayNode> {
+class RollNode : public ArrayOutputMixin<EqualityMixin<ArrayNode, RollNode>> {
  public:
     /// Construct a RollNode.
     ///
@@ -391,6 +421,9 @@ class RollNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Array::diff()
     std::span<const Update> diff(const State& state) const override;
 
+    /// @copydoc Node::equal_to()
+    bool equal_to(const RollNode& rhs) const override;
+
     /// @copydoc Node::initialize_state()
     void initialize_state(State& state) const override;
 
@@ -428,6 +461,9 @@ class RollNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Array::size_diff()
     ssize_t size_diff(const State& state) const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     // Rotate the given array by shift in-place.
     static void rotate_(std::span<double> array, ssize_t shift);
@@ -457,7 +493,7 @@ class RollNode : public ArrayOutputMixin<ArrayNode> {
     const SizeInfo sizeinfo_;
 };
 
-class SizeNode : public ScalarOutputMixin<ArrayNode, true> {
+class SizeNode : public ScalarOutputMixin<EqualityMixin<ArrayNode>, true> {
  public:
     explicit SizeNode(ArrayNode* node_ptr);
 
@@ -474,6 +510,9 @@ class SizeNode : public ScalarOutputMixin<ArrayNode, true> {
 
     void propagate(State& state) const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     // we could dynamically cast each time, but it's easier to just keep separate
     // pointer to the "array" part of the predecessor
@@ -483,7 +522,7 @@ class SizeNode : public ScalarOutputMixin<ArrayNode, true> {
 };
 
 // Compute the transpose of predecessor
-class TransposeNode : public ArrayNode {
+class TransposeNode : public EqualityMixin<ArrayNode> {
  public:
     TransposeNode(ArrayNode* array_ptr);
 
@@ -537,6 +576,9 @@ class TransposeNode : public ArrayNode {
 
     void revert(State&) const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     const Array* array_ptr_;
 
diff --git a/dwave/optimization/include/dwave-optimization/nodes/naryop.hpp b/dwave/optimization/include/dwave-optimization/nodes/naryop.hpp
index 22ebded2..ff9d50f6 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/naryop.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/naryop.hpp
@@ -30,7 +30,7 @@
 namespace dwave::optimization {
 
 template <class BinaryOp>
-class NaryOpNode : public ArrayOutputMixin<ArrayNode> {
+class NaryOpNode : public ArrayOutputMixin<EqualityMixin<ArrayNode>> {
  public:
     // Need at least one node to start with to determine the shape
     explicit NaryOpNode(ArrayNode* node_ptr);
@@ -69,6 +69,9 @@ class NaryOpNode : public ArrayOutputMixin<ArrayNode> {
         return operands_;
     }
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     BinaryOp op;
 
diff --git a/dwave/optimization/include/dwave-optimization/nodes/quadratic_model.hpp b/dwave/optimization/include/dwave-optimization/nodes/quadratic_model.hpp
index 9a2f2d9f..d8d2dbb6 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/quadratic_model.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/quadratic_model.hpp
@@ -78,7 +78,7 @@ class QuadraticModel {
     friend class QuadraticModelNode;
 };
 
-class QuadraticModelNode : public ScalarOutputMixin<ArrayNode> {
+class QuadraticModelNode : public ScalarOutputMixin<EqualityMixin<ArrayNode>> {
  public:
     QuadraticModelNode(ArrayNode* state_node_ptr, QuadraticModel&& quadratic_model);
 
diff --git a/dwave/optimization/include/dwave-optimization/nodes/reduce.hpp b/dwave/optimization/include/dwave-optimization/nodes/reduce.hpp
index 9037edc4..6e379e92 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/reduce.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/reduce.hpp
@@ -29,7 +29,7 @@
 namespace dwave::optimization {
 
 template <class BinaryOp>
-class ReduceNode : public ArrayOutputMixin<ArrayNode> {
+class ReduceNode : public ArrayOutputMixin<EqualityMixin<ArrayNode, ReduceNode<BinaryOp>>> {
  public:
     ReduceNode(ArrayNode* array_ptr);
 
@@ -56,6 +56,9 @@ class ReduceNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Array::diff()
     std::span<const Update> diff(const State& state) const override;
 
+    /// @copydoc Node::equal_to()
+    bool equal_to(const ReduceNode& rhs) const override;
+
     /// @copydoc Node::initialize_state()
     void initialize_state(State& state) const override;
 
@@ -70,11 +73,11 @@ class ReduceNode : public ArrayOutputMixin<ArrayNode> {
 
     // The predecessor of the reduction, as an Array*.
     std::span<Array* const> operands() {
-        assert(predecessors().size() == 1);
+        assert(this->predecessors().size() == 1);
         return std::span<Array* const, 1>(&array_ptr_, 1);
     }
     std::span<const Array* const> operands() const {
-        assert(predecessors().size() == 1);
+        assert(this->predecessors().size() == 1);
         return std::span<const Array* const, 1>(&array_ptr_, 1);
     }
 
@@ -85,11 +88,11 @@ class ReduceNode : public ArrayOutputMixin<ArrayNode> {
     void revert(State& state) const override;
 
     /// @copydoc Array::shape()
-    using ArrayOutputMixin::shape;
+    using ArrayOutputMixin<EqualityMixin<ArrayNode, ReduceNode<BinaryOp>>>::shape;
     std::span<const ssize_t> shape(const State& state) const override;
 
     /// @copydoc Array::size()
-    using ArrayOutputMixin::size;
+    using ArrayOutputMixin<EqualityMixin<ArrayNode, ReduceNode<BinaryOp>>>::size;
     ssize_t size(const State& state) const override;
 
     /// @copydoc Array::sizeinfo()
@@ -102,6 +105,9 @@ class ReduceNode : public ArrayOutputMixin<ArrayNode> {
     /// Otherwise uses the first element in the reduction.
     const std::optional<double> initial;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     // Perform a reduction over the reduction space associated with the given
     // index. The return type is determined by the reduction_type, see
@@ -110,7 +116,7 @@ class ReduceNode : public ArrayOutputMixin<ArrayNode> {
 
     BinaryOp op;
 
-    Array* const array_ptr_;
+    Array* array_ptr_;
 
     // The axes we're reducing in a sorted, unique vector.
     // An empty vector means we're reducing over everything
diff --git a/dwave/optimization/include/dwave-optimization/nodes/set_routines.hpp b/dwave/optimization/include/dwave-optimization/nodes/set_routines.hpp
index 0ce7def6..862f9619 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/set_routines.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/set_routines.hpp
@@ -22,7 +22,7 @@
 
 namespace dwave::optimization {
 
-class IsInNode : public ArrayOutputMixin<ArrayNode> {
+class IsInNode : public ArrayOutputMixin<EqualityMixin<ArrayNode>> {
  public:
     IsInNode(ArrayNode* element_ptr, ArrayNode* test_elements_ptr);
 
@@ -69,6 +69,9 @@ class IsInNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Array::sizeinfo()
     SizeInfo sizeinfo() const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     // these are redundant, but convenient
     const Array* element_ptr_;
diff --git a/dwave/optimization/include/dwave-optimization/nodes/softmax.hpp b/dwave/optimization/include/dwave-optimization/nodes/softmax.hpp
index 26b66dd5..c2560d27 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/softmax.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/softmax.hpp
@@ -24,7 +24,7 @@
 
 namespace dwave::optimization {
 
-class SoftMaxNode : public ArrayOutputMixin<ArrayNode> {
+class SoftMaxNode : public ArrayOutputMixin<EqualityMixin<ArrayNode>> {
  public:
     SoftMaxNode(ArrayNode* arr_ptr);
 
@@ -71,6 +71,9 @@ class SoftMaxNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Array::sizeinfo()
     SizeInfo sizeinfo() const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     // these are redundant, but convenient
     const Array* arr_ptr_;
diff --git a/dwave/optimization/include/dwave-optimization/nodes/sorting.hpp b/dwave/optimization/include/dwave-optimization/nodes/sorting.hpp
index 55ea5abd..6effdae4 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/sorting.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/sorting.hpp
@@ -23,7 +23,7 @@
 
 namespace dwave::optimization {
 
-class ArgSortNode : public ArrayOutputMixin<ArrayNode> {
+class ArgSortNode : public ArrayOutputMixin<EqualityMixin<ArrayNode>> {
  public:
     ArgSortNode(ArrayNode* arr_ptr);
 
@@ -70,6 +70,9 @@ class ArgSortNode : public ArrayOutputMixin<ArrayNode> {
     /// @copydoc Array::sizeinfo()
     SizeInfo sizeinfo() const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     // these are redundant, but convenient
     const Array* arr_ptr_;
diff --git a/dwave/optimization/include/dwave-optimization/nodes/statistics.hpp b/dwave/optimization/include/dwave-optimization/nodes/statistics.hpp
index a6c3b407..0283c6e1 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/statistics.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/statistics.hpp
@@ -20,7 +20,7 @@
 
 namespace dwave::optimization {
 
-class MeanNode : public ScalarOutputMixin<ArrayNode, true> {
+class MeanNode : public ScalarOutputMixin<EqualityMixin<ArrayNode>, true> {
  public:
     MeanNode(ArrayNode* arr_ptr);
 
@@ -39,6 +39,9 @@ class MeanNode : public ScalarOutputMixin<ArrayNode, true> {
     /// @copydoc Node::propagate()
     void propagate(State& state) const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
  private:
     // these are redundant, but convenient
     const Array* arr_ptr_;
diff --git a/dwave/optimization/include/dwave-optimization/nodes/testing.hpp b/dwave/optimization/include/dwave-optimization/nodes/testing.hpp
index 53e08a61..51fb0234 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/testing.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/testing.hpp
@@ -23,7 +23,7 @@
 
 namespace dwave::optimization {
 
-class ArrayValidationNode : public Node {
+class ArrayValidationNode : public EqualityMixin<Node> {
  public:
     explicit ArrayValidationNode(ArrayNode* node_ptr);
 
@@ -33,6 +33,11 @@ class ArrayValidationNode : public Node {
     void propagate(State& state) const override;
     void revert(State& state) const override;
 
+ protected:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override {
+        assert(false and "ArrayValidationNode cannot have its predecessor replaced");
+    }
+
  private:
     const ArrayNode* array_ptr;
 
diff --git a/dwave/optimization/include/dwave-optimization/nodes/unaryop.hpp b/dwave/optimization/include/dwave-optimization/nodes/unaryop.hpp
index 6f633ae5..7e279533 100644
--- a/dwave/optimization/include/dwave-optimization/nodes/unaryop.hpp
+++ b/dwave/optimization/include/dwave-optimization/nodes/unaryop.hpp
@@ -35,6 +35,9 @@ class UnaryOpNode : public ArrayOutputMixin<ArrayNode> {
     double const* buff(const State& state) const override;
     std::span<const Update> diff(const State& state) const override;
 
+    bool equal_to(const Node& rhs) const override;
+    bool equal_to(const UnaryOpNode& rhs) const;
+
     /// @copydoc Array::integral()
     bool integral() const override;
 
@@ -57,21 +60,23 @@ class UnaryOpNode : public ArrayOutputMixin<ArrayNode> {
     void propagate(State& state) const override;
 
     // The predecessor of the operation, as an Array*.
-    std::span<Array* const> operands() {
+    std::span<ArrayNode* const> operands() {
         assert(predecessors().size() == 1);
-        return std::span<Array* const, 1>(&array_ptr_, 1);
+        return std::span<ArrayNode* const, 1>(&array_ptr_, 1);
     }
-    std::span<const Array* const> operands() const {
+    std::span<const ArrayNode* const> operands() const {
         assert(predecessors().size() == 1);
-        return std::span<const Array* const, 1>(&array_ptr_, 1);
+        return std::span<const ArrayNode* const, 1>(&array_ptr_, 1);
     }
 
  private:
+    void replace_predecessor_(ssize_t index, Node* node_ptr) override;
+
     UnaryOp op;
 
-    // There are redundant, because we could dynamic_cast each time from
+    // This is redundant, because we could dynamic_cast each time from
     // predecessors(), but this is more performant
-    Array* const array_ptr_;
+    ArrayNode* array_ptr_;
 
     const ValuesInfo values_info_;
     const SizeInfo sizeinfo_;
diff --git a/dwave/optimization/libcpp/graph.pxd b/dwave/optimization/libcpp/graph.pxd
index c4cade7f..4cf32453 100644
--- a/dwave/optimization/libcpp/graph.pxd
+++ b/dwave/optimization/libcpp/graph.pxd
@@ -64,9 +64,11 @@ cdef extern from "dwave-optimization/graph.hpp" namespace "dwave::optimization"
         void recursive_initialize(State&, Node*) except+
         @staticmethod
         void recursive_reset(State&, Node*) except+
+        Py_ssize_t remove_redundant_nodes(bool) except+
+        Py_ssize_t remove_redundant_nodes(bool, double) except+
+        Py_ssize_t remove_unused_nodes(bool) except+
         void reset_topological_sort()
         void set_objective(ArrayNode*) except+
         void add_constraint(ArrayNode*) except+
         void topological_sort()
         bool topologically_sorted() const
-        Py_ssize_t remove_unused_nodes()
diff --git a/dwave/optimization/libcpp/nodes/lambda_.pxd b/dwave/optimization/libcpp/nodes/lambda_.pxd
index 0cf2d53d..08d4ec08 100644
--- a/dwave/optimization/libcpp/nodes/lambda_.pxd
+++ b/dwave/optimization/libcpp/nodes/lambda_.pxd
@@ -22,4 +22,4 @@ cdef extern from "dwave-optimization/nodes/lambda.hpp" namespace "dwave::optimiz
     cdef cppclass AccumulateZipNode(ArrayNode):
         ctypedef variant["ArrayNode*", double] array_or_double
         shared_ptr[Graph] expression_ptr()
-        const array_or_double initial
+        const array_or_double& initial() const
diff --git a/dwave/optimization/src/graph.cpp b/dwave/optimization/src/graph.cpp
index 733e94ae..8b53fdad 100644
--- a/dwave/optimization/src/graph.cpp
+++ b/dwave/optimization/src/graph.cpp
@@ -15,6 +15,7 @@
 #include "dwave-optimization/graph.hpp"
 
 #include <algorithm>
+#include <chrono>
 #include <deque>
 #include <ranges>
 #include <sstream>
@@ -241,6 +242,175 @@ void Graph::recursive_reset(State& state, const Node* ptr) {
     }
 }
 
+ssize_t Graph::remove_redundant_nodes(bool ignore_listeners, double time_limit_s) {
+    if (topologically_sorted_) throw std::logic_error("cannot remove nodes from a locked model");
+
+    // This function can get quite expensive, so we have a timeout we'll try to respect.
+    const auto stop_time =
+        std::chrono::steady_clock::now() + std::chrono::duration<double>(time_limit_s);
+    auto out_of_time = [&stop_time]() -> bool {
+        return std::chrono::steady_clock::now() >= stop_time;
+    };
+
+    // We need to know how many nodes we started with in order to know how many we removed
+    const ssize_t num_nodes = this->num_nodes();
+
+    // We want the nodes in topological order, and we will respect that order
+    // while we're transferring successors.
+    // We will, however, abuse the topological_index_ to store information
+    topological_sort();
+
+    // We'll set the topological index to one of these values to track what we've
+    // already done.
+    // They are arbitrary values, but we steer away from -1 because that's used to
+    // indicate unsorted
+    constexpr ssize_t seen = -2;  // already checked for duplicates
+    constexpr ssize_t drop = -3;  // marked for deletion
+
+    // The actions that we always want to do before returning, whether it's because
+    // we hit the time limit or because we have no more work to do.
+    auto cleanup = [&]() -> ssize_t {
+        // we should only ever have swapped the objective_ptr
+        assert(objective_ptr_ == nullptr or objective_ptr_->topological_index() != drop);
+
+        // Whether a node was dropped
+        auto dropped = [&drop, &ignore_listeners](const auto& ptr) {
+            if (not ignore_listeners and ptr->num_listeners() > 0) return false;
+            assert(ptr->topological_index_ != drop or ptr->successors_.empty());
+            return ptr->topological_index() == drop;
+        };
+
+        // Go through our "special" nodes and drop anything
+        assert(std::ranges::none_of(decisions_, dropped));  // should never be dropped
+        assert(std::ranges::none_of(inputs_, dropped));  // should never be dropped
+        std::erase_if(constants_, dropped);
+
+        std::erase_if(constraints_, dropped);
+        assert(objective_ptr_ == nullptr or objective_ptr_->topological_index_ != drop);
+
+        // This is the step that actually deallocates the node
+        for (auto& uptr : nodes_) {
+            if (not dropped(uptr)) continue;
+
+            // Remove the node from its predecessor's successor vectors.
+            // This very leaves the node in an invalid state, very briefly
+            for (auto* pred_ptr : uptr->predecessors_) {
+                [[maybe_unused]] ssize_t num_removed = pred_ptr->remove_successor_(uptr.get());
+                assert(num_removed > 0);
+            }
+
+            // And then reset the pointer, thereby dellocating the node
+            uptr.reset();
+
+        }
+        // Finally, remove the nullptrs from the nodelist
+        std::erase_if(nodes_, [](const auto& uptr) { return not uptr; });
+
+        // Reset the topological index of everything that's left so we clean up
+        // everything
+        reset_topological_sort();
+
+        // Finally, report the number of nodes we dropped
+        return num_nodes - this->num_nodes();
+    };
+
+    // Given two equal nodes, transfer successors and the objective marker
+    // from `from_ptr` to `to_ptr`. This does not fix the constraints, we handle
+    // that as part of `cleanup()`.
+    auto transfer = [&](Node* from_ptr, Node* to_ptr) -> void {
+        assert(from_ptr->topological_index_ > to_ptr->topological_index_);
+
+        // Transfer the successors
+        to_ptr->take_successors(*from_ptr);
+
+        // Mark from for dropping
+        from_ptr->topological_index_ = drop;
+
+        // We also want to fix the objective if relevant
+        if (objective_ptr_ != nullptr and static_cast<Node*>(objective_ptr_) == from_ptr) {
+            objective_ptr_ = dynamic_cast<ArrayNode*>(to_ptr);
+            assert(objective_ptr_ != nullptr);
+            assert(objective_ptr_->size() == 1);
+        }
+    };
+
+    // Ok, all that setup done, now let's start checking for redundancy.
+
+    // Decisions are never redundant, so we skip over them
+
+    // Nor are inputs, so likewise we skip them
+
+    // The first set of nodes we're worried about are the constants - the only
+    // class of root node that can have redundancy
+    for (ssize_t i = 0, num_constants = constants_.size(); i < num_constants; ++i) {
+        if (constants_[i]->topological_index_ == drop) continue;
+
+        for (ssize_t j = i + 1; j < num_constants; ++j) {
+            if (constants_[j]->topological_index_ == drop) continue;
+
+            // the topological indices of our constants should not yet be touched
+            // and because they are added in order, they should always be ascending
+            assert(constants_[i]->topological_index_ >= 0);
+            assert(constants_[j]->topological_index_ >= 0);
+            assert(constants_[i]->topological_index_ < constants_[j]->topological_index_);
+
+            // Check against or current time limit before doing the potentially
+            // expensive (O(buffer_size)) equality check.
+            if (out_of_time()) return cleanup();
+
+            // nothing to do if they are not equal
+            if (not constants_[i]->equal_to(*constants_[j])) continue;
+
+            transfer(constants_[j], constants_[i]);
+        }
+
+        constants_[i]->topological_index_ = seen;
+    }
+
+    // For each node, we check pairwise among all of its successors.
+    // This is because for nodes to be equal, they *must* have the same set of
+    // predecessors.
+    for (auto& uptr : nodes_) {
+        auto& successors = uptr->successors();
+
+        for (ssize_t i = 0, num_successors = successors.size(); i < num_successors; ++i) {
+            // We've already seen or dropped this node
+            if (successors[i]->topological_index_ < 0) continue;
+
+            for (ssize_t j = i + 1; j < num_successors; ++j) {
+                // We've already seen or dropped this node
+                if (successors[j]->topological_index_ < 0) continue;
+
+                // A node cannot be redundant with itself
+                if (successors[i] == successors[j]) continue;
+
+                // Check against or current time limit before doing the potentially
+                // expensive (O(num_nodes^2)) equality check.
+                if (out_of_time()) return cleanup();
+
+                // If lhs != rhs there's nothing to do, so keep looking
+                if (not successors[i]->equal_to(*successors[j])) continue;
+
+                // We have a redundant node!
+
+                // We want to transfer to the node with the lower topological order
+                if (successors[i]->topological_index_ < successors[j]->topological_index_) {
+                    transfer(successors[j], successors[i]);
+                } else {
+                    transfer(successors[i], successors[j]);
+                    break;  // stop comparing i to other nodes because we dropped it
+                }
+            }
+
+            // Ok, we've checked i against everything, so if we didn't drop it
+            // we can mark it as seen
+            if (successors[i]->topological_index_ >= 0) successors[i]->topological_index_ = seen;
+        }
+    }
+
+    return cleanup();
+}
+
 ssize_t Graph::remove_unused_nodes(bool ignore_listeners) {
     if (topologically_sorted_) throw std::logic_error("cannot remove nodes from a locked model");
 
@@ -523,6 +693,22 @@ std::string Node::repr() const {
 
 std::string Node::str() const { return classname(); }
 
+void Node::take_successors(Node& from) {
+    assert(this != &from and "a node cannot take successors from itself");
+
+    for (const auto& sv : from.successors_) {
+        sv.ptr->replace_predecessor_(sv.index, this);
+        this->successors_.emplace_back(sv);
+    }
+    from.successors_.clear();
+}
+
+void Node::replace_predecessor_(ssize_t previous_index, Node* node_ptr) {
+    assert(0 <= previous_index);
+    assert(static_cast<size_t>(previous_index) < predecessors_.size());
+    predecessors_[previous_index] = node_ptr;
+}
+
 std::ostream& operator<<(std::ostream& os, const Node& node) { return os << node.repr(); }
 
 [[noreturn]] void DecisionNode::update(State& state, int index) const {
diff --git a/dwave/optimization/src/nodes/binaryop.cpp b/dwave/optimization/src/nodes/binaryop.cpp
index 2241cc92..d7029f6d 100644
--- a/dwave/optimization/src/nodes/binaryop.cpp
+++ b/dwave/optimization/src/nodes/binaryop.cpp
@@ -187,7 +187,7 @@ bool calculate_integral(const Array* lhs_ptr, const Array* rhs_ptr) {
 
 template <class BinaryOp>
 BinaryOpNode<BinaryOp>::BinaryOpNode(ArrayNode* a_ptr, ArrayNode* b_ptr) :
-    ArrayOutputMixin(broadcast_shapes(a_ptr->shape(), b_ptr->shape())),
+    ArrayOutputMixin<EqualityMixin<ArrayNode, BinaryOpNode<BinaryOp>>>(broadcast_shapes(a_ptr->shape(), b_ptr->shape())),
     operands_({a_ptr, b_ptr}),
     values_info_(
         calculate_values_minmax<BinaryOp>(operands_[0], operands_[1]),
@@ -200,17 +200,58 @@ BinaryOpNode<BinaryOp>::BinaryOpNode(ArrayNode* a_ptr, ArrayNode* b_ptr) :
 
 template <class BinaryOp>
 double const* BinaryOpNode<BinaryOp>::buff(const State& state) const {
-    return data_ptr_<ArrayNodeStateData>(state)->buff();
+    return this->template data_ptr_<ArrayNodeStateData>(state)->buff();
 }
 
 template <class BinaryOp>
 std::span<const Update> BinaryOpNode<BinaryOp>::diff(const State& state) const {
-    return data_ptr_<ArrayNodeStateData>(state)->diff();
+    return this->template data_ptr_<ArrayNodeStateData>(state)->diff();
 }
 
 template <class BinaryOp>
 void BinaryOpNode<BinaryOp>::commit(State& state) const {
-    data_ptr_<ArrayNodeStateData>(state)->commit();
+    this->template data_ptr_<ArrayNodeStateData>(state)->commit();
+}
+
+template <class BinaryOp>
+bool BinaryOpNode<BinaryOp>::equal_to(const BinaryOpNode& rhs) const {
+    // if we're the same type, then we just need to make sure we're operating
+    // on the same operands (subject to whether we're commutative or not.
+    // Once we switch to ufuncs, this gets even easier
+    if (std::ranges::equal(operands_, rhs.operands_)) return true;
+
+    // Once we switch to ufuncs
+    // https://github.com/dwavesystems/dwave-optimization/pull/412
+    // we can use ::communcative. For now we hardcode it.
+
+    // Dev note: This implementation doesn't check the pathalogical case op(x, x).
+    // While it would be mathematically correct, IMO it would be more confusing
+    // than helpful.
+
+    if constexpr (
+        std::same_as<BinaryOpNode, AddNode> or       //
+        std::same_as<BinaryOpNode, AndNode> or       //
+        std::same_as<BinaryOpNode, EqualNode> or     //
+        std::same_as<BinaryOpNode, MaximumNode> or   //
+        std::same_as<BinaryOpNode, MinimumNode> or   //
+        std::same_as<BinaryOpNode, MultiplyNode> or  //
+        std::same_as<BinaryOpNode, OrNode> or        //
+        std::same_as<BinaryOpNode, XorNode>
+    ) {
+        // commutative
+        return std::ranges::equal(operands_, rhs.operands_ | std::views::reverse);
+    } else {
+        // not commutative
+        static_assert(
+            std::same_as<BinaryOpNode, DivideNode> or      //
+            std::same_as<BinaryOpNode, LessEqualNode> or   //
+            std::same_as<BinaryOpNode, ModulusNode> or     //
+            std::same_as<BinaryOpNode, SafeDivideNode> or  //
+            std::same_as<BinaryOpNode, SubtractNode>
+        );
+    }
+
+    return false;
 }
 
 template <class BinaryOp>
@@ -254,7 +295,7 @@ void BinaryOpNode<BinaryOp>::initialize_state(State& state) const {
         unreachable();
     }
 
-    emplace_data_ptr_<ArrayNodeStateData>(state, std::move(values));
+    this->template emplace_data_ptr_<ArrayNodeStateData>(state, std::move(values));
 }
 
 template <class BinaryOp>
@@ -274,7 +315,7 @@ double BinaryOpNode<BinaryOp>::min() const {
 
 template <class BinaryOp>
 void BinaryOpNode<BinaryOp>::propagate(State& state) const {
-    auto ptr = data_ptr_<ArrayNodeStateData>(state);
+    auto ptr = this->template data_ptr_<ArrayNodeStateData>(state);
 
     const Array* lhs_ptr = operands_[0];
     const Array* rhs_ptr = operands_[1];
@@ -387,9 +428,18 @@ void BinaryOpNode<BinaryOp>::propagate(State& state) const {
     if (ptr->diff().size()) Node::propagate(state);
 }
 
+template <class BinaryOp>
+void BinaryOpNode<BinaryOp>::replace_predecessor_(ssize_t previous_index, Node* node_ptr) {
+    Node::replace_predecessor_(previous_index, node_ptr);
+
+    assert(0 <= previous_index and previous_index < 2);
+    operands_[previous_index] = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(operands_[previous_index] != nullptr);
+}
+
 template <class BinaryOp>
 void BinaryOpNode<BinaryOp>::revert(State& state) const {
-    data_ptr_<ArrayNodeStateData>(state)->revert();
+    this->template data_ptr_<ArrayNodeStateData>(state)->revert();
 }
 
 template <class BinaryOp>
@@ -434,7 +484,7 @@ ssize_t BinaryOpNode<BinaryOp>::size(const State& state) const {
 
 template <class BinaryOp>
 ssize_t BinaryOpNode<BinaryOp>::size_diff(const State& state) const {
-    return data_ptr_<ArrayNodeStateData>(state)->size_diff();
+    return this->template data_ptr_<ArrayNodeStateData>(state)->size_diff();
 }
 
 SizeInfo binaryop_calculate_sizeinfo(
diff --git a/dwave/optimization/src/nodes/collections.cpp b/dwave/optimization/src/nodes/collections.cpp
index 0ea4f805..a673f119 100644
--- a/dwave/optimization/src/nodes/collections.cpp
+++ b/dwave/optimization/src/nodes/collections.cpp
@@ -493,6 +493,10 @@ std::span<const Update> DisjointBitSetNode::diff(const State& state) const {
     return pred_data->diffs[set_index_];
 }
 
+bool DisjointBitSetNode::equal_to(const Node& rhs) const {
+    return static_cast<const Node*>(this) == &rhs;
+}
+
 bool DisjointBitSetNode::integral() const { return true; }
 
 double DisjointBitSetNode::min() const { return 0; }
@@ -844,6 +848,10 @@ std::span<const Update> DisjointListNode::diff(const State& state) const {
     return data->all_list_updates[list_index_];
 }
 
+bool DisjointListNode::equal_to(const Node& rhs) const {
+    return static_cast<const Node*>(this) == &rhs;
+}
+
 bool DisjointListNode::integral() const { return true; }
 
 double DisjointListNode::min() const { return 0; }
diff --git a/dwave/optimization/src/nodes/constants.cpp b/dwave/optimization/src/nodes/constants.cpp
index 98f79579..e20de6f0 100644
--- a/dwave/optimization/src/nodes/constants.cpp
+++ b/dwave/optimization/src/nodes/constants.cpp
@@ -65,6 +65,12 @@ ConstantNode::ConstantNode(OwningDataSource&& data_source, const std::span<const
     values_info_(calculate_values_info(std::span<const double>(buffer_ptr_, this->size()))),
     data_source_(std::make_unique<OwningDataSource>(std::move(data_source))) {}
 
+bool ConstantNode::equal_to(const ConstantNode& rhs) const {
+    return this->ndim() == rhs.ndim() and                      // same ndim
+           std::ranges::equal(this->shape(), rhs.shape()) and  // same shape
+           std::ranges::equal(this->data(), rhs.data());      // same content
+}
+
 bool ConstantNode::integral() const { return this->values_info_.integral; }
 
 double ConstantNode::min() const { return this->values_info_.min; }
diff --git a/dwave/optimization/src/nodes/creation.cpp b/dwave/optimization/src/nodes/creation.cpp
index 31860722..d68c32be 100644
--- a/dwave/optimization/src/nodes/creation.cpp
+++ b/dwave/optimization/src/nodes/creation.cpp
@@ -354,6 +354,10 @@ std::span<const Update> ARangeNode::diff(const State& state) const {
     return data_ptr_<ArrayNodeStateData>(state)->diff();
 }
 
+bool ARangeNode::equal_to(const ARangeNode& rhs) const {
+    return start_ == rhs.start_ and stop_ == rhs.stop_ and step_ == rhs.step_;
+}
+
 bool ARangeNode::integral() const { return true; }
 
 void ARangeNode::initialize_state(State& state) const {
@@ -422,6 +426,45 @@ void ARangeNode::propagate(State& state) const {
     if (ptr->diff().size()) Node::propagate(state);
 }
 
+void ARangeNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    assert(0 <= index and static_cast<size_t>(index) < predecessors().size());
+
+    const Array* array_ptr = dynamic_cast<const Array*>(node_ptr);
+    assert(array_ptr != nullptr);
+
+    if (std::holds_alternative<const Array*>(start_)) {
+        if (index == 0) {
+            // found it!
+            start_ = array_ptr;
+            return;
+        }
+
+        // If we haven't found the Array* we want to replace, then "normalize"
+        // the index so that it looks like it would if start wasn't an array.
+        index -= 1;
+    }
+
+    if (std::holds_alternative<const Array*>(stop_)) {
+        if (index == 0) {
+            // found it!
+            stop_ = array_ptr;
+            return;
+        }
+
+        // If we haven't found the Array* we want to replace, then "normalize"
+        // the index so that it looks like it would if stop wasn't an array.
+        index -= 1;
+    }
+
+    // If we're here then step must be an array and index must be 0 (whether
+    // because that's what was passed or because we've changed it).
+    assert(index == 0);
+    assert(std::holds_alternative<const Array*>(step_));
+    step_ = array_ptr;
+}
+
 void ARangeNode::revert(State& state) const { data_ptr_<ArrayNodeStateData>(state)->revert(); }
 
 std::span<const ssize_t> ARangeNode::shape(const State& state) const {
diff --git a/dwave/optimization/src/nodes/flow.cpp b/dwave/optimization/src/nodes/flow.cpp
index ebd464b2..f4b52d74 100644
--- a/dwave/optimization/src/nodes/flow.cpp
+++ b/dwave/optimization/src/nodes/flow.cpp
@@ -143,6 +143,19 @@ void ExtractNode::propagate(State& state) const {
     node_data->assign(std::move(new_values), count);
 }
 
+void ExtractNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    if (index == 0) {
+        condition_ptr_ = dynamic_cast<ArrayNode*>(node_ptr);
+        assert(condition_ptr_ != nullptr);
+    } else {
+        assert(index == 1);
+        arr_ptr_ = dynamic_cast<ArrayNode*>(node_ptr);
+        assert(condition_ptr_ != nullptr);
+    }
+}
+
 void ExtractNode::revert(State& state) const { data_ptr_<ArrayNodeStateData>(state)->revert(); }
 
 std::span<const ssize_t> ExtractNode::shape(const State& state) const {
@@ -364,6 +377,22 @@ void WhereNode::propagate(State& state) const {
     }
 }
 
+void WhereNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    if (index == 0) {
+        condition_ptr_ = dynamic_cast<ArrayNode*>(node_ptr);
+        assert(condition_ptr_ != nullptr);
+    } else if (index == 1) {
+        x_ptr_ = dynamic_cast<ArrayNode*>(node_ptr);
+        assert(x_ptr_ != nullptr);
+    } else {
+        assert(index == 2);
+        y_ptr_ = dynamic_cast<ArrayNode*>(node_ptr);
+        assert(y_ptr_ != nullptr);
+    }
+}
+
 void WhereNode::revert(State& state) const { data_ptr_<WhereNodeData>(state)->revert(); }
 
 std::span<const ssize_t> WhereNode::shape(const State& state) const {
diff --git a/dwave/optimization/src/nodes/indexing.cpp b/dwave/optimization/src/nodes/indexing.cpp
index c3eb6470..776a0eff 100644
--- a/dwave/optimization/src/nodes/indexing.cpp
+++ b/dwave/optimization/src/nodes/indexing.cpp
@@ -17,8 +17,11 @@
 #include <algorithm>
 #include <ranges>
 #include <unordered_set>
+#include <variant>
 
 #include "_state.hpp"
+#include "dwave-optimization/common.hpp"
+#include "dwave-optimization/graph.hpp"
 #include "dwave-optimization/nodes/constants.hpp"
 #include "dwave-optimization/state.hpp"
 
@@ -527,6 +530,10 @@ std::span<const Update> AdvancedIndexingNode::diff(const State& state) const {
     return data_ptr_<AdvancedIndexingNodeData>(state)->diff;
 }
 
+bool AdvancedIndexingNode::equal_to(const AdvancedIndexingNode& rhs) const {
+    return array_ptr_ == rhs.array_ptr_ and std::ranges::equal(indices_, rhs.indices_);
+}
+
 void AdvancedIndexingNode::fill_subspace(
     State& state,
     ssize_t array_offset,
@@ -743,6 +750,39 @@ void AdvancedIndexingNode::commit(State& state) const {
     data_ptr_<AdvancedIndexingNodeData>(state)->commit();
 }
 
+void AdvancedIndexingNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    assert(index >= 0);
+
+    // First check if we're replacing the array. If not, we change index to refer
+    // to the index of the indexing arrays
+    if (index == 0) {
+        // replace the array
+        array_ptr_ = dynamic_cast<ArrayNode*>(node_ptr);
+        assert(array_ptr_ != nullptr);
+        return;
+    }
+    index -= 1;
+
+    //
+    for (array_or_slice& indexer : indices_) {
+        if (std::holds_alternative<Slice>(indexer)) continue;  // carry on
+
+        assert(std::holds_alternative<ArrayNode*>(indexer));
+        if (index == 0) {
+            // found the indexer we want to replace
+            indexer = dynamic_cast<ArrayNode*>(node_ptr);
+            assert(std::get<ArrayNode*>(indexer) != nullptr);
+            return;
+        }
+        index -= 1;
+    }
+
+    assert(false and "should not be able to get here");
+    unreachable();
+}
+
 void AdvancedIndexingNode::revert(State& state) const {
     data_ptr_<AdvancedIndexingNodeData>(state)->revert();
 
@@ -1461,6 +1501,26 @@ std::span<const Update> BasicIndexingNode::diff(const State& state) const {
     return data_ptr_<BasicIndexingNodeData>(state)->diff;
 }
 
+bool BasicIndexingNode::equal_to(const BasicIndexingNode& rhs) const {
+    if (array_ptr_ == rhs.array_ptr_ and              //
+        ndim_ == rhs.ndim_ and                        //
+        std::ranges::equal(shape(), rhs.shape()) and  //
+        std::ranges::equal(strides(), rhs.strides())) {
+        if (size_ < 0) {
+            // Dynamic basic indexing, so we need to inspect the first slice.
+            // Note that we don't normalize the first slice in the constructor
+            // so we can get false negatives here. E.g. x[:10000], x[:10001] in
+            // some cases could evaluate to equal.
+            return axis0_slice_ == rhs.axis0_slice_;
+        } else {
+            // Otherwise we inspect the start_
+            return start_ == rhs.start_;
+        }
+    }
+
+    return false;
+}
+
 std::vector<BasicIndexingNode::slice_or_int> BasicIndexingNode::infer_indices() const {
     std::vector<slice_or_int> indices;
 
@@ -1889,6 +1949,19 @@ void BasicIndexingNode::propagate(State& state) const {
     if (diff.size()) Node::propagate(state);
 }
 
+void BasicIndexingNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    assert(index == 0);
+
+    ArrayNode* array_ptr = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr != nullptr);
+    assert(std::ranges::equal(this->array_ptr_->shape(), array_ptr->shape()));
+    assert(std::ranges::equal(this->array_ptr_->strides(), array_ptr->strides()));
+
+    this->array_ptr_ = array_ptr;
+}
+
 void BasicIndexingNode::revert(State& state) const {
     auto node_data = data_ptr_<BasicIndexingNodeData>(state);
     node_data->diff.clear();
@@ -1981,6 +2054,10 @@ std::span<const Update> PermutationNode::diff(const State& state) const {
     return data_ptr_<IndexingNodeData>(state)->diff;
 }
 
+bool PermutationNode::equal_to(const PermutationNode& rhs) const {
+    return array_ptr_ == rhs.array_ptr_ and order_ptr_ == rhs.order_ptr_;
+}
+
 bool PermutationNode::integral() const { return values_info_.integral; }
 
 double PermutationNode::max() const { return values_info_.max; }
@@ -2088,6 +2165,19 @@ void PermutationNode::propagate(State& state) const {
     if (updates.size()) Node::propagate(state);
 }
 
+void PermutationNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    if (index == 0) {
+        array_ptr_ = dynamic_cast<ArrayNode*>(node_ptr);
+        assert(array_ptr_ != nullptr);
+    } else {
+        assert(index == 1);
+        order_ptr_ = dynamic_cast<ArrayNode*>(node_ptr);
+        assert(order_ptr_ != nullptr);
+    }
+}
+
 void PermutationNode::revert(State& state) const { data_ptr_<IndexingNodeData>(state)->revert(); }
 
 }  // namespace dwave::optimization
diff --git a/dwave/optimization/src/nodes/inputs.cpp b/dwave/optimization/src/nodes/inputs.cpp
index ba0967be..63f0b0d9 100644
--- a/dwave/optimization/src/nodes/inputs.cpp
+++ b/dwave/optimization/src/nodes/inputs.cpp
@@ -99,6 +99,10 @@ void InputNode::initialize_state(State& state, std::span<const double> data) con
     emplace_data_ptr_<ArrayNodeStateData>(state, std::move(copy));
 }
 
+void InputNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    assert(false and "InputNode never has any predecessors");
+}
+
 void InputNode::revert(State& state) const noexcept {
     data_ptr_<ArrayNodeStateData>(state)->revert();
 }
diff --git a/dwave/optimization/src/nodes/interpolation.cpp b/dwave/optimization/src/nodes/interpolation.cpp
index 45bc9f0f..0f3cca56 100644
--- a/dwave/optimization/src/nodes/interpolation.cpp
+++ b/dwave/optimization/src/nodes/interpolation.cpp
@@ -127,6 +127,11 @@ void BSplineNode::commit(State& state) const {
 std::span<const Update> BSplineNode::diff(const State& state) const {
     return data_ptr_<ArrayNodeStateData>(state)->diff();
 }
+
+bool BSplineNode::equal_to(const BSplineNode& rhs) const {
+    return array_ptr_ == rhs.array_ptr_ and k_ == rhs.k_ and t_ == rhs.t_ and c_ == rhs.c_;
+}
+
 bool BSplineNode::integral() const { return false; }
 
 void BSplineNode::revert(State& state) const {
@@ -157,4 +162,12 @@ void BSplineNode::propagate(State& state) const {
     }
 }
 
+void BSplineNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    assert(index == 0);
+    array_ptr_ = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr_ != nullptr);
+}
+
 }  // namespace dwave::optimization
diff --git a/dwave/optimization/src/nodes/lambda.cpp b/dwave/optimization/src/nodes/lambda.cpp
index bdb0bd4d..84976b51 100644
--- a/dwave/optimization/src/nodes/lambda.cpp
+++ b/dwave/optimization/src/nodes/lambda.cpp
@@ -16,6 +16,7 @@
 
 #include "_state.hpp"
 #include "dwave-optimization/array.hpp"
+#include "dwave-optimization/graph.hpp"
 #include "dwave-optimization/nodes/inputs.hpp"
 #include "dwave-optimization/state.hpp"
 
@@ -44,15 +45,15 @@ AccumulateZipNode::AccumulateZipNode(
     array_or_double initial
 ) :
     ArrayOutputMixin(operands.empty() ? std::span<ssize_t, 0>() : operands[0]->shape()),
-    initial(initial),
     expression_ptr_(std::move(expression_ptr)),
+    initial_(initial),
     operands_(operands),
     sizeinfo_(operands.empty() ? SizeInfo(0) : operands_[0]->sizeinfo()) {
-    check(*expression_ptr_, operands, initial);
+    check(*expression_ptr_, operands, initial_);
 
-    if (std::holds_alternative<ArrayNode*>(initial)) {
+    if (std::holds_alternative<ArrayNode*>(initial_)) {
         // was checked in check() method
-        add_predecessor_(std::get<ArrayNode*>(initial));
+        add_predecessor_(std::get<ArrayNode*>(initial_));
     }
     for (const auto& op : operands_) {
         add_predecessor_(op);
@@ -242,10 +243,10 @@ double AccumulateZipNode::evaluate_expression(State& register_) const {
 }
 
 double AccumulateZipNode::get_initial_value(const State& state) const {
-    if (std::holds_alternative<double>(initial)) {
-        return std::get<double>(initial);
+    if (std::holds_alternative<double>(initial_)) {
+        return std::get<double>(initial_);
     } else {
-        return std::get<ArrayNode*>(initial)->view(state)[0];
+        return std::get<ArrayNode*>(initial_)->view(state)[0];
     }
 }
 
@@ -301,6 +302,16 @@ void AccumulateZipNode::initialize_state(State& state) const {
     );
 }
 
+bool AccumulateZipNode::equal_to(const AccumulateZipNode& rhs) const {
+    // note that we don't have a notion of Graph equality, so we instead
+    // just check whether we have the same underlying shared ptr
+    return (
+        expression_ptr_ == rhs.expression_ptr_ and  //
+        operands_ == rhs.operands_ and              //
+        initial_ == rhs.initial_
+    );
+}
+
 bool AccumulateZipNode::integral() const { return expression_ptr_->objective()->integral(); }
 
 double AccumulateZipNode::max() const { return expression_ptr_->objective()->max(); }
@@ -354,6 +365,24 @@ const InputNode* const AccumulateZipNode::accumulate_input() const {
     return expression_ptr_->inputs()[0];
 }
 
+void AccumulateZipNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    if (std::holds_alternative<ArrayNode*>(initial_)) {
+        if (index == 0) {
+            // we're replacing the initial array
+            initial_ = dynamic_cast<ArrayNode*>(node_ptr);
+            assert(std::get<ArrayNode*>(initial_) != nullptr);
+            return;
+        }
+
+        index -= 1;
+    }
+
+    operands_[index] = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(operands_[index] != nullptr);
+}
+
 void AccumulateZipNode::revert(State& state) const {
     data_ptr_<AccumulateZipNodeData>(state)->revert();
 }
diff --git a/dwave/optimization/src/nodes/linear_algebra.cpp b/dwave/optimization/src/nodes/linear_algebra.cpp
index 267b0793..94f85b43 100644
--- a/dwave/optimization/src/nodes/linear_algebra.cpp
+++ b/dwave/optimization/src/nodes/linear_algebra.cpp
@@ -15,15 +15,16 @@
 #include "dwave-optimization/nodes/linear_algebra.hpp"
 
 #include <algorithm>
-#include <array>
 #include <ranges>
+#include <vector>
+
+#include "dwave-optimization/graph.hpp"
 
 #if __has_include(<openblas_config.h>) and __has_include(<cblas.h>)
 #include <cblas.h>
 #define HAS_BLAS_
 #endif
 
-#include "../functional_.hpp"
 #include "_state.hpp"
 #include "dwave-optimization/array.hpp"
 #include "dwave-optimization/state.hpp"
@@ -381,8 +382,7 @@ class MatrixMultiplyNodeData : public ArrayNodeStateData {
 
 MatrixMultiplyNode::MatrixMultiplyNode(ArrayNode* x_ptr, ArrayNode* y_ptr) :
     ArrayOutputMixin(output_shape(x_ptr, y_ptr)),
-    x_ptr_(x_ptr),
-    y_ptr_(y_ptr),
+    operands_{x_ptr, y_ptr},
     sizeinfo_(get_sizeinfo(x_ptr, y_ptr)),
     values_info_(get_values_info(x_ptr, y_ptr)) {
     add_predecessor_(x_ptr);
@@ -421,15 +421,18 @@ void MatrixMultiplyNode::matmul_(
     // to handle x/y with more than 2 dimensions. We do all that handling in terms of
     // "leaps", i.e. the number of pointer increments to apply in each dimension.
 
-    // const ssize_t x_penultimate_axis_size = get_axis_size(x_ptr_->shape(state), -2, true);
+    const ArrayNode* const x_ptr = operands_[0];
+    const ArrayNode* const y_ptr = operands_[1];
+
+    // const ssize_t x_penultimate_axis_size = get_axis_size(x_ptr->shape(state), -2, true);
     const ssize_t leading_subspace_size =
-        get_leading_subspace_size(x_ptr_->shape(state), y_ptr_->shape(state));
+        get_leading_subspace_size(x_ptr->shape(state), y_ptr->shape(state));
 
-    const ssize_t x_leading_leap = get_leading_leap(x_ptr_->shape(state));
-    const ssize_t y_leading_leap = get_leading_leap(y_ptr_->shape(state));
+    const ssize_t x_leading_leap = get_leading_leap(x_ptr->shape(state));
+    const ssize_t y_leading_leap = get_leading_leap(y_ptr->shape(state));
     const ssize_t out_leading_leap = [&]() -> ssize_t {
-        if (x_ptr_->ndim() >= 2 and y_ptr_->ndim() >= 2) return get_leading_leap(out_shape);
-        if (x_ptr_->ndim() == 1 and y_ptr_->ndim() == 1) return 0;
+        if (x_ptr->ndim() >= 2 and y_ptr->ndim() >= 2) return get_leading_leap(out_shape);
+        if (x_ptr->ndim() == 1 and y_ptr->ndim() == 1) return 0;
         return out_shape.back();
     }();
 
@@ -443,9 +446,9 @@ void MatrixMultiplyNode::matmul_(
         return vals[static_cast<ssize_t>(vals.size()) + index];
     };
 
-    const ssize_t m = (x_ptr_->ndim() > 1) ? get(x_ptr_->shape(state), -2) : 1;
-    const ssize_t n = (y_ptr_->ndim() > 1) ? get(y_ptr_->shape(state), -1) : 1;
-    const ssize_t k = get(x_ptr_->shape(state), -1);
+    const ssize_t m = (x_ptr->ndim() > 1) ? get(x_ptr->shape(state), -2) : 1;
+    const ssize_t n = (y_ptr->ndim() > 1) ? get(y_ptr->shape(state), -1) : 1;
+    const ssize_t k = get(x_ptr->shape(state), -1);
 
     // We also need the stride information about x/y/out
 
@@ -454,14 +457,14 @@ void MatrixMultiplyNode::matmul_(
             return std::array<ssize_t, 2>{k * static_cast<ssize_t>(sizeof(double)), strides[0]};
         }
         return std::array<ssize_t, 2>{get(strides, -2), get(strides, -1)};
-    }(x_ptr_->strides());
+    }(x_ptr->strides());
 
     std::array<ssize_t, 2> y_matmul_strides = [&get](std::span<const ssize_t> strides) {
         if (strides.size() == 1) {
             return std::array<ssize_t, 2>{strides[0], sizeof(double)};
         }
         return std::array<ssize_t, 2>{get(strides, -2), get(strides, -1)};
-    }(y_ptr_->strides());
+    }(y_ptr->strides());
 
     std::array<ssize_t, 2> out_matmul_strides{
         n * static_cast<ssize_t>(sizeof(double)), sizeof(double)
@@ -472,8 +475,8 @@ void MatrixMultiplyNode::matmul_(
         // and checking the strides of both x/y, we use ArrayIterators to
         // get us to the correct subspace, and then use the strides of the
         // predecessors to iterate through the last one or two dimensions.
-        const double* const x_data = &x_ptr_->view(state).begin()[w * x_leading_leap];
-        const double* const y_data = &y_ptr_->view(state).begin()[w * y_leading_leap];
+        const double* const x_data = &x_ptr->view(state).begin()[w * x_leading_leap];
+        const double* const y_data = &y_ptr->view(state).begin()[w * y_leading_leap];
         double* const out_data = out.data() + w * out_leading_leap;
 
         gemm(
@@ -494,7 +497,7 @@ void MatrixMultiplyNode::initialize_state(State& state) const {
     ssize_t start_size = this->size();
     std::vector<ssize_t> shape(this->shape().begin(), this->shape().end());
     if (this->dynamic()) {
-        shape[0] = x_ptr_->shape(state)[0];
+        shape[0] = operands_[0]->shape(state)[0];
         start_size = Array::shape_to_size(shape);
     }
 
@@ -523,12 +526,12 @@ double MatrixMultiplyNode::min() const { return values_info_.min; }
 
 void MatrixMultiplyNode::update_shape_(State& state) const {
     if (this->dynamic()) {
-        data_ptr_<MatrixMultiplyNodeData>(state)->shape[0] = x_ptr_->shape(state)[0];
+        data_ptr_<MatrixMultiplyNodeData>(state)->shape[0] = operands_[0]->shape(state)[0];
     }
 }
 
 void MatrixMultiplyNode::propagate(State& state) const {
-    if (x_ptr_->diff(state).size() == 0 and y_ptr_->diff(state).size() == 0) return;
+    if (operands_[0]->diff(state).size() == 0 and operands_[1]->diff(state).size() == 0) return;
 
     auto data = data_ptr_<MatrixMultiplyNodeData>(state);
 
@@ -541,6 +544,21 @@ void MatrixMultiplyNode::propagate(State& state) const {
     data->assign(data->output);
 }
 
+void MatrixMultiplyNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    const ArrayNode* array_ptr = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr != nullptr);
+
+    if (index == 0) {
+        assert(std::ranges::equal(operands_[0]->shape(), array_ptr->shape()));
+        operands_[0] = array_ptr;
+    } else {
+        assert(std::ranges::equal(operands_[1]->shape(), array_ptr->shape()));
+        operands_[1] = array_ptr;
+    }
+}
+
 void MatrixMultiplyNode::revert(State& state) const {
     auto data = data_ptr_<MatrixMultiplyNodeData>(state);
     data->revert();
diff --git a/dwave/optimization/src/nodes/lp.cpp b/dwave/optimization/src/nodes/lp.cpp
index b7ea1714..ddf281f4 100644
--- a/dwave/optimization/src/nodes/lp.cpp
+++ b/dwave/optimization/src/nodes/lp.cpp
@@ -18,6 +18,8 @@
 
 #include "../simplex.hpp"
 #include "_state.hpp"
+#include "dwave-optimization/common.hpp"
+#include "dwave-optimization/graph.hpp"
 
 namespace dwave::optimization {
 
@@ -60,6 +62,12 @@ void LinearProgramFeasibleNode::propagate(State& state) const {
     set_state(state, lp_ptr_->feasible(state));
 }
 
+void LinearProgramFeasibleNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+    lp_ptr_ = dynamic_cast<LinearProgramNodeBase*>(node_ptr);
+    assert(lp_ptr_ != nullptr);
+}
+
 void LinearProgramNodeBase::check_input_arguments(
     const ArrayNode* c_ptr,
     const ArrayNode* b_lb_ptr,
@@ -207,6 +215,19 @@ void LinearProgramNode::commit(State& state) const {};
 
 bool LinearProgramNode::deterministic_state() const { return false; }
 
+bool LinearProgramNode::equal_to(const LinearProgramNode& rhs) const {
+    return (
+        c_ptr_ == rhs.c_ptr_ and        //
+        b_lb_ptr_ == rhs.b_lb_ptr_ and  //
+        A_ptr_ == rhs.A_ptr_ and        //
+        b_ub_ptr_ == rhs.b_ub_ptr_ and  //
+        A_eq_ptr_ == rhs.A_eq_ptr_ and  //
+        b_eq_ptr_ == rhs.b_eq_ptr_ and  //
+        lb_ptr_ == rhs.lb_ptr_ and      //
+        ub_ptr_ == rhs.ub_ptr_          //
+    );
+}
+
 bool LinearProgramNode::feasible(const State& state) const {
     return data_ptr_<LinearProgramNodeData>(state)->result.feasible();
 }
@@ -340,6 +361,41 @@ void LinearProgramNode::propagate(State& state) const {
     Node::propagate(state);
 }
 
+void LinearProgramNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    const ArrayNode* array_ptr = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr != nullptr);
+
+    auto check_and_replace = [&array_ptr, &index](const ArrayNode*& to_replace) -> bool {
+        if (to_replace == nullptr) return false;  // nothing to replace and no need to re-index
+        if (index-- > 0) return false;  // found a node, but it's not the one we're replacing
+
+        // sanity check the size. There are other things we could check, but this is easy
+        // and simple
+        assert(std::ranges::equal(to_replace->shape(), array_ptr->shape()));
+        to_replace = array_ptr;
+
+        return true;  // we found the one we wanted to replace.
+    };
+
+    assert(c_ptr_ != nullptr);  // never null
+    if (check_and_replace(c_ptr_)) return;
+
+    if (check_and_replace(b_lb_ptr_)) return;
+    if (check_and_replace(A_ptr_)) return;
+    if (check_and_replace(b_ub_ptr_)) return;
+
+    if (check_and_replace(A_eq_ptr_)) return;
+    if (check_and_replace(b_eq_ptr_)) return;
+
+    if (check_and_replace(lb_ptr_)) return;
+    if (check_and_replace(ub_ptr_)) return;
+
+    unreachable();
+    assert(false and "should never get here");
+}
+
 void LinearProgramNode::revert(State& state) const {
     // Nothing to do on revert as all changes are tracked by successor nodes
 }
@@ -373,6 +429,12 @@ void LinearProgramObjectiveValueNode::propagate(State& state) const {
     // We could consider setting ourselves to max()
 }
 
+void LinearProgramObjectiveValueNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+    lp_ptr_ = dynamic_cast<LinearProgramNodeBase*>(node_ptr);
+    assert(lp_ptr_ != nullptr);
+}
+
 LinearProgramSolutionNode::LinearProgramSolutionNode(LinearProgramNodeBase* lp_ptr) :
     ArrayOutputMixin(lp_ptr->variables_shape()), lp_ptr_(lp_ptr) {
     add_predecessor_(lp_ptr);
@@ -439,6 +501,12 @@ void LinearProgramSolutionNode::propagate(State& state) const {
     }
 }
 
+void LinearProgramSolutionNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+    lp_ptr_ = dynamic_cast<LinearProgramNodeBase*>(node_ptr);
+    assert(lp_ptr_ != nullptr);
+}
+
 void LinearProgramSolutionNode::revert(State& state) const {
     data_ptr_<ArrayNodeStateData>(state)->revert();
 }
diff --git a/dwave/optimization/src/nodes/manipulation.cpp b/dwave/optimization/src/nodes/manipulation.cpp
index 7260c123..57223f87 100644
--- a/dwave/optimization/src/nodes/manipulation.cpp
+++ b/dwave/optimization/src/nodes/manipulation.cpp
@@ -206,6 +206,10 @@ std::span<const Update> BroadcastToNode::diff(const State& state) const {
     return data_ptr_<BroadcastToNodeData>(state)->diff;
 }
 
+bool BroadcastToNode::equal_to(const BroadcastToNode& rhs) const {
+    return array_ptr_ == rhs.array_ptr_ and std::ranges::equal(shape(), rhs.shape());
+}
+
 void BroadcastToNode::initialize_state(State& state) const {
     std::vector<ssize_t> offsets = diff_offsets(array_ptr_->shape(), this->shape());
 
@@ -355,6 +359,16 @@ ssize_t BroadcastToNode::convert_predecessor_index_(ssize_t index) const {
     return flat_index;
 }
 
+void BroadcastToNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    assert(index == 0);  // we should only ever have one predecessor
+    ArrayNode* array_ptr = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(std::ranges::equal(array_ptr_->shape(), array_ptr->shape()));
+    assert(array_ptr != nullptr);
+    array_ptr_ = array_ptr;
+}
+
 void BroadcastToNode::revert(State& state) const {
     data_ptr_<BroadcastToNodeData>(state)->revert();
 }
@@ -430,6 +444,10 @@ std::span<const Update> ConcatenateNode::diff(const State& state) const {
     return data_ptr_<ArrayNodeStateData>(state)->diff();
 }
 
+bool ConcatenateNode::equal_to(const ConcatenateNode& rhs) const {
+    return axis_ == rhs.axis_ and std::ranges::equal(array_ptrs_, rhs.array_ptrs_);
+}
+
 void ConcatenateNode::initialize_state(State& state) const {
     std::vector<double> values;
     values.resize(size());
@@ -539,6 +557,17 @@ void ConcatenateNode::propagate(State& state) const {
     }
 }
 
+void ConcatenateNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    ArrayNode* array_ptr = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr != nullptr);
+    assert(std::ranges::equal(array_ptrs_[index]->shape(), array_ptr->shape()));
+
+    assert(0 <= index and static_cast<size_t>(index) < array_ptrs_.size());
+    array_ptrs_[index] = array_ptr;
+}
+
 void ConcatenateNode::revert(State& state) const { data_ptr_<ArrayNodeStateData>(state)->revert(); }
 
 CopyNode::CopyNode(ArrayNode* array_ptr) :
@@ -570,6 +599,17 @@ void CopyNode::propagate(State& state) const {
     data_ptr_<ArrayNodeStateData>(state)->update(array_ptr_->diff(state));
 }
 
+void CopyNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    ArrayNode* array_ptr = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr != nullptr);
+    assert(std::ranges::equal(array_ptr_->shape(), array_ptr->shape()));
+
+    assert(index == 0);
+    array_ptr_ = array_ptr;
+}
+
 void CopyNode::revert(State& state) const { data_ptr_<ArrayNodeStateData>(state)->revert(); }
 
 std::span<const ssize_t> CopyNode::shape(const State& state) const {
@@ -860,6 +900,25 @@ void PutNode::propagate(State& state) const {
     if (ptr->diff().size()) Node::propagate(state);
 }
 
+void PutNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    ArrayNode* array_ptr = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr != nullptr);
+
+    if (index == 0) {
+        assert(std::ranges::equal(array_ptr_->shape(), array_ptr->shape()));
+        array_ptr_ = array_ptr;
+    } else if (index == 1) {
+        assert(std::ranges::equal(indices_ptr_->shape(), array_ptr->shape()));
+        indices_ptr_ = array_ptr;
+    } else {
+        assert(index == 2);
+        assert(std::ranges::equal(values_ptr_->shape(), array_ptr->shape()));
+        values_ptr_ = array_ptr;
+    }
+}
+
 void PutNode::revert(State& state) const { return data_ptr_<PutNodeState>(state)->revert(); }
 
 // Reshape allows one shape dimension to be -1. In that case the size is inferred.
@@ -1002,6 +1061,10 @@ std::span<const Update> ReshapeNode::diff(const State& state) const {
     return array_ptr_->diff(state);
 }
 
+bool ReshapeNode::equal_to(const ReshapeNode& rhs) const {
+    return array_ptr_ == rhs.array_ptr_ and std::ranges::equal(shape(), rhs.shape());
+}
+
 void ReshapeNode::initialize_state(State& state) const {
     if (!this->dynamic()) return Node::initialize_state(state);  // stateless
     emplace_data_ptr_<DynamicReshapeNodeData>(state, this->shape(), array_ptr_->size(state));
@@ -1018,6 +1081,17 @@ void ReshapeNode::propagate(State& state) const {
     data_ptr_<DynamicReshapeNodeData>(state)->set_size(array_ptr_->size(state));
 }
 
+void ReshapeNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    ArrayNode* array_ptr = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr != nullptr);
+    assert(std::ranges::equal(array_ptr_->shape(), array_ptr->shape()));
+
+    assert(index == 0);
+    array_ptr_ = array_ptr;
+}
+
 void ReshapeNode::revert(State& state) const {
     if (!this->dynamic()) return;  // stateless
     return data_ptr_<DynamicReshapeNodeData>(state)->revert();
@@ -1094,6 +1168,14 @@ std::span<const Update> ResizeNode::diff(const State& state) const {
     return data_ptr_<ArrayNodeStateData>(state)->diff();
 }
 
+bool ResizeNode::equal_to(const ResizeNode& rhs) const {
+    return (
+        array_ptr_ == rhs.array_ptr_ and              //
+        std::ranges::equal(shape(), rhs.shape()) and  //
+        fill_value_ == rhs.fill_value_
+    );
+}
+
 void ResizeNode::initialize_state(State& state) const {
     const ssize_t size = this->size();  // the desired size of our state
     assert(size >= 0);                  // we're never dynamic
@@ -1147,6 +1229,17 @@ void ResizeNode::propagate(State& state) const {
     if (data_ptr->diff().size()) Node::propagate(state);
 }
 
+void ResizeNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    ArrayNode* array_ptr = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr != nullptr);
+    assert(std::ranges::equal(array_ptr_->shape(), array_ptr->shape()));
+
+    assert(index == 0);
+    array_ptr_ = array_ptr;
+}
+
 void ResizeNode::revert(State& state) const {
     return data_ptr_<ArrayNodeStateData>(state)->revert();
 }
@@ -1242,6 +1335,14 @@ std::span<const Update> RollNode::diff(const State& state) const {
     return data_ptr_<ArrayNodeStateData>(state)->diff();
 }
 
+bool RollNode::equal_to(const RollNode& rhs) const {
+    return (
+        array_ptr_ == rhs.array_ptr_ and  //
+        shift_ == rhs.shift_ and          //
+        std::ranges::equal(axis_, rhs.axis_)
+    );
+}
+
 void RollNode::initialize_state(State& state) const {
     // Get the predecessor's state as a vector that will become our state.
     std::vector<double> buffer(array_ptr_->begin(state), array_ptr_->end(state));
@@ -1358,6 +1459,23 @@ void RollNode::propagate(State& state) const {
     }
 }
 
+void RollNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    ArrayNode* array_ptr = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr != nullptr);
+
+    if (index == 0) {
+        assert(std::ranges::equal(array_ptr_->shape(), array_ptr->shape()));
+        array_ptr_ = array_ptr;
+    } else {
+        assert(index == 1);
+        assert(std::holds_alternative<const Array*>(shift_));
+        assert(std::ranges::equal(std::get<const Array*>(shift_)->shape(), array_ptr->shape()));
+        shift_ = array_ptr;
+    }
+}
+
 void RollNode::revert(State& state) const { data_ptr_<ArrayNodeStateData>(state)->revert(); }
 
 // Act on the given array *in place*.
@@ -1500,6 +1618,14 @@ double SizeNode::max() const { return minmax_.second; }
 
 void SizeNode::propagate(State& state) const { set_state(state, array_ptr_->size(state)); }
 
+void SizeNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    assert(index == 0);
+    array_ptr_ = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr_ != nullptr);
+}
+
 // TransposeNode **************************************************************
 
 ArrayNode* TransposeNode::predeccesor_check_(ArrayNode* array_ptr) const {
@@ -1692,4 +1818,16 @@ void TransposeNode::revert(State& state) const {
     }  // otherwise, stateless
 }
 
+void TransposeNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    ArrayNode* array_ptr = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr != nullptr);
+    assert(std::ranges::equal(array_ptr_->shape(), array_ptr->shape()));
+    assert(std::ranges::equal(array_ptr_->strides(), array_ptr->strides()));
+
+    assert(index == 0);
+    array_ptr_ = array_ptr;
+}
+
 }  // namespace dwave::optimization
diff --git a/dwave/optimization/src/nodes/naryop.cpp b/dwave/optimization/src/nodes/naryop.cpp
index 420d9ebb..76377d49 100644
--- a/dwave/optimization/src/nodes/naryop.cpp
+++ b/dwave/optimization/src/nodes/naryop.cpp
@@ -315,6 +315,18 @@ void NaryOpNode<BinaryOp>::propagate(State& state) const {
     }
 }
 
+template <class BinaryOp>
+void NaryOpNode<BinaryOp>::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    assert(0 <= index and static_cast<size_t>(index) < operands_.size());
+    ArrayNode* array_ptr = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr != nullptr);
+
+    assert(std::ranges::equal(operands_[index]->shape(), array_ptr->shape()));
+    operands_[index] = array_ptr;
+}
+
 template class NaryOpNode<functional::max<double>>;
 template class NaryOpNode<functional::min<double>>;
 template class NaryOpNode<std::multiplies<double>>;
diff --git a/dwave/optimization/src/nodes/reduce.cpp b/dwave/optimization/src/nodes/reduce.cpp
index 487c5ee7..8b1794aa 100644
--- a/dwave/optimization/src/nodes/reduce.cpp
+++ b/dwave/optimization/src/nodes/reduce.cpp
@@ -801,13 +801,13 @@ ReduceNode<BinaryOp>::ReduceNode(
     std::span<const ssize_t> axes,
     std::optional<double> initial
 ) :
-    ArrayOutputMixin(reduce_shape(array_ptr, axes)),
+    ArrayOutputMixin<EqualityMixin<ArrayNode, ReduceNode<BinaryOp>>>(reduce_shape(array_ptr, axes)),
     initial(initial),
     array_ptr_(array_ptr),
     axes_(normalize_axes(array_ptr, axes)),
     values_info_(values_info<BinaryOp>(array_ptr_, axes_, initial)),
     sizeinfo_(reducenode_calculate_sizeinfo(this, array_ptr_, axes_)) {
-    add_predecessor_(array_ptr);
+    this->add_predecessor_(array_ptr);
 }
 
 template <class BinaryOp>
@@ -820,17 +820,26 @@ ReduceNode<BinaryOp>::ReduceNode(
 
 template <class BinaryOp>
 double const* ReduceNode<BinaryOp>::buff(const State& state) const {
-    return data_ptr_<ReduceNodeData<BinaryOp>>(state)->buff();
+    return this->template data_ptr_<ReduceNodeData<BinaryOp>>(state)->buff();
 }
 
 template <class BinaryOp>
 void ReduceNode<BinaryOp>::commit(State& state) const {
-    return data_ptr_<ReduceNodeData<BinaryOp>>(state)->commit();
+    return this->template data_ptr_<ReduceNodeData<BinaryOp>>(state)->commit();
 }
 
 template <class BinaryOp>
 std::span<const Update> ReduceNode<BinaryOp>::diff(const State& state) const {
-    return data_ptr_<ReduceNodeData<BinaryOp>>(state)->diff();
+    return this->template data_ptr_<ReduceNodeData<BinaryOp>>(state)->diff();
+}
+
+template <class BinaryOp>
+bool ReduceNode<BinaryOp>::equal_to(const ReduceNode& rhs) const {
+    return (
+        array_ptr_ == rhs.array_ptr_ and  //
+        initial == rhs.initial and        //
+        std::ranges::equal(axes_, rhs.axes_)
+    );
 }
 
 template <class BinaryOp>
@@ -846,12 +855,14 @@ void ReduceNode<BinaryOp>::initialize_state(State& state) const {
             reductions.emplace_back(reduce_(state, index));
         }
 
-        emplace_data_ptr_<ReduceNodeData<BinaryOp>>(state, std::move(reductions), this->shape());
+        this->template emplace_data_ptr_<ReduceNodeData<BinaryOp>>(
+            state, std::move(reductions), this->shape()
+        );
     } else {
         for (ssize_t index = 0; index < this->size(); ++index) {
             reductions.emplace_back(reduce_(state, index));
         }
-        emplace_data_ptr_<ReduceNodeData<BinaryOp>>(state, std::move(reductions));
+        this->template emplace_data_ptr_<ReduceNodeData<BinaryOp>>(state, std::move(reductions));
     }
 }
 
@@ -944,7 +955,7 @@ ssize_t ReduceNode<BinaryOp>::convert_predecessor_index_(ssize_t index) const {
 
 template <class BinaryOp>
 void ReduceNode<BinaryOp>::propagate(State& state) const {
-    auto* const state_ptr = data_ptr_<ReduceNodeData<BinaryOp>>(state);
+    auto* const state_ptr = this->template data_ptr_<ReduceNodeData<BinaryOp>>(state);
 
     // We are reducing over all axes, so this is nice and simple
     if (axes_.empty() or axes_.size() == static_cast<std::size_t>(array_ptr_->ndim())) {
@@ -1066,27 +1077,39 @@ auto ReduceNode<BinaryOp>::reduce_(const State& state, const ssize_t index) cons
     return ufunc.reduce(std::ranges::subrange(it, std::default_sentinel), initial);
 }
 
+template <class BinaryOp>
+void ReduceNode<BinaryOp>::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    ArrayNode* array_ptr = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr != nullptr);
+    assert(std::ranges::equal(array_ptr_->shape(), array_ptr->shape()));
+
+    assert(index == 0);
+    array_ptr_ = array_ptr;
+}
+
 template <class BinaryOp>
 void ReduceNode<BinaryOp>::revert(State& state) const {
-    return data_ptr_<ReduceNodeData<BinaryOp>>(state)->revert();
+    return this->template data_ptr_<ReduceNodeData<BinaryOp>>(state)->revert();
 }
 
 template <class BinaryOp>
 std::span<const ssize_t> ReduceNode<BinaryOp>::shape(const State& state) const {
     if (ssize_t size = this->size(); size >= 0) return this->shape();  // if we're not dynamic
-    return data_ptr_<ReduceNodeData<BinaryOp>>(state)->shape();
+    return this->template data_ptr_<ReduceNodeData<BinaryOp>>(state)->shape();
 }
 
 template <class BinaryOp>
 ssize_t ReduceNode<BinaryOp>::size(const State& state) const {
     if (ssize_t size = this->size(); size >= 0) return size;  // if we're not dynamic
-    return data_ptr_<ReduceNodeData<BinaryOp>>(state)->size();
+    return this->template data_ptr_<ReduceNodeData<BinaryOp>>(state)->size();
 }
 
 template <class BinaryOp>
 ssize_t ReduceNode<BinaryOp>::size_diff(const State& state) const {
     if (ssize_t size = this->size(); size >= 0) return 0;  // if we're not dynamic
-    return data_ptr_<ReduceNodeData<BinaryOp>>(state)->size_diff();
+    return this->template data_ptr_<ReduceNodeData<BinaryOp>>(state)->size_diff();
 }
 
 template class ReduceNode<functional::max<double>>;
diff --git a/dwave/optimization/src/nodes/set_routines.cpp b/dwave/optimization/src/nodes/set_routines.cpp
index a15301a4..8251f76b 100644
--- a/dwave/optimization/src/nodes/set_routines.cpp
+++ b/dwave/optimization/src/nodes/set_routines.cpp
@@ -230,6 +230,20 @@ void IsInNode::propagate(State& state) const {
     assert(set_data_is_correct(state, element_ptr_, *node_data));
 }
 
+void IsInNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    ArrayNode* array_ptr = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr != nullptr);
+
+    if (index == 0) {
+        element_ptr_ = array_ptr;
+    } else {
+        assert(index == 1);
+        test_elements_ptr_ = array_ptr;
+    }
+}
+
 void IsInNode::revert(State& state) const {
     IsInNodeData* node_data = data_ptr_<IsInNodeData>(state);
 
diff --git a/dwave/optimization/src/nodes/softmax.cpp b/dwave/optimization/src/nodes/softmax.cpp
index dce5136a..467a413a 100644
--- a/dwave/optimization/src/nodes/softmax.cpp
+++ b/dwave/optimization/src/nodes/softmax.cpp
@@ -139,6 +139,17 @@ void SoftMaxNode::propagate(State& state) const {
     node_data->prior_denominator = prior_denominator;
 }
 
+void SoftMaxNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    ArrayNode* array_ptr = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr != nullptr);
+    assert(std::ranges::equal(arr_ptr_->shape(), array_ptr->shape()));
+
+    assert(index == 0);
+    arr_ptr_ = array_ptr;
+}
+
 void SoftMaxNode::revert(State& state) const {
     auto node_data = data_ptr_<SoftMaxNodeStateData>(state);
     node_data->revert();
diff --git a/dwave/optimization/src/nodes/sorting.cpp b/dwave/optimization/src/nodes/sorting.cpp
index 28131562..c194a708 100644
--- a/dwave/optimization/src/nodes/sorting.cpp
+++ b/dwave/optimization/src/nodes/sorting.cpp
@@ -112,6 +112,17 @@ void ArgSortNode::propagate(State& state) const {
     );
 }
 
+void ArgSortNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    ArrayNode* array_ptr = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr != nullptr);
+    assert(std::ranges::equal(arr_ptr_->shape(), array_ptr->shape()));
+
+    assert(index == 0);
+    arr_ptr_ = array_ptr;
+}
+
 void ArgSortNode::revert(State& state) const {
     auto node_data = data_ptr_<ArgSortNodeData>(state);
 
diff --git a/dwave/optimization/src/nodes/statistics.cpp b/dwave/optimization/src/nodes/statistics.cpp
index 341fdc21..7157447a 100644
--- a/dwave/optimization/src/nodes/statistics.cpp
+++ b/dwave/optimization/src/nodes/statistics.cpp
@@ -38,7 +38,7 @@ std::pair<double, double> calculate_values_minmax_(const Array* arr_ptr) {
 }
 
 MeanNode::MeanNode(ArrayNode* arr_ptr) :
-    ScalarOutputMixin<ArrayNode, true>(),
+    ScalarOutputMixin<EqualityMixin<ArrayNode>, true>(),
     arr_ptr_(arr_ptr),
     minmax_(calculate_values_minmax_(arr_ptr_)) {
     add_predecessor_(arr_ptr);
@@ -95,4 +95,16 @@ void MeanNode::propagate(State& state) const {
     }
     set_state(state, sum / static_cast<double>(state_size));
 }
+
+void MeanNode::replace_predecessor_(ssize_t index, Node* node_ptr) {
+    Node::replace_predecessor_(index, node_ptr);
+
+    ArrayNode* array_ptr = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr != nullptr);
+    assert(std::ranges::equal(arr_ptr_->shape(), array_ptr->shape()));
+
+    assert(index == 0);
+    arr_ptr_ = array_ptr;
+}
+
 }  // namespace dwave::optimization
diff --git a/dwave/optimization/src/nodes/unaryop.cpp b/dwave/optimization/src/nodes/unaryop.cpp
index 5aaae206..6cfd3980 100644
--- a/dwave/optimization/src/nodes/unaryop.cpp
+++ b/dwave/optimization/src/nodes/unaryop.cpp
@@ -185,6 +185,19 @@ std::span<const Update> UnaryOpNode<UnaryOp>::diff(const State& state) const {
     return data_ptr_<ArrayNodeStateData>(state)->diff();
 }
 
+template <class UnaryOp>
+bool UnaryOpNode<UnaryOp>::equal_to(const Node& rhs) const {
+    const auto* rhs_ptr = dynamic_cast<const UnaryOpNode*>(&rhs);
+    if (rhs_ptr == nullptr) return false;  // not same type so not equal
+    return this->equal_to(*rhs_ptr);
+}
+
+template <class UnaryOp>
+bool UnaryOpNode<UnaryOp>::equal_to(const UnaryOpNode& rhs) const {
+    // If we're the same type, then we just need to check we have the same predecessor
+    return this->array_ptr_ == rhs.array_ptr_;
+}
+
 template <class UnaryOp>
 void UnaryOpNode<UnaryOp>::initialize_state(State& state) const {
     std::vector<double> values;
@@ -232,6 +245,13 @@ void UnaryOpNode<UnaryOp>::propagate(State& state) const {
     if (node_data->diff().size()) Node::propagate(state);
 }
 
+template <class UnaryOp>
+void UnaryOpNode<UnaryOp>::replace_predecessor_(ssize_t previous_index, Node* node_ptr) {
+    Node::replace_predecessor_(previous_index, node_ptr);
+    array_ptr_ = dynamic_cast<ArrayNode*>(node_ptr);
+    assert(array_ptr_ != nullptr);
+}
+
 template <class UnaryOp>
 void UnaryOpNode<UnaryOp>::revert(State& state) const {
     data_ptr_<ArrayNodeStateData>(state)->revert();
diff --git a/dwave/optimization/symbols/accumulate_zip.pyx b/dwave/optimization/symbols/accumulate_zip.pyx
index f9d80139..0319a8a3 100644
--- a/dwave/optimization/symbols/accumulate_zip.pyx
+++ b/dwave/optimization/symbols/accumulate_zip.pyx
@@ -191,10 +191,10 @@ cdef class AccumulateZip(ArraySymbol):
         return AccumulateZip(expression, operands, initial=initial)
 
     def initial(self):
-        if holds_alternative["ArrayNode*"](self.ptr.initial):
-            return symbol_from_ptr(self.model, get["ArrayNode*"](self.ptr.initial))
+        if holds_alternative["ArrayNode*"](self.ptr.initial()):
+            return symbol_from_ptr(self.model, get["ArrayNode*"](self.ptr.initial()))
         else:
-            return get[double](self.ptr.initial)
+            return get[double](self.ptr.initial())
 
     def _into_zipfile(self, zf, directory):
         """Store a AccumulateZip symbol as a compressed file.
@@ -224,13 +224,13 @@ cdef class AccumulateZip(ArraySymbol):
 
         # Now save information about the initial state
         encoder = json.JSONEncoder(separators=(',', ':'))
-        if holds_alternative["ArrayNode*"](self.ptr.initial):
+        if holds_alternative["ArrayNode*"](self.ptr.initial()):
             initial_info = dict(
                 type="node",
-                value=get["ArrayNode*"](self.ptr.initial).topological_index()
+                value=get["ArrayNode*"](self.ptr.initial()).topological_index()
             )
         else:
-            initial_info = dict(type="double", value=get[double](self.ptr.initial))
+            initial_info = dict(type="double", value=get[double](self.ptr.initial()))
         zf.writestr(directory + "initial.json", encoder.encode(initial_info))
 
     cdef AccumulateZipNode* ptr
diff --git a/releasenotes/notes/feature-redundant-node-removal-9e6335836d0f5603.yaml b/releasenotes/notes/feature-redundant-node-removal-9e6335836d0f5603.yaml
new file mode 100644
index 00000000..2582a1ea
--- /dev/null
+++ b/releasenotes/notes/feature-redundant-node-removal-9e6335836d0f5603.yaml
@@ -0,0 +1,8 @@
+---
+features:
+  - Add ``MatrixMultiplyNode::operands()`` method.
+  - Add ``Graph::remove_redundant_nodes()`` method. See `#563 <https://github.com/dwavesystems/dwave-optimization/issues/563>`_.
+  - Add ``Model.remove_redundant_symbols()`` method. See `#563 <https://github.com/dwavesystems/dwave-optimization/issues/563>`_.
+  - Add ``Node::take_successors()`` method.
+upgrade:
+  - Remove ``AccumulateZipNode::initial`` attribute and replace it with ``AccumulateZipNode::initial()`` method.
diff --git a/tests/cpp/nodes/indexing/test_advanced.cpp b/tests/cpp/nodes/indexing/test_advanced.cpp
index c01db7bd..f7609ffb 100644
--- a/tests/cpp/nodes/indexing/test_advanced.cpp
+++ b/tests/cpp/nodes/indexing/test_advanced.cpp
@@ -1758,6 +1758,68 @@ TEST_CASE("AdvancedIndexingNode") {
             );
         }
     }
+
+    SECTION("equality") {
+        std::vector<double> values = {0, 1, 2, 3, 4, 5, 6, 7};
+        auto* arr0_ptr = graph.emplace_node<ConstantNode>(values, std::array<ssize_t, 3>{2, 2, 2});
+        auto* arr1_ptr = graph.emplace_node<ConstantNode>(values, std::array<ssize_t, 3>{2, 2, 2});
+
+        auto* index0_ptr = graph.emplace_node<IntegerNode>(std::array<ssize_t, 0>{}, 0, 1);
+        auto* index1_ptr = graph.emplace_node<IntegerNode>(std::array<ssize_t, 0>{}, 0, 1);
+        auto* index2_ptr = graph.emplace_node<IntegerNode>(std::array<ssize_t, 0>{}, 0, 1);
+
+        Node* a_ptr =
+            graph.emplace_node<AdvancedIndexingNode>(arr0_ptr, index0_ptr, Slice(), index1_ptr);
+        Node* b_ptr =
+            graph.emplace_node<AdvancedIndexingNode>(arr0_ptr, index0_ptr, Slice(), index1_ptr);
+        Node* c_ptr =
+            graph.emplace_node<AdvancedIndexingNode>(arr0_ptr, index0_ptr, index1_ptr, Slice());
+        Node* d_ptr =
+            graph.emplace_node<AdvancedIndexingNode>(arr1_ptr, index0_ptr, Slice(), index1_ptr);
+        Node* e_ptr =
+            graph.emplace_node<AdvancedIndexingNode>(arr0_ptr, index0_ptr, Slice(), index2_ptr);
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+
+        CHECK(not a_ptr->equal_to(*arr0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+        CHECK(not a_ptr->equal_to(*d_ptr));
+        CHECK(not a_ptr->equal_to(*e_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        std::vector<double> values0 = {0, 1, 2, 3, 4, 5, 6, 7};
+        auto* arr0_ptr = graph.emplace_node<ConstantNode>(values0, std::array<ssize_t, 3>{2, 2, 2});
+        std::vector<double> values1 = {8, 9, 10, 11, 12, 13, 14, 15};
+        auto* arr1_ptr = graph.emplace_node<ConstantNode>(values1, std::array<ssize_t, 3>{2, 2, 2});
+
+        auto* index0_ptr = graph.emplace_node<IntegerNode>(std::array<ssize_t, 0>{}, 0, 1);
+        auto* index1_ptr = graph.emplace_node<IntegerNode>(std::array<ssize_t, 0>{}, 0, 1);
+        auto* index2_ptr = graph.emplace_node<IntegerNode>(std::array<ssize_t, 0>{}, 0, 1);
+        auto* index3_ptr = graph.emplace_node<IntegerNode>(std::array<ssize_t, 0>{}, 0, 1);
+
+        auto* adv_ptr =
+            graph.emplace_node<AdvancedIndexingNode>(arr0_ptr, index0_ptr, Slice(), index1_ptr);
+
+        arr1_ptr->take_successors(*arr0_ptr);
+        index2_ptr->take_successors(*index0_ptr);
+        index3_ptr->take_successors(*index1_ptr);
+
+        CHECK_THAT(
+            adv_ptr->predecessors(),
+            RangeEquals(std::array<Node*, 3>{arr1_ptr, index2_ptr, index3_ptr})
+        );
+
+        auto state = graph.empty_state();
+        index0_ptr->initialize_state(state, {0});
+        index1_ptr->initialize_state(state, {0});
+        index2_ptr->initialize_state(state, {1});
+        index3_ptr->initialize_state(state, {1});
+        graph.initialize_state(state);
+
+        CHECK_THAT(adv_ptr->view(state), RangeEquals({13, 15}));
+    }
 }
 
 }  // namespace dwave::optimization
diff --git a/tests/cpp/nodes/indexing/test_basic.cpp b/tests/cpp/nodes/indexing/test_basic.cpp
index ae8d9f13..f4b9178d 100644
--- a/tests/cpp/nodes/indexing/test_basic.cpp
+++ b/tests/cpp/nodes/indexing/test_basic.cpp
@@ -2458,6 +2458,61 @@ TEST_CASE("BasicIndexingNode") {
         }
     }
 
+    SECTION("equality") {
+        auto* x0_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{5, 5});
+        auto* x1_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{5, 5});
+
+        Node* a_ptr = graph.emplace_node<BasicIndexingNode>(x0_ptr, 0, 1);
+        Node* b_ptr = graph.emplace_node<BasicIndexingNode>(x0_ptr, 0, 1);
+        Node* c_ptr = graph.emplace_node<BasicIndexingNode>(x1_ptr, 0, 1);
+        Node* d_ptr = graph.emplace_node<BasicIndexingNode>(x0_ptr, 1, 1);
+        Node* e_ptr = graph.emplace_node<BasicIndexingNode>(x0_ptr, Slice(0, 3), Slice(1, 3));
+        Node* f_ptr = graph.emplace_node<BasicIndexingNode>(x0_ptr, Slice(0, 3), Slice(1, 3));
+        Node* g_ptr = graph.emplace_node<BasicIndexingNode>(x0_ptr, Slice(0, 2), Slice(1, 3));
+        Node* h_ptr = graph.emplace_node<BasicIndexingNode>(x0_ptr, Slice(0, 2), Slice(1, 10));
+        Node* i_ptr = graph.emplace_node<BasicIndexingNode>(x0_ptr, Slice(0, 2), Slice(1, 100));
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+        CHECK(e_ptr->equal_to(*f_ptr));
+        CHECK(h_ptr->equal_to(*i_ptr));
+
+        CHECK(not a_ptr->equal_to(*x0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+        CHECK(not a_ptr->equal_to(*d_ptr));
+        CHECK(not e_ptr->equal_to(*g_ptr));
+
+        auto* y0_ptr = graph.emplace_node<SetNode>(5);
+
+        Node* j_ptr = graph.emplace_node<BasicIndexingNode>(y0_ptr, Slice(-3, 1000));
+        Node* k_ptr = graph.emplace_node<BasicIndexingNode>(y0_ptr, Slice(-2, 1000));
+
+        CHECK(not j_ptr->equal_to(*k_ptr));
+
+        // dev note: We don't currently support this case but could in the future, see
+        // note in .cpp file
+        // Node* l_ptr = graph.emplace_node<BasicIndexingNode>(y0_ptr, Slice(-3, 10000));
+        // CHECK(j_ptr->equal_to(*l_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto* x0_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{2, 2});
+        auto* x1_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{2, 2});
+
+        auto* basic_ptr = graph.emplace_node<BasicIndexingNode>(x0_ptr, 0, 1);
+
+        x1_ptr->take_successors(*x0_ptr);
+
+        CHECK_THAT(basic_ptr->predecessors(), RangeEquals({x1_ptr}));
+
+        auto state = graph.empty_state();
+        x0_ptr->initialize_state(state, {0, 1, 2, 3});
+        x1_ptr->initialize_state(state, {4, 5, 6, 7});
+        graph.initialize_state(state);
+
+        CHECK_THAT(basic_ptr->view(state), RangeEquals({5}));
+    }
+
     // todo: slicing on multidimensional dynamic array, when we have one to test...
 }
 
diff --git a/tests/cpp/nodes/indexing/test_permutation.cpp b/tests/cpp/nodes/indexing/test_permutation.cpp
index 62dec181..9c19c938 100644
--- a/tests/cpp/nodes/indexing/test_permutation.cpp
+++ b/tests/cpp/nodes/indexing/test_permutation.cpp
@@ -16,6 +16,7 @@
 #include <catch2/generators/catch_generators.hpp>
 #include <catch2/matchers/catch_matchers_all.hpp>
 
+#include "catch2/matchers/catch_matchers.hpp"
 #include "dwave-optimization/graph.hpp"
 #include "dwave-optimization/nodes/collections.hpp"
 #include "dwave-optimization/nodes/constants.hpp"
@@ -174,6 +175,54 @@ TEST_CASE("PermutationNode") {
             }
         }
     }
+
+    SECTION("equality") {
+        std::vector<double> values = {0, 1, 2, 3};
+        auto* arr0_ptr = graph.emplace_node<ConstantNode>(values, std::array<ssize_t, 2>{2, 2});
+        auto* arr1_ptr = graph.emplace_node<ConstantNode>(values, std::array<ssize_t, 2>{2, 2});
+
+        auto* order0_ptr = graph.emplace_node<ListNode>(2);
+        auto* order1_ptr = graph.emplace_node<ListNode>(2);
+
+        Node* a_ptr = graph.emplace_node<PermutationNode>(arr0_ptr, order0_ptr);
+        Node* b_ptr = graph.emplace_node<PermutationNode>(arr0_ptr, order0_ptr);
+        Node* c_ptr = graph.emplace_node<PermutationNode>(arr1_ptr, order0_ptr);
+        Node* d_ptr = graph.emplace_node<PermutationNode>(arr0_ptr, order1_ptr);
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+        CHECK(not a_ptr->equal_to(*arr0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+        CHECK(not a_ptr->equal_to(*d_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        std::vector<double> values0 = {0, 1, 2, 3};
+        auto* arr0_ptr = graph.emplace_node<ConstantNode>(values0, std::array<ssize_t, 2>{2, 2});
+        std::vector<double> values1 = {5, 6, 7, 8};
+        auto* arr1_ptr = graph.emplace_node<ConstantNode>(values1, std::array<ssize_t, 2>{2, 2});
+
+        auto* order0_ptr = graph.emplace_node<ListNode>(2);
+        auto* order1_ptr = graph.emplace_node<ListNode>(2);
+
+        auto* permutation_ptr = graph.emplace_node<PermutationNode>(arr0_ptr, order0_ptr);
+
+        arr1_ptr->take_successors(*arr0_ptr);
+        CHECK_THAT(
+            permutation_ptr->predecessors(), RangeEquals(std::array<Node*, 2>{arr1_ptr, order0_ptr})
+        );
+        order1_ptr->take_successors(*order0_ptr);
+        CHECK_THAT(
+            permutation_ptr->predecessors(), RangeEquals(std::array<Node*, 2>{arr1_ptr, order1_ptr})
+        );
+
+        auto state = graph.empty_state();
+        order0_ptr->initialize_state(state, {0, 1});
+        order1_ptr->initialize_state(state, {1, 0});
+        graph.initialize_state(state);
+
+        CHECK_THAT(permutation_ptr->view(state), RangeEquals({8, 7, 6, 5}));
+    }
 }
 
 }  // namespace dwave::optimization
diff --git a/tests/cpp/nodes/test_binaryop.cpp b/tests/cpp/nodes/test_binaryop.cpp
index d0e59c92..fcd9f98c 100644
--- a/tests/cpp/nodes/test_binaryop.cpp
+++ b/tests/cpp/nodes/test_binaryop.cpp
@@ -383,6 +383,70 @@ TEMPLATE_TEST_CASE(
             }
         }
     }
+
+    SECTION("equality") {
+        auto* x_ptr = graph.emplace_node<IntegerNode>();
+        auto* y_ptr = graph.emplace_node<IntegerNode>();
+        auto* z_ptr = graph.emplace_node<IntegerNode>();
+
+        Node* a_ptr = graph.emplace_node<BinaryOpNode<TestType>>(x_ptr, y_ptr);
+        Node* b_ptr = graph.emplace_node<BinaryOpNode<TestType>>(x_ptr, y_ptr);
+        Node* c_ptr = graph.emplace_node<BinaryOpNode<TestType>>(y_ptr, x_ptr);  // commutative
+        Node* d_ptr = graph.emplace_node<BinaryOpNode<TestType>>(x_ptr, z_ptr);  // not equal
+
+        CHECK(a_ptr->equal_to(*b_ptr));
+        CHECK(b_ptr->equal_to(*a_ptr));
+
+        CHECK(not a_ptr->equal_to(*d_ptr));
+        CHECK(not d_ptr->equal_to(*a_ptr));
+
+        // Once we switch to ufuncs
+        // https://github.com/dwavesystems/dwave-optimization/pull/412
+        // we can use ::communcative. For now we hardcode it.
+
+        if constexpr (
+            std::same_as<BinaryOpNode<TestType>, AddNode> or
+            std::same_as<BinaryOpNode<TestType>, AndNode> or
+            std::same_as<BinaryOpNode<TestType>, EqualNode> or
+            std::same_as<BinaryOpNode<TestType>, MaximumNode> or
+            std::same_as<BinaryOpNode<TestType>, MinimumNode> or
+            std::same_as<BinaryOpNode<TestType>, MultiplyNode> or
+            std::same_as<BinaryOpNode<TestType>, OrNode> or
+            std::same_as<BinaryOpNode<TestType>, XorNode>
+        ) {
+            // commutative
+            CHECK(a_ptr->equal_to(*c_ptr));
+            CHECK(c_ptr->equal_to(*a_ptr));
+        } else {
+            // not commutative
+            static_assert(
+                std::same_as<BinaryOpNode<TestType>, DivideNode> or
+                std::same_as<BinaryOpNode<TestType>, LessEqualNode> or
+                std::same_as<BinaryOpNode<TestType>, ModulusNode> or
+                std::same_as<BinaryOpNode<TestType>, SafeDivideNode> or
+                std::same_as<BinaryOpNode<TestType>, SubtractNode>
+            );
+
+            CHECK(not a_ptr->equal_to(*c_ptr));
+            CHECK(not c_ptr->equal_to(*a_ptr));
+        }
+    }
+
+    SECTION("predecessor replacement") {
+        auto* x_ptr = graph.emplace_node<IntegerNode>();
+        auto* y_ptr = graph.emplace_node<IntegerNode>();
+        auto* z_ptr = graph.emplace_node<IntegerNode>();
+
+        auto* a_ptr = graph.emplace_node<BinaryOpNode<TestType>>(x_ptr, y_ptr);
+
+        CHECK_THAT(a_ptr->predecessors(), RangeEquals({x_ptr, y_ptr}));
+        CHECK_THAT(a_ptr->operands(), RangeEquals({x_ptr, y_ptr}));
+
+        z_ptr->take_successors(*y_ptr);
+
+        CHECK_THAT(a_ptr->predecessors(), RangeEquals({x_ptr, z_ptr}));
+        CHECK_THAT(a_ptr->operands(), RangeEquals({x_ptr, z_ptr}));
+    }
 }
 
 TEST_CASE("BinaryOpNode - LessEqualNode") {
diff --git a/tests/cpp/nodes/test_collections.cpp b/tests/cpp/nodes/test_collections.cpp
index 4090de04..b6e45baf 100644
--- a/tests/cpp/nodes/test_collections.cpp
+++ b/tests/cpp/nodes/test_collections.cpp
@@ -49,6 +49,17 @@ TEST_CASE("DisjointBitSetsNode") {
                 graph.emplace_node<ArrayValidationNode>(sets.at(i));
             }
 
+            THEN("node equality works as expected") {
+                CHECK(ptr->equal_to(*ptr));
+
+                CHECK(sets[0]->equal_to(*sets[0]));
+                CHECK(not sets[0]->equal_to(*sets[1]));
+                CHECK(static_cast<Node*>(sets[0])->equal_to(*sets[0]));
+                CHECK(not static_cast<Node*>(sets[0])->equal_to(*sets[1]));
+                CHECK(sets[0]->equal_to(*static_cast<Node*>(sets[0])));
+                CHECK(not sets[0]->equal_to(*static_cast<Node*>(sets[1])));
+            }
+
             THEN("We shouldn't be able to add any more successors") {
                 CHECK_THROWS(graph.emplace_node<DisjointBitSetNode>(ptr));
             }
diff --git a/tests/cpp/nodes/test_creation.cpp b/tests/cpp/nodes/test_creation.cpp
index e0db6ae9..754fb096 100644
--- a/tests/cpp/nodes/test_creation.cpp
+++ b/tests/cpp/nodes/test_creation.cpp
@@ -16,6 +16,7 @@
 #include <catch2/matchers/catch_matchers_all.hpp>
 
 #include "dwave-optimization/nodes/collections.hpp"
+#include "dwave-optimization/nodes/constants.hpp"
 #include "dwave-optimization/nodes/creation.hpp"
 #include "dwave-optimization/nodes/manipulation.hpp"
 #include "dwave-optimization/nodes/numbers.hpp"
@@ -566,6 +567,55 @@ TEST_CASE("ARangeNode") {
             CHECK_THAT(arange_ptr->view(state), RangeEquals({17, 14, 11}));
         }
     }
+
+    SECTION("equality") {
+        auto graph = Graph();
+
+        auto* x_ptr = graph.emplace_node<IntegerNode>();
+        auto* y_ptr = graph.emplace_node<IntegerNode>();
+
+        Node* a_ptr = graph.emplace_node<ARangeNode>(0, 10, 1);
+        Node* b_ptr = graph.emplace_node<ARangeNode>(0, 10, 1);
+        Node* c_ptr = graph.emplace_node<ARangeNode>(1, 10, 2);
+
+        CHECK(a_ptr->equal_to(*b_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+
+        Node* d_ptr = graph.emplace_node<ARangeNode>(0, x_ptr);
+        Node* e_ptr = graph.emplace_node<ARangeNode>(0, x_ptr);
+        Node* f_ptr = graph.emplace_node<ARangeNode>(0, y_ptr);
+
+        CHECK(d_ptr->equal_to(*e_ptr));
+        CHECK(not d_ptr->equal_to(*b_ptr));
+        CHECK(not d_ptr->equal_to(*f_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto graph = Graph();
+
+        auto* x_ptr = graph.emplace_node<ConstantNode>(1);
+        auto* y_ptr = graph.emplace_node<ConstantNode>(1);
+
+        SECTION("range(x, x, x); y.take_sucessors(x)") {
+            auto* arange_ptr = graph.emplace_node<ARangeNode>(x_ptr, x_ptr, x_ptr);
+            y_ptr->take_successors(*x_ptr);
+
+            CHECK_THAT(arange_ptr->predecessors(), RangeEquals({y_ptr, y_ptr, y_ptr}));
+            CHECK(std::get<const Array*>(arange_ptr->start()) == y_ptr);
+            CHECK(std::get<const Array*>(arange_ptr->stop()) == y_ptr);
+            CHECK(std::get<const Array*>(arange_ptr->step()) == y_ptr);
+        }
+
+        SECTION("range(x, 1, x); y.take_sucessors(x)") {
+            auto* arange_ptr = graph.emplace_node<ARangeNode>(x_ptr, 1, x_ptr);
+            y_ptr->take_successors(*x_ptr);
+
+            CHECK_THAT(arange_ptr->predecessors(), RangeEquals({y_ptr, y_ptr}));
+            CHECK(std::get<const Array*>(arange_ptr->start()) == y_ptr);
+            CHECK(std::get<ssize_t>(arange_ptr->stop()) == 1);
+            CHECK(std::get<const Array*>(arange_ptr->step()) == y_ptr);
+        }
+    }
 }
 
 }  // namespace dwave::optimization
diff --git a/tests/cpp/nodes/test_flow.cpp b/tests/cpp/nodes/test_flow.cpp
index 37a98395..bebe3d75 100644
--- a/tests/cpp/nodes/test_flow.cpp
+++ b/tests/cpp/nodes/test_flow.cpp
@@ -201,6 +201,49 @@ TEST_CASE("ExtractNode") {
             }
         }
     }
+
+    SECTION("equality") {
+        auto* c0_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, 0, 5);
+        auto* x0_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, -10, 10);
+
+        auto* c1_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, 0, 5);
+        auto* x1_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, -10, 10);
+
+        Node* a_ptr = graph.emplace_node<ExtractNode>(c0_ptr, x0_ptr);
+        Node* b_ptr = graph.emplace_node<ExtractNode>(c0_ptr, x0_ptr);
+        Node* c_ptr = graph.emplace_node<ExtractNode>(c1_ptr, x0_ptr);
+        Node* d_ptr = graph.emplace_node<ExtractNode>(c0_ptr, x1_ptr);
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+        CHECK(not a_ptr->equal_to(*x0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+        CHECK(not a_ptr->equal_to(*d_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto* c0_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, 0, 5);
+        auto* x0_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, -10, 10);
+
+        auto* c1_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, 0, 5);
+        auto* x1_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, -10, 10);
+
+        auto* extract_ptr = graph.emplace_node<ExtractNode>(c0_ptr, x0_ptr);
+
+        c1_ptr->take_successors(*c0_ptr);
+        x1_ptr->take_successors(*x0_ptr);
+
+        CHECK_THAT(extract_ptr->predecessors(), RangeEquals({c1_ptr, x1_ptr}));
+
+        auto state = graph.empty_state();
+        c0_ptr->initialize_state(state, {0});
+        x0_ptr->initialize_state(state, {1});
+        c1_ptr->initialize_state(state, {1});
+        x1_ptr->initialize_state(state, {5});
+        graph.initialize_state(state);
+
+        CHECK_THAT(extract_ptr->view(state), RangeEquals({5}));
+    }
 }
 
 TEST_CASE("WhereNode") {
@@ -555,6 +598,58 @@ TEST_CASE("WhereNode") {
             CHECK_THROWS_AS(WhereNode(condition_ptr, x_ptr, y_ptr), std::invalid_argument);
         }
     }
+
+    SECTION("equality") {
+        auto* c0_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, 0, 5);
+        auto* x0_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, -10, 10);
+        auto* y0_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, 0, 10);
+
+        auto* c1_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, 0, 5);
+        auto* x1_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, -10, 10);
+        auto* y1_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, 0, 10);
+
+        Node* a_ptr = graph.emplace_node<WhereNode>(c0_ptr, x0_ptr, y0_ptr);
+        Node* b_ptr = graph.emplace_node<WhereNode>(c0_ptr, x0_ptr, y0_ptr);
+        Node* c_ptr = graph.emplace_node<WhereNode>(c1_ptr, x0_ptr, y0_ptr);
+        Node* d_ptr = graph.emplace_node<WhereNode>(c0_ptr, x1_ptr, y0_ptr);
+        Node* e_ptr = graph.emplace_node<WhereNode>(c0_ptr, x0_ptr, y1_ptr);
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+        CHECK(not a_ptr->equal_to(*x0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+        CHECK(not a_ptr->equal_to(*d_ptr));
+        CHECK(not a_ptr->equal_to(*e_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto* c0_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, 0, 5);
+        auto* x0_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, -10, 10);
+        auto* y0_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, 0, 10);
+
+        auto* c1_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, 0, 5);
+        auto* x1_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, -10, 10);
+        auto* y1_ptr = graph.emplace_node<IntegerNode>(std::vector<ssize_t>{}, 0, 10);
+
+        auto* where_ptr = graph.emplace_node<WhereNode>(c0_ptr, x0_ptr, y0_ptr);
+
+        c1_ptr->take_successors(*c0_ptr);
+        x1_ptr->take_successors(*x0_ptr);
+        y1_ptr->take_successors(*y0_ptr);
+
+        CHECK_THAT(where_ptr->predecessors(), RangeEquals({c1_ptr, x1_ptr, y1_ptr}));
+
+        auto state = graph.empty_state();
+        c0_ptr->initialize_state(state, {0});
+        x0_ptr->initialize_state(state, {1});
+        y0_ptr->initialize_state(state, {2});
+        c1_ptr->initialize_state(state, {1});
+        x1_ptr->initialize_state(state, {3});
+        y1_ptr->initialize_state(state, {4});
+        graph.initialize_state(state);
+
+        CHECK_THAT(where_ptr->view(state), RangeEquals({3}));
+    }
 }
 
 }  // namespace dwave::optimization
diff --git a/tests/cpp/nodes/test_interpolation.cpp b/tests/cpp/nodes/test_interpolation.cpp
index 51768d4b..f07da232 100644
--- a/tests/cpp/nodes/test_interpolation.cpp
+++ b/tests/cpp/nodes/test_interpolation.cpp
@@ -13,12 +13,16 @@
 //    limitations under the License.
 
 #include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
 
+#include "catch2/matchers/catch_matchers_range_equals.hpp"
 #include "dwave-optimization/graph.hpp"
 #include "dwave-optimization/nodes/constants.hpp"
 #include "dwave-optimization/nodes/interpolation.hpp"
 #include "dwave-optimization/nodes/numbers.hpp"
 
+using Catch::Matchers::RangeEquals;
+
 namespace dwave::optimization {
 
 TEST_CASE("BSpline") {
@@ -114,6 +118,50 @@ TEST_CASE("BSpline") {
             }
         }
     }
+
+    SECTION("equality") {
+        auto graph = Graph();
+
+        auto* x0_ptr = graph.emplace_node<ConstantNode>(std::vector{2.5});
+        auto* x1_ptr = graph.emplace_node<ConstantNode>(std::vector{3.5});
+
+        int k = 2;
+        std::vector<double> t0 = {0, 1, 2, 3, 4, 5, 6};
+        std::vector<double> t1 = {1, 1, 2, 3, 4, 5, 6};
+        std::vector<double> c0 = {-1, 2, 0, -1};
+        std::vector<double> c1 = {1, 2, 0, -1};
+
+        Node* a_ptr = graph.emplace_node<BSplineNode>(x0_ptr, k, t0, c0);
+        Node* b_ptr = graph.emplace_node<BSplineNode>(x0_ptr, k, t0, c0);
+        Node* c_ptr = graph.emplace_node<BSplineNode>(x1_ptr, k, t1, c0);
+        Node* d_ptr = graph.emplace_node<BSplineNode>(x0_ptr, k, t0, c1);
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+        CHECK(not a_ptr->equal_to(*x0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+        CHECK(not a_ptr->equal_to(*d_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto graph = Graph();
+
+        auto* x0_ptr = graph.emplace_node<ConstantNode>(std::vector{2.5});
+        auto* x1_ptr = graph.emplace_node<ConstantNode>(std::vector{3.5});
+
+        int k = 2;
+        std::vector<double> t = {0, 1, 2, 3, 4, 5, 6};
+        std::vector<double> c = {-1, 2, 0, -1};
+
+        auto* bspline_ptr = graph.emplace_node<BSplineNode>(x0_ptr, k, t, c);
+
+        x1_ptr->take_successors(*x0_ptr);
+
+        CHECK_THAT(bspline_ptr->predecessors(), RangeEquals({x1_ptr}));
+
+        auto state = graph.initialize_state();
+        CHECK_THAT(bspline_ptr->view(state), RangeEquals({0.125}));
+    }
 }
 
 }  // namespace dwave::optimization
diff --git a/tests/cpp/nodes/test_lambda.cpp b/tests/cpp/nodes/test_lambda.cpp
index 0382bdb1..83641e9c 100644
--- a/tests/cpp/nodes/test_lambda.cpp
+++ b/tests/cpp/nodes/test_lambda.cpp
@@ -12,6 +12,9 @@
 //    See the License for the specific language governing permissions and
 //    limitations under the License.
 
+#include <memory>
+#include <vector>
+
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_all.hpp>
 
@@ -454,6 +457,105 @@ TEST_CASE("AccumulateZipNode") {
             }
         }
     }
+
+    SECTION("equality") {
+        auto* i0_ptr = graph.emplace_node<ConstantNode>(std::vector<double>{0, 1, 2, 2});
+        auto* i1_ptr = graph.emplace_node<ConstantNode>(std::vector<double>{0, 1, 2, 2});
+
+        auto* j0_ptr = graph.emplace_node<ConstantNode>(std::vector<double>{1, 2, 4, 3});
+        auto* j1_ptr = graph.emplace_node<ConstantNode>(std::vector<double>{1, 2, 4, 3});
+
+        // x0 * x1 + x2
+        auto make_expression = []() -> std::shared_ptr<Graph> {
+            auto expression = Graph();
+
+            std::vector<InputNode*> inputs = {
+                expression.emplace_node<InputNode>(InputNode::unbounded_scalar()),
+                expression.emplace_node<InputNode>(InputNode::unbounded_scalar()),
+                expression.emplace_node<InputNode>(InputNode::unbounded_scalar())
+            };
+            auto output_ptr = expression.emplace_node<AddNode>(
+                expression.emplace_node<MultiplyNode>(inputs[1], inputs[2]), inputs[0]
+            );
+            expression.set_objective(output_ptr);
+            expression.topological_sort();
+
+            return std::make_shared<Graph>(std::move(expression));
+        };
+
+        auto expr0_ptr = make_expression();
+        auto expr1_ptr = make_expression();
+
+        Node* a_ptr = graph.emplace_node<AccumulateZipNode>(
+            expr0_ptr, std::vector<ArrayNode*>{i0_ptr, j0_ptr}, 5.0
+        );
+        Node* b_ptr = graph.emplace_node<AccumulateZipNode>(
+            expr0_ptr, std::vector<ArrayNode*>{i0_ptr, j0_ptr}, 5.0
+        );
+        Node* c_ptr = graph.emplace_node<AccumulateZipNode>(
+            expr1_ptr, std::vector<ArrayNode*>{i0_ptr, j0_ptr}, 5.0
+        );
+        Node* d_ptr = graph.emplace_node<AccumulateZipNode>(
+            expr0_ptr, std::vector<ArrayNode*>{i1_ptr, j0_ptr}, 5.0
+        );
+        Node* e_ptr = graph.emplace_node<AccumulateZipNode>(
+            expr0_ptr, std::vector<ArrayNode*>{i0_ptr, j1_ptr}, 5.0
+        );
+        Node* f_ptr = graph.emplace_node<AccumulateZipNode>(
+            expr0_ptr, std::vector<ArrayNode*>{i0_ptr, j0_ptr}, 4.0
+        );
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+        CHECK(not a_ptr->equal_to(*i0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+        CHECK(not a_ptr->equal_to(*d_ptr));
+        CHECK(not a_ptr->equal_to(*e_ptr));
+        CHECK(not a_ptr->equal_to(*f_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto* i0_ptr = graph.emplace_node<ConstantNode>(std::vector<double>{0, 1, 2, 2});
+        auto* i1_ptr = graph.emplace_node<ConstantNode>(std::vector<double>{1, 2, 2, 0});
+
+        auto* j0_ptr = graph.emplace_node<ConstantNode>(std::vector<double>{1, 2, 4, 3});
+        auto* j1_ptr = graph.emplace_node<ConstantNode>(std::vector<double>{2, 4, 3, 1});
+
+        auto* init0_ptr = graph.emplace_node<ConstantNode>(5);
+        auto* init1_ptr = graph.emplace_node<ConstantNode>(6);
+
+        // x0 * x1 + x2
+        auto make_expression = []() -> std::shared_ptr<Graph> {
+            auto expression = Graph();
+
+            std::vector<InputNode*> inputs = {
+                expression.emplace_node<InputNode>(InputNode::unbounded_scalar()),
+                expression.emplace_node<InputNode>(InputNode::unbounded_scalar()),
+                expression.emplace_node<InputNode>(InputNode::unbounded_scalar())
+            };
+            auto output_ptr = expression.emplace_node<AddNode>(
+                expression.emplace_node<MultiplyNode>(inputs[1], inputs[2]), inputs[0]
+            );
+            expression.set_objective(output_ptr);
+            expression.topological_sort();
+
+            return std::make_shared<Graph>(std::move(expression));
+        };
+
+        auto* acc_ptr = graph.emplace_node<AccumulateZipNode>(
+            make_expression(), std::vector<ArrayNode*>{i0_ptr, j0_ptr}, init0_ptr
+        );
+
+        i1_ptr->take_successors(*i0_ptr);
+        j1_ptr->take_successors(*j0_ptr);
+        init1_ptr->take_successors(*init0_ptr);
+
+        CHECK_THAT(acc_ptr->predecessors(), RangeEquals({init1_ptr, i1_ptr, j1_ptr}));
+
+        auto state = graph.initialize_state();
+
+        CHECK_THAT(acc_ptr->view(state), RangeEquals({8, 16, 22, 22}));
+    }
 }
 
 }  // namespace dwave::optimization
diff --git a/tests/cpp/nodes/test_linear_algebra.cpp b/tests/cpp/nodes/test_linear_algebra.cpp
index 9eef3894..2c9ea065 100644
--- a/tests/cpp/nodes/test_linear_algebra.cpp
+++ b/tests/cpp/nodes/test_linear_algebra.cpp
@@ -547,6 +547,60 @@ TEST_CASE("MatrixMultiplyNode") {
             }
         }
     }
+
+    SECTION("equality") {
+        auto* x0_ptr = graph.emplace_node<ConstantNode>(
+            std::vector{1, 2, 3, 4}, std::vector<ssize_t>{2, 2}
+        );
+        auto* x1_ptr = graph.emplace_node<ConstantNode>(
+            std::vector{6, 5, 4, 3}, std::vector<ssize_t>{2, 2}
+        );
+
+        auto* y0_ptr = graph.emplace_node<ConstantNode>(
+            std::vector{7, 8, 9, 10}, std::vector<ssize_t>{2, 2}
+        );
+        auto* y1_ptr = graph.emplace_node<ConstantNode>(
+            std::vector{12, 11, 10, 9}, std::vector<ssize_t>{2, 2}
+        );
+
+        Node* a_ptr = graph.emplace_node<MatrixMultiplyNode>(x0_ptr, y0_ptr);
+        Node* b_ptr = graph.emplace_node<MatrixMultiplyNode>(x0_ptr, y0_ptr);
+        Node* c_ptr = graph.emplace_node<MatrixMultiplyNode>(x1_ptr, y0_ptr);
+        Node* d_ptr = graph.emplace_node<MatrixMultiplyNode>(x0_ptr, y1_ptr);
+        Node* e_ptr = graph.emplace_node<MatrixMultiplyNode>(y0_ptr, x0_ptr);
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+
+        CHECK(not a_ptr->equal_to(*x0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+        CHECK(not a_ptr->equal_to(*d_ptr));
+        CHECK(not a_ptr->equal_to(*e_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto* x0_ptr = graph.emplace_node<ConstantNode>(
+            std::vector{1, 2, 3, 4, 5, 6}, std::vector<ssize_t>{2, 3}
+        );
+        auto* x1_ptr = graph.emplace_node<ConstantNode>(
+            std::vector{6, 5, 4, 3, 2, 1}, std::vector<ssize_t>{2, 3}
+        );
+
+        auto* y0_ptr = graph.emplace_node<ConstantNode>(
+            std::vector{7, 8, 9, 10, 11, 12}, std::vector<ssize_t>{3, 2}
+        );
+        auto* y1_ptr = graph.emplace_node<ConstantNode>(
+            std::vector{12, 11, 10, 9, 8, 7}, std::vector<ssize_t>{3, 2}
+        );
+
+        auto* matmul_ptr = graph.emplace_node<MatrixMultiplyNode>(x0_ptr, y0_ptr);
+
+        x1_ptr->take_successors(*x0_ptr);
+        y1_ptr->take_successors(*y0_ptr);
+
+        CHECK_THAT(matmul_ptr->predecessors(), RangeEquals({x1_ptr, y1_ptr}));
+        CHECK_THAT(matmul_ptr->operands(), RangeEquals({x1_ptr, y1_ptr}));
+    }
 }
 
 }  // namespace dwave::optimization
diff --git a/tests/cpp/nodes/test_lp.cpp b/tests/cpp/nodes/test_lp.cpp
index 08e4020b..f28850c8 100644
--- a/tests/cpp/nodes/test_lp.cpp
+++ b/tests/cpp/nodes/test_lp.cpp
@@ -16,8 +16,12 @@
 
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <catch2/matchers/catch_matchers_range_equals.hpp>
 
-#include "dwave-optimization/nodes.hpp"
+#include "dwave-optimization/nodes/constants.hpp"
+#include "dwave-optimization/nodes/lp.hpp"
+#include "dwave-optimization/nodes/numbers.hpp"
+#include "dwave-optimization/nodes/testing.hpp"
 
 using namespace Catch::Matchers;
 
@@ -738,6 +742,164 @@ TEST_CASE("LinearProgramNode") {
             }
         }
     }
+
+    SECTION("equality") {
+        auto graph = Graph();
+
+        const ssize_t num_variables = 3;
+        const ssize_t num_eq = 2;
+        const ssize_t num_ineq = 2;
+
+        const ssize_t min = IntegerNode::minimum_lower_bound;
+        const ssize_t max = IntegerNode::maximum_upper_bound;
+
+        auto* c0 = graph.emplace_node<IntegerNode>(std::vector{num_variables}, min, max);
+        auto* c1 = graph.emplace_node<IntegerNode>(std::vector{num_variables}, min, max);
+
+        auto* b_lb0 = graph.emplace_node<IntegerNode>(std::vector{num_ineq}, min, max);
+        auto* b_lb1 = graph.emplace_node<IntegerNode>(std::vector{num_ineq}, min, max);
+        auto* A0 = graph.emplace_node<IntegerNode>(std::vector{num_ineq, num_variables}, min, max);
+        auto* A1 = graph.emplace_node<IntegerNode>(std::vector{num_ineq, num_variables}, min, max);
+        auto* b_ub0 = graph.emplace_node<IntegerNode>(std::vector{num_ineq}, min, max);
+        auto* b_ub1 = graph.emplace_node<IntegerNode>(std::vector{num_ineq}, min, max);
+
+        auto* A_eq0 = graph.emplace_node<IntegerNode>(std::vector{num_eq, num_variables}, min, max);
+        auto* A_eq1 = graph.emplace_node<IntegerNode>(std::vector{num_eq, num_variables}, min, max);
+        auto* b_eq0 = graph.emplace_node<IntegerNode>(std::vector{num_eq}, min, max);
+        auto* b_eq1 = graph.emplace_node<IntegerNode>(std::vector{num_eq}, min, max);
+
+        auto* lb0 = graph.emplace_node<IntegerNode>(std::vector{num_variables}, min, max);
+        auto* lb1 = graph.emplace_node<IntegerNode>(std::vector{num_variables}, min, max);
+        auto* ub0 = graph.emplace_node<IntegerNode>(std::vector{num_variables}, min, max);
+        auto* ub1 = graph.emplace_node<IntegerNode>(std::vector{num_variables}, min, max);
+
+        Node* lp0 =
+            graph.emplace_node<LinearProgramNode>(c0, b_lb0, A0, b_ub0, A_eq0, b_eq0, lb0, ub0);
+        Node* lp1 =
+            graph.emplace_node<LinearProgramNode>(c0, b_lb0, A0, b_ub0, A_eq0, b_eq0, lb0, ub0);
+        Node* lp2 =
+            graph.emplace_node<LinearProgramNode>(c1, b_lb0, A0, b_ub0, A_eq0, b_eq0, lb0, ub0);
+        Node* lp3 =
+            graph.emplace_node<LinearProgramNode>(c0, b_lb1, A0, b_ub0, A_eq0, b_eq0, lb0, ub0);
+        Node* lp4 =
+            graph.emplace_node<LinearProgramNode>(c0, b_lb0, A1, b_ub0, A_eq0, b_eq0, lb0, ub0);
+        Node* lp5 =
+            graph.emplace_node<LinearProgramNode>(c0, b_lb0, A0, b_ub1, A_eq0, b_eq0, lb0, ub0);
+        Node* lp6 =
+            graph.emplace_node<LinearProgramNode>(c0, b_lb0, A0, b_ub0, A_eq1, b_eq0, lb0, ub0);
+        Node* lp7 =
+            graph.emplace_node<LinearProgramNode>(c0, b_lb0, A0, b_ub0, A_eq0, b_eq1, lb0, ub0);
+        Node* lp8 =
+            graph.emplace_node<LinearProgramNode>(c0, b_lb0, A0, b_ub0, A_eq0, b_eq0, lb1, ub0);
+        Node* lp9 =
+            graph.emplace_node<LinearProgramNode>(c0, b_lb0, A0, b_ub0, A_eq0, b_eq0, lb0, ub1);
+
+        CHECK(lp0->equal_to(*lp0));
+        CHECK(lp0->equal_to(*lp1));
+
+        CHECK(not lp0->equal_to(*c0));
+        CHECK(not lp0->equal_to(*lp2));
+        CHECK(not lp0->equal_to(*lp3));
+        CHECK(not lp0->equal_to(*lp4));
+        CHECK(not lp0->equal_to(*lp5));
+        CHECK(not lp0->equal_to(*lp6));
+        CHECK(not lp0->equal_to(*lp7));
+        CHECK(not lp0->equal_to(*lp8));
+        CHECK(not lp0->equal_to(*lp9));
+
+        Node* feas0 =
+            graph.emplace_node<LinearProgramFeasibleNode>(dynamic_cast<LinearProgramNode*>(lp0));
+        Node* feas1 =
+            graph.emplace_node<LinearProgramFeasibleNode>(dynamic_cast<LinearProgramNode*>(lp0));
+        Node* feas2 =
+            graph.emplace_node<LinearProgramFeasibleNode>(dynamic_cast<LinearProgramNode*>(lp1));
+
+        CHECK(feas0->equal_to(*feas0));
+        CHECK(feas0->equal_to(*feas1));
+        CHECK(not feas0->equal_to(*feas2));
+        CHECK(not feas0->equal_to(*c0));
+
+        Node* obj0 = graph.emplace_node<LinearProgramObjectiveValueNode>(
+            dynamic_cast<LinearProgramNode*>(lp0)
+        );
+        Node* obj1 = graph.emplace_node<LinearProgramObjectiveValueNode>(
+            dynamic_cast<LinearProgramNode*>(lp0)
+        );
+        Node* obj2 = graph.emplace_node<LinearProgramObjectiveValueNode>(
+            dynamic_cast<LinearProgramNode*>(lp1)
+        );
+
+        CHECK(obj0->equal_to(*obj0));
+        CHECK(obj0->equal_to(*obj1));
+        CHECK(not obj0->equal_to(*obj2));
+        CHECK(not obj0->equal_to(*c0));
+
+        Node* sol0 =
+            graph.emplace_node<LinearProgramSolutionNode>(dynamic_cast<LinearProgramNode*>(lp0));
+        Node* sol1 =
+            graph.emplace_node<LinearProgramSolutionNode>(dynamic_cast<LinearProgramNode*>(lp0));
+        Node* sol2 =
+            graph.emplace_node<LinearProgramSolutionNode>(dynamic_cast<LinearProgramNode*>(lp1));
+
+        CHECK(sol0->equal_to(*sol0));
+        CHECK(sol0->equal_to(*sol1));
+        CHECK(not sol0->equal_to(*sol2));
+        CHECK(not sol0->equal_to(*c0));
+    }
+
+    SECTION("predecessor replacement") {
+        auto graph = Graph();
+
+        const ssize_t num_variables = 3;
+        const ssize_t num_eq = 2;
+        const ssize_t num_ineq = 2;
+
+        const ssize_t min = IntegerNode::minimum_lower_bound;
+        const ssize_t max = IntegerNode::maximum_upper_bound;
+
+        auto* c0 = graph.emplace_node<IntegerNode>(std::vector{num_variables}, min, max);
+        auto* c1 = graph.emplace_node<IntegerNode>(std::vector{num_variables}, min, max);
+
+        auto* b_lb0 = graph.emplace_node<IntegerNode>(std::vector{num_ineq}, min, max);
+        auto* b_lb1 = graph.emplace_node<IntegerNode>(std::vector{num_ineq}, min, max);
+        auto* A0 = graph.emplace_node<IntegerNode>(std::vector{num_ineq, num_variables}, min, max);
+        auto* A1 = graph.emplace_node<IntegerNode>(std::vector{num_ineq, num_variables}, min, max);
+
+        auto* A_eq0 = graph.emplace_node<IntegerNode>(std::vector{num_eq, num_variables}, min, max);
+        auto* A_eq1 = graph.emplace_node<IntegerNode>(std::vector{num_eq, num_variables}, min, max);
+        auto* b_eq0 = graph.emplace_node<IntegerNode>(std::vector{num_eq}, min, max);
+        auto* b_eq1 = graph.emplace_node<IntegerNode>(std::vector{num_eq}, min, max);
+
+        auto* ub0 = graph.emplace_node<IntegerNode>(std::vector{num_variables}, min, max);
+        auto* ub1 = graph.emplace_node<IntegerNode>(std::vector{num_variables}, min, max);
+
+        auto* lp0 = graph.emplace_node<LinearProgramNode>(
+            c0, b_lb0, A0, nullptr, A_eq0, b_eq0, nullptr, ub0
+        );
+        auto* lp1 = graph.emplace_node<LinearProgramNode>(
+            c0, b_lb0, A0, nullptr, A_eq0, b_eq0, nullptr, ub0
+        );
+
+        auto* feas = graph.emplace_node<LinearProgramFeasibleNode>(lp0);
+        auto* obj = graph.emplace_node<LinearProgramObjectiveValueNode>(lp0);
+        auto* sol = graph.emplace_node<LinearProgramSolutionNode>(lp0);
+
+        lp1->take_successors(*lp0);
+
+        CHECK_THAT(feas->predecessors(), RangeEquals({lp1}));
+        CHECK_THAT(obj->predecessors(), RangeEquals({lp1}));
+        CHECK_THAT(sol->predecessors(), RangeEquals({lp1}));
+
+        c1->take_successors(*c0);
+        b_lb1->take_successors(*b_lb0);
+        A1->take_successors(*A0);
+        A_eq1->take_successors(*A_eq0);
+        b_eq1->take_successors(*b_eq0);
+        ub1->take_successors(*ub0);
+
+        CHECK_THAT(lp0->predecessors(), RangeEquals({c1, b_lb1, A1, A_eq1, b_eq1, ub1}));
+        CHECK_THAT(lp1->predecessors(), RangeEquals({c1, b_lb1, A1, A_eq1, b_eq1, ub1}));
+    }
 }
 
 }  // namespace dwave::optimization
diff --git a/tests/cpp/nodes/test_manipulation.cpp b/tests/cpp/nodes/test_manipulation.cpp
index 71d2bc1b..6071763c 100644
--- a/tests/cpp/nodes/test_manipulation.cpp
+++ b/tests/cpp/nodes/test_manipulation.cpp
@@ -337,6 +337,37 @@ TEST_CASE("BroadcastToNode") {
             }
         }
     }
+
+    SECTION("equality") {
+        auto* arr0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2});
+        auto* arr1_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2});
+
+        Node* a_ptr = graph.emplace_node<BroadcastToNode>(arr0_ptr, std::vector<ssize_t>{2, 3});
+        Node* b_ptr = graph.emplace_node<BroadcastToNode>(arr0_ptr, std::vector<ssize_t>{2, 3});
+        Node* c_ptr = graph.emplace_node<BroadcastToNode>(arr1_ptr, std::vector<ssize_t>{2, 3});
+        Node* d_ptr = graph.emplace_node<BroadcastToNode>(arr0_ptr, std::vector<ssize_t>{5, 3});
+    
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+
+        CHECK(not a_ptr->equal_to(*arr0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+        CHECK(not a_ptr->equal_to(*d_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto* arr0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2});
+        auto* arr1_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 4});
+
+        auto* b_ptr = graph.emplace_node<BroadcastToNode>(arr0_ptr, std::vector<ssize_t>{2, 3});
+
+        arr1_ptr->take_successors(*arr0_ptr);
+
+        CHECK_THAT(b_ptr->predecessors(), RangeEquals({arr1_ptr}));
+
+        auto state = graph.initialize_state();
+        CHECK_THAT(b_ptr->view(state), RangeEquals({0, 1, 4, 0, 1, 4}));
+    }
 }
 
 TEST_CASE("ConcatenateNode") {
@@ -664,6 +695,61 @@ TEST_CASE("ConcatenateNode") {
             THEN("The concatenated node is not integral") { CHECK(c.integral() == false); }
         }
     }
+
+    SECTION("equality") {
+        auto graph = Graph();
+
+        auto* arr0_ptr = graph.emplace_node<ConstantNode>(
+            std::vector<double>{1, 2.2, 3}, std::vector<ssize_t>{1, 3}
+        );
+        auto* arr1_ptr = graph.emplace_node<ConstantNode>(
+            std::vector<double>{1, 2, 3}, std::vector<ssize_t>{1, 3}
+        );
+        auto* arr2_ptr = graph.emplace_node<ConstantNode>(
+            std::vector<double>{1, 2, 3}, std::vector<ssize_t>{1, 3}
+        );
+
+        Node* a_ptr =
+            graph.emplace_node<ConcatenateNode>(std::vector<ArrayNode*>{arr0_ptr, arr1_ptr}, 0);
+        Node* b_ptr =
+            graph.emplace_node<ConcatenateNode>(std::vector<ArrayNode*>{arr0_ptr, arr1_ptr}, 0);
+        Node* c_ptr =
+            graph.emplace_node<ConcatenateNode>(std::vector<ArrayNode*>{arr0_ptr, arr2_ptr}, 0);
+        Node* d_ptr =
+            graph.emplace_node<ConcatenateNode>(std::vector<ArrayNode*>{arr0_ptr, arr1_ptr}, 1);
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+
+        CHECK(not a_ptr->equal_to(*arr0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+        CHECK(not a_ptr->equal_to(*d_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto graph = Graph();
+
+        auto* arr0_ptr = graph.emplace_node<ConstantNode>(
+            std::vector<double>{1, 2.2, 3}, std::vector<ssize_t>{1, 3}
+        );
+        auto* arr1_ptr = graph.emplace_node<ConstantNode>(
+            std::vector<double>{1, 2, 3}, std::vector<ssize_t>{1, 3}
+        );
+        auto* arr2_ptr = graph.emplace_node<ConstantNode>(
+            std::vector<double>{4, 5, 6}, std::vector<ssize_t>{1, 3}
+        );
+
+        auto* concat_ptr =
+            graph.emplace_node<ConcatenateNode>(std::vector<ArrayNode*>{arr0_ptr, arr1_ptr}, 0);
+
+        arr2_ptr->take_successors(*arr1_ptr);
+        arr1_ptr->take_successors(*arr0_ptr);
+
+        CHECK_THAT(concat_ptr->predecessors(), RangeEquals({arr1_ptr, arr2_ptr}));
+
+        auto state = graph.initialize_state();
+        CHECK_THAT(concat_ptr->view(state), RangeEquals({1, 2, 3, 4, 5, 6}));
+    }
 }
 
 TEST_CASE("CopyNode") {
@@ -768,6 +854,39 @@ TEST_CASE("CopyNode") {
             }
         }
     }
+
+    SECTION("equality") {
+        auto graph = Graph();
+
+        auto* x0_ptr = graph.emplace_node<SetNode>(6);
+        auto* x1_ptr = graph.emplace_node<SetNode>(6);
+        
+        Node* a_ptr = graph.emplace_node<CopyNode>(x0_ptr);
+        Node* b_ptr = graph.emplace_node<CopyNode>(x0_ptr);
+        Node* c_ptr = graph.emplace_node<CopyNode>(x1_ptr);
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+
+        CHECK(not a_ptr->equal_to(*x0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto graph = Graph();
+
+        auto* x0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2});
+        auto* x1_ptr = graph.emplace_node<ConstantNode>(std::vector{3, 4, 5});
+        
+        auto* copy_ptr = graph.emplace_node<CopyNode>(x0_ptr);
+
+        x1_ptr->take_successors(*x0_ptr);
+
+        CHECK_THAT(copy_ptr->predecessors(), RangeEquals({x1_ptr}));
+
+        auto state = graph.initialize_state();
+        CHECK_THAT(copy_ptr->view(state), RangeEquals({3, 4, 5}));
+    }
 }
 
 TEST_CASE("PutNode") {
@@ -1057,6 +1176,57 @@ TEST_CASE("PutNode") {
             CHECK_THAT(put_ptr->mask(state), RangeEquals({2, 0, 0, 0}));
         }
     }
+
+    SECTION("equality") {
+        auto graph = Graph();
+
+        auto* a0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3, 4});
+        auto* a1_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3, 4});
+
+        auto* i0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 2});
+        auto* i1_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1});
+    
+        auto* v0_ptr = graph.emplace_node<ConstantNode>(std::vector{-44, -55});
+        auto* v1_ptr = graph.emplace_node<ConstantNode>(std::vector{44, 55});
+
+        Node* a_ptr = graph.emplace_node<PutNode>(a0_ptr, i0_ptr, v0_ptr);
+        Node* b_ptr = graph.emplace_node<PutNode>(a0_ptr, i0_ptr, v0_ptr);
+        Node* c_ptr = graph.emplace_node<PutNode>(a1_ptr, i0_ptr, v0_ptr);
+        Node* d_ptr = graph.emplace_node<PutNode>(a0_ptr, i1_ptr, v0_ptr);
+        Node* e_ptr = graph.emplace_node<PutNode>(a0_ptr, i0_ptr, v1_ptr);
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+
+        CHECK(not a_ptr->equal_to(*a0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+        CHECK(not a_ptr->equal_to(*d_ptr));
+        CHECK(not a_ptr->equal_to(*e_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto graph = Graph();
+
+        auto* a0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3, 4});
+        auto* a1_ptr = graph.emplace_node<ConstantNode>(std::vector{5, 6, 7, 8, 9});
+
+        auto* i0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 2});
+        auto* i1_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1});
+    
+        auto* v0_ptr = graph.emplace_node<ConstantNode>(std::vector{-44, -55});
+        auto* v1_ptr = graph.emplace_node<ConstantNode>(std::vector{44, 55});
+
+        auto* put_ptr = graph.emplace_node<PutNode>(a0_ptr, i0_ptr, v0_ptr);
+
+        a1_ptr->take_successors(*a0_ptr);
+        i1_ptr->take_successors(*i0_ptr);
+        v1_ptr->take_successors(*v0_ptr);
+
+        CHECK_THAT(put_ptr->predecessors(), RangeEquals({a1_ptr, i1_ptr, v1_ptr}));
+
+        auto state = graph.initialize_state();
+        CHECK_THAT(put_ptr->view(state), RangeEquals({44, 55, 7, 8, 9}));
+    }
 }
 
 TEST_CASE("ReshapeNode") {
@@ -1208,6 +1378,41 @@ TEST_CASE("ReshapeNode") {
 
         CHECK(graph.num_nodes() == 1);  // no side effects
     }
+
+    SECTION("equality") {
+        auto graph = Graph();
+
+        auto* arr0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3});
+        auto* arr1_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3});
+
+        Node* a_ptr = graph.emplace_node<ReshapeNode>(arr0_ptr, std::array<ssize_t, 2>{2, 2});
+        Node* b_ptr = graph.emplace_node<ReshapeNode>(arr0_ptr, std::array<ssize_t, 2>{2, 2});
+        Node* c_ptr = graph.emplace_node<ReshapeNode>(arr1_ptr, std::array<ssize_t, 2>{2, 2});
+        Node* d_ptr = graph.emplace_node<ReshapeNode>(arr0_ptr, std::array<ssize_t, 2>{4, 1});
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+
+        CHECK(not a_ptr->equal_to(*arr0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+        CHECK(not a_ptr->equal_to(*d_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto graph = Graph();
+
+        auto* arr0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3});
+        auto* arr1_ptr = graph.emplace_node<ConstantNode>(std::vector{2, 3, 4, 5});
+
+        auto* reshape_ptr = graph.emplace_node<ReshapeNode>(arr0_ptr, std::array<ssize_t, 2>{2, 2});
+    
+        arr1_ptr->take_successors(*arr0_ptr);
+
+        CHECK_THAT(reshape_ptr->predecessors(), RangeEquals({arr1_ptr}));
+
+        auto state = graph.initialize_state();
+        CHECK_THAT(reshape_ptr->view(state), RangeEquals({2, 3, 4, 5}));
+    }
 }
 
 TEST_CASE("ResizeNode") {
@@ -1297,6 +1502,45 @@ TEST_CASE("ResizeNode") {
         auto state = graph.initialize_state();
         CHECK_THAT(resize_ptr->view(state), RangeEquals({0, 1, 2, 3, 4, -1}));
     }
+
+    SECTION("equality") {
+        auto graph = Graph();
+
+        auto* set0_ptr = graph.emplace_node<SetNode>(10);
+        auto* set1_ptr = graph.emplace_node<SetNode>(10);
+
+        Node* a_ptr = graph.emplace_node<ResizeNode>(set0_ptr, std::array{5}, 15.1);
+        Node* b_ptr = graph.emplace_node<ResizeNode>(set0_ptr, std::array{5}, 15.1);
+        Node* c_ptr = graph.emplace_node<ResizeNode>(set1_ptr, std::array{5}, 15.1);
+        Node* d_ptr = graph.emplace_node<ResizeNode>(set0_ptr, std::array{6}, 15.1);
+        Node* e_ptr = graph.emplace_node<ResizeNode>(set0_ptr, std::array{5}, -15.1);
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+
+        CHECK(not a_ptr->equal_to(*set0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+        CHECK(not a_ptr->equal_to(*d_ptr));
+        CHECK(not a_ptr->equal_to(*e_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto graph = Graph();
+
+        auto* set0_ptr = graph.emplace_node<SetNode>(10);
+        auto* set1_ptr = graph.emplace_node<SetNode>(10);
+
+        auto* resize_ptr = graph.emplace_node<ResizeNode>(set0_ptr, std::array{5}, 15.1);
+
+        set1_ptr->take_successors(*set0_ptr);
+
+        CHECK_THAT(resize_ptr->predecessors(), RangeEquals({set1_ptr}));
+
+        auto state = graph.empty_state();
+        set1_ptr->initialize_state(state, {0, 1, 2});
+        graph.initialize_state(state);
+        CHECK_THAT(resize_ptr->view(state), RangeEquals(std::vector<double>{0, 1, 2, 15.1, 15.1}));
+    }
 }
 
 TEST_CASE("RollNode") {
@@ -1495,6 +1739,52 @@ TEST_CASE("RollNode") {
             }
         }
     }
+
+    SECTION("equality") {
+        auto graph = Graph();
+
+        auto* arr0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3}, std::array<ssize_t, 2>{2, 2});
+        auto* arr1_ptr = graph.emplace_node<ConstantNode>(std::vector{4, 5, 6, 7}, std::array<ssize_t, 2>{2, 2});
+
+        auto* shift0_ptr = graph.emplace_node<ConstantNode>(2);
+        auto* shift1_ptr = graph.emplace_node<ConstantNode>(3);
+
+        Node* a_ptr = graph.emplace_node<RollNode>(arr0_ptr, shift0_ptr, std::vector<ssize_t>{1});
+        Node* b_ptr = graph.emplace_node<RollNode>(arr0_ptr, shift0_ptr, std::vector<ssize_t>{1});
+        Node* c_ptr = graph.emplace_node<RollNode>(arr1_ptr, shift0_ptr, std::vector<ssize_t>{1});
+        Node* d_ptr = graph.emplace_node<RollNode>(arr0_ptr, shift1_ptr, std::vector<ssize_t>{1});
+        Node* e_ptr = graph.emplace_node<RollNode>(arr0_ptr, shift0_ptr, std::vector<ssize_t>{0});
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+
+        CHECK(not a_ptr->equal_to(*arr0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+        CHECK(not a_ptr->equal_to(*d_ptr));
+        CHECK(not a_ptr->equal_to(*e_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto graph = Graph();
+
+        auto* arr0_ptr =
+            graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3}, std::array<ssize_t, 2>{2, 2});
+        auto* arr1_ptr =
+            graph.emplace_node<ConstantNode>(std::vector{4, 5, 6, 7}, std::array<ssize_t, 2>{2, 2});
+
+        auto* shift0_ptr = graph.emplace_node<ConstantNode>(2);
+        auto* shift1_ptr = graph.emplace_node<ConstantNode>(3);
+
+        auto* roll_ptr = graph.emplace_node<RollNode>(arr0_ptr, shift0_ptr);
+
+        arr1_ptr->take_successors(*arr0_ptr);
+        shift1_ptr->take_successors(*shift0_ptr);
+
+        CHECK_THAT(roll_ptr->predecessors(), RangeEquals({arr1_ptr, shift1_ptr}));
+
+        auto state = graph.initialize_state();
+        CHECK_THAT(roll_ptr->view(state), RangeEquals(std::vector<double>{5, 6, 7, 4}));
+    }
 }
 
 TEST_CASE("SizeNode") {
@@ -1623,6 +1913,25 @@ TEST_CASE("SizeNode") {
             }
         }
     }
+
+    SECTION("predecessor replacement") {
+        auto graph = Graph();
+
+        auto* set0_ptr = graph.emplace_node<SetNode>(5);
+        auto* set1_ptr = graph.emplace_node<SetNode>(5);
+
+        auto* size_ptr = graph.emplace_node<SizeNode>(set0_ptr);
+
+        set1_ptr->take_successors(*set0_ptr);
+
+        CHECK_THAT(size_ptr->predecessors(), RangeEquals({set1_ptr}));
+
+        auto state = graph.empty_state();
+        set0_ptr->initialize_state(state, {0});
+        set1_ptr->initialize_state(state, {0, 1, 2, 3});
+        graph.initialize_state(state);
+        CHECK_THAT(size_ptr->view(state), RangeEquals({4}));
+    }
 }
 
 TEST_CASE("TransposeNode") {
@@ -1993,6 +2302,24 @@ TEST_CASE("TransposeNode") {
             }
         }
     }
+
+    SECTION("predecessor replacement") {
+        auto graph = Graph();
+
+        auto* arr0_ptr =
+            graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3}, std::array<ssize_t, 2>{2, 2});
+        auto* arr1_ptr =
+            graph.emplace_node<ConstantNode>(std::vector{4, 5, 6, 7}, std::array<ssize_t, 2>{2, 2});
+
+        auto* transpose_ptr = graph.emplace_node<TransposeNode>(arr0_ptr);
+
+        arr1_ptr->take_successors(*arr0_ptr);
+
+        CHECK_THAT(transpose_ptr->predecessors(), RangeEquals({arr1_ptr}));
+
+        auto state = graph.initialize_state();
+        CHECK_THAT(transpose_ptr->view(state), RangeEquals({4, 6, 5, 7}));
+    }
 }
 
 }  // namespace dwave::optimization
diff --git a/tests/cpp/nodes/test_naryop.cpp b/tests/cpp/nodes/test_naryop.cpp
index 9d3d2fc2..1aee9753 100644
--- a/tests/cpp/nodes/test_naryop.cpp
+++ b/tests/cpp/nodes/test_naryop.cpp
@@ -14,6 +14,7 @@
 
 #include <catch2/catch_template_test_macros.hpp>
 #include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_range_equals.hpp>
 
 #include "dwave-optimization/graph.hpp"
 #include "dwave-optimization/nodes/binaryop.hpp"
@@ -24,6 +25,8 @@
 #include "dwave-optimization/nodes/numbers.hpp"
 #include "dwave-optimization/nodes/testing.hpp"
 
+using Catch::Matchers::RangeEquals;
+
 namespace dwave::optimization {
 
 TEMPLATE_TEST_CASE(
@@ -186,6 +189,19 @@ TEMPLATE_TEST_CASE(
                 }
             }
         }
+
+        THEN("We can replace predecessors") {
+            auto* a0_ptr = graph.emplace_node<ListNode>(4);
+            auto* b0_ptr = graph.emplace_node<ListNode>(4);
+            auto* c0_ptr = graph.emplace_node<ListNode>(4);
+
+            a0_ptr->take_successors(*a_ptr);
+            b0_ptr->take_successors(*b_ptr);
+            c0_ptr->take_successors(*c_ptr);
+
+            CHECK_THAT(p_ptr->predecessors(), RangeEquals({a0_ptr, b0_ptr, c0_ptr}));
+            CHECK_THAT(p_ptr->operands(), RangeEquals({a0_ptr, b0_ptr, c0_ptr}));
+        }
     }
 
     GIVEN("A vector<Node*> of constants") {
diff --git a/tests/cpp/nodes/test_reduce.cpp b/tests/cpp/nodes/test_reduce.cpp
index 794e687e..31c28d5e 100644
--- a/tests/cpp/nodes/test_reduce.cpp
+++ b/tests/cpp/nodes/test_reduce.cpp
@@ -261,6 +261,44 @@ TEMPLATE_TEST_CASE(
             }
         }
     }
+
+    SECTION("equality") {
+        auto graph = Graph();
+
+        auto* c0_ptr =
+            graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3}, std::vector<ssize_t>{2, 2});
+        auto* c1_ptr =
+            graph.emplace_node<ConstantNode>(std::vector{4, 5, 6, 7}, std::vector<ssize_t>{2, 2});
+
+        Node* a_ptr = graph.emplace_node<ReduceNode<TestType>>(c0_ptr, std::vector<ssize_t>{}, 1.5);
+        Node* b_ptr = graph.emplace_node<ReduceNode<TestType>>(c0_ptr, std::vector<ssize_t>{}, 1.5);
+        Node* c_ptr = graph.emplace_node<ReduceNode<TestType>>(c1_ptr, std::vector<ssize_t>{}, 1.5);
+        Node* d_ptr = graph.emplace_node<ReduceNode<TestType>>(c0_ptr, std::vector<ssize_t>{0}, 1.5);
+        Node* e_ptr = graph.emplace_node<ReduceNode<TestType>>(c0_ptr, std::vector<ssize_t>{}, 2.5);
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+        CHECK(not a_ptr->equal_to(*c0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+        CHECK(not a_ptr->equal_to(*d_ptr));
+        CHECK(not a_ptr->equal_to(*e_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto graph = Graph();
+
+        auto* c0_ptr =
+            graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3}, std::vector<ssize_t>{2, 2});
+        auto* c1_ptr =
+            graph.emplace_node<ConstantNode>(std::vector{4, 5, 6, 7}, std::vector<ssize_t>{2, 2});
+
+        auto* reduce_ptr = graph.emplace_node<ReduceNode<TestType>>(c0_ptr, std::vector<ssize_t>{0});
+
+        c1_ptr->take_successors(*c0_ptr);
+
+        CHECK_THAT(reduce_ptr->predecessors(), RangeEquals({c1_ptr}));
+        CHECK_THAT(reduce_ptr->operands(), RangeEquals({c1_ptr}));
+    }
 }
 
 TEST_CASE("AllNode/AnyNode") {
diff --git a/tests/cpp/nodes/test_set_routines.cpp b/tests/cpp/nodes/test_set_routines.cpp
index 6770605f..de79611f 100644
--- a/tests/cpp/nodes/test_set_routines.cpp
+++ b/tests/cpp/nodes/test_set_routines.cpp
@@ -265,5 +265,23 @@ TEST_CASE("IsInNode") {
                 }
             }
         }
+
+    SECTION("predecessor replacement") {
+        auto* e0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 6, 7});
+        auto* e1_ptr = graph.emplace_node<ConstantNode>(std::vector{2, 3, 4, 5});
+
+        auto* t0_ptr = graph.emplace_node<ConstantNode>(std::vector{1, 2});
+        auto* t1_ptr = graph.emplace_node<ConstantNode>(std::vector{6, 5});
+
+        auto* isin_ptr = graph.emplace_node<IsInNode>(e0_ptr, t0_ptr);
+
+        e1_ptr->take_successors(*e0_ptr);
+        t1_ptr->take_successors(*t0_ptr);
+
+        CHECK_THAT(isin_ptr->predecessors(), RangeEquals({e1_ptr, t1_ptr}));
+
+        auto state = graph.initialize_state();
+        CHECK_THAT(isin_ptr->view(state), RangeEquals({0, 0, 0, 1}));
     }
+}
 }  // namespace dwave::optimization
diff --git a/tests/cpp/nodes/test_softmax.cpp b/tests/cpp/nodes/test_softmax.cpp
index 848065f9..92595f0d 100644
--- a/tests/cpp/nodes/test_softmax.cpp
+++ b/tests/cpp/nodes/test_softmax.cpp
@@ -14,6 +14,7 @@
 
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <catch2/matchers/catch_matchers_range_equals.hpp>
 
 #include "dwave-optimization/nodes/constants.hpp"
 #include "dwave-optimization/nodes/numbers.hpp"
@@ -22,6 +23,7 @@
 
 namespace dwave::optimization {
 
+using Catch::Matchers::RangeEquals;
 using Catch::Matchers::WithinRel;
 
 TEST_CASE("SoftMaxNode") {
@@ -236,5 +238,37 @@ TEST_CASE("SoftMaxNode") {
             }
         }
     }
+
+    SECTION("equality") {
+        auto graph = Graph();
+
+        auto* c0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3});
+        auto* c1_ptr = graph.emplace_node<ConstantNode>(std::vector{4, 5, 6, 7});
+
+        Node* a_ptr = graph.emplace_node<SoftMaxNode>(c0_ptr);
+        Node* b_ptr = graph.emplace_node<SoftMaxNode>(c0_ptr);
+        Node* c_ptr = graph.emplace_node<SoftMaxNode>(c1_ptr);
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+        CHECK(not a_ptr->equal_to(*c0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto graph = Graph();
+
+        auto* c0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3});
+        auto* c1_ptr = graph.emplace_node<ConstantNode>(std::vector{4, 5, 6, 7});
+
+        auto* softmax_ptr = graph.emplace_node<SoftMaxNode>(c0_ptr);
+
+        c1_ptr->take_successors(*c0_ptr);
+
+        CHECK_THAT(softmax_ptr->predecessors(), RangeEquals({c1_ptr}));
+
+        auto state = graph.initialize_state();
+        CHECK_THAT(softmax_ptr->view(state)[0], WithinRel(0.03205860328008499, 1e-9));
+    }
 }
 }  // namespace dwave::optimization
diff --git a/tests/cpp/nodes/test_sorting.cpp b/tests/cpp/nodes/test_sorting.cpp
index df231077..393b7583 100644
--- a/tests/cpp/nodes/test_sorting.cpp
+++ b/tests/cpp/nodes/test_sorting.cpp
@@ -14,7 +14,7 @@
 
 #include <catch2/benchmark/catch_benchmark_all.hpp>
 #include <catch2/catch_test_macros.hpp>
-#include <catch2/matchers/catch_matchers_all.hpp>
+#include <catch2/matchers/catch_matchers_range_equals.hpp>
 
 #include "dwave-optimization/nodes/collections.hpp"
 #include "dwave-optimization/nodes/constants.hpp"
@@ -151,6 +151,38 @@ TEST_CASE("ArgSortNode") {
             }
         }
     }
+
+    SECTION("equality") {
+        auto graph = Graph();
+
+        auto* c0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3});
+        auto* c1_ptr = graph.emplace_node<ConstantNode>(std::vector{4, 5, 6, 7});
+
+        Node* a_ptr = graph.emplace_node<ArgSortNode>(c0_ptr);
+        Node* b_ptr = graph.emplace_node<ArgSortNode>(c0_ptr);
+        Node* c_ptr = graph.emplace_node<ArgSortNode>(c1_ptr);
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+        CHECK(not a_ptr->equal_to(*c0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto graph = Graph();
+
+        auto* c0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3});
+        auto* c1_ptr = graph.emplace_node<ConstantNode>(std::vector{7, 6, 5, 4});
+
+        auto* argsort_ptr = graph.emplace_node<ArgSortNode>(c0_ptr);
+
+        c1_ptr->take_successors(*c0_ptr);
+
+        CHECK_THAT(argsort_ptr->predecessors(), RangeEquals({c1_ptr}));
+
+        auto state = graph.initialize_state();
+        CHECK_THAT(argsort_ptr->view(state), RangeEquals({3, 2, 1, 0}));
+    }
 }
 
 }  // namespace dwave::optimization
diff --git a/tests/cpp/nodes/test_statistics.cpp b/tests/cpp/nodes/test_statistics.cpp
index 2efbe896..faeffc60 100644
--- a/tests/cpp/nodes/test_statistics.cpp
+++ b/tests/cpp/nodes/test_statistics.cpp
@@ -21,6 +21,8 @@
 #include "dwave-optimization/nodes/statistics.hpp"
 #include "dwave-optimization/nodes/testing.hpp"
 
+using Catch::Matchers::RangeEquals;
+
 namespace dwave::optimization {
 
 TEST_CASE("MeanNode") {
@@ -155,5 +157,38 @@ TEST_CASE("MeanNode") {
             }
         }
     }
+
+    SECTION("equality") {
+        auto graph = Graph();
+
+        auto* c0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3});
+        auto* c1_ptr = graph.emplace_node<ConstantNode>(std::vector{4, 5, 6, 7});
+
+        Node* a_ptr = graph.emplace_node<MeanNode>(c0_ptr);
+        Node* b_ptr = graph.emplace_node<MeanNode>(c0_ptr);
+        Node* c_ptr = graph.emplace_node<MeanNode>(c1_ptr);
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+        CHECK(not a_ptr->equal_to(*c0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto graph = Graph();
+
+        auto* c0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3});
+        auto* c1_ptr = graph.emplace_node<ConstantNode>(std::vector{7, 6, 5, 4});
+
+        auto* argsort_ptr = graph.emplace_node<MeanNode>(c0_ptr);
+
+        c1_ptr->take_successors(*c0_ptr);
+
+        CHECK_THAT(argsort_ptr->predecessors(), RangeEquals({c1_ptr}));
+
+        auto state = graph.initialize_state();
+        CHECK_THAT(argsort_ptr->view(state), RangeEquals({5.5}));
+    }
+
 }
 }  // namespace dwave::optimization
diff --git a/tests/cpp/nodes/test_unaryop.cpp b/tests/cpp/nodes/test_unaryop.cpp
index e4f85304..6788c149 100644
--- a/tests/cpp/nodes/test_unaryop.cpp
+++ b/tests/cpp/nodes/test_unaryop.cpp
@@ -244,6 +244,37 @@ TEMPLATE_TEST_CASE(
             }
         }
     }
+
+    SECTION("equality") {
+        auto graph = Graph();
+
+        auto* c0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3});
+        auto* c1_ptr = graph.emplace_node<ConstantNode>(std::vector{4, 5, 6, 7});
+
+        Node* a_ptr = graph.emplace_node<UnaryOpNode<TestType>>(c0_ptr);
+        Node* b_ptr = graph.emplace_node<UnaryOpNode<TestType>>(c0_ptr);
+        Node* c_ptr = graph.emplace_node<UnaryOpNode<TestType>>(c1_ptr);
+
+        CHECK(a_ptr->equal_to(*a_ptr));
+        CHECK(a_ptr->equal_to(*b_ptr));
+        CHECK(not a_ptr->equal_to(*c0_ptr));
+        CHECK(not a_ptr->equal_to(*c_ptr));
+    }
+
+    SECTION("predecessor replacement") {
+        auto graph = Graph();
+
+        auto* c0_ptr = graph.emplace_node<ConstantNode>(std::vector{0, 1, 2, 3});
+        auto* c1_ptr = graph.emplace_node<ConstantNode>(std::vector{7, 6, 5, 4});
+
+        auto* argsort_ptr = graph.emplace_node<UnaryOpNode<TestType>>(c0_ptr);
+
+        c1_ptr->take_successors(*c0_ptr);
+
+        CHECK_THAT(argsort_ptr->predecessors(), RangeEquals({c1_ptr}));
+        CHECK_THAT(argsort_ptr->operands(), RangeEquals({c1_ptr}));
+    }
+
 }
 
 TEST_CASE("UnaryOpNode - AbsoluteNode") {
diff --git a/tests/cpp/test_graph.cpp b/tests/cpp/test_graph.cpp
index cfad9aa2..9ae62875 100644
--- a/tests/cpp/test_graph.cpp
+++ b/tests/cpp/test_graph.cpp
@@ -13,11 +13,13 @@
 //    limitations under the License.
 
 #include <catch2/catch_test_macros.hpp>
-#include <catch2/matchers/catch_matchers_all.hpp>
+#include <catch2/matchers/catch_matchers_range_equals.hpp>
 
 #include "dwave-optimization/graph.hpp"
 #include "dwave-optimization/nodes.hpp"
 
+using Catch::Matchers::RangeEquals;
+
 namespace dwave::optimization {
 
 TEST_CASE("Topological Sort", "[topological_sort]") {
@@ -243,7 +245,7 @@ TEST_CASE("Graph::commit(), Graph::descendants(), Graph::propagate(), and Graph:
 
     SECTION("Find descendants") {
         auto descendants = graph.descendants({x_ptr});
-        CHECK_THAT(descendants, Catch::Matchers::RangeEquals(std::vector<Node*>{x_ptr, z_ptr}));
+        CHECK_THAT(descendants, RangeEquals(std::vector<Node*>{x_ptr, z_ptr}));
     }
     SECTION("Propagate all") {
         CHECK(x_ptr->view(state).front() == 0);
@@ -316,6 +318,106 @@ TEST_CASE("Graph::objective()") {
     }
 }
 
+TEST_CASE("Graph::remove_redundant_nodes()") {
+    GIVEN("A model with two redundant nodes") {
+        auto graph = Graph();
+
+        auto* x_ptr = graph.emplace_node<BinaryNode>();
+        auto* y_ptr = graph.emplace_node<BinaryNode>();
+
+        auto* left_x_plus_y = graph.emplace_node<AddNode>(x_ptr, y_ptr);
+        auto* right_x_plus_y = graph.emplace_node<AddNode>(x_ptr, y_ptr);
+
+        CHECK(graph.num_nodes() == 4);
+
+        THEN("remove_redundant_nodes() removes one") {
+            CHECK(graph.remove_redundant_nodes() == 1);
+
+            CHECK(graph.num_nodes() == 3);
+
+            // we kept the one that's earlier in the topological order
+            CHECK(graph.nodes()[2].get() == left_x_plus_y);
+        }
+
+        AND_GIVEN("Two more redundant nodes decended from them") {
+            auto* negative_left_x_plus_y = graph.emplace_node<NegativeNode>(left_x_plus_y);
+            auto* negative_right_x_plus_y = graph.emplace_node<NegativeNode>(right_x_plus_y);
+
+            THEN("remove_redundant_nodes() removes 2") {
+                CHECK(graph.remove_redundant_nodes() == 2);
+
+                CHECK(graph.num_nodes() == 4);
+
+                // we kept the ones earlier in the topological order
+                CHECK(graph.nodes()[2].get() == left_x_plus_y);
+                CHECK(graph.nodes()[3].get() == negative_left_x_plus_y);
+            }
+
+            WHEN("we have a listener on one of the redundant nodes") {
+                auto listener_ptr = negative_right_x_plus_y->expired_ptr();
+
+                THEN("by default we don't remove the listened to node") {
+                    CHECK(graph.remove_redundant_nodes() == 1);
+                    CHECK(graph.num_nodes() == 5);
+
+                    // we kept the ones earlier in the topological order
+                    CHECK(graph.nodes()[2].get() == left_x_plus_y);
+                    CHECK(graph.nodes()[3].get() == negative_left_x_plus_y);
+                    CHECK(graph.nodes()[4].get() == negative_right_x_plus_y);
+                }
+
+                THEN("we can force the removal of the listened-to node") {
+                    CHECK(graph.remove_redundant_nodes(true) == 2);
+
+                    CHECK(graph.num_nodes() == 4);
+
+                    // we kept the ones earlier in the topological order
+                    CHECK(graph.nodes()[2].get() == left_x_plus_y);
+                    CHECK(graph.nodes()[3].get() == negative_left_x_plus_y);
+                }
+            }
+        }
+
+        AND_GIVEN("that the later redundant node is used in the objective") {
+            graph.set_objective(right_x_plus_y);
+
+            CHECK(graph.remove_redundant_nodes() == 1);
+
+            CHECK(graph.nodes()[2].get() == left_x_plus_y);
+            CHECK(graph.objective() == left_x_plus_y);  // objective was updated
+        }
+
+        AND_GIVEN("two more redundant logical nodes") {
+            auto* logical_left_x_plus_y = graph.emplace_node<LogicalNode>(left_x_plus_y);
+            auto* logical_right_x_plus_y = graph.emplace_node<LogicalNode>(right_x_plus_y);
+
+            WHEN("both are registered as constraints") {
+                graph.add_constraint(logical_left_x_plus_y);
+                graph.add_constraint(logical_right_x_plus_y);
+
+                CHECK(graph.remove_redundant_nodes() == 2);
+                CHECK(graph.constraints().size() == 1);
+                CHECK_THAT(graph.constraints(), RangeEquals({logical_left_x_plus_y}));
+            }
+        }
+    }
+
+    GIVEN("A model with two redundant constants used in the same unary op") {
+        auto graph = Graph();
+
+        auto* a = graph.emplace_node<ConstantNode>(5);
+        auto* negative_a = graph.emplace_node<NegativeNode>(a);
+        auto* b = graph.emplace_node<ConstantNode>(5);
+        graph.emplace_node<NegativeNode>(b);
+
+        THEN("remove_redundant_nodes() removes the redundant path") {
+            CHECK(graph.remove_redundant_nodes() == 2);
+            CHECK_THAT(graph.constants(), RangeEquals({a}));
+            CHECK_THAT(negative_a->predecessors(), RangeEquals({a}));
+        }
+    }
+}
+
 TEST_CASE("Graph::remove_unused_nodes()") {
     GIVEN("A single integer variable") {
         auto graph = Graph();
diff --git a/tests/test_model.py b/tests/test_model.py
index 41de0f34..f798781f 100644
--- a/tests/test_model.py
+++ b/tests/test_model.py
@@ -326,6 +326,22 @@ def test_objective(self):
         with self.assertRaises(TypeError):
             model.objective = "hello"
 
+    def test_remove_redundant_symbols(self):
+        model = Model()
+
+        x = model.binary()
+        y = model.binary()
+
+        first = x + y
+        second = y + x
+
+        self.assertEqual(model.remove_redundant_symbols(), 0)  # both are named symbols
+
+        del second
+        self.assertEqual(model.remove_redundant_symbols(), 1)
+
+        self.assertEqual(list(model.iter_symbols())[2].id(), first.id())
+
     def test_remove_unused_symbols(self):
         with self.subTest("all unused"):
             model = Model()