diff --git a/src/hidimstat/adaptive_permutation_threshold.py b/src/hidimstat/adaptive_permutation_threshold.py
index 933fb9e..7132f5d 100644
--- a/src/hidimstat/adaptive_permutation_threshold.py
+++ b/src/hidimstat/adaptive_permutation_threshold.py
@@ -37,9 +37,8 @@ def ada_svr(X, y, rcond=1e-3):
     """
 
     X = np.asarray(X)
-    n_samples, n_features = X.shape
 
-    K = _manual_inversion(np.dot(X, X.T), rcond=rcond)
+    K = np.linalg.pinv(np.dot(X, X.T), rcond=rcond)
     sum_K = np.sum(K)
 
     L = -np.outer(np.sum(K, axis=0), np.sum(K, axis=1)) / sum_K
@@ -50,25 +49,3 @@ def ada_svr(X, y, rcond=1e-3):
     scale = np.sqrt(np.sum(C**2, axis=1))
 
     return beta_hat, scale
-
-
-def _manual_inversion(X, rcond=1e-3, full_rank=False):
-    "Inverting taking care of low eigenvalues to increase numerical stability"
-
-    X = np.asarray(X)
-    n_samples, n_features = X.shape
-
-    if n_samples != n_features:
-        raise ValueError("The matrix is not a square matrix")
-
-    U, s, V = np.linalg.svd(X, full_matrices=False)
-    rank = np.sum(s > rcond * s.max())
-    s_inv = np.zeros(np.size(s))
-    s_inv[:rank] = 1 / s[:rank]
-
-    if full_rank:
-        s_inv[rank:] = 1 / (rcond * s.max())
-
-    X_inv = np.linalg.multi_dot([U, np.diag(s_inv), V])
-
-    return X_inv
diff --git a/test/test_adaptive_permutation_threshold.py b/test/test_adaptive_permutation_threshold.py
index 7886253..6fca754 100644
--- a/test/test_adaptive_permutation_threshold.py
+++ b/test/test_adaptive_permutation_threshold.py
@@ -11,32 +11,33 @@
 
 
 def test_ada_svr():
-    """Testing the procedure on a simulation with no structure and a support
+    """
+    Testing the procedure on a simulation with no structure and a support
     of size 1. Computing one-sided p-values, we want a low p-value
-    for the first feature and p-values close to 0.5 for the others."""
+    for the first feature and p-values close to 0.5 for the others.
+    """
 
+    # Parameters for the generation of data
     n_samples, n_features = 20, 50
-    support_size = 1
-    sigma = 0.1
-    rho = 0.0
+    support_size = 4
 
-    X_init, y, beta, noise = multivariate_1D_simulation(
+    X_init, y, _, _ = multivariate_1D_simulation(
         n_samples=n_samples,
         n_features=n_features,
         support_size=support_size,
-        sigma=sigma,
-        rho=rho,
+        sigma=0.1,
         shuffle=False,
-        seed=3,
+        seed=42,
     )
 
-    y = y - np.mean(y)
-    X_init = X_init - np.mean(X_init, axis=0)
-
+    # Run the procedure
     beta_hat, scale_hat = ada_svr(X_init, y)
 
+    # Compute p-values
     pval, pval_corr, _, _ = pval_from_scale(beta_hat, scale_hat)
 
+    # Check that the p-values are close to 0.5 for the features not in the support
+    # and close to 0 for the feature in the support
     expected = 0.5 * np.ones(n_features)
     expected[:support_size] = 0.0