Compute cross-correlation matrices without matrix inversion

equinor · Oct 16, 2023 · b117757 · b117757
1 parent 20cfcde
commit b117757
Showing 1 changed file with 8 additions and 10 deletions.
diff --git a/src/ert/analysis/_es_update.py b/src/ert/analysis/_es_update.py
@@ -457,9 +457,6 @@ def analysis_ES(
         noise = rng.standard_normal(size=(num_obs, ensemble_size))
 
         if module.localization():
-            Y_prime = S - S.mean(axis=1, keepdims=True)
-            C_YY = Y_prime @ Y_prime.T / (ensemble_size - 1)
-            Sigma_Y = np.diag(np.sqrt(np.diag(C_YY)))
             batch_size: int = 1000
             correlation_threshold = module.localization_correlation_threshold(
                 ensemble_size
@@ -476,17 +473,18 @@ def analysis_ES(
                 batches = _split_by_batchsize(np.arange(0, num_params), batch_size)
                 for param_batch_idx in tqdm(batches):
                     X_local = temp_storage[parameter.name][param_batch_idx, :]
-                    A = X_local - X_local.mean(axis=1, keepdims=True)
-                    C_AA = A @ A.T / (ensemble_size - 1)
+
+                    # Standard deviations
+                    Sigma_Y = np.std(S, axis=1, ddof=1)
+                    Sigma_A = np.std(X_local, axis=1, ddof=1)
 
                     # State-measurement covariance matrix
+                    Y_prime = S - S.mean(axis=1, keepdims=True)
+                    A = X_local - X_local.mean(axis=1, keepdims=True)
                     C_AY = A @ Y_prime.T / (ensemble_size - 1)
-                    Sigma_A = np.diag(np.sqrt(np.diag(C_AA)))
 
-                    # State-measurement correlation matrix
-                    c_AY = np.abs(
-                        np.linalg.inv(Sigma_A) @ C_AY @ np.linalg.inv(Sigma_Y)
-                    )
+                    # Absolute values of the correlation matrix
+                    c_AY = np.abs((C_AY / Sigma_Y.reshape(1, -1)) / Sigma_A.reshape(-1, 1))
                     c_bool = c_AY > correlation_threshold
                     # Some parameters might be significantly correlated
                     # to the exact same responses,