diff --git a/mlx/compile.cpp b/mlx/compile.cpp
index d110404190..e6bb2953b5 100644
--- a/mlx/compile.cpp
+++ b/mlx/compile.cpp
@@ -650,11 +650,16 @@ void compile_simplify(
     return pa.is_equivalent(pb);
   };
 
-  // Merge scalars
+  std::unordered_map<uintptr_t, array> output_map;
+  for (auto& o : outputs) {
+    output_map.insert({o.id(), o});
+  }
+
+  // Merge scalars, except outputs
   std::vector<array> new_tape;
   for (auto& arr : tape) {
     // Check if we can merge scalars
-    if (is_scalar(arr)) {
+    if (is_scalar(arr) && output_map.find(arr.id()) == output_map.end()) {
       auto scalar = scalars.find(get_scalar_rep(arr));
       if (scalar->second.id() != arr.id()) {
         merge(scalar->second, arr, parents_map);
@@ -668,10 +673,6 @@ void compile_simplify(
 
   // Remove no-ops
   {
-    std::unordered_map<uintptr_t, array> output_map;
-    for (auto& o : outputs) {
-      output_map.insert({o.id(), o});
-    }
     for (auto& arr : tape) {
       if (!arr.has_primitive() || !is_noop(arr.primitive())) {
         new_tape.push_back(std::move(arr));
@@ -693,11 +694,6 @@ void compile_simplify(
     tape_order.insert({tape[i].id(), i});
   }
 
-  std::unordered_set<uintptr_t> output_set;
-  for (auto& o : outputs) {
-    output_set.insert(o.id());
-  }
-
   // Multi-pass merge only keeping non-orphaned arrays in the tape
   for (int pass = 0; pass < passes; ++pass) {
     for (auto& arr : tape) {
@@ -717,7 +713,7 @@ void compile_simplify(
             auto& src = parents->second[src_idx].first;
             auto& dst = parents->second[dst_idx].first;
             if (src.id() != dst.id() && array_equivalent(src, dst) &&
-                output_set.find(src.id()) == output_set.end()) {
+                output_map.find(src.id()) == output_map.end()) {
               merge(dst, src, parents_map);
               mask[src_idx] = true;
             }
@@ -783,7 +779,7 @@ void compile_simplify(
           }
           return false;
         } else {
-          return output_set.find(a.id()) == output_set.end();
+          return output_map.find(a.id()) == output_map.end();
         }
       };
       bool discard = maybe_merge_parents(arr);
diff --git a/python/tests/test_compile.py b/python/tests/test_compile.py
index 9f1ced1867..7436e58cd2 100644
--- a/python/tests/test_compile.py
+++ b/python/tests/test_compile.py
@@ -765,6 +765,24 @@ def fn(x):
         self.assertEqual(cfn(x4).shape, fn(x4).shape)
         self.assertEqual(cfn(x8).shape, fn(x8).shape)
 
+    def test_compile_scalar_constant_output(self):
+        for shapeless in [False, True]:
+            fun = mx.compile(lambda x: (x + 1.0, mx.array(1.0)), shapeless=shapeless)
+            out, c = fun(mx.array(3.0))
+            self.assertEqual(out.item(), 4.0)
+            self.assertEqual(c.item(), 1.0)
+
+        s = mx.array(7.0)
+        out, c = mx.compile(lambda x: (x + 7.0, s))(mx.array([1.0, 2.0]))
+        self.assertEqual(out.tolist(), [8.0, 9.0])
+        self.assertEqual(c.item(), 7.0)
+
+        # The gradient is the seed constant of the backward pass
+        fun = mx.compile(mx.value_and_grad(lambda x: x + 1.0))
+        value, grad = fun(mx.array(2.0))
+        self.assertEqual(value.item(), 3.0)
+        self.assertEqual(grad.item(), 1.0)
+
     def test_compile_with_constant(self):
         # Test float
         @partial(mx.compile)