Merge pull request statsmodels#5198 from bashtage/fix-corr-nearest-sparse

BUG/TST/CLN: Fix corr nearest factor

Author: josef-pkt

statsmodels/stats/correlation_tools.py

@@ -154,8 +154,7 @@ def corr_clipped(corr, threshold=1e-15):
 
 def cov_nearest(cov, method='clipped', threshold=1e-15, n_fact=100,
                 return_all=False):
-
-    '''
+    """
     Find the nearest covariance matrix that is postive (semi-) definite
 
     This leaves the diagonal, i.e. the variance, unchanged
@@ -169,7 +168,7 @@ def cov_nearest(cov, method='clipped', threshold=1e-15, n_fact=100,
         used.if "nearest", then ``corr_nearest`` is used
     threshold : float
         clipping threshold for smallest eigen value, see Notes
-    nfact : int or float
+    n_fact : int or float
         factor to determine the maximum number of iterations in
         ``corr_nearest``. See its doc string
     return_all : bool
@@ -205,14 +204,13 @@ def cov_nearest(cov, method='clipped', threshold=1e-15, n_fact=100,
     --------
     corr_nearest
     corr_clipped
-
-    '''
+    """
 
     from statsmodels.stats.moment_helpers import cov2corr, corr2cov
     cov_, std_ = cov2corr(cov, return_std=True)
     if method == 'clipped':
         corr_ = corr_clipped(cov_, threshold=threshold)
-    elif method == 'nearest':
+    else:  # method == 'nearest'
         corr_ = corr_nearest(cov_, threshold=threshold, n_fact=n_fact)
 
     cov_ = corr2cov(corr_, std_)
@@ -604,10 +602,10 @@ def corr_nearest_factor(corr, rank, ctol=1e-6, lam_min=1e-30,
 
     References
     ----------
-    R Borsdof, N Higham, M Raydan (2010).  Computing a nearest
-    correlation matrix with factor structure. SIAM J Matrix Anal
-    Appl, 31:5, 2603-2622.
-    http://eprints.ma.man.ac.uk/1523/01/covered/MIMS_ep2009_87.pdf
+    .. [*] R Borsdof, N Higham, M Raydan (2010).  Computing a nearest
+       correlation matrix with factor structure. SIAM J Matrix Anal Appl,
+       31:5, 2603-2622.
+       http://eprints.ma.man.ac.uk/1523/01/covered/MIMS_ep2009_87.pdf
 
     Examples
     --------
@@ -681,7 +679,8 @@ def func(X):
                 ir += bs
             return fval
 
-    rslt = _spg_optim(func, grad, X, _project_correlation_factors)
+    rslt = _spg_optim(func, grad, X, _project_correlation_factors, ctol=ctol,
+                      lam_min=lam_min, lam_max=lam_max, maxiter=maxiter)
     root = rslt.params
     diag = 1 - (root**2).sum(1)
     soln = FactoredPSDMatrix(diag, root)

statsmodels/stats/tests/test_corrpsd.py

@@ -87,6 +87,7 @@ class Holder(object):
      0.7445177382760834, 0.7766852947049376, 0.8107916469252828, 1.01
     ]).reshape(6,6, order='F')
 
+
 def test_corr_psd():
     # test positive definite matrix is unchanged
     x = np.array([[1, -0.2, -0.9], [-0.2, 1, -0.2], [-0.9, -0.2, 1]])
@@ -122,7 +123,6 @@ def test_nearest(self):
         #print evals[0] / 1e-7 - 1
         assert_allclose(evals[0], 1e-7, rtol=1e-6)
 
-
     def test_clipped(self):
         x = self.x
         res_r = self.res
@@ -190,37 +190,30 @@ def setup_class(cls):
         cls.x = x
         cls.res = cov1_r
 
-def test_corrpsd_threshold():
+
+@pytest.mark.parametrize('threshold', [0, 1e-15, 1e-10, 1e-6])
+def test_corrpsd_threshold(threshold):
     x = np.array([[1, -0.9, -0.9], [-0.9, 1, -0.9], [-0.9, -0.9, 1]])
 
-    #print np.linalg.eigvalsh(x)
-    for threshold in [0, 1e-15, 1e-10, 1e-6]:
+    y = corr_nearest(x, n_fact=100, threshold=threshold)
+    evals = np.linalg.eigvalsh(y)
+    assert_allclose(evals[0], threshold, rtol=1e-6, atol=1e-15)
 
-        y = corr_nearest(x, n_fact=100, threshold=threshold)
-        evals = np.linalg.eigvalsh(y)
-        #print 'evals', evals, threshold
-        assert_allclose(evals[0], threshold, rtol=1e-6, atol=1e-15)
+    y = corr_clipped(x, threshold=threshold)
+    evals = np.linalg.eigvalsh(y)
+    assert_allclose(evals[0], threshold, rtol=0.25, atol=1e-15)
 
-        y = corr_clipped(x, threshold=threshold)
-        evals = np.linalg.eigvalsh(y)
-        #print 'evals', evals, threshold
-        assert_allclose(evals[0], threshold, rtol=0.25, atol=1e-15)
+    y = cov_nearest(x, method='nearest', n_fact=100, threshold=threshold)
+    evals = np.linalg.eigvalsh(y)
+    assert_allclose(evals[0], threshold, rtol=1e-6, atol=1e-15)
 
-        y = cov_nearest(x, method='nearest', n_fact=100, threshold=threshold)
-        evals = np.linalg.eigvalsh(y)
-        #print 'evals', evals, threshold
-        #print evals[0] / threshold - 1
-        assert_allclose(evals[0], threshold, rtol=1e-6, atol=1e-15)
+    y = cov_nearest(x, n_fact=100, threshold=threshold)
+    evals = np.linalg.eigvalsh(y)
+    assert_allclose(evals[0], threshold, rtol=0.25, atol=1e-15)
 
-        y = cov_nearest(x, n_fact=100, threshold=threshold)
-        evals = np.linalg.eigvalsh(y)
-        #print 'evals', evals, threshold
-        #print evals[0] / threshold - 1
-        assert_allclose(evals[0], threshold, rtol=0.25, atol=1e-15)
 
 class Test_Factor(object):
 
-
     def test_corr_nearest_factor_arrpack(self):
 
         # regression results for svds call
@@ -283,74 +276,67 @@ def test_corr_nearest_factor_arrpack(self):
         assert_allclose(u, dsign * u2, rtol=1e-6, atol=1e-14)
         assert_allclose(s, s2, rtol=1e-6)
 
-
-    def test_corr_nearest_factor(self):
+    @pytest.mark.parametrize('dm', [1, 2])
+    def test_corr_nearest_factor(self, dm):
 
         objvals = [np.array([6241.8, 6241.8, 579.4, 264.6, 264.3]),
                    np.array([2104.9, 2104.9, 710.5, 266.3, 286.1])]
 
         d = 100
 
-        for dm in 1,2:
-
-            # Construct a test matrix with exact factor structure
-            X = np.zeros((d,dm), dtype=np.float64)
-            x = np.linspace(0, 2*np.pi, d)
-            np.random.seed(10)
-            for j in range(dm):
-                X[:,j] = np.sin(x*(j+1)) + 1e-10 * np.random.randn(d)
-
-            _project_correlation_factors(X)
-            assert np.isfinite(X).all()
-            X *= 0.7
-            mat = np.dot(X, X.T)
-            np.fill_diagonal(mat, 1.)
-
-            # Try to recover the structure
-            rslt = corr_nearest_factor(mat, dm)
-            err_msg = 'rank=%d, niter=%d' % (dm, len(rslt.objective_values))
-            assert_allclose(rslt.objective_values[:5], objvals[dm - 1],
-                            rtol=0.5, err_msg=err_msg)
-            assert_equal(rslt.Converged, True, err_msg=err_msg)
-            mat1 = rslt.corr.to_matrix()
-            assert_allclose(mat, mat1, rtol=0.25, atol=1e-3, err_msg=err_msg)
-
+        # Construct a test matrix with exact factor structure
+        X = np.zeros((d, dm), dtype=np.float64)
+        x = np.linspace(0, 2 * np.pi, d)
+        np.random.seed(10)
+        for j in range(dm):
+            X[:, j] = np.sin(x * (j + 1)) + 1e-10 * np.random.randn(d)
 
-    # Test that we get the same result if the input is dense or sparse
-    def test_corr_nearest_factor_sparse(self):
+        _project_correlation_factors(X)
+        assert np.isfinite(X).all()
+        X *= 0.7
+        mat = np.dot(X, X.T)
+        np.fill_diagonal(mat, 1.)
 
+        # Try to recover the structure
+        rslt = corr_nearest_factor(mat, dm, maxiter=10000)
+        err_msg = 'rank=%d, niter=%d' % (dm, len(rslt.objective_values))
+        assert_allclose(rslt.objective_values[:5], objvals[dm - 1],
+                        rtol=0.5, err_msg=err_msg)
+        assert_equal(rslt.Converged, True, err_msg=err_msg)
+        mat1 = rslt.corr.to_matrix()
+        assert_allclose(mat, mat1, rtol=0.25, atol=1e-3, err_msg=err_msg)
+
+    @pytest.mark.parametrize('dm', [1, 2])
+    def test_corr_nearest_factor_sparse(self, dm):
+        # Test that result is the same if the input is dense or sparse
         d = 100
 
-        for dm in 1,2:
-
-            # Generate a test matrix of factors
-            X = np.zeros((d,dm), dtype=np.float64)
-            x = np.linspace(0, 2*np.pi, d)
-            np.random.seed(10)
-            for j in range(dm):
-                X[:,j] = np.sin(x*(j+1)) + 1e-10 * np.random.randn(d)
-
-            # Get the correlation matrix
-            _project_correlation_factors(X)
-            X *= 0.7
-            mat = np.dot(X, X.T)
-            np.fill_diagonal(mat, 1)
+        # Generate a test matrix of factors
+        X = np.zeros((d, dm), dtype=np.float64)
+        x = np.linspace(0, 2 * np.pi, d)
+        np.random.seed(10)
+        for j in range(dm):
+            X[:, j] = np.sin(x * (j + 1)) + 1e-10 * np.random.randn(d)
 
-            # Threshold it
-            mat *= (np.abs(mat) >= 0.35)
-            smat = sparse.csr_matrix(mat)
+        # Get the correlation matrix
+        _project_correlation_factors(X)
+        X *= 0.7
+        mat = np.dot(X, X.T)
+        np.fill_diagonal(mat, 1)
 
-            rslt = corr_nearest_factor(smat, dm)
-            assert_equal(rslt.Converged, True)
-            mat_dense = rslt.corr.to_matrix()
+        # Threshold it
+        mat *= (np.abs(mat) >= 0.35)
+        smat = sparse.csr_matrix(mat)
 
-            rslt = corr_nearest_factor(smat, dm)
-            assert_equal(rslt.Converged, True)
-            mat_sparse = rslt.corr.to_matrix()
+        rslt = corr_nearest_factor(mat, dm, maxiter=10000)
+        assert rslt.Converged is True
+        mat_dense = rslt.corr.to_matrix()
 
-            assert_allclose(mat_dense, mat_sparse, rtol=0.25,
-                            atol=1e-3)
+        rslt = corr_nearest_factor(smat, dm, maxiter=10000)
+        assert rslt.Converged is True
+        mat_sparse = rslt.corr.to_matrix()
 
+        assert_allclose(mat_dense, mat_sparse, rtol=.25, atol=1e-3)
 
     # Test on a quadratic function.
     def test_spg_optim(self):
@@ -419,52 +405,49 @@ def test_solve(self):
         sr2 = np.linalg.solve(mat, rhs)
         assert_almost_equal(sr1, sr2)
 
-    def test_cov_nearest_factor_homog(self):
+    @pytest.mark.parametrize('dm', [1, 2])
+    def test_cov_nearest_factor_homog(self, dm):
 
         d = 100
 
-        for dm in 1,2:
-
-            # Construct a test matrix with exact factor structure
-            X = np.zeros((d,dm), dtype=np.float64)
-            x = np.linspace(0, 2*np.pi, d)
-            for j in range(dm):
-                X[:,j] = np.sin(x*(j+1))
-            mat = np.dot(X, X.T)
-            np.fill_diagonal(mat, np.diag(mat) + 3.1)
-
-            # Try to recover the structure
-            rslt = cov_nearest_factor_homog(mat, dm)
-            mat1 = rslt.to_matrix()
+        # Construct a test matrix with exact factor structure
+        X = np.zeros((d,dm), dtype=np.float64)
+        x = np.linspace(0, 2*np.pi, d)
+        for j in range(dm):
+            X[:,j] = np.sin(x*(j+1))
+        mat = np.dot(X, X.T)
+        np.fill_diagonal(mat, np.diag(mat) + 3.1)
 
-            assert_allclose(mat, mat1, rtol=0.25, atol=1e-3)
+        # Try to recover the structure
+        rslt = cov_nearest_factor_homog(mat, dm)
+        mat1 = rslt.to_matrix()
 
+        assert_allclose(mat, mat1, rtol=0.25, atol=1e-3)
 
-    # Check that dense and sparse inputs give the same result
-    def test_cov_nearest_factor_homog_sparse(self):
+    @pytest.mark.parametrize('dm', [1, 2])
+    def test_cov_nearest_factor_homog_sparse(self, dm):
+        # Check that dense and sparse inputs give the same result
 
         d = 100
 
-        for dm in 1,2:
-
-            # Construct a test matrix with exact factor structure
-            X = np.zeros((d,dm), dtype=np.float64)
-            x = np.linspace(0, 2*np.pi, d)
-            for j in range(dm):
-                X[:,j] = np.sin(x*(j+1))
-            mat = np.dot(X, X.T)
-            np.fill_diagonal(mat, np.diag(mat) + 3.1)
+        # Construct a test matrix with exact factor structure
+        X = np.zeros((d,dm), dtype=np.float64)
+        x = np.linspace(0, 2*np.pi, d)
+        for j in range(dm):
+            X[:,j] = np.sin(x*(j+1))
+        mat = np.dot(X, X.T)
+        np.fill_diagonal(mat, np.diag(mat) + 3.1)
 
-            # Fit to dense
-            rslt = cov_nearest_factor_homog(mat, dm)
-            mat1 = rslt.to_matrix()
+        # Fit to dense
+        rslt = cov_nearest_factor_homog(mat, dm)
+        mat1 = rslt.to_matrix()
 
-            # Fit to sparse
-            smat = sparse.csr_matrix(mat)
-            rslt = cov_nearest_factor_homog(smat, dm)
-            mat2 = rslt.to_matrix()
+        # Fit to sparse
+        smat = sparse.csr_matrix(mat)
+        rslt = cov_nearest_factor_homog(smat, dm)
+        mat2 = rslt.to_matrix()
 
-            assert_allclose(mat1, mat2, rtol=0.25, atol=1e-3)
+        assert_allclose(mat1, mat2, rtol=0.25, atol=1e-3)
 
     def test_corr_thresholded(self):