From b11775763171473d190052c044b654fc2d0033c2 Mon Sep 17 00:00:00 2001 From: Blunde1 Date: Mon, 16 Oct 2023 14:49:52 +0200 Subject: [PATCH] Compute cross-correlation matrices without matrix inversion --- src/ert/analysis/_es_update.py | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/ert/analysis/_es_update.py b/src/ert/analysis/_es_update.py index ef9cc282414..b9c7ad5e1c1 100644 --- 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,