diff --git a/test/property_tests/test_optimizers_complexity.py b/test/property_tests/test_optimizers_complexity.py
index c893ee645..1bbe42aab 100644
--- a/test/property_tests/test_optimizers_complexity.py
+++ b/test/property_tests/test_optimizers_complexity.py
@@ -3,6 +3,8 @@
 from hypothesis import given
 from hypothesis import settings
 from hypothesis.strategies import integers
+from numpy.random import randint
+from numpy.random import seed
 from sklearn.datasets import make_regression
 from sklearn.linear_model import ElasticNet
 from sklearn.linear_model import Lasso
@@ -20,7 +22,7 @@
     n_informative=integers(min_value=3, max_value=9),
     random_state=integers(min_value=0, max_value=2 ** 32 - 1),
 )
-@settings(max_examples=20)
+@settings(max_examples=20, deadline=None)
 def test_complexity(n_samples, n_features, n_informative, random_state):
     """Behaviour test for complexity.
 
@@ -28,31 +30,39 @@ def test_complexity(n_samples, n_features, n_informative, random_state):
     """
     assume(n_informative < n_features)
 
-    x, y = make_regression(
-        n_samples, n_features, n_informative, 1, 0, noise=0.1, random_state=random_state
-    )
-    y = y.reshape(-1, 1)
+    # Average complexity over multiple datasets
+    n_datasets = 5
+    complexities = [0] * 7
 
-    opt_kwargs = dict(fit_intercept=True, normalize=False)
-    optimizers = [
-        SR3(thresholder="l0", threshold=0.1, **opt_kwargs),
-        SR3(thresholder="l1", threshold=0.1, **opt_kwargs),
-        Lasso(**opt_kwargs),
-        STLSQ(**opt_kwargs),
-        ElasticNet(**opt_kwargs),
-        Ridge(**opt_kwargs),
-        LinearRegression(**opt_kwargs),
-    ]
+    seed(random_state)
+    for rs in randint(low=0, high=2 ** 32 - 1, size=n_datasets):
 
-    optimizers = [SINDyOptimizer(o, unbias=True) for o in optimizers]
+        x, y = make_regression(
+            n_samples, n_features, n_informative, 1, 0, noise=0.1, random_state=rs
+        )
+        y = y.reshape(-1, 1)
 
-    for opt in optimizers:
-        opt.fit(x, y)
+        opt_kwargs = dict(fit_intercept=True, normalize=False)
+        optimizers = [
+            SR3(thresholder="l0", threshold=0.1, **opt_kwargs),
+            SR3(thresholder="l1", threshold=0.1, **opt_kwargs),
+            Lasso(**opt_kwargs),
+            STLSQ(**opt_kwargs),
+            ElasticNet(**opt_kwargs),
+            Ridge(**opt_kwargs),
+            LinearRegression(**opt_kwargs),
+        ]
 
-    for less_complex, more_complex in zip(optimizers, optimizers[1:]):
+        optimizers = [SINDyOptimizer(o, unbias=True) for o in optimizers]
+
+        for k, opt in enumerate(optimizers):
+            opt.fit(x, y)
+            complexities[k] += opt.complexity
+
+    for less_complex, more_complex in zip(complexities, complexities[1:]):
         # relax the condition to account for
         # noise and non-normalized threshold parameters
-        assert less_complex.complexity <= more_complex.complexity + 2
+        assert less_complex <= more_complex + 2
 
 
 @pytest.mark.parametrize(