diff --git a/pytensor/link/mlx/dispatch/subtensor.py b/pytensor/link/mlx/dispatch/subtensor.py
index 6fded559df..482695b235 100644
--- a/pytensor/link/mlx/dispatch/subtensor.py
+++ b/pytensor/link/mlx/dispatch/subtensor.py
@@ -26,8 +26,27 @@ def subtensor(x, *ilists):
 
 @mlx_funcify.register(AdvancedSubtensor)
 def mlx_funcify_AdvancedSubtensor(op, node, **kwargs):
-    def advanced_subtensor(x, *ilists):
-        indices = indices_from_subtensor(ilists, op.idx_list)
+    # MLX slices reject array-typed bounds, so every flat index position a slice
+    # uses as a bound is coerced to a Python int; array-valued indices are left
+    # alone.
+    bound_positions = {
+        bound
+        for idx_entry in op.idx_list
+        if isinstance(idx_entry, slice)
+        for bound in (idx_entry.start, idx_entry.stop, idx_entry.step)
+        if bound is not None
+    }
+
+    def advanced_subtensor(
+        x, *ilists, idx_list=op.idx_list, bound_positions=bound_positions
+    ):
+        indices = indices_from_subtensor(
+            tuple(
+                int(index) if position in bound_positions else index
+                for position, index in enumerate(ilists)
+            ),
+            idx_list,
+        )
         if len(indices) == 1:
             indices = indices[0]
 
@@ -96,9 +115,40 @@ def mlx_fn(x, indices, y):
         def mlx_fn(x, indices, y):
             return x.at[indices].add(y)
 
-    def advancedincsubtensor(x, y, *ilist, mlx_fn=mlx_fn):
+    # MLX slices reject array-typed bounds, so every flat index position a slice
+    # uses as a bound is coerced to a Python int; array-valued indices are left
+    # alone.
+    bound_positions = {
+        bound
+        for idx_entry in op.idx_list
+        if isinstance(idx_entry, slice)
+        for bound in (idx_entry.start, idx_entry.stop, idx_entry.step)
+        if bound is not None
+    }
+
+    def advancedincsubtensor(
+        x,
+        y,
+        *ilist,
+        mlx_fn=mlx_fn,
+        idx_list=op.idx_list,
+        bound_positions=bound_positions,
+    ):
         op._check_runtime_broadcast_of_vector_index(node, x, y, ilist[0])
 
-        return mlx_fn(x, ilist, y)
+        # Slices and plain integers live in `idx_list`, not in `ilist`, and have
+        # to be spliced back in; without them `x[:, idx] = y` would scatter
+        # along the leading axis instead.
+        indices = indices_from_subtensor(
+            tuple(
+                int(index) if position in bound_positions else index
+                for position, index in enumerate(ilist)
+            ),
+            idx_list,
+        )
+        if len(indices) == 1:
+            indices = indices[0]
+
+        return mlx_fn(x, indices, y)
 
     return advancedincsubtensor
diff --git a/pytensor/link/mlx/dispatch/tensor_basic.py b/pytensor/link/mlx/dispatch/tensor_basic.py
index 3cdc47323f..450808352d 100644
--- a/pytensor/link/mlx/dispatch/tensor_basic.py
+++ b/pytensor/link/mlx/dispatch/tensor_basic.py
@@ -11,6 +11,7 @@
     Eye,
     Join,
     MakeVector,
+    Nonzero,
     ScalarFromTensor,
     Split,
     TensorFromScalar,
@@ -206,6 +207,23 @@ def arange(*_args):
     return arange
 
 
+@mlx_funcify.register(Nonzero)
+def mlx_funcify_Nonzero(op, node, **kwargs):
+    ndim = node.inputs[0].type.ndim
+
+    def nonzero(a):
+        # How many entries are selected is data-dependent, so `a` has to be
+        # counted on the host. `MLXLinker` keeps graphs containing `Nonzero` out
+        # of `mx.compile`, which forbids the evaluation this needs.
+        indices = np.nonzero(np.asarray(a))
+        if ndim == 1:
+            return mx.array(indices[0])
+
+        return [mx.array(index) for index in indices]
+
+    return nonzero
+
+
 def _extract_static_dims(shape_inputs):
     static_dims = []
     for dim in shape_inputs:
diff --git a/pytensor/link/mlx/linker.py b/pytensor/link/mlx/linker.py
index 9d662308cf..1f54dcdaea 100644
--- a/pytensor/link/mlx/linker.py
+++ b/pytensor/link/mlx/linker.py
@@ -1,6 +1,41 @@
+import warnings
+
 from pytensor.link.basic import JITLinker
 
 
+UNCOMPILABLE_SHAPE_WARNING = (
+    "This graph contains an `Op` whose output shape depends on its input "
+    "values, which `mx.compile` cannot trace, so the whole graph will run "
+    "uncompiled. To get compilation back, give the graph static shapes: a "
+    "boolean mask that is constant becomes integer indices, and a mask only "
+    "known at runtime can often be replaced by shape-preserving arithmetic, "
+    "e.g. `(x * mask).sum()` rather than `x[mask].sum()`."
+)
+
+
+def _has_data_dependent_shape(fgraph):
+    """Whether any `Op` in ``fgraph`` has an output shape that depends on data.
+
+    `mx.compile` traces a graph once per input *shape*, so an output whose shape
+    is only known from the values cannot be traced at all.
+    """
+    # Imported here because `pytensor.compile` imports this module while
+    # `pytensor.tensor` is still initializing.
+    from pytensor.graph.op import HasInnerGraph
+    from pytensor.tensor.basic import Nonzero
+
+    for node in fgraph.apply_nodes:
+        if isinstance(node.op, Nonzero):
+            return True
+
+        if isinstance(node.op, HasInnerGraph):
+            inner_fgraph = getattr(node.op, "fgraph", None)
+            if inner_fgraph is not None and _has_data_dependent_shape(inner_fgraph):
+                return True
+
+    return False
+
+
 class MLXLinker(JITLinker):
     """A `Linker` that JIT-compiles NumPy-based operations using Apple's MLX."""
 
@@ -44,7 +79,12 @@ def jit_compile(self, fn):
 
         from pytensor.link.mlx.dispatch import mlx_typify
 
-        if not self.use_compile:
+        use_compile = self.use_compile
+        if use_compile and _has_data_dependent_shape(self.fgraph):
+            warnings.warn(UNCOMPILABLE_SHAPE_WARNING, UserWarning, stacklevel=2)
+            use_compile = False
+
+        if not use_compile:
             # Skip compilation and just return the function with MLX typification
             def fn_no_compile(*inputs):
                 return fn(*(mlx_typify(inp) for inp in inputs))
diff --git a/pytensor/tensor/rewriting/jax.py b/pytensor/tensor/rewriting/jax.py
index 6e0ee1c1d6..2b4c42d4b0 100644
--- a/pytensor/tensor/rewriting/jax.py
+++ b/pytensor/tensor/rewriting/jax.py
@@ -48,7 +48,10 @@ def boolean_indexing_set_or_inc(fgraph, node):
     "jax_boolean_indexing_set_or_inc",
     dfs_rewriter(boolean_indexing_set_or_inc),
     "jax",
-    position=100,
+    "mlx",
+    # Must run before `specialize`: `bool_idx_to_nonzero` lives there and rewrites
+    # the mask into integer indices, leaving nothing for this to match.
+    position=1.95,
 )
 
 
@@ -96,7 +99,12 @@ def boolean_indexing_sum(fgraph, node):
 
 
 optdb.register(
-    "jax_boolean_indexing_sum", dfs_rewriter(boolean_indexing_sum), "jax", position=100
+    "jax_boolean_indexing_sum",
+    dfs_rewriter(boolean_indexing_sum),
+    "jax",
+    "mlx",
+    # Before `specialize`, for the same reason as `boolean_indexing_set_or_inc`.
+    position=1.95,
 )
 
 
diff --git a/pytensor/tensor/rewriting/subtensor.py b/pytensor/tensor/rewriting/subtensor.py
index 189375a2eb..ebee7970cb 100644
--- a/pytensor/tensor/rewriting/subtensor.py
+++ b/pytensor/tensor/rewriting/subtensor.py
@@ -2947,6 +2947,7 @@ def bool_idx_to_nonzero(fgraph, node):
     bool_idx_to_nonzero.__name__,
     bool_idx_to_nonzero,
     "numba",
+    "mlx",
     "shape_unsafe",  # It can mask invalid mask sizes
     use_db_name_as_tag=False,  # Not included if only "specialize" is requested
 )
diff --git a/tests/link/mlx/test_subtensor.py b/tests/link/mlx/test_subtensor.py
index c5bd707c35..9461b65bef 100644
--- a/tests/link/mlx/test_subtensor.py
+++ b/tests/link/mlx/test_subtensor.py
@@ -7,11 +7,16 @@
 from pytensor.compile.mode import Mode
 from pytensor.tensor import subtensor as pt_subtensor
 from pytensor.tensor import tensor
+from pytensor.tensor.basic import Nonzero
 from tests.link.mlx.test_basic import compare_mlx_and_py, mlx_mode, py_mode
 
 
 mx = pytest.importorskip("mlx.core")
 
+# `bool_idx_to_nonzero`, which converts boolean masks into the integer indices MLX
+# supports, lives in `specialize`, which `mlx_mode` does not include.
+bool_idx_mode = mlx_mode.including("specialize")
+
 
 def test_mlx_Subtensor_basic():
     """Test basic subtensor operations with constant indices."""
@@ -79,9 +84,6 @@ def test_mlx_AdvancedSubtensor():
     compare_mlx_and_py([x_pt], [out_pt], [x_np])
 
 
-@pytest.mark.xfail(
-    raises=ValueError, reason="MLX does not support boolean indexing yet"
-)
 def test_mlx_AdvancedSubtensor_boolean():
     """Test advanced subtensor operations with boolean indexing."""
     shape = (3, 4, 5)
@@ -92,7 +94,7 @@ def test_mlx_AdvancedSubtensor_boolean():
     bool_mask = np.array([True, False, True])
     out_pt = x_pt[bool_mask]
     assert isinstance(out_pt.owner.op, pt_subtensor.AdvancedSubtensor)
-    compare_mlx_and_py([x_pt], [out_pt], [x_np])
+    compare_mlx_and_py([x_pt], [out_pt], [x_np], mlx_mode=bool_idx_mode)
 
 
 def test_mlx_IncSubtensor_set():
@@ -347,3 +349,186 @@ def test_mlx_AdvancedIncSubtensor_ignore_duplicates():
     assert out.owner.op.ignore_duplicates
 
     compare_mlx_and_py([x], [out], [np.zeros(3, dtype=np.float32)])
+
+
+@pytest.mark.parametrize("set_instead_of_inc", [True, False], ids=["set", "inc"])
+@pytest.mark.parametrize(
+    "index, y_shape",
+    [
+        ((slice(None), np.array([0, 1, 2])), (6, 3)),
+        ((np.array([0, 1, 2]), slice(None)), (3, 6)),
+        ((slice(None), np.array([0, 0, 1])), (6, 3)),
+        ((1, np.array([0, 1, 2])), (3,)),
+    ],
+    ids=["slice-first", "slice-last", "duplicate-columns", "int-and-array"],
+)
+def test_mlx_AdvancedIncSubtensor_mixed_indices(index, y_shape, set_instead_of_inc):
+    """Slices and integers in ``idx_list`` must be spliced back into the scatter.
+
+    Without them the update lands on the leading axis instead of the indexed one.
+    """
+    x = pt.matrix("x", shape=(6, 6), dtype="float32")
+    y = tensor("y", shape=y_shape, dtype="float32")
+
+    inc_fn = (
+        pt_subtensor.set_subtensor if set_instead_of_inc else pt_subtensor.inc_subtensor
+    )
+    out = inc_fn(x[index], y)
+    assert isinstance(out.owner.op, pt_subtensor.AdvancedIncSubtensor)
+    assert out.owner.op.set_instead_of_inc == set_instead_of_inc
+
+    x_np = np.zeros((6, 6), dtype="float32")
+    y_np = np.arange(np.prod(y_shape), dtype="float32").reshape(y_shape)
+    compare_mlx_and_py([x, y], [out], [x_np, y_np])
+
+
+def test_mlx_AdvancedIncSubtensor_batched_diagonal_grad():
+    """The adjoint of ``x[..., i, i]`` scatters through an ellipsis of slices."""
+    x = pt.tensor("x", shape=(5, 3, 3), dtype="float32")
+    idx = pt.arange(3)
+    grad = pytensor.grad(x[..., idx, idx].sum(), x)
+
+    compare_mlx_and_py([x], [grad], [np.zeros((5, 3, 3), dtype="float32")])
+
+
+@pytest.mark.parametrize(
+    "index",
+    [
+        np.array([True, False, True, False, True, False]),
+        (slice(None), np.array([True, False, True, False, True, False])),
+        (np.array([True, False, True, False, True, False]), slice(1, 4)),
+        np.zeros((6,), dtype=bool),
+        np.arange(36).reshape(6, 6) % 3 == 0,
+    ],
+    ids=["rows", "columns", "rows-and-slice", "all-false", "2d-mask"],
+)
+def test_mlx_AdvancedSubtensor_boolean_mask(index):
+    """MLX has no boolean gather; `bool_idx_to_nonzero` supplies integer indices."""
+    x = pt.matrix("x", shape=(6, 6), dtype="float32")
+    out = x[index]
+    assert isinstance(out.owner.op, pt_subtensor.AdvancedSubtensor)
+
+    x_np = np.arange(36, dtype="float32").reshape(6, 6)
+    compare_mlx_and_py([x], [out], [x_np], mlx_mode=bool_idx_mode)
+
+
+@pytest.mark.parametrize("set_instead_of_inc", [True, False], ids=["set", "inc"])
+@pytest.mark.parametrize(
+    "index, y_shape",
+    [
+        (np.array([True, False, True, False, True, False]), (3, 6)),
+        (np.arange(36).reshape(6, 6) % 3 == 0, (12,)),
+    ],
+    ids=["rows", "2d-mask"],
+)
+def test_mlx_AdvancedIncSubtensor_boolean_mask(index, y_shape, set_instead_of_inc):
+    """MLX has no boolean scatter; `bool_idx_to_nonzero` supplies integer indices."""
+    x = pt.matrix("x", shape=(6, 6), dtype="float32")
+    y = tensor("y", shape=y_shape, dtype="float32")
+
+    inc_fn = (
+        pt_subtensor.set_subtensor if set_instead_of_inc else pt_subtensor.inc_subtensor
+    )
+    out = inc_fn(x[index], y)
+    assert isinstance(out.owner.op, pt_subtensor.AdvancedIncSubtensor)
+
+    x_np = np.arange(36, dtype="float32").reshape(6, 6)
+    y_np = np.arange(np.prod(y_shape), dtype="float32").reshape(y_shape)
+    compare_mlx_and_py([x, y], [out], [x_np, y_np], mlx_mode=bool_idx_mode)
+
+
+@pytest.mark.parametrize(
+    "index",
+    [
+        (np.array([0, 2, 4]), slice(1, 4)),
+        (np.array([0, 2, 4]), slice(None, None, 2)),
+        (slice(1, 4), np.array([0, 2, 4])),
+    ],
+    ids=["stop-bound", "step-bound", "slice-first"],
+)
+def test_mlx_AdvancedSubtensor_slice_bounds(index):
+    """Slice bounds arrive as MLX scalars and must be coerced to Python ints."""
+    x = pt.matrix("x", shape=(6, 6), dtype="float32")
+    out = x[index]
+    assert isinstance(out.owner.op, pt_subtensor.AdvancedSubtensor)
+
+    compare_mlx_and_py([x], [out], [np.arange(36, dtype="float32").reshape(6, 6)])
+
+
+@pytest.mark.parametrize(
+    "index_fn, y_shape",
+    [
+        (lambda x, mask: x[mask], (3, 6)),
+        (lambda x, mask: x[:, mask], (6, 3)),
+        (lambda x, mask: x[mask, 1:4], (3, 3)),
+    ],
+    ids=["rows", "columns", "rows-and-slice"],
+)
+def test_mlx_runtime_boolean_mask(index_fn, y_shape):
+    """A mask only known at runtime needs `Nonzero`, so the linker runs it eagerly."""
+    x = pt.matrix("x", shape=(6, 6), dtype="float32")
+    mask = pt.vector("mask", shape=(6,), dtype="bool")
+    y = tensor("y", shape=y_shape, dtype="float32")
+
+    x_np = np.arange(36, dtype="float32").reshape(6, 6)
+    mask_np = np.array([True, False, True, False, True, False])
+    y_np = np.arange(np.prod(y_shape), dtype="float32").reshape(y_shape)
+
+    compare_mlx_and_py(
+        [x, mask],
+        [index_fn(x, mask)],
+        [x_np, mask_np],
+        mlx_mode=bool_idx_mode,
+    )
+    compare_mlx_and_py(
+        [x, mask, y],
+        [pt_subtensor.set_subtensor(index_fn(x, mask), y)],
+        [x_np, mask_np, y_np],
+        mlx_mode=bool_idx_mode,
+    )
+
+
+def test_mlx_runtime_boolean_mask_grad():
+    """`Nonzero` has to survive being differentiated through and scattered back."""
+    x = pt.matrix("x", shape=(6, 6), dtype="float32")
+    mask = pt.matrix("mask", shape=(6, 6), dtype="bool")
+    grad = pytensor.grad((x[mask] ** 2).sum(), x)
+
+    compare_mlx_and_py(
+        [x, mask],
+        [grad],
+        [
+            np.arange(36, dtype="float32").reshape(6, 6),
+            np.arange(36).reshape(6, 6) % 3 == 0,
+        ],
+        mlx_mode=bool_idx_mode,
+    )
+
+
+@pytest.mark.filterwarnings("error")  # The eager fallback warns; the rewrite avoids it
+def test_mlx_boolean_mask_sum_reexpressible():
+    """`boolean_indexing_sum` rewrites a masked sum into `where`, which compiles."""
+    x = pt.matrix("x", shape=(5, 5), dtype="float32")
+    out = x[x < 0].sum()
+
+    x_np = np.arange(-12, 13, dtype="float32").reshape(5, 5)
+    fn, _ = compare_mlx_and_py([x], [out], [x_np], mlx_mode=bool_idx_mode)
+
+    assert not any(isinstance(node.op, Nonzero) for node in fn.maker.fgraph.apply_nodes)
+
+
+@pytest.mark.parametrize("set_instead_of_inc", [True, False], ids=["set", "inc"])
+@pytest.mark.filterwarnings("error")  # The eager fallback warns; the rewrite avoids it
+def test_mlx_boolean_mask_scalar_update_reexpressible(set_instead_of_inc):
+    """`boolean_indexing_set_or_inc` rewrites a scalar masked update into `where`."""
+    x = pt.matrix("x", shape=(4, 5), dtype="float32")
+    inc_fn = (
+        pt_subtensor.set_subtensor if set_instead_of_inc else pt_subtensor.inc_subtensor
+    )
+    out = inc_fn(x[x > 0], np.float32(0.0 if set_instead_of_inc else 1.0))
+    assert isinstance(out.owner.op, pt_subtensor.AdvancedIncSubtensor)
+
+    x_np = np.arange(-10, 10, dtype="float32").reshape(4, 5)
+    fn, _ = compare_mlx_and_py([x], [out], [x_np], mlx_mode=bool_idx_mode)
+
+    assert not any(isinstance(node.op, Nonzero) for node in fn.maker.fgraph.apply_nodes)
diff --git a/tests/link/mlx/test_tensor_basic.py b/tests/link/mlx/test_tensor_basic.py
index 8ab3e07542..a086841172 100644
--- a/tests/link/mlx/test_tensor_basic.py
+++ b/tests/link/mlx/test_tensor_basic.py
@@ -164,3 +164,49 @@ def test_arange():
     out = arange(1, 10, 2)
 
     compare_mlx_and_py([], [out], [])
+
+
+@pytest.mark.parametrize("ndim", [1, 2], ids=["vector", "matrix"])
+def test_nonzero(ndim):
+    shape = (6,) * ndim
+    x = pt.tensor("x", shape=shape, dtype="float32")
+    out = pt.stack(pt.nonzero(x)) if ndim > 1 else pt.nonzero(x)[0]
+
+    x_np = (np.arange(np.prod(shape), dtype="float32") % 3).reshape(shape)
+    compare_mlx_and_py([x], [out], [x_np], mlx_mode=mlx_mode_no_compile)
+
+
+def test_nonzero_skips_compilation():
+    """A data-dependent output shape cannot be traced, so the graph runs eagerly."""
+    x = pt.vector("x", shape=(6,), dtype="float32")
+    x_np = np.array([0, 1, 0, 2, 0, 3], dtype="float32")
+
+    with pytest.warns(UserWarning, match="output shape depends on its input values"):
+        compiled_f = pytensor.function([x], pt.nonzero(x)[0], mode=compile_mode)
+
+    np.testing.assert_array_equal(compiled_f(x_np), [1, 3, 5])
+
+
+def test_nonzero_inside_inner_graph_skips_compilation():
+    """The check has to see into inner graphs, where an `Op` can hide a `Nonzero`."""
+    from pytensor.compile.builders import OpFromGraph
+
+    a = pt.vector("a", shape=(6,), dtype="float32")
+    x = pt.vector("x", shape=(6,), dtype="float32")
+    out = OpFromGraph([a], [pt.nonzero(a)[0]], inline=False)(x)
+
+    with pytest.warns(UserWarning, match="output shape depends on its input values"):
+        f = pytensor.function([x], out, mode=compile_mode)
+
+    assert any(isinstance(node.op, OpFromGraph) for node in f.maker.fgraph.apply_nodes)
+    np.testing.assert_array_equal(
+        f(np.array([0, 1, 0, 2, 0, 3], dtype="float32")), [1, 3, 5]
+    )
+
+
+@pytest.mark.filterwarnings("error")  # Raise if we did not expect to skip compilation
+def test_compiles_without_data_dependent_shape():
+    x = pt.vector("x", shape=(6,), dtype="float32")
+    f = pytensor.function([x], x * 2, mode=compile_mode)
+
+    np.testing.assert_allclose(f(np.ones(6, dtype="float32")), 2.0)