diff --git a/geometric_kernels/_logging.py b/geometric_kernels/_logging.py
index d7b73c89..624a4ebb 100644
--- a/geometric_kernels/_logging.py
+++ b/geometric_kernels/_logging.py
@@ -1,4 +1,4 @@
-""" Setup logging """
+"""Setup logging"""
 
 import logging
 
diff --git a/geometric_kernels/feature_maps/deterministic.py b/geometric_kernels/feature_maps/deterministic.py
index 39875d79..bbad9c1a 100644
--- a/geometric_kernels/feature_maps/deterministic.py
+++ b/geometric_kernels/feature_maps/deterministic.py
@@ -43,15 +43,28 @@ def __init__(
         self.num_levels = num_levels
 
         if repeated_eigenvalues_laplacian is None:
-            assert eigenfunctions is None
+            if eigenfunctions is not None:
+                raise ValueError(
+                    "You must either provide both `repeated_eigenvalues_laplacian` and `eigenfunctions` or none of the two."
+                )
             repeated_eigenvalues_laplacian = self.space.get_repeated_eigenvalues(
                 self.num_levels
             )
             eigenfunctions = self.space.get_eigenfunctions(self.num_levels)
         else:
-            assert eigenfunctions is not None
-            assert repeated_eigenvalues_laplacian.shape == (num_levels, 1)
-            assert eigenfunctions.num_levels == num_levels
+            if eigenfunctions is None:
+                raise ValueError(
+                    "You must either provide both `repeated_eigenvalues_laplacian` and `eigenfunctions` or none of the two."
+                )
+            if repeated_eigenvalues_laplacian.shape != (num_levels, 1):
+                raise ValueError(
+                    f"Expected `repeated_eigenvalues_laplacian` to have shape [num_levels={num_levels}, 1] but got {B.shape(repeated_eigenvalues_laplacian)}."
+                )
+            if eigenfunctions.num_levels != num_levels:
+                raise ValueError(
+                    f"`num_levels` must coincide with `num_levels` in the provided `eigenfunctions`,"
+                    f"but `num_levels`={num_levels} and `eigenfunctions.num_levels`={eigenfunctions.num_levels}"
+                )
 
         self._repeated_eigenvalues = repeated_eigenvalues_laplacian
         self._eigenfunctions = eigenfunctions
diff --git a/geometric_kernels/feature_maps/probability_densities.py b/geometric_kernels/feature_maps/probability_densities.py
index 1e472ccc..97c210a7 100644
--- a/geometric_kernels/feature_maps/probability_densities.py
+++ b/geometric_kernels/feature_maps/probability_densities.py
@@ -22,7 +22,13 @@
     eigvalsh,
     from_numpy,
 )
-from geometric_kernels.utils.utils import ordered_pairwise_differences
+from geometric_kernels.utils.utils import (
+    _check_1_vector,
+    _check_field_in_params,
+    _check_matrix,
+    _check_rank_1_array,
+    ordered_pairwise_differences,
+)
 
 
 def student_t_sample(
@@ -74,13 +80,17 @@ def student_t_sample(
         samples of type `dtype`, and `key` is the updated random key for `jax`,
         or the similar random state (generator) for any other backend.
     """
-    assert B.shape(df) == (1,), "df must be a 1-vector."
+    _check_1_vector(df, "df")
 
-    n = int(B.length(loc))
+    _check_rank_1_array(loc, "loc")
+    _check_matrix(shape, "shape")
 
-    assert B.shape(loc) == (n,), "loc must be a 1-dim vector"
-    assert B.shape(shape) == (n, n), "shape must be a matrix"
+    n = B.shape(loc)[0]
 
+    if tuple(B.shape(shape)) != (n, n):
+        raise ValueError(
+            f"`Expected `shape` matrix to have shape [{n}, {n}], but got {B.shape(shape)}."
+        )
     shape_sqrt = B.chol(shape)
     dtype = dtype or dtype_double(key)
     key, z = B.randn(key, dtype, *size, n)
@@ -140,10 +150,11 @@ def base_density_sample(
         of samples, and `key` is the updated random key for `jax`, or the
         similar random state (generator) for any other backend.
     """
-    assert "lengthscale" in params
-    assert params["lengthscale"].shape == (1,)
-    assert "nu" in params
-    assert params["nu"].shape == (1,)
+    _check_field_in_params(params, "lengthscale")
+    _check_1_vector(params["lengthscale"], 'params["lengthscale"]')
+
+    _check_field_in_params(params, "nu")
+    _check_1_vector(params["nu"], 'params["nu"]')
 
     nu = params["nu"]
     L = params["lengthscale"]
@@ -232,7 +243,8 @@ def _alphas(n: int) -> B.Numeric:
     .. todo::
         Update proposition numbers when the paper gets published.
     """
-    assert n >= 2
+    if n < 2:
+        raise ValueError("Dimension of the hyperbolic space `n` must be >= 2.")
     x, j = symbols("x, j")
     if (n % 2) == 0:
         m = n // 2
@@ -269,9 +281,10 @@ def _sample_mixture_heat(
     .. todo::
         Update proposition numbers when the paper gets published.
     """
-    assert B.rank(alpha) == 1
+    _check_rank_1_array(alpha, "alpha")
     m = B.shape(alpha)[0] - 1
-    assert m >= 0
+    if m < 0:
+        raise ValueError("The mixture must contain at least 1 component.")
     dtype = B.dtype(lengthscale)
     js = B.range(dtype, 0, m + 1)
 
@@ -332,9 +345,10 @@ def _sample_mixture_matern(
     .. todo::
         Update proposition numbers when the paper gets published.
     """
-    assert B.rank(alpha) == 1
+    _check_rank_1_array(alpha, "alpha")
     m = B.shape(alpha)[0] - 1
-    assert m >= 0
+    if m < 0:
+        raise ValueError("The mixture must contain at least 1 component.")
     dtype = B.dtype(lengthscale)
     js = B.range(dtype, 0, m + 1)
     if shifted_laplacian:
@@ -397,10 +411,11 @@ def hyperbolic_density_sample(
         samples, and `key` is the updated random key for `jax`, or the similar
         random state (generator) for any other backend.
     """
-    assert "lengthscale" in params
-    assert params["lengthscale"].shape == (1,)
-    assert "nu" in params
-    assert params["nu"].shape == (1,)
+    _check_field_in_params(params, "lengthscale")
+    _check_1_vector(params["lengthscale"], 'params["lengthscale"]')
+
+    _check_field_in_params(params, "nu")
+    _check_1_vector(params["nu"], 'params["nu"]')
 
     nu = params["nu"]
     L = params["lengthscale"]
@@ -477,10 +492,11 @@ def spd_density_sample(
         samples, and `key` is the updated random key for `jax`, or the similar
         random state (generator) for any other backend.
     """
-    assert "lengthscale" in params
-    assert params["lengthscale"].shape == (1,)
-    assert "nu" in params
-    assert params["nu"].shape == (1,)
+    _check_field_in_params(params, "nu")
+    _check_1_vector(params["nu"], 'params["nu"]')
+
+    _check_field_in_params(params, "lengthscale")
+    _check_1_vector(params["lengthscale"], 'params["lengthscale"]')
 
     nu = params["nu"]
     L = params["lengthscale"]
@@ -514,7 +530,6 @@ def spd_density_sample(
         diffp = B.pi * B.abs(diffp)
         logprod = B.sum(B.log(B.tanh(diffp)), axis=-1)
         prod = B.exp(logprod)
-        assert B.all(prod > 0)
 
         # accept with probability `prod`
         key, u = B.rand(key, B.dtype(L), 1)
diff --git a/geometric_kernels/kernels/feature_map.py b/geometric_kernels/kernels/feature_map.py
index b224be43..f8a506fa 100644
--- a/geometric_kernels/kernels/feature_map.py
+++ b/geometric_kernels/kernels/feature_map.py
@@ -11,7 +11,11 @@
 from geometric_kernels.feature_maps import FeatureMap
 from geometric_kernels.kernels.base import BaseGeometricKernel
 from geometric_kernels.spaces.base import Space
-from geometric_kernels.utils.utils import make_deterministic
+from geometric_kernels.utils.utils import (
+    _check_1_vector,
+    _check_field_in_params,
+    make_deterministic,
+)
 
 
 class MaternFeatureMapKernel(BaseGeometricKernel):
@@ -108,10 +112,11 @@ def K(
         X2: Optional[B.Numeric] = None,
         **kwargs,
     ):
-        assert "lengthscale" in params
-        assert params["lengthscale"].shape == (1,)
-        assert "nu" in params
-        assert params["nu"].shape == (1,)
+        _check_field_in_params(params, "lengthscale")
+        _check_1_vector(params["lengthscale"], 'params["lengthscale"]')
+
+        _check_field_in_params(params, "nu")
+        _check_1_vector(params["nu"], 'params["nu"]')
 
         _, features_X = self.feature_map(
             X, params, normalize=self.normalize, **kwargs
@@ -127,10 +132,11 @@ def K(
         return feature_product
 
     def K_diag(self, params: Dict[str, B.Numeric], X: B.Numeric, **kwargs):
-        assert "lengthscale" in params
-        assert params["lengthscale"].shape == (1,)
-        assert "nu" in params
-        assert params["nu"].shape == (1,)
+        _check_field_in_params(params, "lengthscale")
+        _check_1_vector(params["lengthscale"], 'params["lengthscale"]')
+
+        _check_field_in_params(params, "nu")
+        _check_1_vector(params["nu"], 'params["nu"]')
 
         _, features_X = self.feature_map(
             X, params, normalize=self.normalize, **kwargs
diff --git a/geometric_kernels/kernels/hodge_compositional.py b/geometric_kernels/kernels/hodge_compositional.py
index 1180c67d..2d60a2df 100644
--- a/geometric_kernels/kernels/hodge_compositional.py
+++ b/geometric_kernels/kernels/hodge_compositional.py
@@ -11,6 +11,7 @@
 from geometric_kernels.kernels.base import BaseGeometricKernel
 from geometric_kernels.kernels.karhunen_loeve import MaternKarhunenLoeveKernel
 from geometric_kernels.spaces import HodgeDiscreteSpectrumSpace
+from geometric_kernels.utils.utils import _check_1_vector, _check_field_in_params
 
 
 class MaternHodgeCompositionalKernel(BaseGeometricKernel):
@@ -129,13 +130,9 @@ def K(
         inputs, or batches of matrices of inputs, depending on the space.
         """
 
-        assert all(
-            key in params for key in ["harmonic", "gradient", "curl"]
-        ), "MaternHodgeCompositionalKernel's parameters must contain keys 'harmonic', 'gradient', 'curl'."
-        assert all(
-            B.shape(params[key]["logit"]) == (1,)
-            for key in ["harmonic", "gradient", "curl"]
-        ), "The 'logit' parameters of MaternHodgeCompositionalKernel must have shape (1,)."
+        for key in ("harmonic", "gradient", "curl"):
+            _check_field_in_params(params, key)
+            _check_1_vector(params[key]["logit"], f'params["{key}"]["logit"]')
 
         # Copy the parameters to avoid modifying the original dict.
         params = {key: params[key].copy() for key in ["harmonic", "gradient", "curl"]}
@@ -162,13 +159,9 @@ def K_diag(
         diagonal.
         """
 
-        assert all(
-            key in params for key in ["harmonic", "gradient", "curl"]
-        ), "MaternHodgeCompositionalKernel's parameters must contain keys 'harmonic', 'gradient', 'curl'."
-        assert all(
-            B.shape(params[key]["logit"]) == (1,)
-            for key in ["harmonic", "gradient", "curl"]
-        ), "The 'logit' parameters of MaternHodgeCompositionalKernel must have shape (1,)."
+        for key in ("harmonic", "gradient", "curl"):
+            _check_field_in_params(params, key)
+            _check_1_vector(params[key]["logit"], f'params["{key}"]["logit"]')
 
         # Copy the parameters to avoid modifying the original dict.
         params = {key: params[key].copy() for key in ["harmonic", "gradient", "curl"]}
diff --git a/geometric_kernels/kernels/karhunen_loeve.py b/geometric_kernels/kernels/karhunen_loeve.py
index 48f38e96..0aa3625d 100644
--- a/geometric_kernels/kernels/karhunen_loeve.py
+++ b/geometric_kernels/kernels/karhunen_loeve.py
@@ -11,6 +11,7 @@
 from geometric_kernels.lab_extras import from_numpy, is_complex
 from geometric_kernels.spaces import DiscreteSpectrumSpace
 from geometric_kernels.spaces.eigenfunctions import Eigenfunctions
+from geometric_kernels.utils.utils import _check_1_vector, _check_field_in_params
 
 
 class MaternKarhunenLoeveKernel(BaseGeometricKernel):
@@ -73,13 +74,26 @@ def __init__(
         self.num_levels = num_levels  # in code referred to as `L`.
 
         if eigenvalues_laplacian is None:
-            assert eigenfunctions is None
+            if eigenfunctions is not None:
+                raise ValueError(
+                    "You must either provide both `eigenfunctions` and `eigenvalues_laplacian`, or none of the two."
+                )
             eigenvalues_laplacian = self.space.get_eigenvalues(self.num_levels)
             eigenfunctions = self.space.get_eigenfunctions(self.num_levels)
         else:
-            assert eigenfunctions is not None
-            assert eigenvalues_laplacian.shape == (num_levels, 1)
-            assert eigenfunctions.num_levels == num_levels
+            if eigenfunctions is None:
+                raise ValueError(
+                    "You must either provide both `eigenfunctions` and `eigenvalues_laplacian`, or none of the two."
+                )
+            if eigenvalues_laplacian.shape != (num_levels, 1):
+                raise ValueError(
+                    f"Expected `eigenvalues_laplacian` to have shape [num_levels={num_levels}, 1] but got {eigenvalues_laplacian.shape}"
+                )
+            if eigenfunctions.num_levels != num_levels:
+                raise ValueError(
+                    f"`num_levels` must coincide with `num_levels` in the provided `eigenfunctions`,"
+                    f"but `num_levels`={num_levels} and `eigenfunctions.num_levels`={eigenfunctions.num_levels}"
+                )
 
         self._eigenvalues_laplacian = eigenvalues_laplacian
         self._eigenfunctions = eigenfunctions
@@ -134,8 +148,8 @@ def spectrum(
         :return:
             The spectrum of the Matérn kernel.
         """
-        assert lengthscale.shape == (1,)
-        assert nu.shape == (1,)
+        _check_1_vector(lengthscale, "lengthscale")
+        _check_1_vector(nu, "nu")
 
         # Note: 1.0 in safe_nu can be replaced by any finite positive value
         safe_nu = B.where(nu == np.inf, B.cast(B.dtype(lengthscale), np.r_[1.0]), nu)
@@ -180,10 +194,11 @@ def eigenvalues(self, params: Dict[str, B.Numeric]) -> B.Numeric:
         :return:
             An [L, 1]-shaped array.
         """
-        assert "lengthscale" in params
-        assert params["lengthscale"].shape == (1,)
-        assert "nu" in params
-        assert params["nu"].shape == (1,)
+        _check_field_in_params(params, "lengthscale")
+        _check_1_vector(params["lengthscale"], 'params["lengthscale"]')
+
+        _check_field_in_params(params, "nu")
+        _check_1_vector(params["nu"], 'params["nu"]')
 
         spectral_values = self.spectrum(
             self.eigenvalues_laplacian,
@@ -210,10 +225,11 @@ def eigenvalues(self, params: Dict[str, B.Numeric]) -> B.Numeric:
     def K(
         self, params: Dict[str, B.Numeric], X: B.Numeric, X2: Optional[B.Numeric] = None, **kwargs  # type: ignore
     ) -> B.Numeric:
-        assert "lengthscale" in params
-        assert params["lengthscale"].shape == (1,)
-        assert "nu" in params
-        assert params["nu"].shape == (1,)
+        _check_field_in_params(params, "lengthscale")
+        _check_1_vector(params["lengthscale"], 'params["lengthscale"]')
+
+        _check_field_in_params(params, "nu")
+        _check_1_vector(params["nu"], 'params["nu"]')
 
         weights = B.cast(B.dtype(params["nu"]), self.eigenvalues(params))  # [L, 1]
         Phi = self.eigenfunctions
@@ -224,10 +240,11 @@ def K(
             return K
 
     def K_diag(self, params: Dict[str, B.Numeric], X: B.Numeric, **kwargs) -> B.Numeric:
-        assert "lengthscale" in params
-        assert params["lengthscale"].shape == (1,)
-        assert "nu" in params
-        assert params["nu"].shape == (1,)
+        _check_field_in_params(params, "lengthscale")
+        _check_1_vector(params["lengthscale"], 'params["lengthscale"]')
+
+        _check_field_in_params(params, "nu")
+        _check_1_vector(params["nu"], 'params["nu"]')
 
         weights = B.cast(B.dtype(params["nu"]), self.eigenvalues(params))  # [L, 1]
         Phi = self.eigenfunctions
diff --git a/geometric_kernels/kernels/product.py b/geometric_kernels/kernels/product.py
index 69dcdde6..85ce7d30 100644
--- a/geometric_kernels/kernels/product.py
+++ b/geometric_kernels/kernels/product.py
@@ -63,7 +63,10 @@ def __init__(
         self.spaces: List[Space] = []
         for kernel in self.kernels:
             # Make sure there is no product kernel in the list of kernels.
-            assert isinstance(kernel.space, Space)
+            if not isinstance(kernel.space, Space):  # as opposed to List[Space]
+                raise NotImplementedError(
+                    "One of the provided kernels is a product kernel itself."
+                )
             self.spaces.append(kernel.space)
         self.element_shapes = [space.element_shape for space in self.spaces]
         self.element_dtypes = [space.element_dtype for space in self.spaces]
@@ -77,9 +80,16 @@ def __init__(
                 self.dimension_indices.append(inds[i : i + dim])
                 i += dim
         else:
-            assert len(dimension_indices) == len(self.kernels)
+            if len(dimension_indices) != len(self.kernels):
+                raise ValueError(
+                    f"`dimension_indices` must correspond to `kernels`, but got {len(kernels)} kernels and {len(dimension_indices)} dimension indices."
+                )
             for idx_list in dimension_indices:
-                assert all(idx >= 0 for idx in idx_list)
+                for idx in idx_list:
+                    if idx < 0:
+                        raise ValueError(
+                            "Expected all `dimension_indices` to be non-negative."
+                        )
 
             self.dimension_indices = dimension_indices
 
@@ -99,10 +109,19 @@ def init_params(self) -> Dict[str, B.NPNumeric]:
         nu_list: List[B.NPNumeric] = []
         lengthscale_list: List[B.NPNumeric] = []
 
-        for kernel in self.kernels:
+        for kernel_idx, kernel in enumerate(self.kernels):
             cur_params = kernel.init_params()
-            assert cur_params["lengthscale"].shape == (1,)
-            assert cur_params["nu"].shape == (1,)
+            if B.shape(cur_params["lengthscale"]) != (1,):
+                raise ValueError(
+                    f"All kernels' `lengthscale`s must be have shape [1,], but {kernel_idx}th kernel "
+                    f"({kernel}) violates this with shape {B.shape(cur_params['lengthscale'])}."
+                )
+            if B.shape(cur_params["nu"]) != (1,):
+                raise ValueError(
+                    f"All kernels' `nu`s must be have [1,], but {kernel_idx}th kernel "
+                    f"({kernel}) violates this with shape {B.shape(cur_params['nu'])}."
+                )
+
             nu_list.append(cur_params["nu"])
             lengthscale_list.append(cur_params["lengthscale"])
 
diff --git a/geometric_kernels/lab_extras/numpy/sparse_extras.py b/geometric_kernels/lab_extras/numpy/sparse_extras.py
index e4a938ce..a500aa5b 100644
--- a/geometric_kernels/lab_extras/numpy/sparse_extras.py
+++ b/geometric_kernels/lab_extras/numpy/sparse_extras.py
@@ -90,5 +90,6 @@ def pinv(a: Union[SparseArray]):
 B.shape.register(lambda a: a.shape, _SparseArray)
 B.sqrt.register(lambda a: a.sqrt(), _SparseArray)
 B.any.register(lambda a: bool((a == True).sum()), _SparseArray)  # noqa
+B.rank.register(lambda a: a.ndim, _SparseArray)
 
 B.linear_algebra.pinv.register(pinv, _SparseArray)
diff --git a/geometric_kernels/spaces/circle.py b/geometric_kernels/spaces/circle.py
index d80d4d51..40f1b3b0 100644
--- a/geometric_kernels/spaces/circle.py
+++ b/geometric_kernels/spaces/circle.py
@@ -28,7 +28,8 @@ class SinCosEigenfunctions(EigenfunctionsWithAdditionTheorem):
     """
 
     def __init__(self, num_levels: int):
-        assert num_levels >= 1
+        if num_levels < 1:
+            raise ValueError("`num_levels` must be a positive integer.")
 
         self._num_eigenfunctions = num_levels * 2 - 1
         self._num_levels = num_levels
diff --git a/geometric_kernels/spaces/graph.py b/geometric_kernels/spaces/graph.py
index bfdf37f5..fed8e88d 100644
--- a/geometric_kernels/spaces/graph.py
+++ b/geometric_kernels/spaces/graph.py
@@ -18,6 +18,7 @@
     Eigenfunctions,
     EigenfunctionsFromEigenvectors,
 )
+from geometric_kernels.utils.utils import _check_matrix
 
 
 class Graph(DiscreteSpectrumSpace):
@@ -66,15 +67,16 @@ def __str__(self):
         return f"Graph({self.num_vertices}, {'normalized' if self._normalized else 'unnormalized'})"
 
     @staticmethod
-    def _checks(adjacency):
+    def _checks(adjacency_matrix):
         """
-        Checks if `adjacency` is a square symmetric matrix.
+        Checks if `adjacency_matrix` is a square symmetric matrix.
         """
-        assert (
-            len(adjacency.shape) == 2 and adjacency.shape[0] == adjacency.shape[1]
-        ), "Matrix is not square."
+        _check_matrix(adjacency_matrix, "adjacency_matrix")
+        if B.shape(adjacency_matrix)[0] != B.shape(adjacency_matrix)[1]:
+            raise ValueError("`adjacency_matrix` must be a square matrix.")
 
-        assert not B.any(adjacency != B.T(adjacency)), "Adjacency is not symmetric"
+        if B.any(adjacency_matrix != B.T(adjacency_matrix)):
+            raise ValueError("`adjacency_matrix` must be a symmetric matrix.")
 
     @property
     def dimension(self) -> int:
@@ -118,9 +120,10 @@ def get_eigensystem(self, num):
         :return:
             A tuple of eigenvectors [n, num], eigenvalues [num, 1].
         """
-        assert (
-            num <= self.num_vertices
-        ), "Number of eigenpairs cannot exceed the number of vertices"
+        if num > self.num_vertices:
+            raise ValueError(
+                "Number of eigenpairs cannot exceed the number of vertices."
+            )
         if num not in self.cache:
             evals, evecs = eigenpairs(self._laplacian, num)
 
diff --git a/geometric_kernels/spaces/graph_edges.py b/geometric_kernels/spaces/graph_edges.py
index b712b47b..ae603115 100644
--- a/geometric_kernels/spaces/graph_edges.py
+++ b/geometric_kernels/spaces/graph_edges.py
@@ -21,6 +21,7 @@
     Eigenfunctions,
     EigenfunctionsFromEigenvectors,
 )
+from geometric_kernels.utils.utils import _check_matrix, _check_rank_1_array
 
 
 class GraphEdges(HodgeDiscreteSpectrumSpace):
@@ -217,7 +218,7 @@ def _compute_index(num_nodes: int, oriented_edges: B.Numeric) -> csr_matrix:
             result[oriented_edges[i, 1], oriented_edges[i, 0]] = -i - 1
         return result.tocsr()
 
-    def _checks_oriented_edges(
+    def _checks_oriented_edges(  # NOQA: C901
         self, oriented_edges: B.Numeric, num_nodes: int, comprehensive: bool = False
     ):
         """
@@ -229,43 +230,42 @@ def _checks_oriented_edges(
         :param comprehensive:
             If True, perform more extensive checks.
         """
+        _check_matrix(oriented_edges, "oriented_edges")
 
-        assert (
-            B.rank(oriented_edges) == 2
-        ), "The oriented_edges array must be 2-dimensional."
-
-        assert B.shape(oriented_edges)[1] == 2, "oriented_edges must have shape (*, 2)."
-
-        assert B.dtype(oriented_edges) == int_like(
-            oriented_edges
-        ), "The oriented_edges must be an array of integers."
-        assert B.all(
-            oriented_edges >= 0
-        ), "The oriented_edges array must contain only non-negative values."
-        assert B.all(
-            oriented_edges < self.num_nodes
-        ), "The values in the oriented_edges array must be < self.num_nodes."
-        assert B.all(
-            oriented_edges[:, 0] - oriented_edges[:, 1] != 0
-        ), "Loops are not allowed."
-
-        if comprehensive:
-            num_edges = oriented_edges.shape[0]
-
-            for i in range(num_edges):
-                for j in range(i + 1, num_edges):
-                    assert B.any(
-                        oriented_edges[i, :] != oriented_edges[j, :]
-                    ), "The oriented_edges array must not contain duplicate edges."
-                    assert B.any(
-                        oriented_edges[i, :] != oriented_edges[j, ::-1]
-                    ), "The oriented_edges array must not contain duplicate edges."
-
-            assert set(range(self.num_nodes)) == set(
-                B.to_numpy(B.flatten(oriented_edges))
-            ), "The oriented_edges array must contain all nodes."
-
-    def _checks_oriented_triangles(
+        if B.shape(oriented_edges)[1] != 2:
+            raise ValueError("`oriented_edges` must have shape (*, 2).")
+
+        if B.dtype(oriented_edges) != int_like(oriented_edges):
+            raise ValueError("`oriented_edges` must be an array of integers.")
+        if B.any(oriented_edges < 0):
+            raise ValueError("`oriented_edges` must contain only non-negative values.")
+        if B.any(oriented_edges >= self.num_nodes):
+            raise ValueError(
+                "The values in the `oriented_edges` array must be less than `self.num_nodes.`"
+            )
+        if B.any(oriented_edges[:, 0] - oriented_edges[:, 1] == 0):
+            raise ValueError("Loops are not allowed.")
+
+        if not comprehensive:
+            return
+
+        num_edges = oriented_edges.shape[0]
+
+        for i in range(num_edges):
+            for j in range(i + 1, num_edges):
+                if B.all(oriented_edges[i, :] == oriented_edges[j, :]):
+                    raise ValueError(
+                        "`oriented_edges` must not contain duplicate edges."
+                    )
+                if B.all(oriented_edges[i, :] == oriented_edges[j, ::-1]):
+                    raise ValueError(
+                        "`oriented_edges` must not contain duplicate edges."
+                    )
+
+        if set(range(self.num_nodes)) != set(B.to_numpy(B.flatten(oriented_edges))):
+            raise ValueError("`oriented_edges` must contain all nodes.")
+
+    def _checks_oriented_triangles(  # NOQA: C901
         self, oriented_triangles: B.Numeric, comprehensive=False
     ):
         """
@@ -277,48 +277,43 @@ def _checks_oriented_triangles(
         :param comprehensive:
             If True, perform more extensive checks.
         """
+        _check_matrix(oriented_triangles, "oriented_triangles")
+        if B.shape(oriented_triangles)[1] != 3:
+            raise ValueError("`oriented_triangles` must have shape (*, 3).")
+
+        if B.dtype(oriented_triangles) != int_like(oriented_triangles):
+            raise ValueError("`oriented_triangles` must be an array of integers.")
+        if B.any(B.abs(oriented_triangles) < 1):
+            raise ValueError("`oriented_triangles` must contain only non-zero values.")
+        if B.any(B.abs(oriented_triangles) > self.num_edges):
+            raise ValueError(
+                "Absolute values in `oriented_triangles` array must be less than or equal to `self.num_edges`."
+            )
 
-        assert (
-            B.rank(oriented_triangles) == 2
-        ), "The oriented_triangles array must be 2-dimensional."
-
-        assert (
-            B.shape(oriented_triangles)[1] == 3
-        ), "oriented_triangles must have shape (*, 3)."
-
-        assert B.dtype(oriented_triangles) == int_like(
-            oriented_triangles
-        ), "The oriented_triangles must be an array of integers."
-        assert B.all(
-            B.abs(oriented_triangles) >= 1
-        ), "The oriented_triangles array must contain only non-zero values."
-        assert B.all(
-            B.abs(oriented_triangles) <= self.num_edges
-        ), "The absolute values in the oriented_triangles array must be <= self.num_edges."
-
-        assert B.all(
-            B.abs(oriented_triangles) < self.num_edges
-        ), "The absolute values in the oriented_triangles array must be less than self.num_edges."
-        assert (
-            B.all(
-                B.abs(oriented_triangles[:, 0]) - B.abs(oriented_triangles[:, 1]) != 0
+        if (
+            B.any(
+                B.abs(oriented_triangles[:, 0]) - B.abs(oriented_triangles[:, 1]) == 0
             )
-            or B.all(
-                B.abs(oriented_triangles[:, 0]) - B.abs(oriented_triangles[:, 2]) != 0
+            and B.any(
+                B.abs(oriented_triangles[:, 0]) - B.abs(oriented_triangles[:, 2]) == 0
             )
-            or B.all(
-                B.abs(oriented_triangles[:, 1]) - B.abs(oriented_triangles[:, 2]) != 0
+            and B.any(
+                B.abs(oriented_triangles[:, 1]) - B.abs(oriented_triangles[:, 2]) == 0
             )
-        ), "Triangles must consist of 3 different edges."
+        ):
+            raise ValueError("Triangles must consist of 3 different edges.")
 
-        if comprehensive:
-            num_triangles = oriented_triangles.shape[0]
+        if not comprehensive:
+            return
 
-            for i in range(num_triangles):
-                for j in range(i + 1, num_triangles):
-                    assert B.any(
-                        oriented_triangles[i, :] != oriented_triangles[j, :]
-                    ), "The oriented_triangles array must not contain duplicate triangles."
+        num_triangles = oriented_triangles.shape[0]
+
+        for i in range(num_triangles):
+            for j in range(i + 1, num_triangles):
+                if B.all(oriented_triangles[i, :] == oriented_triangles[j, :]):
+                    raise ValueError(
+                        "The oriented_triangles array must not contain duplicate triangles."
+                    )
 
     def _checks_compatible(
         self,
@@ -331,44 +326,42 @@ def _checks_compatible(
             The oriented triangles array.
         """
 
-        assert B.dtype(self.oriented_edges) == B.dtype(
-            oriented_triangles
-        ), "The oriented_edges and oriented_triangles arrays must have the same dtype."
+        if B.dtype(self.oriented_edges) != B.dtype(oriented_triangles):
+            raise ValueError(
+                "`oriented_edges` and `oriented_triangles` must have the same dtype."
+            )
 
         num_triangles = oriented_triangles.shape[0]
         for t in range(num_triangles):
             resolved_edges = self.resolve_edges(oriented_triangles[t, :])
-            assert (
-                resolved_edges[0, 1] == resolved_edges[1, 0]
-            ), "The edges in the triangle must be connected."
-            assert (
-                resolved_edges[1, 1] == resolved_edges[2, 0]
-            ), "The edges in the triangle must be connected."
-            assert (
-                resolved_edges[2, 1] == resolved_edges[0, 0]
-            ), "The edges in the triangle must be connected."
+            if resolved_edges[0, 1] != resolved_edges[1, 0]:
+                raise ValueError("The edges in the triangle must be connected.")
+            if resolved_edges[1, 1] != resolved_edges[2, 0]:
+                raise ValueError("The edges in the triangle must be connected.")
+            if resolved_edges[2, 1] != resolved_edges[0, 0]:
+                raise ValueError("The edges in the triangle must be connected.")
 
     def _check_index(self, index: csr_matrix):
         edges = []
         for e in range(1, self.oriented_edges.shape[0] + 1):
             i, j = self.oriented_edges[e - 1, :]
-            assert (
-                index[i, j] == e
-            ), "The index matrix must be compatible with oriented_edges."
-            assert (
-                index[j, i] == -e
-            ), "The index matrix must be compatible with oriented_edges."
+            if index[i, j] != e:
+                raise ValueError("`index` must be compatible with `oriented_edges`.")
+            if index[j, i] != -e:
+                raise ValueError("`index` must be compatible with `oriented_edges`.")
             edges.append((min(i, j), max(i, j)))
 
         for i in range(self.num_nodes):
             for j in range(i + 1, self.num_nodes):
                 if (i, j) not in edges:
-                    assert (
-                        index[i, j] == 0
-                    ), "The index matrix must be compatible with oriented_edges."
-                    assert (
-                        index[j, i] == 0
-                    ), "The index matrix must be compatible with oriented_edges."
+                    if index[i, j] != 0:
+                        raise ValueError(
+                            "`index` must be compatible with `oriented_edges`."
+                        )
+                    if index[j, i] != 0:
+                        raise ValueError(
+                            "`index` must be compatible with `oriented_edges`."
+                        )
 
     def resolve_edges(self, es: B.Int) -> B.Int:
         r"""
@@ -382,8 +375,9 @@ def resolve_edges(self, es: B.Int) -> B.Int:
             A 2-dimensional array `result` such that `result[e, :]` is `[i, j]`
             where \|e\| = (i, j) if e > 0 and \|e\| = (j, i) if e < 0.
         """
-        assert B.rank(es) == 1
-        assert B.all(B.abs(es) >= 1) and B.all(B.abs(es) <= self.num_edges)
+        _check_rank_1_array(es, "es")
+        if not (B.all(B.abs(es) >= 1) and B.all(B.abs(es) <= self.num_edges)):
+            raise ValueError("`abs(es)` must lie in the interval [1, `num_edges`].")
 
         result = self.oriented_edges[B.abs(es) - 1]
         result = B.where(B.expand_dims(es > 0, axis=-1), result, result[:, ::-1])
@@ -402,8 +396,9 @@ def resolve_triangles(self, ts: B.Int) -> B.Int:
             where i = e1[0], j = e2[0], k = e3[0], and e1, e2, e3 are the
             oriented edges constituting the triangle `t`.
         """
-        assert B.rank(ts) == 1
-        assert B.all(B.abs(ts) >= 0) and B.all(B.abs(ts) < self.num_triangles)
+        _check_rank_1_array(ts, "ts")
+        if not (B.all(B.abs(ts) >= 0) and B.all(B.abs(ts) < self.num_triangles)):
+            raise ValueError("`abs(ts)` must lie in the interval [1, `num_edges`].")
 
         edge_indices = B.flatten(
             self.oriented_triangles[ts]
@@ -453,14 +448,15 @@ def from_adjacency(  # noqa: C901
                 f"The adjacency matrix must be a numpy array or a scipy sparse matrix not {type(adjacency_matrix)}. Use `type_reference` to specify the backend."
             )
 
-        if len(index.shape) != 2:
-            raise ValueError("Adjacency matrix must be a square matrix.")
+        _check_matrix(index, "adjacency_matrix")
+        if B.shape(index)[0] != B.shape(index)[1]:
+            raise ValueError("`adjacency_matrix` must be a square matrix.")
         if (abs(index - index.T) > 1e-10).nnz != 0:
-            raise ValueError("Adjacency matrix must be symmetric.")
+            raise ValueError("`adjacency_matrix` must be symmetric.")
         if (index.diagonal() != 0).any():
-            raise ValueError("Adjacency matrix must have zeros on the diagonal.")
+            raise ValueError("`adjacency_matrix` must have zeros on the diagonal.")
         if np.sum(index.data == 1) + np.sum(index.data == 0) != len(index.data):
-            raise ValueError("Adjacency matrix can only contain zeros and ones.")
+            raise ValueError("`adjacency_matrix` can only contain zeros and ones.")
 
         number_of_nodes = index.shape[0]
         number_of_edges = np.sum(index.data) // 2
@@ -477,9 +473,11 @@ def from_adjacency(  # noqa: C901
                     index[i, j] = cur_edge_ind
                     index[j, i] = -index[i, j]
                     cur_edge_ind += 1
-        assert (
-            cur_edge_ind == number_of_edges + 1
-        )  # double check that we have the right number of edges
+        if cur_edge_ind != number_of_edges + 1:
+            # double check that we have the right number of edges
+            raise RuntimeError(
+                "This should have never happened, please report a bug at https://github.com/geometric-kernels/GeometricKernels/issues."
+            )
         oriented_edges = B.cast(dtype_integer(type_reference), oriented_edges)
 
         if triangles is None:
diff --git a/geometric_kernels/spaces/hypercube_graph.py b/geometric_kernels/spaces/hypercube_graph.py
index 2034fead..f81f41ab 100644
--- a/geometric_kernels/spaces/hypercube_graph.py
+++ b/geometric_kernels/spaces/hypercube_graph.py
@@ -45,7 +45,8 @@ class WalshFunctions(EigenfunctionsWithAdditionTheorem):
     """
 
     def __init__(self, dim: int, num_levels: int) -> None:
-        assert num_levels <= dim + 1, "The number of levels should be at most dim+1."
+        if num_levels > dim + 1:
+            raise ValueError("The number of levels should be at most `dim`+1.")
         self.dim = dim
         self._num_levels = num_levels
         self._num_eigenfunctions: Optional[int] = None  # To be computed when needed.
diff --git a/geometric_kernels/spaces/mesh.py b/geometric_kernels/spaces/mesh.py
index 122806f7..b761b1b6 100644
--- a/geometric_kernels/spaces/mesh.py
+++ b/geometric_kernels/spaces/mesh.py
@@ -60,7 +60,10 @@ class Mesh(DiscreteSpectrumSpace):
 
     def __init__(self, vertices: np.ndarray, faces: np.ndarray):
         self._vertices = vertices
-        assert self._vertices.shape[1] == 3  # make sure we all is in R^3.
+        if B.shape(self._vertices)[1] != 3:
+            # make sure we are in R^3.
+            raise ValueError("The last dimension (axis) of `_vertices` must be 3.")
+
         self._faces = faces
         self._eigenvalues = None
         self._eigenfunctions = None
diff --git a/geometric_kernels/spaces/product.py b/geometric_kernels/spaces/product.py
index c2b7de1c..14810deb 100644
--- a/geometric_kernels/spaces/product.py
+++ b/geometric_kernels/spaces/product.py
@@ -215,7 +215,13 @@ def __init__(
             self.eigenindicies, self.nums_per_level
         )
 
-        assert self.eigenindicies.shape[-1] == len(self.eigenfunctions)
+        if self.eigenindicies.shape[-1] != len(self.eigenfunctions):
+            raise ValueError(
+                "Expected to have S `eigenfunctions` and `eigenindicies` of shape [L, S], "
+                "where S is the number of spaces and L is the number of levels, "
+                f"but got S1={len(self.eigenfunctions)} eigenfunctions and "
+                f"the shape of `eigenindicies` is {self.eigenindicies.shape}, which is incompatible."
+            )
 
     def __call__(self, X: B.Numeric, **kwargs) -> B.Numeric:
         """
@@ -416,9 +422,10 @@ def __init__(
         num_levels_per_space: Optional[int] = None,
     ):
         for space in spaces:
-            assert isinstance(
-                space, DiscreteSpectrumSpace
-            ), "One of the spaces is not an instance of DiscreteSpectrumSpace."
+            if not isinstance(space, DiscreteSpectrumSpace):
+                raise ValueError(
+                    "One of the spaces is not an instance of DiscreteSpectrumSpace."
+                )
 
         self.factor_spaces = spaces  # List of length S
         self.num_levels = num_levels
@@ -430,9 +437,10 @@ def __init__(
 
         if num_levels_per_space is None:
             num_levels_per_space = num_levels
-        assert num_levels <= num_levels_per_space ** len(
-            spaces
-        ), "Cannot have more levels than there are possible combinations"
+        if num_levels > num_levels_per_space ** len(spaces):
+            raise ValueError(
+                "Cannot have more levels than there are possible combinations."
+            )
 
         # prefetch the eigenvalues of the subspaces
         factor_space_eigenvalues = B.stack(
@@ -493,7 +501,10 @@ def get_eigenfunctions(self, num: int) -> Eigenfunctions:
             Number of levels. Cannot be larger than the `num_levels` parameter
             of the constructor.
         """
-        assert num <= self.num_levels
+        if num > self.num_levels:
+            raise ValueError(
+                "`num` cannot be larger than the `num_levels` provided in the constructor."
+            )
 
         max_level = int(self.factor_space_eigenindices[:num, :].max() + 1)
 
@@ -518,7 +529,10 @@ def get_eigenvalues(self, num: int) -> B.Numeric:
         :return:
             (num, 1)-shaped array containing the eigenvalues.
         """
-        assert num <= self.num_levels
+        if num > self.num_levels:
+            raise ValueError(
+                "`num` cannot be larger than the `num_levels` provided in the constructor."
+            )
 
         return self._eigenvalues[:num, None]
 
@@ -535,7 +549,10 @@ def get_repeated_eigenvalues(self, num: int) -> B.Numeric:
             (J, 1)-shaped array containing the repeated eigenvalues,`J is
             the resulting number of the repeated eigenvalues.
         """
-        assert num <= self.num_levels
+        if num > self.num_levels:
+            raise ValueError(
+                "`num` cannot be larger than the `num_levels` provided in the constructor."
+            )
 
         eigenfunctions = self.get_eigenfunctions(num)
         eigenvalues = self._eigenvalues[:num]
diff --git a/geometric_kernels/utils/kernel_formulas/euclidean.py b/geometric_kernels/utils/kernel_formulas/euclidean.py
index f76f347d..0a3d762d 100644
--- a/geometric_kernels/utils/kernel_formulas/euclidean.py
+++ b/geometric_kernels/utils/kernel_formulas/euclidean.py
@@ -27,7 +27,8 @@ def euclidean_matern_12_kernel(
         The kernel values evaluated at `r`, an array of shape [...].
     """
 
-    assert B.all(r >= 0.0)
+    if not B.all(r >= 0.0):
+        raise ValueError("Distances must be non-negative.")
 
     return B.exp(-r / lengthscale)
 
@@ -49,7 +50,8 @@ def euclidean_matern_32_kernel(
         The kernel values evaluated at `r`, an array of shape [...].
     """
 
-    assert B.all(r >= 0.0)
+    if not B.all(r >= 0.0):
+        raise ValueError("Distances must be non-negative.")
 
     sqrt3 = sqrt(3.0)
     r = r / lengthscale
@@ -73,7 +75,8 @@ def euclidean_matern_52_kernel(
         The kernel values evaluated at `r`, an array of shape [...].
     """
 
-    assert B.all(r >= 0.0)
+    if not B.all(r >= 0.0):
+        raise ValueError("Distances must be non-negative.")
 
     sqrt5 = sqrt(5.0)
     r = r / lengthscale
@@ -97,7 +100,8 @@ def euclidean_rbf_kernel(
         The kernel values evaluated at `r`, an array of shape [...].
     """
 
-    assert B.all(r >= 0.0)
+    if not B.all(r >= 0.0):
+        raise ValueError("Distances must be non-negative.")
 
     r = r / lengthscale
     return B.exp(-0.5 * r**2)
diff --git a/geometric_kernels/utils/kernel_formulas/hypercube_graph.py b/geometric_kernels/utils/kernel_formulas/hypercube_graph.py
index 254d7f23..783aca08 100644
--- a/geometric_kernels/utils/kernel_formulas/hypercube_graph.py
+++ b/geometric_kernels/utils/kernel_formulas/hypercube_graph.py
@@ -10,7 +10,11 @@
 from beartype.typing import Optional
 
 from geometric_kernels.lab_extras import float_like
-from geometric_kernels.utils.utils import hamming_distance
+from geometric_kernels.utils.utils import (
+    _check_1_vector,
+    _check_matrix,
+    hamming_distance,
+)
 
 
 def hypercube_graph_heat_kernel(
@@ -36,9 +40,12 @@ def hypercube_graph_heat_kernel(
     if X2 is None:
         X2 = X
 
-    assert lengthscale.shape == (1,)
-    assert X.ndim == 2 and X2.ndim == 2
-    assert X.shape[-1] == X2.shape[-1]
+    _check_1_vector(lengthscale, "lengthscale")
+    _check_matrix(X, "X")
+    _check_matrix(X2, "X2")
+
+    if X.shape[-1] != X2.shape[-1]:
+        raise ValueError("`X` and `X2` must live in a same-dimensional space.")
 
     if normalized_laplacian:
         d = X.shape[-1]
diff --git a/geometric_kernels/utils/kernel_formulas/spd.py b/geometric_kernels/utils/kernel_formulas/spd.py
index 842348c1..e628cc93 100644
--- a/geometric_kernels/utils/kernel_formulas/spd.py
+++ b/geometric_kernels/utils/kernel_formulas/spd.py
@@ -43,8 +43,10 @@ def _spd_heat_kernel_2x2_base(
     if x2 is None:
         x2 = x
 
-    assert x.shape == (2, 2)
-    assert x2.shape == (2, 2)
+    if B.shape(x) != (2, 2):
+        raise ValueError("`x` must have shape [2, 2].")
+    if x2.shape != (2, 2):
+        raise ValueError("`x2` must have shape [2, 2].")
 
     cl_1 = np.linalg.cholesky(x)
     cl_2 = np.linalg.cholesky(x2)
@@ -53,7 +55,8 @@ def _spd_heat_kernel_2x2_base(
     # Note: singular values that np.linalg.svd outputs are sorted, the following
     # code relies on this fact.
     H1, H2 = np.log(singular_values[0]), np.log(singular_values[1])
-    assert H1 >= H2
+    if H1 < H2:
+        raise RuntimeError("Expected `np.linalg.svd` to return sorted eigenvalues.")
 
     r_H_sq = H1 * H1 + H2 * H2
     alpha = H1 - H2
diff --git a/geometric_kernels/utils/manifold_utils.py b/geometric_kernels/utils/manifold_utils.py
index 42280cac..108e24f0 100644
--- a/geometric_kernels/utils/manifold_utils.py
+++ b/geometric_kernels/utils/manifold_utils.py
@@ -1,4 +1,4 @@
-""" Utilities for dealing with manifolds.  """
+"""Utilities for dealing with manifolds."""
 
 import lab as B
 import numpy as np
@@ -21,9 +21,11 @@ def minkowski_inner_product(vector_a: B.Numeric, vector_b: B.Numeric) -> B.Numer
     :return:
         An [...,]-shaped array of inner products.
     """
-    assert vector_a.shape == vector_b.shape
+    if B.shape(vector_a) != B.shape(vector_b):
+        raise ValueError("`vector_a` and `vector_b` must have the same shapes.")
     n = vector_a.shape[-1] - 1
-    assert n > 0
+    if n == 0:
+        raise ValueError("Must have at least 1 point.")
     diagonal = from_numpy(vector_a, [-1.0] + [1.0] * n)  # (n+1)
     diagonal = B.cast(B.dtype(vector_a), diagonal)
     return B.einsum("...i,...i->...", diagonal * vector_a, vector_b)
@@ -161,6 +163,7 @@ def tangent_onb(manifold, x):
     projected_onb_eigvals = projected_onb_eigvals[ambient_dim - manifold_dim :]
     projected_onb_eigvecs = projected_onb_eigvecs[:, ambient_dim - manifold_dim :]
 
-    assert np.all(np.isclose(projected_onb_eigvals, 1.0))
+    if not np.all(np.isclose(projected_onb_eigvals, 1.0)):
+        raise RuntimeError("Expected `projected_onb_eigvals` to be close to 1")
 
     return projected_onb_eigvecs
diff --git a/geometric_kernels/utils/product.py b/geometric_kernels/utils/product.py
index 58baded3..cfc182db 100644
--- a/geometric_kernels/utils/product.py
+++ b/geometric_kernels/utils/product.py
@@ -1,9 +1,10 @@
-""" Utilities for dealing with product spaces and product kernels.  """
+"""Utilities for dealing with product spaces and product kernels."""
 
 import lab as B
 from beartype.typing import Dict, List
 
 from geometric_kernels.lab_extras import smart_cast
+from geometric_kernels.utils.utils import _check_rank_1_array
 
 
 def params_to_params_list(
@@ -20,13 +21,20 @@ def params_to_params_list(
     :param params:
         Parameters of the product kernel.
     """
-    assert params["lengthscale"].shape == params["nu"].shape
-    assert len(params["nu"].shape) == 1
+    if B.shape(params["lengthscale"]) != B.shape(params["nu"]):
+        raise ValueError(
+            'Shape mismatch between `params["lengthscale"]` and `params["nu"].`'
+        )
+
+    _check_rank_1_array(params["nu"], 'params["nu"]')
 
     if params["nu"].shape[0] == 1:
         return [params] * number_of_factors
 
-    assert params["nu"].shape[0] == number_of_factors
+    if B.shape(params["nu"])[0] != number_of_factors:
+        raise ValueError(
+            "Shapes of the kernel parameters `lengthscale`, `nu` must be [`number_of_factors`]."
+        )
 
     list_of_params: List[Dict[str, B.Numeric]] = []
     for i in range(number_of_factors):
diff --git a/geometric_kernels/utils/special_functions.py b/geometric_kernels/utils/special_functions.py
index b222191c..41189737 100644
--- a/geometric_kernels/utils/special_functions.py
+++ b/geometric_kernels/utils/special_functions.py
@@ -11,6 +11,7 @@
     int_like,
     take_along_axis,
 )
+from geometric_kernels.utils.utils import _check_matrix
 
 
 def walsh_function(d: int, combination: List[int], x: B.Bool) -> B.Float:
@@ -35,8 +36,9 @@ def walsh_function(d: int, combination: List[int], x: B.Bool) -> B.Float:
         batch. An array of shape [N].
 
     """
-    assert x.ndim == 2
-    assert x.shape[-1] == d
+    _check_matrix(x, "x")
+    if B.shape(x)[-1] != d:
+        raise ValueError("`x` must live in `d`-dimensional space.")
 
     indices = B.cast(int_like(x), from_numpy(x, combination))[None, :]
 
@@ -91,9 +93,12 @@ def kravchuk_normalized(
     :return:
         $G_{d, j, m}/G_{d, j, 0}$ where $G_{d, j, m}$ is the Kravchuk polynomial.
     """
-    assert d > 0
-    assert 0 <= j and j <= d
-    assert B.all(0 <= m) and B.all(m <= d)
+    if d <= 0:
+        raise ValueError("`d` must be positive.")
+    if not (0 <= j and j <= d):
+        raise ValueError("`j` must lie in the interval [0, d].")
+    if not (B.all(0 <= m) and B.all(m <= d)):
+        raise ValueError("`m` must lie in the interval [0, d].")
 
     m = B.cast(B.dtype_float(m), m)
 
diff --git a/geometric_kernels/utils/utils.py b/geometric_kernels/utils/utils.py
index fc0b8d3d..452ca8ef 100644
--- a/geometric_kernels/utils/utils.py
+++ b/geometric_kernels/utils/utils.py
@@ -240,7 +240,7 @@ def partition_dominance_cone(partition: Tuple[int, ...]) -> Set[Tuple[int, ...]]
 
 
 def partition_dominance_or_subpartition_cone(
-    partition: Tuple[int, ...]
+    partition: Tuple[int, ...],
 ) -> Set[Tuple[int, ...]]:
     """
     Calculates subpartitions and partitions dominated by a given one and having
@@ -326,7 +326,8 @@ def log_binomial(n: B.Int, k: B.Int) -> B.Float:
     :return:
         The logarithm of the binomial coefficient binom(n, k).
     """
-    assert B.all(0 <= k <= n)
+    if not B.all(0 <= k <= n):
+        raise ValueError("Incorrect parameters of the binomial coefficient.")
 
     return B.loggamma(n + 1) - B.loggamma(k + 1) - B.loggamma(n - k + 1)
 
@@ -356,3 +357,39 @@ def binary_vectors_and_subsets(d: int):
             i += 1
 
     return x, combs
+
+
+def _check_field_in_params(params, field):
+    """
+    Raise an error if `params` does not contain a `field`.
+    """
+    if field not in params:
+        raise ValueError(f"`params` must contain `{field}`.")
+
+
+def _check_1_vector(x, desc):
+    """
+    Raise an error if `x` is not a vector of shape [1,].
+    """
+    if B.shape(x) != (1,):
+        raise ValueError(
+            f"`{desc}` must have shape `[1,]`, but has shape {B.shape(x)}."
+        )
+
+
+def _check_rank_1_array(x, desc):
+    """
+    Raise an error if `x` is not a rank-1 array.
+    """
+    if B.rank(x) != 1:
+        raise ValueError(
+            f"`{desc}` must have 1 dimension (`ndim` == 1), but has shape {B.shape(x)}."
+        )
+
+
+def _check_matrix(x, desc):
+    """
+    Raise an error if `x` is not a matrix.
+    """
+    if B.rank(x) != 2:
+        raise ValueError(f"`{desc}` must be a matrix, but has shape {B.shape(x)}.")