diff --git a/autoarray/inversion/mesh/mesh_geometry/delaunay.py b/autoarray/inversion/mesh/mesh_geometry/delaunay.py
index eadbbfbd6..093fa6172 100644
--- a/autoarray/inversion/mesh/mesh_geometry/delaunay.py
+++ b/autoarray/inversion/mesh/mesh_geometry/delaunay.py
@@ -1,5 +1,4 @@
 import numpy as np
-import scipy.spatial
 from typing import Tuple
 
 from autonerves import cached_property
@@ -147,6 +146,7 @@ def neighbors(self) -> Neighbors:
         The neighbors of a Voronoi mesh are computed using the `ridge_points` attribute of the scipy `Voronoi`
         object, as described in the method `voronoi_neighbors_from`.
         """
+        import scipy.spatial
 
         delaunay = scipy.spatial.Delaunay(self.mesh_grid_xy)
 
diff --git a/autoarray/operators/coarse_interp_util.py b/autoarray/operators/coarse_interp_util.py
index e9477babc..58c9b3a60 100644
--- a/autoarray/operators/coarse_interp_util.py
+++ b/autoarray/operators/coarse_interp_util.py
@@ -17,7 +17,6 @@
 """
 
 import numpy as np
-from scipy.sparse import csr_matrix
 
 from autoarray import exc
 from autoarray import numba_util
@@ -242,6 +241,8 @@ def coarse_interp_matrix_from(
     A ``scipy.sparse.csr_matrix`` of shape
     [n_unmasked_fine_pixels, n_unmasked_coarse_pixels].
     """
+    from scipy.sparse import csr_matrix
+
     mask_itp_box = np.asarray(mask_itp_box)
     if np.count_nonzero(~mask_itp_box) == 0:
         raise exc.MeshException(
diff --git a/autoarray/operators/derivative_util.py b/autoarray/operators/derivative_util.py
index 62bf96471..b32606949 100644
--- a/autoarray/operators/derivative_util.py
+++ b/autoarray/operators/derivative_util.py
@@ -27,7 +27,6 @@
 """
 
 import numpy as np
-from scipy.sparse import csr_matrix
 
 from autoarray import exc
 from autoarray import numba_util
@@ -297,9 +296,11 @@ def derivative_1st_operators_from(mask, pixel_scale: float = 1.0):
     -------
     The (Hy, Hx) operators as ``scipy.sparse.csr_matrix``.
     """
+    from scipy.sparse import csr_matrix
+
     mask, diff_types = _diff_types_of_cleaned_mask_from(mask)
-    rows_hx, cols_hx, data_hx, rows_hy, cols_hy, data_hy = (
-        derivative_1st_triplets_from(mask, diff_types, dpix=pixel_scale)
+    rows_hx, cols_hx, data_hx, rows_hy, cols_hy, data_hy = derivative_1st_triplets_from(
+        mask, diff_types, dpix=pixel_scale
     )
 
     n_unmasked = np.count_nonzero(~mask)
@@ -443,6 +444,8 @@ def derivative_2nd_operators_from(mask, pixel_scale: float = 1.0):
     -------
     The (Hyy, Hxx) operators as ``scipy.sparse.csr_matrix``.
     """
+    from scipy.sparse import csr_matrix
+
     mask, diff_types = _diff_types_of_cleaned_mask_from(mask)
     rows_hxx, cols_hxx, data_hxx, rows_hyy, cols_hyy, data_hyy = (
         derivative_2nd_triplets_from(mask, diff_types, dpix=pixel_scale)
@@ -574,6 +577,8 @@ def forward_difference_operators_from(
     -------
     The (Hy, Hx) operators as ``scipy.sparse.csr_matrix``.
     """
+    from scipy.sparse import csr_matrix
+
     if max_order not in (1, 2, 3, 4):
         raise ValueError(f"max_order must be in 1..4, got {max_order}")