diff --git a/.github/workflows/build-and-test.yaml b/.github/workflows/build-and-test.yaml
index 6092b78..aafc0c6 100644
--- a/.github/workflows/build-and-test.yaml
+++ b/.github/workflows/build-and-test.yaml
@@ -19,7 +19,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.8", "3.9", "3.10"]
+        python-version: ["3.8", "3.9", "3.10", "3.11"]
         os: [ubuntu-20.04, windows-2022]
         openmp-on: [true, false]
         exclude:
diff --git a/README.md b/README.md
index e51fe12..9fee4af 100644
--- a/README.md
+++ b/README.md
@@ -136,11 +136,11 @@ Execute an example script for measuring time.
 OMP_NUM_THREADS=1 python time_measurement.py
 
 # Results [s]
-# ICP(Open3D):  0.02030642901081592
-# CPD:  3.6435861150093842
-# SVR:  0.5795929960149806
-# GMMTree:  0.34479290700983256
-# FilterReg:  0.039795294986106455
+# ICP(Open3D):  0.02212289097951725
+# CPD:  0.1632105699973181
+# SVR:  0.1983455469890032
+# GMMTree:  0.27672973900916986
+# FilterReg:  0.015339541976572946
 ```
 
 ## Citing
diff --git a/probreg/cpd.py b/probreg/cpd.py
index 479505c..67a6967 100644
--- a/probreg/cpd.py
+++ b/probreg/cpd.py
@@ -8,6 +8,7 @@
 import numpy as np
 import open3d as o3
 import six
+from scipy.spatial import distance as scipy_distance
 
 from . import math_utils as mu
 from . import transformation as tf
@@ -44,14 +45,17 @@ def __init__(self, source: Optional[np.ndarray] = None, use_cuda: bool = False)
         self._callbacks = []
         if use_cuda:
             import cupy as cp
+            from cupyx.scipy.spatial import distance as cupy_distance
 
             from . import cupy_utils
 
             self.xp = cp
+            self.distance_module = cupy_distance
             self.cupy_utils = cupy_utils
             self._squared_kernel_sum = cupy_utils.squared_kernel_sum
         else:
             self.xp = np
+            self.distance_module = scipy_distance
             self._squared_kernel_sum = mu.squared_kernel_sum
 
     def set_source(self, source: np.ndarray) -> None:
@@ -67,7 +71,8 @@ def _initialize(self, target: np.ndarray) -> MstepResult:
     def expectation_step(self, t_source: np.ndarray, target: np.ndarray, sigma2: float, w: float = 0.0) -> EstepResult:
         """Expectation step for CPD"""
         assert t_source.ndim == 2 and target.ndim == 2, "source and target must have 2 dimensions."
-        pmat = self.xp.stack([self.xp.sum(self.xp.square(target - ts), axis=1) for ts in t_source])
+        pmat = self.distance_module.cdist(t_source, target, "sqeuclidean")
+        # pmat = self.xp.stack([self.xp.sum(self.xp.square(target - ts), axis=1) for ts in t_source])
         pmat = self.xp.exp(-pmat / (2.0 * sigma2))
 
         c = (2.0 * np.pi * sigma2) ** (t_source.shape[1] * 0.5)
diff --git a/probreg/version.py b/probreg/version.py
index d7b30e1..8879c6c 100644
--- a/probreg/version.py
+++ b/probreg/version.py
@@ -1 +1 @@
-__version__ = "0.3.6"
+__version__ = "0.3.7"
diff --git a/pyproject.toml b/pyproject.toml
index 942d813..1ad3390 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,7 +1,7 @@
 [tool.poetry]
 name = "probreg"
-version = "0.3.6"
-description = ""
+version = "0.3.7"
+description = "Probablistic point cloud resitration algorithms"
 authors = ["nekanat <nekanat.stock@gmail.com>"]
 license = "MIT"