Add minres factor analysis (#20)

* Initial commit of minimum residual factor analysis.

* Added minimum residual factor analysis and unittests.

factoranalysis/__init__.py

@@ -1,5 +1,6 @@
 from .pca import *
 from .paf import *
 from .minimum_rank import *
+from .minimum_residual import *
 from .maximum_likelihood import *
 from .rotation import *

factoranalysis/minimum_residual.py

@@ -0,0 +1,63 @@
+import numpy as np
+from scipy.optimize import minimize
+
+from RyStats.factoranalysis import principal_components_analysis as pca
+
+
+def _minres_min_func(unique_variance, correlation_matrix, n_factors, 
+                     correlation_diagonal, indices):
+    """Min function for minimum residual factor analysis"""
+    np.fill_diagonal(correlation_matrix, 
+                     correlation_diagonal - unique_variance)
+
+    eigs, vects = np.linalg.eigh(correlation_matrix)
+
+    eigs = eigs[-n_factors:]
+    vects = vects[:, -n_factors:]
+    vects2 = vects * eigs.reshape(1, -1)
+
+    updated_corr = vects @ vects2.T
+    lower_difference = correlation_matrix[indices] - updated_corr[indices]
+
+    return np.square(lower_difference).sum()
+
+
+def minres_factor_analysis(input_matrix, n_factors, initial_guess=None):
+    """Performs minimum residual factor analysis.
+
+    Minimium Residual factor analysis is equivalent to unweighted least-squares
+    and also equal to Principal Axis Factor if Reduced Matrix remains positive
+    definite.
+
+    Args:
+        input_matrx: Correlation or Covariance Matrix
+        n_factors:  number of factors to extract
+        initial_guess: Guess to seed the search algorithm
+
+    Returns:
+        loadings: unrotated extracted factor loadings
+        eigenvalues: eigenvalues of extracted factor loadings
+        unique_variance: variance unique to each item
+    """
+    working_matrix = input_matrix.copy()
+    diagonal_from_input = np.diag(input_matrix)
+    lower_indices = np.tril_indices_from(input_matrix, k=-1)
+
+    # Initial Guess
+    loads, _, _ = pca(input_matrix, n_factors)
+
+    if initial_guess is None:
+        uvars = np.diag(input_matrix - loads @ loads.T).copy()
+    else:
+        uvars = initial_guess.copy()
+
+    args = (working_matrix, n_factors, diagonal_from_input, lower_indices)
+    bounds = [(0.01, .99 * upper) for upper in diagonal_from_input]
+
+    result = minimize(_minres_min_func, uvars, args, method='SLSQP',
+                      bounds=bounds)
+    unique_variance = result['x']
+
+    loads, eigs, _ = pca(input_matrix - np.diag(unique_variance), n_factors)
+
+    return loads, eigs, unique_variance

factoranalysis/test/test_maximum_likelihood.py

@@ -69,22 +69,5 @@ def test_maximum_likelihood_modern_with_classic(self):
         np.testing.assert_allclose(variance_modern, variance_classic, rtol=1e-4)
         np.testing.assert_allclose(variance_modern, variance_classic2, rtol=1e-4)
 
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
 if __name__ == "__main__":
     unittest.main()

factoranalysis/test/test_minimum_rank.py

@@ -66,10 +66,6 @@ def no_derivative(inverse_half_variance, correlation_cholesky, n_factors):
                                              n_factors)
             np.testing.assert_allclose(result['jac'], derivative_calc[1], atol=1e-5)
 
-
-
-
-
     def test_minimum_zero_eigenvalue(self):
         """Testing Forced Semi-Positive Definite."""
         rng = np.random.default_rng(12473)

factoranalysis/test/test_minimum_residual.py

@@ -0,0 +1,71 @@
+import unittest
+
+import numpy as np
+
+from RyStats.factoranalysis import principal_components_analysis as pca
+from RyStats.factoranalysis import minres_factor_analysis as mfa
+
+from RyStats.common import procrustes_rotation
+
+
+class TestMaximumLikelihood(unittest.TestCase):
+    """Test fixture for minimum residual factor analysis."""
+
+    #TODO: Need algorithm validity test
+
+    def test_minimum_residual_recovery(self):
+        """Testing Minimum Residual Recovery."""
+        rng = np.random.default_rng(49432132341221348721323123324)
+
+        data = rng.uniform(-2, 2, size=(10, 100))
+        unique_var = rng.uniform(0.2, 2, size=10)
+
+        # Create 3 Factor Data
+        cor_matrix = np.cov(data)
+        loadings, eigenvalues, _ = pca(cor_matrix, 3)
+
+        # Add Unique variance
+        cor_matrix2 = loadings @ loadings.T + np.diag(unique_var)
+
+        initial_guess = np.ones((10,)) *.5
+        loadings_paf, _, variance = mfa(cor_matrix2, 3, initial_guess=initial_guess)
+
+        # Remove any rotation
+        rotation = procrustes_rotation(loadings, loadings_paf)
+        updated_loadings = loadings_paf @ rotation
+        updated_eigs = np.square(updated_loadings).sum(0)
+
+        # Did I Recover initial values (upto a rotation)
+        np.testing.assert_allclose(loadings, updated_loadings, atol=1e-3)
+        np.testing.assert_allclose(eigenvalues, updated_eigs, atol=1e-3)
+        np.testing.assert_allclose(unique_var, variance, atol=1e-3)
+
+    def test_minimum_residual_recovery2(self):
+        """Testing Minimum Residual Recovery no initial guess."""
+        rng = np.random.default_rng(498556324111616321324125213)
+
+        data = rng.uniform(-2, 2, size=(10, 100))
+        unique_var = rng.uniform(0.2, 2, size=10)
+
+        # Create 3 Factor Data
+        cor_matrix = np.cov(data)
+        loadings, eigenvalues, _ = pca(cor_matrix, 3)
+
+        # Add Unique variance
+        cor_matrix2 = loadings @ loadings.T + np.diag(unique_var)
+
+        initial_guess = np.ones((10,)) *.5
+        loadings_paf, _, variance = mfa(cor_matrix2, 3)
+
+        # Remove any rotation
+        rotation = procrustes_rotation(loadings, loadings_paf)
+        updated_loadings = loadings_paf @ rotation
+        updated_eigs = np.square(updated_loadings).sum(0)
+
+        # Did I Recover initial values (upto a rotation)
+        np.testing.assert_allclose(loadings, updated_loadings, atol=1e-3)
+        np.testing.assert_allclose(eigenvalues, updated_eigs, atol=1e-3)
+        np.testing.assert_allclose(unique_var, variance, atol=1e-3)        
+
+if __name__ == "__main__":
+    unittest.main()

factoranalysis/test/test_paf.py

@@ -5,6 +5,7 @@
 from RyStats.factoranalysis import principal_components_analysis as pca
 from RyStats.factoranalysis import principal_axis_factor as paf
 
+
 class TestPrincipalAxis(unittest.TestCase):
     """Test fixture for principal components."""
 
@@ -28,6 +29,5 @@ def test_principal_axis_factor(self):
         np.testing.assert_allclose(eigenvalues, eigenvalues2, rtol=1e-4)
         np.testing.assert_allclose(unique_var, variance, rtol=1e-4)
 
-
 if __name__ == "__main__":
     unittest.main()