Nearest correlation

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",
         ]
     )
 

src/fmu/tools/sensitivities/__init__.py

@@ -9,11 +9,13 @@
 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",
     "calc_tornadoinput",
     "DesignMatrix",
     "excel2dict_design",
     "inputdict_to_yaml",
+    "nearcorr",
 ]

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:

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

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}"
+    )