diff --git a/src/operators/reduce-nd.c b/src/operators/reduce-nd.c
index fb34c75..818db69 100644
--- a/src/operators/reduce-nd.c
+++ b/src/operators/reduce-nd.c
@@ -142,12 +142,12 @@
     return xnn_status_unsupported_parameter;
   }
 
-  if (num_reduction_axes > num_input_dims) {
+  if (num_reduction_axes > XNN_MAX_TENSOR_DIMS) {
     xnn_log_error(
         "failed to reshape %s operator with %zu reduction axes: the number of "
-        "reduction axes must not exceed the number of input dimensions %zu",
+        "reduction axes must not exceed %d",
         xnn_operator_type_to_string_v2(reduce_op), num_reduction_axes,
-        num_input_dims);
+        XNN_MAX_TENSOR_DIMS);
     return xnn_status_invalid_parameter;
   }
 
@@ -163,19 +163,6 @@
   assert(num_input_dims <= XNN_MAX_TENSOR_DIMS);
   memcpy(normalized_input_shape, input_shape, num_input_dims * sizeof(size_t));
 
-  for (size_t i = 0; i < num_reduction_axes; i++) {
-    const int64_t signed_num_input_dims = (int64_t)num_input_dims;
-    if (signed_num_input_dims <= reduction_axes[i] ||
-        reduction_axes[i] < -signed_num_input_dims) {
-      xnn_log_error(
-          "failed to reshape %s operator with #%zu reduction axis of %" PRIi64
-          ": the index is out of bounds for a %zuD input shape",
-          xnn_operator_type_to_string_v2(reduce_op), i, reduction_axes[i],
-          num_input_dims);
-      return xnn_status_invalid_parameter;
-    }
-  }
-
   size_t normalized_reduction_axes[XNN_MAX_TENSOR_DIMS];
   assert(num_reduction_axes <= XNN_MAX_TENSOR_DIMS);
   for (int i = 0; i < num_reduction_axes; i++) {
@@ -186,15 +173,20 @@
   qsort(normalized_reduction_axes, num_reduction_axes, sizeof(size_t),
         cmp_value_size_t);
 
-  for (size_t i = 1; i < num_reduction_axes; i++) {
-    if (normalized_reduction_axes[i] <= normalized_reduction_axes[i - 1]) {
-      xnn_log_error(
-          "failed to reshape %s operator with #%zu reduction axis of %" PRIi64
-          ": the reduction axes must be unique",
-          xnn_operator_type_to_string_v2(reduce_op), i, reduction_axes[i]);
-      return xnn_status_invalid_parameter;
+  // Remove duplicate reduction axes.
+  int i = 0;
+  // The array is sorted, we're done if an axis is bigger than num_input_dims.
+  for (int j = 1;
+       j < num_reduction_axes && normalized_reduction_axes[j] < num_input_dims;
+       ++j) {
+    // Shift non-duplicate elements forward.
+    if (normalized_reduction_axes[i] != normalized_reduction_axes[j]) {
+      normalized_reduction_axes[++i] = normalized_reduction_axes[j];
     }
   }
+  num_reduction_axes = i + 1;
+
+  assert(num_reduction_axes <= num_input_dims);
 
   xnn_normalize_reduction(
     &num_reduction_axes, normalized_reduction_axes,
diff --git a/src/subgraph.c b/src/subgraph.c
index 2d785ec..9682395 100644
--- a/src/subgraph.c
+++ b/src/subgraph.c
@@ -2455,11 +2455,24 @@
   return xnn_status_success;
 }
 
+static void convert_static_value_to_fp32(struct xnn_value* value) {
+  assert(xnn_value_is_static(value->allocation_type));
+  if (value->flags & XNN_VALUE_FLAG_NEEDS_CLEANUP) {
+    xnn_release_memory(value->data);
+  }
+  value->data = xnn_allocate_memory(sizeof(float));
+  float data = get_scalar_value_as_float(value);
+  memcpy(value->data, &data, sizeof(float));
+  value->flags |= XNN_VALUE_FLAG_NEEDS_CLEANUP;
+  value->datatype = xnn_datatype_fp32;
+}
+
 // Replace `mul(reduce_sum(x), 1/n)`, `div(reduce_sum(x), n)`  or
 // `mul(reduce_sum_squared(x), 1/n)`, `div(reduce_sum_squared(x), n)`
 // with `reduce_mean(x)` or `reduce_mean_squared(x)`, respectively.
-static enum xnn_status optimize_common_subgraphs_scaled_sum_to_mean(
-    xnn_subgraph_t subgraph, uint32_t node_id, size_t* changes) {
+static enum xnn_status widen_fp16_accumulators(xnn_subgraph_t subgraph,
+                                               uint32_t node_id,
+                                               size_t* changes) {
   struct xnn_node* node = &subgraph->nodes[node_id];
 
   if (node->type != xnn_node_type_binary_elementwise ||
@@ -2468,78 +2481,66 @@
     return xnn_status_success;
   }
 
-  struct xnn_value* reduce_value = &subgraph->values[node->inputs[0]];
+  struct xnn_value* reduced_value = &subgraph->values[node->inputs[0]];
   struct xnn_value* arg_value = &subgraph->values[node->inputs[1]];
   if (xnn_shape_multiply_all_dims(&arg_value->shape) != 1 ||
       !xnn_value_is_static(arg_value->allocation_type)) {
-    if (xnn_shape_multiply_all_dims(&reduce_value->shape) == 1 &&
-        xnn_value_is_static(reduce_value->allocation_type)) {
-      swap_value_pointers(&reduce_value, &arg_value);
+    if (xnn_shape_multiply_all_dims(&reduced_value->shape) == 1 &&
+        xnn_value_is_static(reduced_value->allocation_type)) {
+      swap_value_pointers(&reduced_value, &arg_value);
     } else {
       return xnn_status_success;
     }
   }
 
   // Check that one of the args is a sum or sum2 reduction.
-  if (!(reduce_value->datatype == xnn_datatype_fp16 ||
-        reduce_value->datatype == xnn_datatype_fp32) ||
-      reduce_value->producer == XNN_INVALID_NODE_ID) {
+  if (reduced_value->datatype != xnn_datatype_fp16 ||
+      reduced_value->producer == XNN_INVALID_NODE_ID) {
     return xnn_status_success;
   }
-  struct xnn_node* reduce_node = &subgraph->nodes[reduce_value->producer];
+
+  if (reduced_value->num_consumers > 1 || arg_value->num_consumers > 1) {
+    // Don't rewrite if we might modify an unrelated consumer.
+    return xnn_status_success;
+  }
+
+  struct xnn_node* reduce_node = &subgraph->nodes[reduced_value->producer];
   const enum xnn_node_type reduce_node_type = reduce_node->type;
   if (!(reduce_node_type == xnn_node_type_static_sum ||
         reduce_node_type == xnn_node_type_static_sum_squared)) {
     return xnn_status_success;
   }
 
-  // Check that the other arg is the product of the dimensions of
-  // the reduction axes, or its inverse.
-  struct xnn_value* input_value = &subgraph->values[reduce_node->inputs[0]];
-  const float arg_as_float = get_scalar_value_as_float(arg_value);
-  size_t num_reduced_dims = 1;
-  for (size_t k = 0; k < reduce_node->params.reduce.num_reduction_axes; k++) {
-    num_reduced_dims *= xnn_shape_get_dim(
-        &input_value->shape, reduce_node->params.reduce.reduction_axes[k]);
-  }
-  if (!num_reduced_dims) {
-    return xnn_status_success;
-  }
-  const enum xnn_binary_operator binary_operator = node->binary_operator;
-  float expected_arg = (binary_operator == xnn_binary_multiply)
-                           ? 1.0f / num_reduced_dims
-                           : (float)num_reduced_dims;
-  if (arg_value->datatype == xnn_datatype_fp16) {
-    expected_arg = xnn_float16_to_float(xnn_float16_from_float(expected_arg));
-  }
-  if (arg_as_float != expected_arg) {
-    return xnn_status_success;
-  }
+  // Rewrite the internal values to this subgraph to be fp32.
+  reduced_value->datatype = xnn_datatype_fp32;
+  convert_static_value_to_fp32(arg_value);
 
-  const uint32_t output_id = node->outputs[0];
-  const size_t num_reduction_axes =
-      reduce_node->params.reduce.num_reduction_axes;
-  int64_t reduction_axes[XNN_MAX_TENSOR_DIMS];
-  memcpy(reduction_axes, reduce_node->params.reduce.reduction_axes,
-         num_reduction_axes * sizeof(int64_t));
-  XNN_RETURN_IF_ERROR(xnn_define_static_reduce_v2(
-                          subgraph,
-                          reduce_node_type == xnn_node_type_static_sum
-                              ? xnn_reduce_mean
-                              : xnn_reduce_mean_squared,
-                          num_reduction_axes, reduction_axes, input_value->id,
-                          output_id, reduce_node->flags),
-                      "Failed to create new `Mean` or `Mean Squared` node.");
-  node = move_last_node_to(subgraph, node_id);
+  uint32_t output_id = node->outputs[0];
+  struct xnn_value* output_value = &subgraph->values[output_id];
+  uint32_t output_fp32_id;
+  enum xnn_status status = xnn_define_tensor_value(
+      subgraph, xnn_datatype_fp32, output_value->shape.num_dims,
+      output_value->shape.dim,
+      /*data=*/NULL, XNN_INVALID_VALUE_ID,
+      /*flags=*/0, &output_fp32_id);
+  if (status != xnn_status_success) {
+    return status;
+  }
+  node->outputs[0] = output_fp32_id;
+
+  status = xnn_define_unary(subgraph, xnn_unary_convert, /*params=*/NULL,
+                            output_fp32_id, output_id, /*flags=*/0);
+  if (status != xnn_status_success) {
+    return status;
+  }
+  move_last_node_to(subgraph, node_id);
 
   xnn_log_info(
       "Converted %s[#%u](reduce_sum%s[#%u](v%03u), v%03u) to "
-      "reduce_mean%s[#%u](v%03u).",
-      (binary_operator == xnn_binary_multiply) ? "mul" : "div", node_id,
+      "fp32.",
+      (node->binary_operator == xnn_binary_multiply) ? "mul" : "div", node_id,
       reduce_node_type == xnn_node_type_static_sum_squared ? "_squared" : "",
-      reduce_value->producer, input_value->id, arg_value->id,
-      reduce_node_type == xnn_node_type_static_sum_squared ? "_squared" : "",
-      node_id, input_value->id);
+      reduced_value->producer, reduce_node->inputs[0], arg_value->id);
   (*changes)++;
 
   return xnn_status_success;
@@ -3678,11 +3679,13 @@
         // XNN_RETURN_IF_ERROR(optimize_common_subgraphs_binary_to_const(
         //     subgraph, node_id, changes));
 
-        // Replace `mul(reduce_sum(x), 1/n)`, `div(reduce_sum(x), n)`  or
-        // `mul(reduce_sum_squared(x), 1/n)`, `div(reduce_sum_squared(x), n)`
-        // with `reduce_mean(x)` or `reduce_mean_squared(x)`, respectively.
-        XNN_RETURN_IF_ERROR(optimize_common_subgraphs_scaled_sum_to_mean(
-            subgraph, node_id, changes));
+        // Widen fp16 accumulators for `mul(reduce_sum(x), y)` or
+        // `div(reduce_sum(x), y)`. This is a bit of a hack to make subgraphs
+        // rewritten to be fp16 less likely to overflow. Especially if x is a
+        // squaring operation, and the data is fp16, it is very likely that the
+        // sum will overflow.
+        XNN_RETURN_IF_ERROR(
+            widen_fp16_accumulators(subgraph, node_id, changes));
 
         // Convert min/max operations with a single static value to a unary
         // `clamp` node.
@@ -4298,10 +4301,6 @@
     return xnn_status_unsupported_hardware;
   }
 
-  // Apply some common subgraph optimizations.
-  XNN_RETURN_IF_ERROR(
-      xnn_subgraph_optimize_common_subgraphs(subgraph, optimization_flags));
-
   if ((optimization_flags & XNN_FLAG_FORCE_FP16_INFERENCE) &&
       (!xnn_is_f16_compatible_config(hardware_config))) {
     xnn_log_error(
@@ -4328,6 +4327,10 @@
     }
   }
 
+  // Apply some common subgraph optimizations.
+  XNN_RETURN_IF_ERROR(
+      xnn_subgraph_optimize_common_subgraphs(subgraph, optimization_flags));
+
 #if XNN_ENABLE_SPARSE
   if ((optimization_flags & XNN_FLAG_HINT_SPARSE_INFERENCE) &&
       (xnn_is_chw_compatible_config(hardware_config))) {
diff --git a/test/operators/reduce-nd.cc b/test/operators/reduce-nd.cc
index e79d0c4..ee0834f 100644
--- a/test/operators/reduce-nd.cc
+++ b/test/operators/reduce-nd.cc
@@ -55,14 +55,12 @@
 
   ReduceOperatorTester& reduction_axes(
       std::initializer_list<int64_t> reduction_axes) {
-    assert(reduction_axes.size() <= XNN_MAX_TENSOR_DIMS);
     this->reduction_axes_ = std::vector<int64_t>(reduction_axes);
     return *this;
   }
 
   ReduceOperatorTester& reduction_axes(
       const std::vector<int64_t> reduction_axes) {
-    assert(reduction_axes.size() <= XNN_MAX_TENSOR_DIMS);
     this->reduction_axes_ = reduction_axes;
     return *this;
   }
@@ -550,11 +548,16 @@
     std::vector<int64_t> reduction_axes;
     for (int i = 0; i < param.dims; ++i) {
       if (reduce_dims[i]) {
-        if (param.use_neg_axes) {
-          reduction_axes.push_back(i - param.dims);
-        } else {
-          reduction_axes.push_back(i);
-        }
+        reduction_axes.push_back(i);
+      }
+    }
+
+    // We support negative axes, and we support redundant axes. We can test both
+    // at the same time by just specifying both expressions of the same axis.
+    for (int i = 0;
+         i < param.dims && reduction_axes.size() < XNN_MAX_TENSOR_DIMS; ++i) {
+      if (reduce_dims[i]) {
+        reduction_axes.push_back(i - param.dims);
       }
     }
     return reduction_axes;
diff --git a/test/subgraph/static-reduce.cc b/test/subgraph/static-reduce.cc
index 6b4f322..2dbe64f 100644
--- a/test/subgraph/static-reduce.cc
+++ b/test/subgraph/static-reduce.cc
@@ -4,6 +4,7 @@
 // LICENSE file in the root directory of this source tree.
 
 #include <algorithm>
+#include <bitset>
 #include <cmath>
 #include <cstddef>
 #include <cstdint>
@@ -29,12 +30,11 @@
 namespace xnnpack {
 
 struct Param {
-  using TupleT = std::tuple<xnn_reduce_operator, bool, bool, int>;
+  using TupleT = std::tuple<xnn_reduce_operator, bool, int>;
   explicit Param(TupleT p)
       : reduce_operator(std::get<0>(p)),
         keep_dims(std::get<1>(p)),
-        use_neg_axes(std::get<2>(p)),
-        rank(std::get<3>(p)) {}
+        rank(std::get<2>(p)) {}
 
   std::string Name() const {
     std::stringstream sstr;
@@ -64,33 +64,31 @@
     if (keep_dims) {
       sstr << "_keep_dims";
     }
-    if (use_neg_axes) {
-      sstr << "_use_neg_axes";
-    }
     sstr << "_" << rank;
     return sstr.str();
   }
 
   xnn_reduce_operator reduce_operator;
   bool keep_dims;
-  bool use_neg_axes;
   int rank;
 };
 
-std::vector<int64_t> mask_to_axes(uint32_t mask) {
+std::vector<int64_t> mask_to_axes(uint32_t mask, int32_t rank) {
   std::vector<int64_t> axes;
-  for (uint32_t i = 0; i < XNN_MAX_TENSOR_DIMS; ++i) {
+  for (int32_t i = 0; i < XNN_MAX_TENSOR_DIMS; ++i) {
     if (mask & (1 << i)) {
       axes.push_back(i);
     }
   }
-  return axes;
-}
-
-void negate_axes(int64_t rank, std::vector<int64_t>& axes) {
-  for (int64_t& axis : axes) {
-    axis = axis - rank;
+  // It is legal to specify an axis more than once. To test negative axes,
+  // we can just specify the same axis in both forms.
+  for (int32_t i = 0;
+       i < XNN_MAX_TENSOR_DIMS && axes.size() < XNN_MAX_TENSOR_DIMS; ++i) {
+    if (mask & (1 << i)) {
+      axes.push_back(i - rank);
+    }
   }
+  return axes;
 }
 
 template <typename T>
@@ -139,10 +137,9 @@
   ASSERT_EQ(xnn_status_success, xnn_initialize(/*allocator=*/nullptr));
 
   for (uint32_t mask = 1; mask < (1 << p.rank); ++mask) {
-    std::vector<int64_t> reduction_axes = mask_to_axes(mask);
-    if (p.use_neg_axes) {
-      negate_axes(p.rank, reduction_axes);
-    }
+    const int num_reduction_axes =
+        std::bitset<XNN_MAX_TENSOR_DIMS>(mask).count();
+    std::vector<int64_t> reduction_axes = mask_to_axes(mask, p.rank);
 
     xnn_quantization_params input_quantization =
         random_quantization(datatype, rng);
@@ -157,7 +154,7 @@
     // Create a runtime with the reduce op in it.
     SubgraphTester tester(2);
     tester.AddInputTensor(p.rank, datatype, input_quantization, 0)
-        .AddOutputTensor(p.keep_dims ? p.rank : p.rank - reduction_axes.size(),
+        .AddOutputTensor(p.keep_dims ? p.rank : p.rank - num_reduction_axes,
                          datatype, output_quantization, 1)
         .AddReduce(p.reduce_operator, reduction_axes, 0, 1,
                    p.keep_dims ? XNN_FLAG_KEEP_DIMS : 0);
@@ -246,14 +243,12 @@
 
 template <typename T, typename Accum = T>
 void TestSubgraphRewrite(const Param& p) {
-  if (p.reduce_operator != xnn_reduce_sum_squared &&
-      p.reduce_operator != xnn_reduce_mean_squared) {
+  if (p.reduce_operator != xnn_reduce_sum_squared) {
     GTEST_SKIP();
     return;
   }
 
   const size_t rank = p.rank;
-  const bool use_neg_axes = p.use_neg_axes;
   const bool keep_dims = p.keep_dims;
   const xnn_reduce_operator reduce_operator = p.reduce_operator;
   auto reference_op = get_reference_op(reduce_operator);
@@ -265,10 +260,7 @@
 
   for (uint32_t mask = 1; mask < (1 << rank); ++mask) {
     for (uint32_t iter = 0; iter < 10; iter++) {
-      std::vector<int64_t> reduction_axes = mask_to_axes(mask);
-      if (use_neg_axes) {
-        negate_axes(rank, reduction_axes);
-      }
+      std::vector<int64_t> reduction_axes = mask_to_axes(mask, rank);
 
       // Define subgraph
       enum external_value_ids : uint32_t {
@@ -290,7 +282,6 @@
           output_shape.push_back(input_shape[i]);
         }
       }
-      const T inv_n = 1.0 / reduced_elements;
 
       // Generate the input.
       Tensor<T> input(input_shape, xnnpack::XnnExtraBytes);
@@ -315,39 +306,8 @@
       }
 
       // c = reduce_sum(b).
-      switch (reduce_operator) {
-        case xnn_reduce_sum_squared:
-          subgraph.AddReduce(xnn_reduce_sum, reduction_axes, squared_id,
-                             output_id,
-                             /*flags=*/keep_dims ? XNN_FLAG_KEEP_DIMS : 0);
-          break;
-        case xnn_reduce_mean_squared:
-          if (rng() % 2) {
-            subgraph.AddReduce(xnn_reduce_mean, reduction_axes, squared_id,
-                               output_id,
-                               /*flags=*/keep_dims ? XNN_FLAG_KEEP_DIMS : 0);
-          } else {
-            uint32_t sum_squared_id = XNN_INVALID_VALUE_ID;
-            subgraph.AddInternalDynamicTensor(
-                output_shape, xnn_datatype_of<T>(), &sum_squared_id,
-                /*flags=*/0);
-            subgraph.AddReduce(xnn_reduce_sum, reduction_axes, squared_id,
-                               sum_squared_id,
-                               /*flags=*/keep_dims ? XNN_FLAG_KEEP_DIMS : 0);
-            // d = mul(c, inv_n).
-            uint32_t inv_n_id = XNN_INVALID_VALUE_ID;
-            subgraph.AddInternalStaticTensor(
-                /*shape=*/{1}, xnn_datatype_of<T>(), &inv_n_id, &inv_n);
-            if (rng() % 2) {
-              subgraph.AddMultiply(sum_squared_id, inv_n_id, output_id);
-            } else {
-              subgraph.AddMultiply(inv_n_id, sum_squared_id, output_id);
-            }
-          }
-          break;
-        default:
-          XNN_UNREACHABLE;
-      }
+      subgraph.AddReduce(xnn_reduce_sum, reduction_axes, squared_id, output_id,
+                         /*flags=*/keep_dims ? XNN_FLAG_KEEP_DIMS : 0);
 
       // Evaluate once with `XNN_FLAG_NO_OPERATOR_FUSION` enabled to
       // prevent the subgraph replacement.
@@ -450,7 +410,7 @@
 
 auto params = testing::ConvertGenerator<Param::TupleT>(Combine(
     Values(xnn_reduce_sum, xnn_reduce_mean, xnn_reduce_max, xnn_reduce_min),
-    Bool(), Bool(), Range(0, XNN_MAX_TENSOR_DIMS)));
+    Bool(), Range(0, XNN_MAX_TENSOR_DIMS)));
 INSTANTIATE_TEST_SUITE_P(Reduce, ReduceQS8, params,
                          [](auto p) { return p.param.Name(); });
 INSTANTIATE_TEST_SUITE_P(Reduce, ReduceQU8, params,
@@ -459,15 +419,14 @@
 auto params2 = testing::ConvertGenerator<Param::TupleT>(Combine(
     Values(xnn_reduce_sum, xnn_reduce_mean, xnn_reduce_max, xnn_reduce_min,
            xnn_reduce_mean_squared, xnn_reduce_sum_squared),
-    Bool(), Bool(), Range(0, XNN_MAX_TENSOR_DIMS)));
+    Bool(), Range(0, XNN_MAX_TENSOR_DIMS)));
 INSTANTIATE_TEST_SUITE_P(Reduce, ReduceF16, params2,
                          [](auto p) { return p.param.Name(); });
 INSTANTIATE_TEST_SUITE_P(Reduce, ReduceF32, params2,
                          [](auto p) { return p.param.Name(); });
 
-auto params3 = testing::ConvertGenerator<Param::TupleT>(
-    Combine(Values(xnn_reduce_mean_squared, xnn_reduce_sum_squared), Bool(),
-            Bool(), Range(0, XNN_MAX_TENSOR_DIMS)));
+auto params3 = testing::ConvertGenerator<Param::TupleT>(Combine(
+    Values(xnn_reduce_sum_squared), Bool(), Range(0, XNN_MAX_TENSOR_DIMS)));
 INSTANTIATE_TEST_SUITE_P(Reduce, ReduceF16Rewrite, params3,
                          [](auto p) { return p.param.Name(); });
 INSTANTIATE_TEST_SUITE_P(Reduce, ReduceF32Rewrite, params3,