diff --git a/src/fmu/tools/__init__.py b/src/fmu/tools/__init__.py index 01ebff36..4b544cca 100644 --- a/src/fmu/tools/__init__.py +++ b/src/fmu/tools/__init__.py @@ -42,6 +42,7 @@ from fmu.tools.sensitivities import calc_tornadoinput # noqa from fmu.tools.sensitivities import excel2dict_design # noqa from fmu.tools.sensitivities import summarize_design # noqa + from fmu.tools.sensitivities import nearcorr # noqa __all__.extend( [ @@ -50,6 +51,7 @@ "calc_tornadoinput", "excel2dict_design", "summarize_design", + "nearcorr", ] ) diff --git a/src/fmu/tools/sensitivities/__init__.py b/src/fmu/tools/sensitivities/__init__.py index 046b7fdd..6fa3c1c2 100755 --- a/src/fmu/tools/sensitivities/__init__.py +++ b/src/fmu/tools/sensitivities/__init__.py @@ -9,6 +9,7 @@ from fmu.tools.sensitivities._excel2dict import excel2dict_design, inputdict_to_yaml from fmu.tools.sensitivities._tornado_onebyone import calc_tornadoinput from fmu.tools.sensitivities.create_design import DesignMatrix +from fmu.tools.sensitivities.nearest_correlation import nearcorr __all__ = [ "summarize_design", @@ -16,4 +17,5 @@ "DesignMatrix", "excel2dict_design", "inputdict_to_yaml", + "nearcorr", ] diff --git a/src/fmu/tools/sensitivities/design_distributions.py b/src/fmu/tools/sensitivities/design_distributions.py index b9af4974..21aa4484 100644 --- a/src/fmu/tools/sensitivities/design_distributions.py +++ b/src/fmu/tools/sensitivities/design_distributions.py @@ -10,6 +10,8 @@ import pandas as pd import scipy.stats +from fmu.tools.sensitivities.nearest_correlation import nearcorr + def _check_dist_params_normal(dist_params): if len(dist_params) not in [2, 4]: @@ -563,7 +565,7 @@ def make_covariance_matrix(df_correlations, stddevs=None): # Project to nearest symmetric positive definite matrix if not _is_positive_definite(corr_matrix): input_corr_matrix = corr_matrix.copy() - corr_matrix = _nearest_positive_definite(corr_matrix) + corr_matrix = nearcorr(corr_matrix) print("Input correlation matrix: ") with np.printoptions(precision=3, suppress=True): print(input_corr_matrix) @@ -585,35 +587,6 @@ def make_covariance_matrix(df_correlations, stddevs=None): return np.dot(cov_matrix, diag) -def _nearest_positive_definite(a_mat): - """Implementation taken from: - https://stackoverflow.com/questions/43238173/ - python-convert-matrix-to-positive-semi-definite/43244194#43244194 - """ - - b_mat = (a_mat + a_mat.T) / 2 - _, s_mat, v_mat = la.svd(b_mat) - - h_mat = np.dot(v_mat.T, np.dot(np.diag(s_mat), v_mat)) - - a2_mat = (b_mat + h_mat) / 2 - - a3_mat = (a2_mat + a2_mat.T) / 2 - - if _is_positive_definite(a3_mat): - return a3_mat - - spacing = np.spacing(la.norm(a_mat)) - identity = np.eye(a_mat.shape[0]) - kiter = 1 - while not _is_positive_definite(a3_mat): - mineig = np.min(np.real(la.eigvals(a3_mat))) - a3_mat += identity * (-mineig * kiter**2 + spacing) - kiter += 1 - - return a3_mat - - def _is_positive_definite(b_mat): """Returns true when input is positive-definite, via Cholesky""" try: diff --git a/src/fmu/tools/sensitivities/nearest_correlation.py b/src/fmu/tools/sensitivities/nearest_correlation.py new file mode 100644 index 00000000..c2d46f56 --- /dev/null +++ b/src/fmu/tools/sensitivities/nearest_correlation.py @@ -0,0 +1,79 @@ +import numpy as np +from numpy import copy, inf +from numpy.linalg import norm + + +def nearcorr( + A, + tol=None, + max_iterations=100, + weights=None, +): + """ + Finds the nearest correlation matrix to the symmetric matrix A. + Based on https://github.com/mikecroucher/nearest_correlation + + Parameters + ---------- + A: symmetric numpy array + tol: convergence tolerance + max_iterations: maximum number of iterations (default 100) + weights: optional vector defining a diagonal weight matrix diag(W) + + References + ---------- + + https://nhigham.com/2013/02/13/the-nearest-correlation-matrix/ + https://eprints.maths.manchester.ac.uk/232/1/paper3.pdf + + """ + if not np.array_equal(A, A.T): + raise ValueError("A must be symmetric") + + if not tol: + eps = np.spacing(1) + tol = eps * np.shape(A)[0] * np.array([1, 1]) + + if weights is None: + weights = np.ones(np.shape(A)[0]) + + X = copy(A) + Y = copy(A) + ds = np.zeros(np.shape(A)) + rel_diffY = inf + rel_diffX = inf + rel_diffXY = inf + + Whalf = np.sqrt(np.outer(weights, weights)) + + iteration = 0 + while max(rel_diffX, rel_diffY, rel_diffXY) > tol[0]: + iteration += 1 + if iteration > max_iterations: + raise ValueError(f"No convergence after {max_iterations} iterations") + + Xold = copy(X) + R = X - ds + R_wtd = Whalf * R + X = proj_spd(R_wtd) + + X = X / Whalf + ds = X - R + Yold = copy(Y) + Y = copy(X) + np.fill_diagonal(Y, 1) + normY = norm(Y, "fro") + rel_diffX = norm(X - Xold, "fro") / norm(X, "fro") + rel_diffY = norm(Y - Yold, "fro") / normY + rel_diffXY = norm(Y - X, "fro") / normY + + X = copy(Y) + + return X + + +def proj_spd(A): + # NOTE: the input matrix is assumed to be symmetric + d, v = np.linalg.eigh(A) + A = (v * np.maximum(d, 0)).dot(v.T) + return (A + A.T) / 2 diff --git a/tests/sensitivities/test_nearest_correlation.py b/tests/sensitivities/test_nearest_correlation.py new file mode 100644 index 00000000..7c03b44f --- /dev/null +++ b/tests/sensitivities/test_nearest_correlation.py @@ -0,0 +1,80 @@ +import numpy as np +import pytest + +from fmu.tools import nearcorr + +# References +# [1] 'Computing the nearest correlation matrix - +# a problem from finance': Higham, IMA Journal of Numerical Analysis (2002) 22, 329.343 + + +def test_nag_example(): + """Test from NAG Mark 24 documentation for g02aa, originally from [1]""" + A = np.array([[2, -1, 0, 0], [-1, 2, -1, 0], [0, -1, 2, -1], [0, 0, -1, 2]]) + + X = nearcorr(A) + + expected_result = np.array( + [ + [1.0, -0.8084125, 0.1915875, 0.10677505], + [-0.8084125, 1.0, -0.65623269, 0.1915875], + [0.1915875, -0.65623269, 1.0, -0.8084125], + [0.10677505, 0.1915875, -0.8084125, 1.0], + ] + ) + + assert (np.abs((X - expected_result)) < 1e-8).all() + + +def test_higham_example_2002(): + """Example taken from [1]""" + A = np.array([[1, 1, 0], [1, 1, 1], [0, 1, 1]]) + + X = nearcorr(A) + + expected_result = np.array( + [ + [1.0, 0.76068985, 0.15729811], + [0.76068985, 1.0, 0.76068985], + [0.15729811, 0.76068985, 1.0], + ] + ) + + assert (np.abs((X - expected_result)) < 1e-8).all() + + +def test_assert_symmetric(): + """Test that non-symmetric matrix raises ValueError""" + A = np.array([[1, 1, 0], [1, 1, 1], [1, 1, 1]]) + + with pytest.raises(ValueError): + nearcorr(A) + + +def test_exceeded_max_iterations(): + """Test that exceeding max iterations raises ValueError""" + A = np.array([[1, 1, 0], [1, 1, 1], [0, 1, 1]]) + with pytest.raises(ValueError, match="No convergence after 10 iterations"): + nearcorr(A, max_iterations=10) + + +def test_weights_documented_example(): + """Test weights using the exact example from examples: + https://github.com/mikecroucher/nearest_correlation""" + A = np.array([[1, 1, 0], [1, 1, 1], [0, 1, 1]]) + weights = np.array([1, 2, 3]) + X = nearcorr(A, weights=weights) + + expected_result = np.array( + [ + [1.0, 0.66774961, 0.16723692], + [0.66774961, 1.0, 0.84557496], + [0.16723692, 0.84557496, 1.0], + ] + ) + + assert np.allclose(X, expected_result, atol=1e-8), ( + f"Result doesn't match documented example.\n" + f"Expected:\n{expected_result}\n" + f"Got:\n{X}" + )