Repository: OpenSourceEconomics/estimagic

Issue:
Add a robust cholesky decomposition to the utilities
## Background

The standard cholesky decompositions in `numpy.linalg` and `scipy.linalg` only work for positive-definite covariance matrices. However, there are two places where we might want to take a cholesky factor of a covariance matrix that is only positive semi-definite:

- During the reparametrization step (to_internal), if the user provided I covariance matrix that is only semi-definite in the start values
- In user code, when we need a cholesky factor of of a covariance matrix for random number generation and the optimizer tried a semi-definite covariance matrix. We could directly use it in skillmodels and respy and then probably don't need the bounds_distance anymore. 

## Proposed Solution

Add a cholesky function based on an [LDL](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.ldl.html#scipy.linalg.ldl) decomposition. See also [here](https://stackoverflow.com/questions/46100441/cholesky-ldl-decomposition-for-semidefinite-matrices-in-python)

## How much work is it

If we only want to accept square matrices as inputs, we'll just need a few lines. If we want the functions to work with higher dimensional arrays (i.e. shape [..., k, k]) we should wait until [numba-scipy](https://github.com/numba/numba-scipy) supports `scipy.linalg.ldl`





Discussion:

Patch:
diff --git a/estimagic/optimization/reparametrize.py b/estimagic/optimization/reparametrize.py
index db6eded..49b3360 100644
--- a/estimagic/optimization/reparametrize.py
+++ b/estimagic/optimization/reparametrize.py
@@ -8,6 +8,7 @@ from estimagic.optimization.utilities import cov_matrix_to_params
 from estimagic.optimization.utilities import cov_matrix_to_sdcorr_params
 from estimagic.optimization.utilities import cov_params_to_matrix
 from estimagic.optimization.utilities import number_of_triangular_elements_to_dimension
+from estimagic.optimization.utilities import robust_cholesky
 from estimagic.optimization.utilities import sdcorr_params_to_matrix
 
 
@@ -164,7 +165,7 @@ def _covariance_to_internal(params_subset, case, type_, bounds_distance):
         res["lower"] = np.maximum(res["lower"], np.zeros(len(res)))
         assert (res["upper"] >= res["lower"]).all(), "Invalid upper bound for variance."
     elif case == "free":
-        chol = np.linalg.cholesky(cov)
+        chol = robust_cholesky(cov, threshold=1e-8)
         chol_coeffs = chol[np.tril_indices(dim)]
         res["value"] = chol_coeffs
 
diff --git a/estimagic/optimization/utilities.py b/estimagic/optimization/utilities.py
index a55ba59..fde7f76 100644
--- a/estimagic/optimization/utilities.py
+++ b/estimagic/optimization/utilities.py
@@ -1,5 +1,17 @@
+import warnings
+from collections import namedtuple
+
 import numpy as np
 from fuzzywuzzy import process as fw_process
+from scipy.linalg import ldl
+from scipy.linalg import qr
+
+
+def chol_params_to_lower_triangular_matrix(params):
+    dim = number_of_triangular_elements_to_dimension(len(params))
+    mat = np.zeros((dim, dim))
+    mat[np.tril_indices(dim)] = params
+    return mat
 
 
 def cov_params_to_matrix(cov_params):
@@ -13,18 +25,34 @@ def cov_params_to_matrix(cov_params):
         cov (np.array): a covariance matrix
 
     """
-    dim = number_of_triangular_elements_to_dimension(len(cov_params))
-    lower = np.zeros((dim, dim))
-    lower[np.tril_indices(dim)] = cov_params
-    upper = lower.T.copy()
-    upper[np.diag_indices(dim)] = 0
-    cov = lower + upper
+    lower = chol_params_to_lower_triangular_matrix(cov_params)
+    cov = lower + np.tril(lower, k=-1).T
     return cov
 
 
 def cov_matrix_to_params(cov):
-    dim = len(cov)
-    return cov[np.tril_indices(dim)]
+    return cov[np.tril_indices(len(cov))]
+
+
+def sdcorr_params_to_sds_and_corr(sdcorr_params):
+    dim = number_of_triangular_elements_to_dimension(len(sdcorr_params))
+    sds = np.array(sdcorr_params[:dim])
+    corr = np.eye(dim)
+    corr[np.tril_indices(dim, k=-1)] = sdcorr_params[dim:]
+    corr += np.tril(corr, k=-1).T
+    return sds, corr
+
+
+def sds_and_corr_to_cov(sds, corr):
+    diag = np.diag(sds)
+    return diag @ corr @ diag
+
+
+def cov_to_sds_and_corr(cov):
+    sds = np.sqrt(np.diagonal(cov))
+    diag = np.diag(1 / sds)
+    corr = diag @ cov @ diag
+    return sds, corr
 
 
 def sdcorr_params_to_matrix(sdcorr_params):
@@ -40,20 +68,12 @@ def sdcorr_params_to_matrix(sdcorr_params):
         cov (np.array): a covariance matrix
 
     """
-    dim = number_of_triangular_elements_to_dimension(len(sdcorr_params))
-    diag = np.diag(sdcorr_params[:dim])
-    lower = np.zeros((dim, dim))
-    lower[np.tril_indices(dim, k=-1)] = sdcorr_params[dim:]
-    corr = np.eye(dim) + lower + lower.T
-    cov = diag.dot(corr).dot(diag)
-    return cov
+    return sds_and_corr_to_cov(*sdcorr_params_to_sds_and_corr(sdcorr_params))
 
 
 def cov_matrix_to_sdcorr_params(cov):
     dim = len(cov)
-    sds = np.sqrt(np.diagonal(cov))
-    scaling_matrix = np.diag(1 / sds)
-    corr = scaling_matrix.dot(cov).dot(scaling_matrix)
+    sds, corr = cov_to_sds_and_corr(cov)
     correlations = corr[np.tril_indices(dim, k=-1)]
     return np.hstack([sds, correlations])
 
@@ -120,3 +140,123 @@ def propose_algorithms(requested_algo, algos, number=3):
     proposals = [proposal[0] for proposal in proposals_w_probs]
 
     return proposals
+
+
+def robust_cholesky(matrix, threshold=None, warn=True, return_info=False):
+    """Lower triangular cholesky factor of *matrix*.
+
+    Args:
+        matrix (np.array): Square, symmetric and (almost) positive semi-definite matrix
+        threshold (float): Small negative number. Diagonal elements of D from the LDL
+            decomposition between threshold and zero are set to zero. Default is
+            minus machine accuracy.
+        warn (bool): If True, the user is warned about negative elements in D
+        return_info (bool): If True, also return a dictionary with 'method' and
+            'diagonals'. Method can take the values 'np.linalg.cholesky',
+            'LDL cholesky' and 'LDL cholesky with QR decomposition'.
+            'diagonals' are the diagonal entries of D. They can be used to construct
+            penalty terms for non positive definite matrices. Note that diagonals will
+            also be returned if the standard cholesky decomposition did not fail. This
+            allows for smooth penalty terms that start before actually reaching invalid
+            matrices.
+
+    Returns:
+        chol (np.array): Cholesky factor of matrix
+        info (float, optional): see return_info.
+
+    Raises:
+        np.linalg.LinalgError if the diagonal entries of D from the LDL decomposition
+        are smaller than threshold.
+
+    In contrast to a regular cholesky decomposition, this function will also
+    work for matrices that are only positive semi-definite or even indefinite.
+
+    For speed and precision reasons we first try a regular cholesky decomposition.
+    If it fails we switch to more robust methods. For this we first take an LDL
+    decomposition and set diagonal elements in D that are negative but above
+    threshold to 0. Then we use the modified D to construct a candidate cholesky
+    factor. If the candidate is not lower triangular, we use a QR decomposition to
+    make it lower triangular.
+
+
+    """
+    sum_stabilized = 0.0
+
+    try:
+        chol = np.linalg.cholesky(matrix)
+        method = "np.linalg.cholesky"
+        diags = np.diagonal(chol) ** 2
+    except np.linalg.LinAlgError:
+        method = "LDL cholesky"
+        threshold = threshold if threshold is not None else -np.finfo(float).eps
+
+        lu, d, _ = ldl(matrix)
+
+        diags = np.diagonal(d).copy()
+
+        for i in range(len(diags)):
+            if diags[i] >= 0:
+                diags[i] = np.sqrt(diags[i])
+            elif diags[i] > threshold:
+                diags[i] = 0
+                sum_stabilized += diags[i]
+            else:
+                raise np.linalg.LinAlgError(
+                    "Diagonal entry below threshold in D from LDL decomposition."
+                )
+
+        if sum_stabilized < 0 and warn:
+            warnings.warn(
+                "Negative diagonal entry has been set to 0 in robust_cholesky.",
+                category=RuntimeWarning,
+            )
+
+        candidate = lu * diags.reshape(1, len(diags))
+
+        is_triangular = (candidate[np.triu_indices(len(matrix), k=1)] == 0).all()
+
+        if is_triangular:
+            chol = candidate
+            method = "LDL cholesky"
+        else:
+            _, r = qr(candidate.T)
+            chol = r.T
+            method = "LDL cholesky with QR decomposition"
+
+    info = {"method": method, "diagonals": diags}
+
+    res = (chol, info) if return_info else chol
+    return res
+
+
+def namedtuple_from_dict(field_dict):
+    """Filled namedtuple generated from a dictionary.
+
+    Example:
+        >>> namedtuple_from_dict({'a': 1, 'b': 2})
+        NamedTuple(a=1, b=2)
+
+    """
+    return namedtuple("NamedTuple", field_dict)(**field_dict)
+
+
+def namedtuple_from_kwargs(**kwargs):
+    """Filled namedtuple generated from keyword arguments.
+
+    Example:
+        >>> namedtuple_from_kwargs(a=1, b=2)
+        NamedTuple(a=1, b=2)
+
+    """
+    return namedtuple("NamedTuple", kwargs)(**kwargs)
+
+
+def namedtuple_from_iterables(field_names, field_entries):
+    """Filled namedtuple generated from field_names and field_entries.
+
+    Example:
+        >>> namedtuple_from_iterables(field_names=['a', 'b'], field_entries=[1, 2])
+        NamedTuple(a=1, b=2)
+
+    """
+    return namedtuple("NamedTuple", field_names)(*field_entries)
diff --git a/setup.py b/setup.py
index f5d8147..06d9410 100644
--- a/setup.py
+++ b/setup.py
@@ -1,7 +1,7 @@
 from setuptools import find_packages
 from setuptools import setup
 
-version = "0.0.13"
+version = "0.0.14"
 
 setup(
     name="estimagic",


Test patch:
diff --git a/estimagic/tests/optimization/test_utilities.py b/estimagic/tests/optimization/test_utilities.py
index e393cae..fce0d21 100644
--- a/estimagic/tests/optimization/test_utilities.py
+++ b/estimagic/tests/optimization/test_utilities.py
@@ -1,13 +1,26 @@
 import numpy as np
+import pandas as pd
+import pytest
 from numpy.testing import assert_array_almost_equal as aaae
 
+from estimagic.optimization.utilities import chol_params_to_lower_triangular_matrix
 from estimagic.optimization.utilities import cov_matrix_to_params
 from estimagic.optimization.utilities import cov_matrix_to_sdcorr_params
 from estimagic.optimization.utilities import cov_params_to_matrix
+from estimagic.optimization.utilities import cov_to_sds_and_corr
 from estimagic.optimization.utilities import dimension_to_number_of_triangular_elements
 from estimagic.optimization.utilities import index_element_to_string
 from estimagic.optimization.utilities import number_of_triangular_elements_to_dimension
+from estimagic.optimization.utilities import robust_cholesky
 from estimagic.optimization.utilities import sdcorr_params_to_matrix
+from estimagic.optimization.utilities import sdcorr_params_to_sds_and_corr
+from estimagic.optimization.utilities import sds_and_corr_to_cov
+
+
+def test_chol_params_to_lower_triangular_matrix():
+    calculated = chol_params_to_lower_triangular_matrix(pd.Series([1, 2, 3]))
+    expected = np.array([[1, 0], [2, 3]])
+    aaae(calculated, expected)
 
 
 def test_cov_params_to_matrix():
@@ -24,6 +37,15 @@ def test_cov_matrix_to_params():
     aaae(calculated, expected)
 
 
+def test_sdcorr_params_to_sds_and_corr():
+    sdcorr_params = pd.Series([1, 2, 3, 0.1, 0.2, 0.3])
+    exp_corr = np.array([[1, 0.1, 0.2], [0.1, 1, 0.3], [0.2, 0.3, 1]])
+    exp_sds = np.array([1, 2, 3])
+    calc_sds, calc_corr = sdcorr_params_to_sds_and_corr(sdcorr_params)
+    aaae(calc_sds, exp_sds)
+    aaae(calc_corr, exp_corr)
+
+
 def test_sdcorr_params_to_matrix():
     sds = np.sqrt([1, 2, 3])
     corrs = [0.07071068, 0.11547005, 0.08981462]
@@ -42,6 +64,25 @@ def test_cov_matrix_to_sdcorr_params():
     aaae(calculated, expected)
 
 
+def test_sds_and_corr_to_cov():
+    sds = [1, 2, 3]
+    corr = np.ones((3, 3)) * 0.2
+    corr[np.diag_indices(3)] = 1
+    calculated = sds_and_corr_to_cov(sds, corr)
+    expected = np.array([[1.0, 0.4, 0.6], [0.4, 4.0, 1.2], [0.6, 1.2, 9.0]])
+    aaae(calculated, expected)
+
+
+def test_cov_to_sds_and_corr():
+    cov = np.array([[1.0, 0.4, 0.6], [0.4, 4.0, 1.2], [0.6, 1.2, 9.0]])
+    calc_sds, calc_corr = cov_to_sds_and_corr(cov)
+    exp_sds = [1, 2, 3]
+    exp_corr = np.ones((3, 3)) * 0.2
+    exp_corr[np.diag_indices(3)] = 1
+    aaae(calc_sds, exp_sds)
+    aaae(calc_corr, exp_corr)
+
+
 def test_number_of_triangular_elements_to_dimension():
     inputs = [6, 10, 15, 21]
     expected = [3, 4, 5, 6]
@@ -61,3 +102,33 @@ def test_index_element_to_string():
     expected = ["a_b_1", "bla~5~6", "bla"]
     for inp, exp in zip(inputs, expected):
         assert index_element_to_string(*inp) == exp
+
+
+def random_cov(dim, seed):
+    num_elements = int(dim * (dim + 1) / 2)
+    chol = np.zeros((dim, dim))
+    chol[np.tril_indices(dim)] = np.random.uniform(size=num_elements)
+    cov = chol @ chol.T
+    zero_positions = np.random.choice(range(dim), size=int(dim / 5), replace=False)
+    for pos in zero_positions:
+        cov[:, pos] = 0
+        cov[pos] = 0
+    return cov
+
+
+seeds = [58822, 3181, 98855, 44002, 47631, 97741, 10655, 4600, 1151, 58189]
+dims = [8] * 6 + [10, 12, 15, 20]
+
+
+@pytest.mark.parametrize("dim, seed", zip(dims, seeds))
+def test_robust_cholesky_with_zero_variance(dim, seed):
+    cov = random_cov(dim, seed)
+    chol = robust_cholesky(cov)
+    aaae(chol.dot(chol.T), cov)
+    assert (chol[np.triu_indices(len(cov), k=1)] == 0).all()
+
+
+def test_robust_cholesky_with_extreme_cases():
+    for cov in [np.ones((5, 5)), np.zeros((5, 5))]:
+        chol = robust_cholesky(cov)
+        aaae(chol.dot(chol.T), cov)