diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index a4b18cd..b8e61cd 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -7,7 +7,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.8, 3.9]
+        python-version: [3.10, 3.11]
     steps:
       - uses: actions/checkout@v2
       - name: Set up Python ${{ matrix.python-version }}
diff --git a/hvd/newton.py b/hvd/newton.py
index 0dcb8c8..7f3c0fb 100644
--- a/hvd/newton.py
+++ b/hvd/newton.py
@@ -855,12 +855,6 @@ def _compute_netwon_step(self) -> Tuple[np.ndarray, np.ndarray]:
                     newton_step[r, c] = (
                         -1 * np.linalg.lstsq(DR, R_list[r].reshape(-1, 1), rcond=None)[0].ravel()
                     )
-        # heuristic: prevent the step-length to be too large
-        # X = newton_step[:, : self.dim_p]
-        # t = np.max(self.xu - self.xl) / 4
-        # idx = np.linalg.norm(X, axis=1) >= t
-        # X[idx] /= np.linalg.norm(X[idx], axis=1).reshape(-1, 1) / t
-        # newton_step[:, : self.dim_p] = X
         return newton_step, R
 
     def _shift_reference_set(self):
@@ -892,22 +886,7 @@ def _shift_reference_set(self):
         # shift the medoids
         for i, k in enumerate(indices):
             n = self._eta[self.Y_label[k]]
-            # NOTE: a larger initial shift is needed for ZDT6
-            if 11 < 2:
-                if self.iter_count == 0:
-                    m = self.active_indicator._medoids[k]
-                    a = n[1] / n[0]
-                    b = -1
-                    c = m[1] - (n[1] * m[0]) / n[0]
-                    idx = np.argmin(
-                        [abs(a * p[0] + b * p[1] + c) / np.sqrt(a**2 + b**2) for p in self._pareto_front]
-                    )
-                    eps = np.inner(self._pareto_front[idx] - m, n) * 1.03
-                    v = eps * n
-                else:
-                    v = 0.05 * n
-            else:
-                v = 0.05 * n if self.iter_count > 0 else 0.01 * n  # the initial shift is a bit larger
+            v = 0.05 * n if self.iter_count > 0 else 0.01 * n  # the initial shift is a bit larger
             self.active_indicator.shift_medoids(v, k)
 
         if self.iter_count == 0:  # record the initial medoids
@@ -934,6 +913,7 @@ def phi_func(alpha, i):
             self.active_indicator.re_match = False
             R_ = self._compute_R(state)[0][i]
             self.active_indicator.re_match = True
+            self.state.n_jac_evals = state.n_jac_evals
             return np.linalg.norm(R_)
 
         if max_step_size is not None:
diff --git a/hvd/problems/base.py b/hvd/problems/base.py
index e9dc28c..37129ae 100644
--- a/hvd/problems/base.py
+++ b/hvd/problems/base.py
@@ -47,12 +47,12 @@ def __init__(self):
             eq_func = partial(self.__class__._eq_constraint, self)
         if self.n_ieq_constr > 0:
             ieq_func = partial(self.__class__._ieq_constraint, self)
-        self._eq_constraint_jacobian = jacrev(eq_func) if self.n_eq_constr > 0 else None
+        self._eq_constraint_jacobian = jit(jacrev(eq_func)) if self.n_eq_constr > 0 else None
         self._eq_constraint_hessian = hessian(eq_func) if self.n_eq_constr > 0 else None
-        self._ieq_constraint_jacobian = jacrev(ieq_func) if self.n_ieq_constr > 0 else None
+        self._ieq_constraint_jacobian = jit(jacrev(ieq_func)) if self.n_ieq_constr > 0 else None
         self._ieq_constraint_hessian = hessian(ieq_func) if self.n_ieq_constr > 0 else None
 
-    def _eq_constraint(self, x: np.ndarray) -> np.ndarray:
+    def eq_constraint(self, x: np.ndarray) -> np.ndarray:
         return np.array(self._eq_constraint(x))
 
     def ieq_constraint(self, x: np.ndarray) -> np.ndarray:
diff --git a/hvd/utils.py b/hvd/utils.py
index a6cd667..37c17ca 100644
--- a/hvd/utils.py
+++ b/hvd/utils.py
@@ -86,7 +86,7 @@ def precondition_hessian(H: np.ndarray) -> np.ndarray:
         v = np.min(np.diag(H))
         tau = 0 if v > 0 else -v + beta
         I = np.eye(H.shape[0])
-        for _ in range(35):
+        for _ in range(45):
             try:
                 L = cholesky(H + tau * I, lower=True)
                 break
diff --git a/requirements.txt b/requirements.txt
index 26eb00c..6a4a754 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,40 +1,25 @@
 about-time==4.2.1
 alive-progress==3.1.5
-asttokens==2.4.1
 autograd==1.6.2
 cma==3.2.2
 contourpy==1.2.0
 cycler==0.12.1
-decorator==5.1.1
 Deprecated==1.2.14
 dill==0.3.7
-executing==2.0.1
-fonttools==4.46.0
+fonttools==4.47.0
 future==0.18.3
 grapheme==0.6.0
-ipdb==0.13.13
-ipython==8.18.1
-jax==0.4.21
-jaxlib==0.4.21
-jedi==0.19.1
-Jinja2==3.1.2
+jax==0.4.23
+jaxlib==0.4.23
 joblib==1.3.2
 kiwisolver==1.4.5
-MarkupSafe==2.1.3
 matplotlib==3.8.2
-matplotlib-inline==0.1.6
 ml-dtypes==0.3.1
 numpy==1.26.2
 opt-einsum==3.3.0
 packaging==23.2
 pandas==2.1.4
-parso==0.8.3
-pexpect==4.9.0
 Pillow==10.1.0
-prompt-toolkit==3.0.41
-ptyprocess==0.7.0
-pure-eval==0.2.2
-Pygments==2.17.2
 pymoo==0.6.1.1
 pyparsing==3.1.1
 python-dateutil==2.8.2
@@ -43,9 +28,6 @@ scikit-learn==1.3.2
 scikit-learn-extra==0.3.0
 scipy==1.11.4
 six==1.16.0
-stack-data==0.6.3
 threadpoolctl==3.2.0
-traitlets==5.14.0
 tzdata==2023.3
-wcwidth==0.2.12
 wrapt==1.16.0
diff --git a/scripts/compute_statistics.py b/scripts/compute_statistics.py
index ae2f630..484d926 100644
--- a/scripts/compute_statistics.py
+++ b/scripts/compute_statistics.py
@@ -12,11 +12,11 @@
     # "NSGA-III",
 ]
 problems = [
-    # "ZDT1",
-    # "ZDT2",
-    # "ZDT3",
+    "ZDT1",
+    "ZDT2",
+    "ZDT3",
     "ZDT4",
-    # "ZDT6",
+    "ZDT6",
     # "DTLZ1",
     # "DTLZ2",
     # "DTLZ3",
@@ -35,7 +35,7 @@
         pvalues.append(pvalue)
         stats.append([stats_func(x), stats_func(y)])
 
-reject, pvals_corrected, _, _ = multipletests(pvalues, alpha=0.05)
+reject, pvals_corrected, _, _ = multipletests(pvalues, alpha=0.05, method="fdr_bh")
 win, tie, loose = 0, 0, 0
 for k, (n, m) in enumerate(stats):
     x, y = n[0], m[0]
@@ -55,7 +55,8 @@
 Problem = np.tile(problems, len(MOEAs)).reshape(-1, 1)
 summary = ["+/$\\approx$/-", "", f"{win}/{tie}/{loose}", ""]
 data = pd.DataFrame(
-    np.vstack([np.hstack([Method, Problem, stats]), summary]), columns=["Method", "Problem", "Newton", "MOEA"]
+    np.vstack([np.hstack([Method, Problem, stats]), summary]),
+    columns=["Method", "Problem", "Newton (iter = 5)", "MOEA (iter = 300 + 5 \\times (4 + 10n)"],
 )
 print(data)
 caption = """Warmstarting the Newton method after 1\,510 iterations of the MOEA, we show the IGD values 
diff --git a/scripts/measure_CPU_time.py b/scripts/measure_CPU_time.py
index 926e53e..c7ec943 100644
--- a/scripts/measure_CPU_time.py
+++ b/scripts/measure_CPU_time.py
@@ -10,7 +10,7 @@
 
 np.random.seed(42)
 
-N = 1e4
+N = 1e5
 res = []
 problems = [ZDT1(), ZDT2(), ZDT3(), ZDT4(), ZDT6()]
 for problem in problems:
@@ -24,7 +24,7 @@
         FE_time.append((t1 - t0) / 1e9)
     res.append(["FE", type(problem).__name__, np.mean(FE_time), np.median(FE_time), np.std(FE_time)])
 
-N = 1e4
+N = 1e5
 for problem in problems:
     f = PymooProblemWithAD(problem)
     FE_time = []
@@ -38,9 +38,9 @@
         # Y = func(x)  # `(N, dim_m)`
         # H = h(x)
         # Jacobians
-        YdX = jac(x)  # `(N, dim_m, dim_d)`
+        # YdX = jac(x)  # `(N, dim_m, dim_d)`
         # Hessians
-        # YdX2 = hessian(x)  # `(N, dim_m, dim_d, dim_d)`
+        YdX2 = hessian(x)  # `(N, dim_m, dim_d, dim_d)`
         # HdX = h_jac(x)
         # HdX2 = h_hessian(x)
         t1 = time.process_time_ns()
diff --git a/tests/test_CFs.py b/tests/test_CFs.py
index a91ed92..98477a3 100644
--- a/tests/test_CFs.py
+++ b/tests/test_CFs.py
@@ -36,7 +36,7 @@ def test_CFs():
     p.objective_jacobian(x)
     p.ieq_constraint(x)
     assert np.all(np.isclose(F, [13.53009123, 17.13769761]))
-    assert np.all(np.isclose(C, [-199.6512955]))
+    assert np.all(np.isclose(C, [-199.64903]))
 
     p = CF4()
     interval = p.xu - p.xl