diff --git a/src/UQpy/utilities/kernels/GaussianKernel.py b/src/UQpy/utilities/kernels/GaussianKernel.py
index b5bf5ec77..28ac42a78 100644
--- a/src/UQpy/utilities/kernels/GaussianKernel.py
+++ b/src/UQpy/utilities/kernels/GaussianKernel.py
@@ -4,6 +4,7 @@
 
 from UQpy.utilities.ValidationTypes import RandomStateType
 from UQpy.utilities.kernels import EuclideanKernel
+from scipy.spatial.distance import pdist
 
 
 class GaussianKernel(EuclideanKernel):
@@ -30,11 +31,8 @@ def kernel_entry(self, xi, xj):
         :param xj: Second point.
         :return: Float representing the kernel entry.
         """
-
-        return np.linalg.norm(xi - xj, "fro") ** 2
-
-    def kernel_function(self, distance_pairs):
-        return np.exp(-sd.squareform(distance_pairs) / (4 * self.epsilon))
+        d = pdist([xi, xj], "sqeuclidean")
+        return np.exp(-d ** 2 / (4*self.epsilon))
 
     def optimize_parameters(self, data: np.ndarray, tolerance: float,
                             n_nearest_neighbors: int,
diff --git a/src/UQpy/utilities/kernels/baseclass/EuclideanKernel.py b/src/UQpy/utilities/kernels/baseclass/EuclideanKernel.py
index 997c9e24c..81d6ae5a1 100644
--- a/src/UQpy/utilities/kernels/baseclass/EuclideanKernel.py
+++ b/src/UQpy/utilities/kernels/baseclass/EuclideanKernel.py
@@ -21,26 +21,18 @@ def calculate_kernel_matrix(self, points: Union[list, NumpyFloatArray]):
         :param points: Coordinates of the points in the Euclidean space
 
         """
-        distance_pairs = None
-        if len(np.shape(points)) == 2:
-            distance_pairs = sd.pdist(points, "sqeuclidean")
+        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]
 
-        elif len(np.shape(points)) == 3:
-            nargs = len(points)
-            indices = range(nargs)
-            pairs = list(itertools.combinations(indices, 2))
-            distance_pairs = []
-            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]
 
-                xi = points[i]
-                xj = points[j]
+            kernel[i, j] = self.kernel_entry(xi, xj)
+            kernel[j, i] = kernel[i, j]
 
-                distance_pairs.append(self.kernel_entry(xi, xj))
-
-        self.kernel_matrix = self.kernel_function(distance_pairs)
-
-    @abstractmethod
-    def kernel_function(self, distance_pairs):
-        pass
+        self.kernel_matrix = kernel