diff --git a/autoarray/dataset/plot/imaging_plots.py b/autoarray/dataset/plot/imaging_plots.py
index 095c97045..53f6499c4 100644
--- a/autoarray/dataset/plot/imaging_plots.py
+++ b/autoarray/dataset/plot/imaging_plots.py
@@ -1,8 +1,7 @@
 from typing import Optional
 
-import matplotlib.pyplot as plt
 
-from autoarray.plot.utils import subplot_save, conf_subplot_figsize, tight_layout
+from autoarray.plot.utils import subplots, subplot_save, conf_subplot_figsize, tight_layout
 
 
 def subplot_imaging_dataset(
@@ -51,7 +50,7 @@ def subplot_imaging_dataset(
 
     from autoarray.plot.array import plot_array
 
-    fig, axes = plt.subplots(3, 3, figsize=conf_subplot_figsize(3, 3))
+    fig, axes = subplots(3, 3, figsize=conf_subplot_figsize(3, 3))
     axes = axes.flatten()
 
     plot_array(
@@ -172,7 +171,7 @@ def subplot_imaging_dataset_list(
     from autoarray.plot.array import plot_array
 
     n = len(dataset_list)
-    fig, axes = plt.subplots(n, 3, figsize=conf_subplot_figsize(n, 3))
+    fig, axes = subplots(n, 3, figsize=conf_subplot_figsize(n, 3))
     if n == 1:
         axes = [axes]
     for i, dataset in enumerate(dataset_list):
diff --git a/autoarray/dataset/plot/interferometer_plots.py b/autoarray/dataset/plot/interferometer_plots.py
index 87c751cd5..359d55b1e 100644
--- a/autoarray/dataset/plot/interferometer_plots.py
+++ b/autoarray/dataset/plot/interferometer_plots.py
@@ -1,12 +1,11 @@
 import numpy as np
 from typing import Optional
 
-import matplotlib.pyplot as plt
 
 from autoarray.plot.array import plot_array
 from autoarray.plot.grid import plot_grid
 from autoarray.plot.yx import plot_yx
-from autoarray.plot.utils import subplot_save, hide_unused_axes, conf_subplot_figsize, tight_layout
+from autoarray.plot.utils import subplots, subplot_save, hide_unused_axes, conf_subplot_figsize, tight_layout
 from autoarray.structures.grids.irregular_2d import Grid2DIrregular
 
 
@@ -39,7 +38,7 @@ def subplot_interferometer_dataset(
     use_log10
         Apply log10 normalisation to image panels.
     """
-    fig, axes = plt.subplots(2, 3, figsize=conf_subplot_figsize(2, 3))
+    fig, axes = subplots(2, 3, figsize=conf_subplot_figsize(2, 3))
     axes = axes.flatten()
 
     plot_grid(dataset.data.in_grid, ax=axes[0], title="Visibilities", xlabel="", ylabel="")
@@ -117,7 +116,7 @@ def subplot_interferometer_dirty_images(
     use_log10
         Apply log10 normalisation.
     """
-    fig, axes = plt.subplots(1, 3, figsize=conf_subplot_figsize(1, 3))
+    fig, axes = subplots(1, 3, figsize=conf_subplot_figsize(1, 3))
 
     plot_array(
         dataset.dirty_image,
diff --git a/autoarray/fit/plot/fit_imaging_plots.py b/autoarray/fit/plot/fit_imaging_plots.py
index 298d6ccd3..ae474e695 100644
--- a/autoarray/fit/plot/fit_imaging_plots.py
+++ b/autoarray/fit/plot/fit_imaging_plots.py
@@ -1,9 +1,8 @@
 from typing import Optional
 
-import matplotlib.pyplot as plt
 
 from autoarray.plot.array import plot_array
-from autoarray.plot.utils import subplot_save, symmetric_vmin_vmax, hide_unused_axes, conf_subplot_figsize, tight_layout
+from autoarray.plot.utils import subplots, subplot_save, symmetric_vmin_vmax, hide_unused_axes, conf_subplot_figsize, tight_layout
 
 
 def subplot_fit_imaging(
@@ -43,7 +42,7 @@ def subplot_fit_imaging(
     grid, positions, lines
         Optional overlays forwarded to every panel.
     """
-    fig, axes = plt.subplots(2, 3, figsize=conf_subplot_figsize(2, 3))
+    fig, axes = subplots(2, 3, figsize=conf_subplot_figsize(2, 3))
     axes = axes.flatten()
 
     plot_array(
diff --git a/autoarray/fit/plot/fit_interferometer_plots.py b/autoarray/fit/plot/fit_interferometer_plots.py
index 15f29fe08..2cc805b95 100644
--- a/autoarray/fit/plot/fit_interferometer_plots.py
+++ b/autoarray/fit/plot/fit_interferometer_plots.py
@@ -1,11 +1,10 @@
 import numpy as np
 from typing import Optional
 
-import matplotlib.pyplot as plt
 
 from autoarray.plot.array import plot_array
 from autoarray.plot.yx import plot_yx
-from autoarray.plot.utils import subplot_save, symmetric_vmin_vmax, hide_unused_axes, conf_subplot_figsize, tight_layout
+from autoarray.plot.utils import subplots, subplot_save, symmetric_vmin_vmax, hide_unused_axes, conf_subplot_figsize, tight_layout
 
 
 def subplot_fit_interferometer(
@@ -40,7 +39,7 @@ def subplot_fit_interferometer(
         Not used here (UV-plane residuals are scatter plots); kept for API
         consistency.
     """
-    fig, axes = plt.subplots(2, 3, figsize=conf_subplot_figsize(2, 3))
+    fig, axes = subplots(2, 3, figsize=conf_subplot_figsize(2, 3))
     axes = axes.flatten()
 
     uv = fit.dataset.uv_distances / 10**3.0
@@ -135,7 +134,7 @@ def subplot_fit_interferometer_dirty_images(
     residuals_symmetric_cmap
         Centre residual colour scale symmetrically around zero.
     """
-    fig, axes = plt.subplots(2, 3, figsize=conf_subplot_figsize(2, 3))
+    fig, axes = subplots(2, 3, figsize=conf_subplot_figsize(2, 3))
     axes = axes.flatten()
 
     plot_array(
diff --git a/autoarray/inversion/inversion/inversion_util.py b/autoarray/inversion/inversion/inversion_util.py
index 7c552c23a..40045ed90 100644
--- a/autoarray/inversion/inversion/inversion_util.py
+++ b/autoarray/inversion/inversion/inversion_util.py
@@ -5,7 +5,6 @@
 from autoarray.settings import Settings
 
 from autoarray import exc
-from autoarray.util.fnnls import fnnls_cholesky
 
 
 def curvature_matrix_diag_via_psf_weighted_noise_from(
@@ -281,6 +280,8 @@ def reconstruction_positive_only_from(
 
     try:
 
+        from autoarray.util.fnnls import fnnls_cholesky
+
         return fnnls_cholesky(
             curvature_reg_matrix,
             (data_vector).T,
diff --git a/autoarray/inversion/mesh/interpolator/delaunay.py b/autoarray/inversion/mesh/interpolator/delaunay.py
index 5d2410370..ea491a2b1 100644
--- a/autoarray/inversion/mesh/interpolator/delaunay.py
+++ b/autoarray/inversion/mesh/interpolator/delaunay.py
@@ -1,6 +1,4 @@
 import numpy as np
-import scipy.spatial
-from scipy.spatial import cKDTree, Delaunay, Voronoi
 
 from autoconf import cached_property
 
@@ -12,6 +10,7 @@
 
 def scipy_delaunay(points_np, query_points_np, areas_factor):
     """Compute Delaunay simplices (simplices_padded) and Voronoi areas in one call."""
+    from scipy.spatial import Delaunay
 
     max_simplices = 2 * points_np.shape[0]
 
@@ -182,6 +181,8 @@ def pix_indexes_for_sub_slim_index_delaunay_from(
     # Case 2: Outside → KDTree NN
     # ---------------------------
     if outside_mask.any():
+        from scipy.spatial import cKDTree
+
         tree = cKDTree(delaunay_points)
         _, idx = tree.query(data_grid[outside_mask], k=1)
         out[outside_mask, 0] = idx.astype(np.int32)
@@ -202,6 +203,7 @@ def scipy_delaunay_matern(points_np, query_points_np):
         typically of shape (Q, 3), where each row gives the indices of the
         Delaunay mesh vertices ("pixels") associated with that query point.
     """
+    from scipy.spatial import Delaunay
 
     max_simplices = 2 * points_np.shape[0]
 
diff --git a/autoarray/inversion/plot/inversion_plots.py b/autoarray/inversion/plot/inversion_plots.py
index e6e55ac53..0187caa4c 100644
--- a/autoarray/inversion/plot/inversion_plots.py
+++ b/autoarray/inversion/plot/inversion_plots.py
@@ -4,12 +4,11 @@
 from pathlib import Path
 from typing import Optional, Union
 
-import matplotlib.pyplot as plt
 from autoconf import conf
 
 from autoarray.inversion.mappers.abstract import Mapper
 from autoarray.plot.array import plot_array
-from autoarray.plot.utils import numpy_grid, numpy_lines, numpy_positions, subplot_save, hide_unused_axes, conf_subplot_figsize, tight_layout
+from autoarray.plot.utils import subplots, numpy_grid, numpy_lines, numpy_positions, subplot_save, hide_unused_axes, conf_subplot_figsize, tight_layout
 from autoarray.inversion.plot.mapper_plots import plot_mapper
 from autoarray.structures.arrays.uniform_2d import Array2D
 
@@ -53,7 +52,7 @@ def subplot_of_mapper(
     """
     mapper = inversion.cls_list_from(cls=Mapper)[mapper_index]
 
-    fig, axes = plt.subplots(3, 4, figsize=conf_subplot_figsize(3, 4))
+    fig, axes = subplots(3, 4, figsize=conf_subplot_figsize(3, 4))
     axes = axes.flatten()
 
     # panel 0: data subtracted
@@ -279,7 +278,7 @@ def subplot_mappings(
     )
     mapper.slim_indexes_for_pix_indexes(pix_indexes=pix_indexes)
 
-    fig, axes = plt.subplots(2, 2, figsize=conf_subplot_figsize(2, 2))
+    fig, axes = subplots(2, 2, figsize=conf_subplot_figsize(2, 2))
     axes = axes.flatten()
 
     # panel 0: data subtracted
diff --git a/autoarray/inversion/plot/mapper_plots.py b/autoarray/inversion/plot/mapper_plots.py
index 088d10c51..6e3d9e6be 100644
--- a/autoarray/inversion/plot/mapper_plots.py
+++ b/autoarray/inversion/plot/mapper_plots.py
@@ -1,11 +1,10 @@
 import logging
 from typing import Optional
 
-import matplotlib.pyplot as plt
 
 from autoarray.plot.array import plot_array
 from autoarray.plot.inversion import plot_inversion_reconstruction
-from autoarray.plot.utils import numpy_grid, numpy_lines, subplot_save, conf_subplot_figsize, tight_layout
+from autoarray.plot.utils import subplots, numpy_grid, numpy_lines, subplot_save, conf_subplot_figsize, tight_layout
 
 logger = logging.getLogger(__name__)
 
@@ -114,7 +113,7 @@ def subplot_image_and_mapper(
     lines
         Lines to overlay on both panels.
     """
-    fig, axes = plt.subplots(1, 2, figsize=conf_subplot_figsize(1, 2))
+    fig, axes = subplots(1, 2, figsize=conf_subplot_figsize(1, 2))
 
     plot_array(
         image,
diff --git a/autoarray/mask/mask_2d_util.py b/autoarray/mask/mask_2d_util.py
index e5d1ccc11..c13b97b11 100644
--- a/autoarray/mask/mask_2d_util.py
+++ b/autoarray/mask/mask_2d_util.py
@@ -1,6 +1,5 @@
 import numpy as np
 import warnings
-from scipy.ndimage import binary_dilation
 from typing import Tuple
 
 from autoarray import exc
@@ -465,7 +464,6 @@ def min_false_distance_to_edge(mask: np.ndarray) -> Tuple[int, int]:
 from typing import Tuple
 
 import numpy as np
-from scipy.ndimage import binary_dilation
 
 
 def required_shape_for_kernel(
@@ -608,6 +606,8 @@ def blurring_mask_2d_from(
     )
 
     # Pixels within kernel footprint of any unmasked pixel
+    from scipy.ndimage import binary_dilation
+
     near_unmasked_padded = binary_dilation(unmasked_padded, structure=structure)
     near_unmasked = near_unmasked_padded[
         pad_y : pad_y + mask_2d.shape[0],
diff --git a/autoarray/operators/convolver.py b/autoarray/operators/convolver.py
index 2a48b57ae..b77057608 100644
--- a/autoarray/operators/convolver.py
+++ b/autoarray/operators/convolver.py
@@ -7,7 +7,6 @@
 
 import numpy as np
 from pathlib import Path
-import scipy
 from typing import Optional, Tuple, Union
 import warnings
 
@@ -116,6 +115,8 @@ class determines how masked real-space data are embedded into a padded array,
         full_shape = tuple(
             s1 + s2 - 1 for s1, s2 in zip(mask_shape, self.kernel.shape_native)
         )
+        import scipy.fft
+
         fft_shape = tuple(scipy.fft.next_fast_len(s, real=True) for s in full_shape)
 
         self.fft_shape = fft_shape
diff --git a/autoarray/plot/array.py b/autoarray/plot/array.py
index 39c0e7793..618c74eb1 100644
--- a/autoarray/plot/array.py
+++ b/autoarray/plot/array.py
@@ -5,11 +5,9 @@
 import os
 from typing import List, Optional, Tuple
 
-import matplotlib.pyplot as plt
 import numpy as np
-from matplotlib.colors import LogNorm, Normalize
-
 from autoarray.plot.utils import (
+    subplots,
     apply_extent,
     apply_labels,
     conf_figsize,
@@ -29,7 +27,7 @@
 
 def plot_array(
     array,
-    ax: Optional[plt.Axes] = None,
+    ax=None,
     # --- spatial metadata -------------------------------------------------------
     extent: Optional[Tuple[float, float, float, float]] = None,
     # --- overlays ---------------------------------------------------------------
@@ -158,7 +156,7 @@ def plot_array(
     owns_figure = ax is None
     if owns_figure:
         figsize = figsize or conf_figsize("figures")
-        fig, ax = plt.subplots(1, 1, figsize=figsize)
+        fig, ax = subplots(1, 1, figsize=figsize)
     else:
         fig = ax.get_figure()
 
@@ -181,8 +179,10 @@ def plot_array(
                 vmax_log = np.nanmax(clipped)
         if not np.isfinite(vmax_log) or vmax_log <= vmin_log:
             vmax_log = vmin_log * 10.0
+        from matplotlib.colors import LogNorm
         norm = LogNorm(vmin=vmin_log, vmax=vmax_log)
     elif vmin is not None or vmax is not None:
+        from matplotlib.colors import Normalize
         norm = Normalize(vmin=vmin, vmax=vmax)
     else:
         norm = None
diff --git a/autoarray/plot/grid.py b/autoarray/plot/grid.py
index 012303447..6a109f61d 100644
--- a/autoarray/plot/grid.py
+++ b/autoarray/plot/grid.py
@@ -6,10 +6,11 @@
 
 from typing import Iterable, List, Optional, Tuple
 
-import matplotlib.pyplot as plt
 import numpy as np
 
 from autoarray.plot.utils import (
+    subplots,
+    get_cmap,
     apply_extent,
     apply_labels,
     conf_figsize,
@@ -20,7 +21,7 @@
 
 def plot_grid(
     grid,
-    ax: Optional[plt.Axes] = None,
+    ax=None,
     # --- errors -----------------------------------------------------------------
     y_errors: Optional[np.ndarray] = None,
     x_errors: Optional[np.ndarray] = None,
@@ -108,13 +109,13 @@ def plot_grid(
     owns_figure = ax is None
     if owns_figure:
         figsize = figsize or conf_figsize("figures")
-        fig, ax = plt.subplots(1, 1, figsize=figsize)
+        fig, ax = subplots(1, 1, figsize=figsize)
     else:
         fig = ax.get_figure()
 
     # --- scatter / errorbar ----------------------------------------------------
     if color_array is not None:
-        cmap = plt.get_cmap(colormap)
+        cmap = get_cmap(colormap)
         colors = cmap((color_array - color_array.min()) / (np.ptp(color_array) or 1))
 
         if y_errors is None and x_errors is None:
diff --git a/autoarray/plot/inversion.py b/autoarray/plot/inversion.py
index fa58b7a69..cd85dcecf 100644
--- a/autoarray/plot/inversion.py
+++ b/autoarray/plot/inversion.py
@@ -6,17 +6,16 @@
 
 from typing import List, Optional, Tuple
 
-import matplotlib.pyplot as plt
 import numpy as np
 from matplotlib.colors import LogNorm, Normalize
 
-from autoarray.plot.utils import apply_extent, apply_labels, conf_figsize, save_figure, _conf_imshow_origin
+from autoarray.plot.utils import subplots, apply_extent, apply_labels, conf_figsize, save_figure, _conf_imshow_origin
 
 
 def plot_inversion_reconstruction(
     pixel_values: np.ndarray,
     mapper,
-    ax: Optional[plt.Axes] = None,
+    ax=None,
     # --- cosmetics --------------------------------------------------------------
     title: str = "Reconstruction",
     xlabel: str = 'x (")',
@@ -84,7 +83,7 @@ def plot_inversion_reconstruction(
     owns_figure = ax is None
     if owns_figure:
         figsize = figsize or conf_figsize("figures")
-        fig, ax = plt.subplots(1, 1, figsize=figsize)
+        fig, ax = subplots(1, 1, figsize=figsize)
     else:
         fig = ax.get_figure()
 
diff --git a/autoarray/plot/utils.py b/autoarray/plot/utils.py
index 8dcebc4b0..232f493b3 100644
--- a/autoarray/plot/utils.py
+++ b/autoarray/plot/utils.py
@@ -6,7 +6,6 @@
 import os
 from typing import List, Optional, Tuple
 
-import matplotlib.pyplot as plt
 import numpy as np
 
 logger = logging.getLogger(__name__)
@@ -15,6 +14,20 @@
 _FAST_PLOTS = os.environ.get("PYAUTO_FAST_PLOTS") == "1"
 
 
+def subplots(*args, **kwargs):
+    """Lazy wrapper around ``plt.subplots`` that defers the matplotlib import."""
+    import matplotlib.pyplot as plt
+
+    return plt.subplots(*args, **kwargs)
+
+
+def get_cmap(name):
+    """Lazy wrapper around ``plt.get_cmap`` that defers the matplotlib import."""
+    import matplotlib.pyplot as plt
+
+    return plt.get_cmap(name)
+
+
 def tight_layout():
     """Call ``plt.tight_layout()`` unless fast-plot mode is active.
 
@@ -24,6 +37,9 @@ def tight_layout():
     """
     if _FAST_PLOTS:
         return
+
+    import matplotlib.pyplot as plt
+
     plt.tight_layout()
 
 
@@ -281,6 +297,8 @@ def set_with_color_values(ax, cmap, color_values, norm=None):
 
 
 def _output_mode_save(fig, filename):
+    import matplotlib.pyplot as plt
+
     """If ``PYAUTOARRAY_OUTPUT_MODE=1``, save *fig* to a numbered file in
     ``./output_mode/<script_name>/`` and return ``True``.  Otherwise return
     ``False`` so the caller can proceed with normal saving.
@@ -335,6 +353,8 @@ def subplot_save(fig, output_path, output_filename, output_format=None):
         File format string, e.g. ``"png"`` or ``"pdf"``.  ``"show"``
         displays the figure interactively.  ``None`` reads from config.
     """
+    import matplotlib.pyplot as plt
+
     if output_format is None:
         output_format = _conf_output_format()
 
@@ -441,7 +461,7 @@ def conf_subplot_figsize(rows: int, cols: int) -> Tuple[int, int]:
 
 
 def apply_labels(
-    ax: plt.Axes,
+    ax,
     title: str = "",
     xlabel: str = "",
     ylabel: str = "",
@@ -474,7 +494,7 @@ def apply_labels(
 
 
 def save_figure(
-    fig: plt.Figure,
+    fig,
     path: str,
     filename: str,
     format: str = None,
@@ -505,6 +525,8 @@ def save_figure(
     dpi
         Resolution in dots per inch.
     """
+    import matplotlib.pyplot as plt
+
     if format is None:
         format = _conf_output_format()
 
@@ -543,7 +565,7 @@ def save_figure(
     plt.close(fig)
 
 
-def plot_visibilities_1d(vis, ax: plt.Axes, title: str = "") -> None:
+def plot_visibilities_1d(vis, ax, title: str = "") -> None:
     """Plot the real and imaginary components of a visibilities array as 1D line plots.
 
     Draws two overlapping lines — one for the real part and one for the
@@ -664,7 +686,7 @@ def _colorbar_tick_labels(tick_values: List[float], cb_unit: Optional[str] = Non
 
 def _apply_colorbar(
     mappable,
-    ax: plt.Axes,
+    ax,
     cb_unit: Optional[str] = None,
     is_subplot: bool = False,
 ) -> None:
@@ -679,6 +701,8 @@ def _apply_colorbar(
         When ``True`` uses ``labelsize_subplot`` from config (default 16) instead of
         the single-figure ``labelsize`` (default 16).
     """
+    import matplotlib.pyplot as plt
+
     tick_values = _colorbar_tick_values(getattr(mappable, "norm", None))
 
     cb = plt.colorbar(
@@ -706,7 +730,7 @@ def _apply_colorbar(
 
 
 def _apply_contours(
-    ax: plt.Axes,
+    ax,
     array: np.ndarray,
     extent,
     use_log10: bool = False,
@@ -879,7 +903,7 @@ def _arcsec_labels(ticks) -> List[str]:
 
 
 def apply_extent(
-    ax: plt.Axes,
+    ax,
     extent: Tuple[float, float, float, float],
 ) -> None:
     """
diff --git a/autoarray/plot/yx.py b/autoarray/plot/yx.py
index e80f54766..cb397daef 100644
--- a/autoarray/plot/yx.py
+++ b/autoarray/plot/yx.py
@@ -6,16 +6,15 @@
 
 from typing import List, Optional, Tuple
 
-import matplotlib.pyplot as plt
 import numpy as np
 
-from autoarray.plot.utils import apply_labels, conf_figsize, save_figure
+from autoarray.plot.utils import subplots, apply_labels, conf_figsize, save_figure
 
 
 def plot_yx(
     y,
     x=None,
-    ax: Optional[plt.Axes] = None,
+    ax=None,
     # --- errors / extras --------------------------------------------------------
     y_errors: Optional[np.ndarray] = None,
     x_errors: Optional[np.ndarray] = None,
@@ -91,7 +90,7 @@ def plot_yx(
     owns_figure = ax is None
     if owns_figure:
         figsize = figsize or conf_figsize("figures")
-        fig, ax = plt.subplots(1, 1, figsize=figsize)
+        fig, ax = subplots(1, 1, figsize=figsize)
     else:
         fig = ax.get_figure()
 
diff --git a/autoarray/util/cholesky_funcs.py b/autoarray/util/cholesky_funcs.py
index bd211eeb5..713237827 100644
--- a/autoarray/util/cholesky_funcs.py
+++ b/autoarray/util/cholesky_funcs.py
@@ -1,5 +1,4 @@
 import numpy as np
-from scipy import linalg
 import math
 import time
 from autoarray import numba_util
@@ -46,6 +45,8 @@ def _cholupdate(U, x):
 
 
 def cholinsert(U, index, x):
+    from scipy import linalg
+
     S = np.insert(np.insert(U, index, 0, axis=0), index, 0, axis=1)
 
     S[:index, index] = S12 = linalg.solve_triangular(
@@ -70,6 +71,8 @@ def cholinsertlast(U, x):
     As in current Cholesky scheme implemented in fnnls, we only use this kind of insertion, so I
         separate it out from the `cholinsert`.
     """
+    from scipy import linalg
+
     index = U.shape[0]
 
     S = np.insert(np.insert(U, index, 0, axis=0), index, 0, axis=1)
diff --git a/autoarray/util/fnnls.py b/autoarray/util/fnnls.py
index 3f49c1f2d..5e95d9b73 100644
--- a/autoarray/util/fnnls.py
+++ b/autoarray/util/fnnls.py
@@ -1,5 +1,4 @@
 import numpy as np
-from scipy import linalg as slg
 
 from autoarray.util.cholesky_funcs import cholinsertlast, choldeleteindexes
 
@@ -27,6 +26,7 @@ def fnnls_cholesky(
     """
     Similar to fnnls, but use solving the lstsq problem by updating Cholesky factorisation.
     """
+    from scipy import linalg as slg
 
     lstsq = lambda A, x: slg.solve(
         A,
@@ -147,6 +147,8 @@ def fix_constraint_cholesky(ZTx, s_chol, d, P, P_inorder, U, tolerance):
     # solve the lstsq problem by cho_solve
 
     if len(P_inorder):
+        from scipy import linalg as slg
+
         # there could be a case where P_inorder is empty.
         s_chol[P_inorder] = slg.cho_solve((U, False), ZTx[P_inorder])