diff --git a/src/UQpy/dimension_reduction/diffusion_maps/DiffusionMaps.py b/src/UQpy/dimension_reduction/diffusion_maps/DiffusionMaps.py
index 5d47354e5..c03c2df3a 100644
--- a/src/UQpy/dimension_reduction/diffusion_maps/DiffusionMaps.py
+++ b/src/UQpy/dimension_reduction/diffusion_maps/DiffusionMaps.py
@@ -16,7 +16,7 @@
 from typing import Annotated, Union
 from beartype.vale import Is
 from UQpy.utilities.ValidationTypes import Numpy2DFloatArray, NumpyFloatArray
-from UQpy.utilities.kernels.baseclass import Kernel
+from UQpy.utilities.kernels.baseclass.Kernel import Kernel
 
 
 class DiffusionMaps:
diff --git a/src/UQpy/dimension_reduction/grassmann_manifold/projections/baseclass/GrassmannProjection.py b/src/UQpy/dimension_reduction/grassmann_manifold/projections/baseclass/GrassmannProjection.py
index b32af3a63..acc3ac6cd 100644
--- a/src/UQpy/dimension_reduction/grassmann_manifold/projections/baseclass/GrassmannProjection.py
+++ b/src/UQpy/dimension_reduction/grassmann_manifold/projections/baseclass/GrassmannProjection.py
@@ -1,7 +1,5 @@
 from abc import ABC, abstractmethod
 
-from UQpy.utilities.kernels.baseclass.Kernel import Kernel
-
 
 class GrassmannProjection(ABC):
     """
diff --git a/src/UQpy/surrogates/gaussian_process/__init__.py b/src/UQpy/surrogates/gaussian_process/__init__.py
index 83b0f56c8..d9439caf9 100644
--- a/src/UQpy/surrogates/gaussian_process/__init__.py
+++ b/src/UQpy/surrogates/gaussian_process/__init__.py
@@ -1,5 +1,4 @@
 from UQpy.surrogates.gaussian_process.GaussianProcessRegression import GaussianProcessRegression
 
-from UQpy.surrogates.gaussian_process.kernels import *
 from UQpy.surrogates.gaussian_process.regression_models import *
 from UQpy.surrogates.gaussian_process.constraints import *
diff --git a/src/UQpy/utilities/kernels/BinetCauchyKernel.py b/src/UQpy/utilities/kernels/BinetCauchyKernel.py
deleted file mode 100644
index 0d10fcb96..000000000
--- a/src/UQpy/utilities/kernels/BinetCauchyKernel.py
+++ /dev/null
@@ -1,25 +0,0 @@
-import numpy as np
-
-from UQpy.utilities.GrassmannPoint import GrassmannPoint
-from UQpy.utilities.kernels import GrassmannianKernel
-
-
-class BinetCauchyKernel(GrassmannianKernel):
-    """
-    A class to calculate the Binet-Cauchy kernel.
-
-    """
-    def apply_method(self, points):
-        points.evaluate_matrix(self, self.calculate_kernel_matrix)
-
-    def kernel_entry(self, xi: GrassmannPoint, xj: GrassmannPoint) -> float:
-        """
-        Compute the Binet-Cauchy kernel entry for two points on the Grassmann manifold.
-
-        :param xi: Orthonormal matrix representing the first point.
-        :param xj: Orthonormal matrix representing the second point.
-
-        """
-        r = np.dot(xi.data.T, xj.data)
-        det = np.linalg.det(r)
-        return det * det
diff --git a/src/UQpy/utilities/kernels/GaussianKernel.py b/src/UQpy/utilities/kernels/GaussianKernel.py
index 3eb0c03be..46326be26 100644
--- a/src/UQpy/utilities/kernels/GaussianKernel.py
+++ b/src/UQpy/utilities/kernels/GaussianKernel.py
@@ -1,10 +1,14 @@
+import itertools
+from typing import Tuple
+
 import numpy as np
 import scipy
-import scipy.spatial.distance as sd
+from scipy.spatial.distance import cdist
 
 from UQpy.utilities.ValidationTypes import RandomStateType, Numpy2DFloatArray
 from UQpy.utilities.kernels import EuclideanKernel
 from scipy.spatial.distance import pdist
+import scipy.spatial.distance as sd
 
 
 class GaussianKernel(EuclideanKernel):
@@ -15,26 +19,26 @@ class GaussianKernel(EuclideanKernel):
         k(x_j, x_i) = \exp[-(x_j - xj)^2/4\epsilon]
 
     """
-    def __init__(self, epsilon: float = 1.0):
+    def __init__(self, kernel_parameter: float = 1.0):
         """
         :param epsilon: Scale parameter of the Gaussian kernel
         """
-        super().__init__()
-        self.epsilon = epsilon
+        super().__init__(kernel_parameter=kernel_parameter)
 
-    def kernel_entry(self, xi: Numpy2DFloatArray, xj: Numpy2DFloatArray):
-        """
-        Given two points, this method computes the Gaussian kernel value between those two points
+    def calculate_kernel_matrix(self, x, s):
+        product = [self.element_wise_operation(point_pair)
+                   for point_pair in list(itertools.product(x, s))]
+        self.kernel_matrix = np.array(product).reshape(len(x), len(s))
+        return self.kernel_matrix
+
+    def element_wise_operation(self, xi_j: Tuple) -> float:
+        xi, xj = xi_j
 
-        :param xi: First point.
-        :param xj: Second point.
-        :return: Float representing the kernel entry.
-        """
         if len(xi.shape) == 1:
             d = pdist(np.array([xi, xj]), "sqeuclidean")
         else:
-            d = np.linalg.norm(xi-xj, 'fro') ** 2
-        return np.exp(-d / (2*self.epsilon**2))
+            d = np.linalg.norm(xi - xj, 'fro') ** 2
+        return np.exp(-d / (2 * self.kernel_parameter ** 2))
 
     def optimize_parameters(self, data: np.ndarray, tolerance: float,
                             n_nearest_neighbors: int,
diff --git a/src/UQpy/utilities/kernels/ProjectionKernel.py b/src/UQpy/utilities/kernels/ProjectionKernel.py
deleted file mode 100644
index c23353baf..000000000
--- a/src/UQpy/utilities/kernels/ProjectionKernel.py
+++ /dev/null
@@ -1,25 +0,0 @@
-import numpy as np
-
-from UQpy.utilities.GrassmannPoint import GrassmannPoint
-from UQpy.utilities.kernels import GrassmannianKernel
-
-
-class ProjectionKernel(GrassmannianKernel):
-    """
-    A class to calculate the Projection kernel
-
-    """
-    def apply_method(self, points):
-        points.evaluate_matrix(self, self.calculate_kernel_matrix)
-
-    def kernel_entry(self, xi: GrassmannPoint, xj: GrassmannPoint) -> float:
-        """
-        Compute the Projection kernel entry for two points on the Grassmann manifold.
-
-        :param xi: Orthonormal matrix representing the first point.
-        :param xj: Orthonormal matrix representing the second point.
-        """
-        r = np.dot(xi.data.T, xj.data)
-        n = np.linalg.norm(r, "fro")
-        kij = n * n
-        return kij
diff --git a/src/UQpy/utilities/kernels/__init__.py b/src/UQpy/utilities/kernels/__init__.py
index a46ed2d4f..6e3877dab 100644
--- a/src/UQpy/utilities/kernels/__init__.py
+++ b/src/UQpy/utilities/kernels/__init__.py
@@ -1,7 +1,7 @@
 from UQpy.utilities.kernels.baseclass import *
 from UQpy.utilities.kernels.euclidean_kernels import *
+from UQpy.utilities.kernels.grassmannian_kernels import *
 
-from UQpy.utilities.kernels.BinetCauchyKernel import BinetCauchyKernel
+from UQpy.utilities.kernels.grassmannian_kernels.BinetCauchyKernel import BinetCauchyKernel
 from UQpy.utilities.kernels.GaussianKernel import GaussianKernel
-from UQpy.utilities.kernels.ProjectionKernel import ProjectionKernel
 
diff --git a/src/UQpy/utilities/kernels/baseclass/EuclideanKernel.py b/src/UQpy/utilities/kernels/baseclass/EuclideanKernel.py
index 1fb0e62a5..a71f7997a 100644
--- a/src/UQpy/utilities/kernels/baseclass/EuclideanKernel.py
+++ b/src/UQpy/utilities/kernels/baseclass/EuclideanKernel.py
@@ -1,39 +1,7 @@
-import itertools
-from abc import ABC, abstractmethod
-from typing import Union
-import scipy.spatial.distance as sd
-import numpy as np
+from abc import ABC
 
-from UQpy.utilities import GrassmannPoint
-from UQpy.utilities.ValidationTypes import NumpyFloatArray, Numpy2DFloatArray
 from UQpy.utilities.kernels.baseclass.Kernel import Kernel
 
 
 class EuclideanKernel(Kernel, ABC):
     """This is a blueprint for Euclidean kernels implemented in the :py:mod:`kernels` module ."""
-
-    def __init__(self):
-        super().__init__()
-
-    def calculate_kernel_matrix(self, points: Numpy2DFloatArray):
-        """
-        Using the kernel-specific :py:meth:`.Kernel.kernel_entry` method, this function assembles the kernel matrix.
-
-        :param points: Set of data points in the Euclidean space
-
-        """
-        nargs = len(points)
-        indices = range(nargs)
-        pairs = list(itertools.combinations_with_replacement(indices, 2))
-        kernel = np.zeros((nargs, nargs))
-        for id_pair in range(np.shape(pairs)[0]):
-            i = pairs[id_pair][0]
-            j = pairs[id_pair][1]
-
-            xi = points[i]
-            xj = points[j]
-
-            kernel[i, j] = self.kernel_entry(xi, xj)
-            kernel[j, i] = kernel[i, j]
-
-        self.kernel_matrix = kernel
diff --git a/src/UQpy/utilities/kernels/baseclass/GrassmannianKernel.py b/src/UQpy/utilities/kernels/baseclass/GrassmannianKernel.py
index ec4249fe9..d2d219196 100644
--- a/src/UQpy/utilities/kernels/baseclass/GrassmannianKernel.py
+++ b/src/UQpy/utilities/kernels/baseclass/GrassmannianKernel.py
@@ -1,5 +1,6 @@
 import itertools
-from abc import ABC
+from abc import ABC, abstractmethod
+from typing import Union, Tuple
 
 import numpy as np
 
@@ -10,36 +11,21 @@
 class GrassmannianKernel(Kernel, ABC):
     """The parent class for Grassmannian kernels implemented in the :py:mod:`kernels` module ."""
 
-    def __init__(self):
-        super().__init__()
-
-    def calculate_kernel_matrix(self, points: list[GrassmannPoint], p: int = None):
+    def __init__(self, kernel_parameter: Union[int, float] = None):
         """
-        Compute the kernel matrix given a list of points on the Grassmann manifold.
-
-        :param points: Points on the Grassmann manifold
-        :param p: Number of independent p-planes of each Grassmann point.
-        :return: :class:`ndarray` The kernel matrix.
+        :param kernel_parameter: Number of independent p-planes of each Grassmann point.
         """
-        nargs = len(points)
-        # Define the pairs of points to compute the entries of the kernel matrix.
-        indices = range(nargs)
-        pairs = list(itertools.combinations_with_replacement(indices, 2))
-
-        # Estimate entries of the kernel matrix.
-        kernel = np.zeros((nargs, nargs))
-        for id_pair in range(np.shape(pairs)[0]):
-            i = pairs[id_pair][0]  # Point i
-            j = pairs[id_pair][1]  # Point j
-            if not p:
-                xi = points[i]
-                xj = points[j]
-            else:
-                xi = GrassmannPoint(points[i].data[:, :p])
-                xj = GrassmannPoint(points[j].data[:, :p])
-
-            # RiemannianDistance.check_rows(xi, xj)
-            kernel[i, j] = self.kernel_entry(xi, xj)
-            kernel[j, i] = kernel[i, j]
-
-        self.kernel_matrix = kernel
+        super().__init__(kernel_parameter)
+
+    def calculate_kernel_matrix(self, x: list[GrassmannPoint], s: list[GrassmannPoint]):
+        p = self.kernel_parameter
+        list1 = [point.data if not p else point.data[:, :p] for point in x]
+        list2 = [point.data if not p else point.data[:, :p] for point in s]
+        product = [self.element_wise_operation(point_pair)
+                   for point_pair in list(itertools.product(list1, list2))]
+        self.kernel_matrix = np.array(product).reshape(len(list1), len(list2))
+        return self.kernel_matrix
+
+    @abstractmethod
+    def element_wise_operation(self, xi_j: Tuple) -> float:
+        pass
diff --git a/src/UQpy/utilities/kernels/baseclass/Kernel.py b/src/UQpy/utilities/kernels/baseclass/Kernel.py
index 5724147cb..4b390b494 100644
--- a/src/UQpy/utilities/kernels/baseclass/Kernel.py
+++ b/src/UQpy/utilities/kernels/baseclass/Kernel.py
@@ -12,6 +12,7 @@ class Kernel(ABC):
 
     def __init__(self, kernel_parameter: Union[int, float]):
         self.__kernel_parameter = kernel_parameter
+        self.kernel_matrix=None
 
     @property
     def kernel_parameter(self):
diff --git a/src/UQpy/utilities/kernels/baseclass/__init__.py b/src/UQpy/utilities/kernels/baseclass/__init__.py
index 1b8fdbf8c..e7ff0f2ce 100644
--- a/src/UQpy/utilities/kernels/baseclass/__init__.py
+++ b/src/UQpy/utilities/kernels/baseclass/__init__.py
@@ -1,3 +1,2 @@
-from UQpy.utilities.kernels.baseclass.DRKernel import Kernel
 from UQpy.utilities.kernels.baseclass.EuclideanKernel import EuclideanKernel
 from UQpy.utilities.kernels.baseclass.GrassmannianKernel import GrassmannianKernel
diff --git a/src/UQpy/utilities/kernels/euclidean_kernels/Matern.py b/src/UQpy/utilities/kernels/euclidean_kernels/Matern.py
index d98d79fe9..e21a72d1a 100644
--- a/src/UQpy/utilities/kernels/euclidean_kernels/Matern.py
+++ b/src/UQpy/utilities/kernels/euclidean_kernels/Matern.py
@@ -1,9 +1,13 @@
-from UQpy.utilities.kernels.baseclass.Kernel import *
+from typing import Union
+
+import numpy as np
+
+from UQpy.utilities.kernels.baseclass.EuclideanKernel import EuclideanKernel
 from scipy.spatial.distance import cdist
 from scipy.special import gamma, kv
 
 
-class Matern(Kernel):
+class Matern(EuclideanKernel):
     def __init__(self, kernel_parameter: Union[int, float] = 1, nu=1.5):
         """
         Matern Kernel is a generalization of Radial Basis Function kernel.
@@ -18,16 +22,17 @@ def calculate_kernel_matrix(self, x, s):
         l = self.kernel_parameter
         stack = cdist(x / l, s / l, metric='euclidean')
         if self.nu == 0.5:
-            return np.exp(-np.abs(stack))
+            self.kernel_matrix = np.exp(-np.abs(stack))
         elif self.nu == 1.5:
-            return (1 + np.sqrt(3) * stack) * np.exp(-np.sqrt(3) * stack)
+            self.kernel_matrix = (1 + np.sqrt(3) * stack) * np.exp(-np.sqrt(3) * stack)
         elif self.nu == 2.5:
-            return (1 + np.sqrt(5) * stack + 5 * (stack ** 2) / 3) * np.exp(-np.sqrt(5) * stack)
+            self.kernel_matrix = (1 + np.sqrt(5) * stack + 5 * (stack ** 2) / 3) * np.exp(-np.sqrt(5) * stack)
         elif self.nu == np.inf:
-            return np.exp(-(stack ** 2) / 2)
+            self.kernel_matrix = np.exp(-(stack ** 2) / 2)
         else:
             stack *= np.sqrt(2 * self.nu)
             tmp = 1 / (gamma(self.nu) * (2 ** (self.nu - 1)))
             tmp1 = stack ** self.nu
             tmp2 = kv(self.nu, stack)
-            return tmp * tmp1 * tmp2
+            self.kernel_matrix = tmp * tmp1 * tmp2
+        return self.kernel_matrix
diff --git a/src/UQpy/utilities/kernels/euclidean_kernels/RBF.py b/src/UQpy/utilities/kernels/euclidean_kernels/RBF.py
index 2a8206055..768e4372b 100644
--- a/src/UQpy/utilities/kernels/euclidean_kernels/RBF.py
+++ b/src/UQpy/utilities/kernels/euclidean_kernels/RBF.py
@@ -1,7 +1,12 @@
-from UQpy.utilities.kernels.baseclass.Kernel import *
+from typing import Union
 
+import numpy as np
 
-class RBF(Kernel):
+from UQpy.utilities.kernels.baseclass.EuclideanKernel import EuclideanKernel
+from UQpy.utilities.kernels.baseclass.Kernel import Kernel
+
+
+class RBF(EuclideanKernel):
     def __init__(self, kernel_parameter: Union[int, float] = 1.0):
         super().__init__(kernel_parameter)
 
@@ -13,4 +18,5 @@ def calculate_kernel_matrix(self, x, s):
         :params s: An array containing input points.
         """
         stack = Kernel.check_samples_and_return_stack(x / self.kernel_parameter, s / self.kernel_parameter)
-        return np.exp(np.sum(-0.5 * (stack ** 2), axis=2))
+        self.kernel_matrix = np.exp(np.sum(-0.5 * (stack ** 2), axis=2))
+        return self.kernel_matrix
diff --git a/src/UQpy/utilities/kernels/grassmannian_kernels/BinetCauchyKernel.py b/src/UQpy/utilities/kernels/grassmannian_kernels/BinetCauchyKernel.py
new file mode 100644
index 000000000..91b9eeec7
--- /dev/null
+++ b/src/UQpy/utilities/kernels/grassmannian_kernels/BinetCauchyKernel.py
@@ -0,0 +1,22 @@
+from typing import Tuple
+
+import numpy as np
+
+from UQpy.utilities.kernels import GrassmannianKernel
+
+
+class BinetCauchyKernel(GrassmannianKernel):
+    """
+    A class to calculate the Binet-Cauchy kernel.
+
+    """
+    def element_wise_operation(self, xi_j: Tuple) -> float:
+        """
+        Compute the Projection kernel entry for a tuple of points on the Grassmann manifold.
+
+        :param xi_j: Tuple of orthonormal matrices representing the grassmann points.
+        """
+        xi, xj = xi_j
+        r = np.dot(xi.T, xj)
+        det = np.linalg.det(r)
+        return det * det
diff --git a/src/UQpy/utilities/kernels/grassmannian_kernels/ProjectionKernel.py b/src/UQpy/utilities/kernels/grassmannian_kernels/ProjectionKernel.py
new file mode 100644
index 000000000..a98fa46b4
--- /dev/null
+++ b/src/UQpy/utilities/kernels/grassmannian_kernels/ProjectionKernel.py
@@ -0,0 +1,25 @@
+from typing import Union, Tuple
+
+import numpy as np
+
+from UQpy.utilities.kernels.baseclass.GrassmannianKernel import GrassmannianKernel
+
+
+class ProjectionKernel(GrassmannianKernel):
+
+    def __init__(self, kernel_parameter: Union[int, float] = None):
+        """
+        :param kernel_parameter: Number of independent p-planes of each Grassmann point.
+        """
+        super().__init__(kernel_parameter)
+
+    def element_wise_operation(self, xi_j: Tuple) -> float:
+        """
+        Compute the Projection kernel entry for a tuple of points on the Grassmann manifold.
+
+        :param xi_j: Tuple of orthonormal matrices representing the grassmann points.
+        """
+        xi, xj = xi_j
+        r = np.dot(xi.T, xj)
+        n = np.linalg.norm(r, "fro")
+        return n * n
diff --git a/src/UQpy/utilities/kernels/grassmannian_kernels/__init__.py b/src/UQpy/utilities/kernels/grassmannian_kernels/__init__.py
new file mode 100644
index 000000000..43b4f509f
--- /dev/null
+++ b/src/UQpy/utilities/kernels/grassmannian_kernels/__init__.py
@@ -0,0 +1,2 @@
+from UQpy.utilities.kernels.grassmannian_kernels.ProjectionKernel import ProjectionKernel
+from UQpy.utilities.kernels.grassmannian_kernels.BinetCauchyKernel import BinetCauchyKernel
diff --git a/tests/unit_tests/dimension_reduction/test_kernel.py b/tests/unit_tests/dimension_reduction/test_kernel.py
index c48f012a2..53434368a 100644
--- a/tests/unit_tests/dimension_reduction/test_kernel.py
+++ b/tests/unit_tests/dimension_reduction/test_kernel.py
@@ -1,9 +1,7 @@
-import sys
-
+from UQpy.utilities import ProjectionKernel
 from UQpy.utilities.GrassmannPoint import GrassmannPoint
 from UQpy.dimension_reduction.grassmann_manifold.projections.SVDProjection import SVDProjection
-from UQpy.utilities.kernels.ProjectionKernel import ProjectionKernel
-from UQpy.utilities.kernels.BinetCauchyKernel import BinetCauchyKernel
+from UQpy.utilities.kernels.grassmannian_kernels.BinetCauchyKernel import BinetCauchyKernel
 from UQpy.utilities.kernels.GaussianKernel import GaussianKernel
 import numpy as np
 
@@ -16,12 +14,11 @@ def test_kernel_projection():
                                   [-0.63305978, 0.51850616]]))
     points = [xi, xj, xk]
     k = ProjectionKernel()
-    k.calculate_kernel_matrix(points)
+    k.calculate_kernel_matrix(points, points)
     kernel = np.matrix.round(k.kernel_matrix, 4)
 
     assert np.allclose(kernel, np.array([[2, 1.0063, 1.2345], [1.0063, 2, 1.0101], [1.2345, 1.0101, 2]]))
 
-
 def test_kernel_binet_cauchy():
     xi = GrassmannPoint(np.array([[-np.sqrt(2) / 2, -np.sqrt(2) / 4], [np.sqrt(2) / 2, -np.sqrt(2) / 4],
                                   [0, -np.sqrt(3) / 2]]))
@@ -30,7 +27,7 @@ def test_kernel_binet_cauchy():
     points = [xi, xj, xk]
 
     kernel = BinetCauchyKernel()
-    kernel.calculate_kernel_matrix(points)
+    kernel.calculate_kernel_matrix(points, points)
     kernel = np.matrix.round(kernel.kernel_matrix, 4)
 
     assert np.allclose(kernel, np.array([[1, 0.0063, 0.2345], [0.0063, 1, 0.0101], [0.2345, 0.0101, 1]]))
@@ -41,13 +38,13 @@ def test_kernel_gaussian_1d():
     xj = np.array([0.2, -1, 3, 5])
     xk = np.array([1, 2, 3, 4])
     points = [xi, xj, xk]
-    gaussian = GaussianKernel(epsilon=2.0)
-    gaussian.calculate_kernel_matrix(points)
+    gaussian = GaussianKernel(kernel_parameter=2.0)
+    gaussian.calculate_kernel_matrix(points, points)
 
     assert np.allclose(np.matrix.round(gaussian.kernel_matrix, 4),
                        np.array([[1., 0.26447726, 1.], [0.26447726, 1., 0.26447726], [1, 0.26447726, 1]]),
                        atol=1e-04)
-    assert np.round(gaussian.epsilon, 4) == 2
+    assert np.round(gaussian.kernel_parameter, 4) == 2
 
 
 def test_kernel_gaussian_2d():
@@ -56,12 +53,12 @@ def test_kernel_gaussian_2d():
     xk = np.array([[-0.69535592, -0.0546034], [-0.34016974, -0.85332868], [-0.63305978, 0.51850616]])
     points = [xi, xj, xk]
     gaussian = GaussianKernel()
-    gaussian.calculate_kernel_matrix(points)
+    gaussian.calculate_kernel_matrix(points, points)
 
     assert np.allclose(np.matrix.round(gaussian.kernel_matrix, 4), np.array([[1., 0.39434829, 0.15306655],
                                                                              [0.39434829, 1., 0.06422136],
                                                                              [0.15306655, 0.06422136, 1.]]), atol=1e-4)
-    assert np.round(gaussian.epsilon, 4) == 1.0
+    assert np.round(gaussian.kernel_parameter, 4) == 1.0
 
 
 sol0 = np.array([[0.61415, 1.03029, 1.02001, 0.57327, 0.79874, 0.73274],
@@ -104,6 +101,6 @@ def test_kernel():
     manifold_projection = SVDProjection(Solutions, p="max")
     kernel = ProjectionKernel()
 
-    kernel.calculate_kernel_matrix(manifold_projection.u)
+    kernel.calculate_kernel_matrix(manifold_projection.u, manifold_projection.u)
 
     assert np.round(kernel.kernel_matrix[0, 1], 8) == 6.0
diff --git a/tests/unit_tests/surrogates/test_gpr.py b/tests/unit_tests/surrogates/test_gpr.py
index e5234aa38..9a8a1be8f 100644
--- a/tests/unit_tests/surrogates/test_gpr.py
+++ b/tests/unit_tests/surrogates/test_gpr.py
@@ -1,18 +1,15 @@
 import pytest
 from beartype.roar import BeartypeCallHintPepParamException
 
+from UQpy.utilities.kernels.euclidean_kernels import RBF, Matern
 from UQpy.utilities.MinimizeOptimizer import MinimizeOptimizer
 from UQpy.utilities.FminCobyla import FminCobyla
 from UQpy.surrogates.gaussian_process.GaussianProcessRegression import GaussianProcessRegression
-# # from UQpy.utilities.strata.Rectangular import Rectangular
-# # from UQpy.sampling.StratifiedSampling import StratifiedSampling
-# # from UQpy.RunModel import RunModel
-# # from UQpy.distributions.collection.Uniform import Uniform
 import numpy as np
 import shutil
 from UQpy.surrogates.gaussian_process.constraints import NonNegative
 from UQpy.surrogates.gaussian_process.regression_models import LinearRegression, ConstantRegression, QuadraticRegression
-from UQpy.surrogates.gaussian_process.kernels import RBF, Matern
+
 
 samples = np.linspace(0, 5, 20).reshape(-1, 1)
 values = np.cos(samples)