diff --git a/src/iterative_ensemble_smoother/experimental.py b/src/iterative_ensemble_smoother/experimental.py
index e5764a12..5a7a5337 100644
--- a/src/iterative_ensemble_smoother/experimental.py
+++ b/src/iterative_ensemble_smoother/experimental.py
@@ -3,6 +3,7 @@
 features of iterative_ensemble_smoother
 """
 import numbers
+import warnings
 from typing import List, Optional, Tuple, Union
 
 import numpy as np
@@ -149,9 +150,16 @@ def _correlation_matrix(self, cov_XY, X, Y):
         # Compute the correlation matrix from the covariance matrix
         corr_XY = (cov_XY / stds_X[:, np.newaxis]) / stds_Y[np.newaxis, :]
 
-        # Perform checks
-        assert corr_XY.max() <= 1
-        assert corr_XY.min() >= -1
+        # Perform checks. There appears to be occasional numerical issues in
+        # the equation. With 2 ensemble members, we get e.g. a max value of
+        # 1.0000000000016778. We allow some leeway and clip the results.
+        eps = 1e-8
+        if not ((corr_XY.max() <= 1 + eps) and (corr_XY.min() >= -1 - eps)):
+            msg = "Cross-correlation matrix has entries not in [-1, 1]."
+            msg += f"The min and max values are: {corr_XY.min()} and {corr_XY.max()}"
+            warnings.warn(msg)
+
+        corr_XY = np.clip(corr_XY, a_min=-1, a_max=1)
         return corr_XY
 
     def assimilate(self, X, Y, D, alpha, correlation_threshold=None, verbose=False):