Add support for positive semidefinite matrices to Iman-Conover

Fallback to SVD when Cholesky fails for semidefinite matrices, typically when correlation matrix contains perfect correlations

Author: dafeda

src/semeio/fmudesign/iman_conover.py

@@ -1,5 +1,5 @@
 """
-An implementation of the Iman-Conover transformation.
+An implementation of the Iman-Conover transformation using hybrid Cholesky/SVD.
 
 Using Iman-Conover with Latin Hybercube sampling
 ------------------------------------------------
@@ -41,16 +41,71 @@ def _is_positive_definite(X: npt.NDArray[Any]) -> bool:
         return False
 
 
+def _is_positive_semidefinite(X: npt.NDArray[np.float64]) -> bool:
+    """Check if matrix is positive semidefinite using eigenvalue decomposition."""
+    try:
+        # PSD matrices must be square and symmetric
+        if X.shape[0] != X.shape[1]:
+            return False
+        if not np.allclose(X, X.T, rtol=1e-12, atol=1e-12):
+            return False
+
+        eigenvals = np.linalg.eigvals(X)
+        tol = 1e-12
+        # Even though the eigenvalues of symmetric real matrices are guaranteed
+        # to be real, they might turn out complex due to numerics.
+        # Hence, the usage of np.real.
+        return bool(np.all(np.real(eigenvals) >= -tol))
+    except np.linalg.LinAlgError:
+        return False
+
+
+def _matrix_sqrt_with_pinv(
+    X: npt.NDArray[np.float64],
+) -> tuple[npt.NDArray[np.floating[Any]], npt.NDArray[np.floating[Any]]]:
+    """Compute matrix square root and its pseudoinverse using Cholesky if possible, otherwise SVD."""
+    if _is_positive_definite(X):
+        # Use Cholesky for positive definite matrices (faster and more stable)
+        L = np.linalg.cholesky(X)
+        # For Cholesky decomposition, inv(L) can be computed efficiently
+        L_inv = sp.linalg.solve_triangular(L, np.eye(L.shape[0]), lower=True)
+        return L, L_inv.T
+    else:
+        # Fall back to SVD for positive semidefinite matrices
+        U, s, _ = np.linalg.svd(X, hermitian=True)
+
+        # Threshold small eigenvalues for numerical stability
+        s_clipped = np.maximum(s, 0)
+        s_sqrt = np.sqrt(s_clipped)
+
+        # Keep all dimensions but zero out small eigenvalues
+        tol = 1e-12
+        mask = s_clipped >= tol
+        s_sqrt[~mask] = 0
+
+        # Square root matrix
+        L = U * s_sqrt
+
+        # Pseudoinverse of square root
+        s_inv = np.zeros_like(s_sqrt)
+        s_inv[mask] = 1.0 / s_sqrt[mask]
+        L_pinv = (U * s_inv).T
+
+        return L, L_pinv
+
+
 class ImanConover:
     def __init__(self, correlation_matrix: npt.NDArray[Any]) -> None:
         """Create an Iman-Conover transform.
 
         Parameters
         ----------
         correlation_matrix : ndarray
-            Target correlation matrix of shape (K, K). The Iman-Conover will
-            try to induce a correlation on the data set X so that corr(X) is
-            as close to `correlation_matrix` as possible.
+            Target correlation matrix of shape (K, K).
+            The Iman-Conover will try to induce a correlation on
+            the data set X so that corr(X) is as close
+            to `correlation_matrix` as possible.
+            Can be positive definite or positive semidefinite.
 
         Notes
         -----
@@ -127,11 +182,13 @@ def __init__(self, correlation_matrix: npt.NDArray[Any]) -> None:
             raise ValueError("Correlation matrix must have 1.0 on diagonal.")
         if not np.allclose(correlation_matrix.T, correlation_matrix):
             raise ValueError("Correlation matrix must be symmetric.")
-        if not _is_positive_definite(correlation_matrix):
-            raise ValueError("Correlation matrix must be positive definite.")
+        if not _is_positive_semidefinite(correlation_matrix):
+            raise ValueError("Correlation matrix must be positive semidefinite.")
 
         self.C = correlation_matrix.copy()
-        self.P = np.linalg.cholesky(self.C)
+
+        # Use hybrid approach: Cholesky if possible, SVD otherwise
+        self.P, _ = _matrix_sqrt_with_pinv(self.C)
 
     def __call__(self, X: npt.NDArray[Any]) -> npt.NDArray[Any]:
         """Transform an input matrix X.
@@ -156,13 +213,13 @@ def __call__(self, X: npt.NDArray[Any]) -> npt.NDArray[Any]:
         if not isinstance(X, np.ndarray):
             raise TypeError("Input argument `X` must be NumPy array.")
         if not X.ndim == 2:
-            raise ValueError("Correlation matrix must be square.")
+            raise ValueError("Input matrix must be 2D.")
 
         N, K = X.shape
 
         if self.P.shape[0] != K:
             msg = f"Shape of `X` ({X.shape}) does not match shape of "
-            msg += f"correlation matrix ({self.P.shape})"
+            msg += f"correlation matrix ({self.C.shape})"
             raise ValueError(msg)
 
         if N <= K:
@@ -173,23 +230,38 @@ def __call__(self, X: npt.NDArray[Any]) -> npt.NDArray[Any]:
         # approximately multivariate normal (but with correlations).
         # The new data has the same rank correlation as the original data.
         ranks = sp.stats.rankdata(X, axis=0) / (N + 1)
-        normal_scores = sp.stats.norm.ppf(ranks)  # + np.random.randn(N, K) * epsilon
+        normal_scores = sp.stats.norm.ppf(ranks)
 
         # STEP TWO - Remove correlations from the transformed data
         empirical_correlation = np.corrcoef(normal_scores, rowvar=False)
-        if not _is_positive_definite(empirical_correlation):
-            msg = "Rank data correlation not positive definite."
-            msg += "There are perfect correlations in the ranked data."
-            msg += "Supply more data (rows in X) or sample differently."
+        if not _is_positive_semidefinite(empirical_correlation):
+            msg = "Rank data correlation not positive semidefinite."
             raise ValueError(msg)
 
-        decorrelation_matrix = np.linalg.cholesky(empirical_correlation)
+        # Fail if input has perfect rank correlations that differ from target
+        # (impossible to change perfect correlations via permutation)
+        off_diagonal = empirical_correlation[~np.eye(K, dtype=bool)]
+        if np.any(
+            np.isclose(np.abs(off_diagonal), 1.0, rtol=1e-10)
+        ) and not np.allclose(empirical_correlation, self.C, rtol=1e-6, atol=1e-6):
+            msg = "Input data has perfect rank correlations that conflict with target correlation structure."
+            raise ValueError(msg)
 
-        # We exploit the fact that Q is lower-triangular and avoid the inverse.
-        # X = N @ inv(Q)^T  =>  X @ Q^T = N  =>  (Q @ X^T)^T = N
-        decorrelated_scores = sp.linalg.solve_triangular(
-            decorrelation_matrix, normal_scores.T, lower=True
-        ).T
+        # Use same hybrid approach for decorrelation
+        decorrelation_matrix, decorrelation_pinv = _matrix_sqrt_with_pinv(
+            empirical_correlation
+        )
+
+        if _is_positive_definite(empirical_correlation):
+            # Use efficient triangular solve for Cholesky case
+            # We exploit the fact that Q is lower-triangular and avoid the inverse.
+            # X = N @ inv(Q)^T  =>  X @ Q^T = N  =>  (Q @ X^T)^T = N
+            decorrelated_scores = sp.linalg.solve_triangular(
+                decorrelation_matrix, normal_scores.T, lower=True
+            ).T
+        else:
+            # Use pseudoinverse for SVD case
+            decorrelated_scores = normal_scores @ decorrelation_pinv
 
         # STEP THREE - Induce correlations in transformed space
         correlated_scores = decorrelated_scores @ self.P.T

tests/fmudesign/test_iman_conover.py

@@ -3,7 +3,11 @@
 import scipy as sp
 from scipy.stats import spearmanr
 
-from semeio.fmudesign.iman_conover import ImanConover, decorrelate
+from semeio.fmudesign.iman_conover import (
+    ImanConover,
+    _matrix_sqrt_with_pinv,
+    decorrelate,
+)
 
 
 @pytest.fixture
@@ -206,9 +210,102 @@ def test_dataset_with_unity_correlation_in_ranks(self):
         desired_corr = np.identity(2)
 
         transform = ImanConover(desired_corr)
-        with pytest.raises(ValueError):
+        with pytest.raises(
+            ValueError,
+            match="Input data has perfect rank correlations that conflict with target correlation structure.",
+        ):
             transform(X)
 
+    def test_perfect_correlation_target(self):
+        """Test that Iman-Conover can transform data to have perfect correlation."""
+        # Use Latin Hypercube sampling for well-distributed, uncorrelated data
+        sampler = sp.stats.qmc.LatinHypercube(d=2, seed=42, scramble=True)
+        X = sampler.random(n=1000)
+
+        # Target: perfect positive correlation
+        target_corr = np.array([[1, 1], [1, 1]])
+
+        ic = ImanConover(target_corr)
+        X_transformed = ic(X)
+
+        # Calculate achieved Spearman correlation
+        spearman_corr, _ = sp.stats.spearmanr(X_transformed[:, 0], X_transformed[:, 1])
+
+        # Should achieve perfect correlation
+        assert np.isclose(spearman_corr, 1.0, atol=0.01), (
+            f"Expected perfect correlation, got {spearman_corr}"
+        )
+
+        # Marginal distributions should be preserved
+        for k in range(X.shape[1]):
+            assert np.allclose(np.sort(X[:, k]), np.sort(X_transformed[:, k]))
+
+        # With perfect correlation, the ranks should be identical
+        ranks_col1 = sp.stats.rankdata(X_transformed[:, 0])
+        ranks_col2 = sp.stats.rankdata(X_transformed[:, 1])
+        assert np.allclose(ranks_col1, ranks_col2)
+
+    def test_that_the_matrix_square_root_has_the_expected_properties_with_perfect_correlations(
+        self,
+    ):
+        X = np.array(
+            [
+                [1.0, 0.0, 0.0],
+                [0.0, 1.0, 1.0],  # Perfect correlation between vars 2 & 3
+                [0.0, 1.0, 1.0],
+            ]
+        )
+        P, P_pinv = _matrix_sqrt_with_pinv(X)
+
+        # Test reconstruction: P @ P.T should equal X
+        reconstructed = P @ P.T
+        reconstruction_error = np.linalg.norm(reconstructed - X, "fro")
+        assert reconstruction_error < 1e-10
+
+        # Test pseudoinverse property: P @ P_pinv should be the projection onto
+        # the column space of P (for rank-deficient matrices)
+        projection = P @ P_pinv
+        projection_error = np.linalg.norm(projection @ P - P, "fro")
+        assert projection_error < 1e-10
+
+    def test_rank_deficient_correlation_matrix(self):
+        """Test handling of rank-deficient correlation matrices."""
+        # Create a 3x3 rank-deficient matrix where X3 = X1 (perfect correlation)
+        target_corr = np.array([[1.0, 0.3, 1.0], [0.3, 1.0, 0.3], [1.0, 0.3, 1.0]])
+
+        # Verify it's actually rank-deficient
+        n_vars = target_corr.shape[0]
+        n_dependencies = 1  # X3 depends on X1
+        expected_rank = n_vars - n_dependencies
+        assert np.linalg.matrix_rank(target_corr) == expected_rank
+
+        sampler = sp.stats.qmc.LatinHypercube(d=3, seed=42, scramble=True)
+        X = sampler.random(n=1000)
+
+        ic = ImanConover(target_corr)
+        X_transformed = ic(X)
+        achieved_corr = sp.stats.spearmanr(X_transformed)[0]
+
+        # Should achieve near-perfect correlation between X1 and X3
+        assert np.isclose(achieved_corr[0, 2], 1.0, atol=0.01), (
+            f"Expected near-perfect correlation, got {achieved_corr[0, 2]:.6f}"
+        )
+
+        # Other correlations should be close to target
+        assert np.isclose(achieved_corr[0, 1], 0.3, atol=0.05), (
+            f"Expected ~0.3 correlation, got {achieved_corr[0, 1]:.6f}"
+        )
+
+        ranks_col1 = sp.stats.rankdata(X_transformed[:, 0])
+        ranks_col3 = sp.stats.rankdata(X_transformed[:, 2])
+        assert np.allclose(ranks_col1, ranks_col3), (
+            "Ranks should be identical for perfectly correlated columns"
+        )
+
+        # Marginal distributions should be preserved
+        for k in range(X.shape[1]):
+            assert np.allclose(np.sort(X[:, k]), np.sort(X_transformed[:, k]))
+
 
 if __name__ == "__main__":
     pytest.main(args=[__file__, "--doctest-modules", "-v", "-l"])