diff --git a/src/sage/graphs/graph.py b/src/sage/graphs/graph.py
index c28b4a4c451..92d0f544581 100644
--- a/src/sage/graphs/graph.py
+++ b/src/sage/graphs/graph.py
@@ -424,6 +424,10 @@
 from sage.misc.rest_index_of_methods import doc_index, gen_thematic_rest_table_index
 from sage.graphs.views import EdgesView
 
+from sage.misc.lazy_import import lazy_import
+from sage.features import PythonModule
+lazy_import('sage.graphs.mcqd', ['mcqd'],
+            feature=PythonModule('sage.graphs.mcqd', spkg='mcqd'))
 
 class Graph(GenericGraph):
     r"""
@@ -2691,10 +2695,11 @@ def treewidth(self, k=None, certificate=False, algorithm=None):
         elif algorithm != "sage":
             raise ValueError("'algorithm' must be equal to 'tdlib', 'sage', or None")
 
-        if algorithm is None and tdlib_found:
-            algorithm = 'tdlib'
-        else:
-            algorithm = 'sage'
+        if algorithm is None:
+            if tdlib_found:
+                algorithm = 'tdlib'
+            else:
+                algorithm = 'sage'
 
         if k is not None and k < 0:
             raise ValueError("k(={}) must be a nonnegative integer".format(k))
@@ -2717,8 +2722,9 @@ def treewidth(self, k=None, certificate=False, algorithm=None):
         # TDLIB
         if algorithm == 'tdlib':
             if not tdlib_found:
-                from sage.misc.package import PackageNotFoundError
-                raise PackageNotFoundError("tdlib")
+                from sage.features import FeatureNotPresentError
+                raise FeatureNotPresentError(PythonModule('sage.graphs.graph_decompositions.tdlib',
+                                                          spkg='tdlib'))
 
             T = tdlib.treedecomposition_exact(g, -1 if k is None else k)
             width = tdlib.get_width(T)
@@ -6410,11 +6416,6 @@ def clique_maximum(self,  algorithm="Cliquer", solver=None, verbose=0):
         elif algorithm == "MILP":
             return self.complement().independent_set(algorithm=algorithm, solver=solver, verbosity=verbose)
         elif algorithm == "mcqd":
-            try:
-                from sage.graphs.mcqd import mcqd
-            except ImportError:
-                from sage.misc.package import PackageNotFoundError
-                raise PackageNotFoundError("mcqd")
             return mcqd(self)
         else:
             raise NotImplementedError("Only 'MILP', 'Cliquer' and 'mcqd' are supported.")
@@ -6513,11 +6514,6 @@ def clique_number(self, algorithm="Cliquer", cliques=None, solver=None, verbose=
         elif algorithm == "MILP":
             return len(self.complement().independent_set(algorithm=algorithm, solver=solver, verbosity=verbose))
         elif algorithm == "mcqd":
-            try:
-                from sage.graphs.mcqd import mcqd
-            except ImportError:
-                from sage.misc.package import PackageNotFoundError
-                raise PackageNotFoundError("mcqd")
             return len(mcqd(self))
         else:
             raise NotImplementedError("Only 'networkx' 'MILP' 'Cliquer' and 'mcqd' are supported.")
diff --git a/src/sage/groups/braid.py b/src/sage/groups/braid.py
index 5e020504800..2f180ee26d0 100644
--- a/src/sage/groups/braid.py
+++ b/src/sage/groups/braid.py
@@ -68,6 +68,7 @@
 from sage.rings.integer import Integer
 from sage.rings.integer_ring import IntegerRing
 from sage.misc.lazy_attribute import lazy_attribute
+from sage.misc.lazy_import import lazy_import
 from sage.misc.cachefunc import cached_method
 from sage.categories.groups import Groups
 from sage.groups.free_group import FreeGroup, is_FreeGroup
@@ -80,7 +81,13 @@
 from sage.groups.artin import FiniteTypeArtinGroup, FiniteTypeArtinGroupElement
 from sage.misc.package import PackageNotFoundError
 from sage.structure.richcmp import richcmp, rich_to_bool
+from sage.features import PythonModule
 
+lazy_import('sage.libs.braiding',
+            ['rightnormalform', 'centralizer', 'supersummitset', 'greatestcommondivisor',
+             'leastcommonmultiple', 'conjugatingbraid', 'ultrasummitset',
+             'thurston_type', 'rigidity', 'sliding_circuits'],
+            feature=PythonModule('sage.libs.braiding', spkg='libbraiding'))
 
 class Braid(FiniteTypeArtinGroupElement):
     """
@@ -1070,10 +1077,6 @@ def right_normal_form(self):
             sage: b.right_normal_form()
             (s1*s0, s0*s2, 1)
         """
-        try:
-            from sage.libs.braiding import rightnormalform
-        except ImportError:
-            raise PackageNotFoundError("libbraiding")
         l = rightnormalform(self)
         B = self.parent()
         return tuple([B(b) for b in l[:-1]] + [B.delta() ** l[-1][0]])
@@ -1090,10 +1093,6 @@ def centralizer(self):
             [s1*s0*s2*s1, s0*s2]
 
         """
-        try:
-            from sage.libs.braiding import centralizer
-        except ImportError:
-            raise PackageNotFoundError("libbraiding")
         l = centralizer(self)
         B = self.parent()
         return [B._element_from_libbraiding(b) for b in l]
@@ -1113,10 +1112,6 @@ def super_summit_set(self):
             s0^-1*s1^-1*s0^-2*s1^-1*s0*s1^3*s0]
 
         """
-        try:
-            from sage.libs.braiding import supersummitset
-        except ImportError:
-            raise PackageNotFoundError("libbraiding")
         l = supersummitset(self)
         B = self.parent()
         return [B._element_from_libbraiding(b) for b in l]
@@ -1139,10 +1134,6 @@ def gcd(self, other):
             sage: c.gcd(b)
             s0^-1*s1^-1*s0^-2*s1^2*s0
         """
-        try:
-            from sage.libs.braiding import greatestcommondivisor
-        except ImportError:
-            raise PackageNotFoundError("libbraiding")
         B = self.parent()
         b = greatestcommondivisor(self, other)
         return B._element_from_libbraiding(b)
@@ -1163,10 +1154,6 @@ def lcm(self, other):
             sage: b.lcm(c)
             (s0*s1)^2*s0
         """
-        try:
-            from sage.libs.braiding import leastcommonmultiple
-        except ImportError:
-            raise PackageNotFoundError("libbraiding")
         B = self.parent()
         b = leastcommonmultiple(self, other)
         return B._element_from_libbraiding(b)
@@ -1193,10 +1180,6 @@ def conjugating_braid(self, other):
             sage: d * a / d == c
             False
         """
-        try:
-            from sage.libs.braiding import conjugatingbraid
-        except ImportError:
-            raise PackageNotFoundError("libbraiding")
         l = conjugatingbraid(self, other)
         if not l:
             return None
@@ -1222,10 +1205,6 @@ def is_conjugated(self, other):
             sage: c.is_conjugated(b)
             False
         """
-        try:
-            from sage.libs.braiding import conjugatingbraid
-        except ImportError:
-            raise PackageNotFoundError("libbraiding")
         l = conjugatingbraid(self, other)
         return bool(l)
 
@@ -1254,10 +1233,6 @@ def ultra_summit_set(self):
             s0^-1*s1^-1*s0^-2*s1^-1*s0^4*s1^2*s0,
             s0^-1*s1^-1*s0^-2*s1^-1*s0^3*s1^2*s0^2]]
         """
-        try:
-            from sage.libs.braiding import ultrasummitset
-        except ImportError:
-            raise PackageNotFoundError("libbraiding")
         uss = ultrasummitset(self)
         B = self.parent()
         return [[B._element_from_libbraiding(i) for i in s] for s in uss]
@@ -1283,10 +1258,6 @@ def thurston_type(self):
             sage: c.thurston_type()
             'periodic'
         """
-        try:
-            from sage.libs.braiding import thurston_type
-        except ImportError:
-            raise PackageNotFoundError("libbraiding")
         return thurston_type(self)
 
     def is_reducible(self):
@@ -1352,10 +1323,6 @@ def rigidity(self):
             sage: b.rigidity()
             0
         """
-        try:
-            from sage.libs.braiding import rigidity
-        except ImportError:
-            raise PackageNotFoundError("libbraiding")
         return Integer(rigidity(self))
 
     def sliding_circuits(self):
@@ -1388,10 +1355,6 @@ def sliding_circuits(self):
             sage: b.sliding_circuits()
             [[s0*s1*s0^2, (s0*s1)^2]]
         """
-        try:
-            from sage.libs.braiding import sliding_circuits
-        except ImportError:
-            raise PackageNotFoundError("libbraiding")
         slc = sliding_circuits(self)
         B = self.parent()
         return [[B._element_from_libbraiding(i) for i in s] for s in slc]
diff --git a/src/sage/matrix/matrix_space.py b/src/sage/matrix/matrix_space.py
index a8fc5e64493..78297f35a76 100644
--- a/src/sage/matrix/matrix_space.py
+++ b/src/sage/matrix/matrix_space.py
@@ -73,6 +73,11 @@
 from sage.categories.fields import Fields
 from sage.categories.enumerated_sets import EnumeratedSets
 
+from sage.misc.lazy_import import lazy_import
+from sage.features import PythonModule
+lazy_import('sage.matrix.matrix_gfpn_dense', ['Matrix_gfpn_dense'],
+            feature=PythonModule('sage.matrix.matrix_gfpn_dense', spkg='meataxe'))
+
 _Rings = Rings()
 _Fields = Fields()
 
@@ -279,13 +284,7 @@ def get_matrix_class(R, nrows, ncols, sparse, implementation):
 
         if implementation == 'meataxe':
             if R.is_field() and R.order() < 256:
-                try:
-                    from . import matrix_gfpn_dense
-                except ImportError:
-                    from sage.misc.package import PackageNotFoundError
-                    raise PackageNotFoundError('meataxe')
-                else:
-                    return matrix_gfpn_dense.Matrix_gfpn_dense
+                return Matrix_gfpn_dense
             raise ValueError("'meataxe' matrix can only deal with finite fields of order < 256")
 
         if implementation == 'numpy':
diff --git a/src/sage/misc/lazy_import.pyx b/src/sage/misc/lazy_import.pyx
index f759d953d98..9de1859fd3f 100644
--- a/src/sage/misc/lazy_import.pyx
+++ b/src/sage/misc/lazy_import.pyx
@@ -68,7 +68,7 @@ import inspect
 from . import sageinspect
 
 from .lazy_import_cache import get_cache_file
-
+from sage.features import FeatureNotPresentError
 
 cdef inline obj(x):
     if type(x) is LazyImport:
@@ -159,8 +159,10 @@ cdef class LazyImport(object):
     cdef _namespace
     cdef bint _at_startup
     cdef _deprecation
+    cdef _feature
 
-    def __init__(self, module, name, as_name=None, at_startup=False, namespace=None, deprecation=None):
+    def __init__(self, module, name, as_name=None, at_startup=False, namespace=None,
+                 deprecation=None, feature=None):
         """
         EXAMPLES::
 
@@ -178,6 +180,7 @@ cdef class LazyImport(object):
         self._namespace = namespace
         self._at_startup = at_startup
         self._deprecation = deprecation
+        self._feature = feature
 
     cdef inline get_object(self):
         """
@@ -217,7 +220,12 @@ cdef class LazyImport(object):
             raise RuntimeError(f"resolving lazy import {self._name} during startup")
         elif self._at_startup and not startup_guard:
             print('Option ``at_startup=True`` for lazy import {0} not needed anymore'.format(self._name))
-        self._object = getattr(__import__(self._module, {}, {}, [self._name]), self._name)
+        try:
+            self._object = getattr(__import__(self._module, {}, {}, [self._name]), self._name)
+        except ImportError as e:
+            if self._feature:
+                raise FeatureNotPresentError(self._feature, reason=f'Importing {self._name} failed: {e}')
+            raise
         name = self._as_name
         if self._deprecation is not None:
             from sage.misc.superseded import deprecation
@@ -953,7 +961,7 @@ cdef class LazyImport(object):
 
 
 def lazy_import(module, names, as_=None, *,
-    at_startup=False, namespace=None, deprecation=None):
+    at_startup=False, namespace=None, deprecation=None, feature=None):
     """
     Create a lazy import object and inject it into the caller's global
     namespace. For the purposes of introspection and calling, this is
@@ -1061,7 +1069,7 @@ def lazy_import(module, names, as_=None, *,
         names[ix:ix+1] = all
         as_[ix:ix+1] = all
     for name, alias in zip(names, as_):
-        namespace[alias] = LazyImport(module, name, alias, at_startup, namespace, deprecation)
+        namespace[alias] = LazyImport(module, name, alias, at_startup, namespace, deprecation, feature)
 
 
 star_imports = None
diff --git a/src/sage/sat/solvers/cryptominisat.py b/src/sage/sat/solvers/cryptominisat.py
index c27ed646162..595eb03224d 100644
--- a/src/sage/sat/solvers/cryptominisat.py
+++ b/src/sage/sat/solvers/cryptominisat.py
@@ -27,6 +27,10 @@
 
 from .satsolver import SatSolver
 
+from sage.misc.lazy_import import lazy_import
+from sage.features import PythonModule
+lazy_import('pycryptosat', ['Solver'],
+            feature=PythonModule('pycryptosat', spkg='cryptominisat'))
 
 class CryptoMiniSat(SatSolver):
     r"""
@@ -65,11 +69,6 @@ def __init__(self, verbosity=0, confl_limit=None, threads=None):
         if confl_limit is None:
             from sys import maxsize
             confl_limit = maxsize
-        try:
-            from pycryptosat import Solver
-        except ImportError:
-            from sage.misc.package import PackageNotFoundError
-            raise PackageNotFoundError("cryptominisat")
         self._solver = Solver(verbose=int(verbosity), confl_limit=int(confl_limit), threads=int(threads))
         self._nvars = 0
         self._clauses = []
diff --git a/src/sage/sat/solvers/picosat.py b/src/sage/sat/solvers/picosat.py
index dea731c0e84..91063de1fb7 100644
--- a/src/sage/sat/solvers/picosat.py
+++ b/src/sage/sat/solvers/picosat.py
@@ -24,6 +24,11 @@
 
 from .satsolver import SatSolver
 
+from sage.misc.lazy_import import lazy_import
+from sage.features import PythonModule
+lazy_import('pycosat', ['solve'],
+            feature=PythonModule('pycosat', spkg='pycosat'))
+
 class PicoSAT(SatSolver):
     r"""
     PicoSAT Solver.
@@ -55,12 +60,7 @@ def __init__(self, verbosity=0, prop_limit=0):
             self._prop_limit = 0
         else:
             self._prop_limit = int(prop_limit)
-        try:
-            import pycosat
-        except ImportError:
-            from sage.misc.package import PackageNotFoundError
-            raise PackageNotFoundError("pycosat")
-        self._solve = pycosat.solve
+        self._solve = solve
         self._nvars = 0
         self._clauses = []