From 51fc20585785a741be222e506daf9119ca75fcec Mon Sep 17 00:00:00 2001
From: szorowi1 <szorowi1@gmail.com>
Date: Tue, 20 Sep 2022 16:37:46 -0400
Subject: [PATCH] random initial values for rotation matrix in GPA rotations

---
 factor_analyzer/rotator.py | 21 ++++++++++++++++++---
 1 file changed, 18 insertions(+), 3 deletions(-)

diff --git a/factor_analyzer/rotator.py b/factor_analyzer/rotator.py
index 097c1e4..5a74595 100644
--- a/factor_analyzer/rotator.py
+++ b/factor_analyzer/rotator.py
@@ -69,6 +69,10 @@ class Rotator(BaseEstimator):
     tol : float, optional
         The convergence threshold. Used for 'varimax' and 'oblique' rotations.
         Defaults to 1e-5.
+    random_state : int, optional
+        The seed for initializing the rotation matrix. Orthogonal matrix
+        generated from 'scipy.stats.ortho_group'. If ``None``, the default
+        is the identity matrix ('np.eye').
 
     Attributes
     ----------
@@ -120,6 +124,7 @@ def __init__(
         delta=0.01,
         max_iter=500,
         tol=1e-5,
+        random_state=None,
     ):
         """Initialize the rotator class."""
         self.method = method
@@ -130,6 +135,7 @@ def __init__(
         self.delta = delta
         self.max_iter = max_iter
         self.tol = tol
+        self.random_state = random_state
 
         self.loadings_ = None
         self.rotation_ = None
@@ -323,7 +329,10 @@ def _oblique(self, loadings, method):
 
         # initialize the rotation matrix
         _, n_cols = loadings.shape
-        rotation_matrix = np.eye(n_cols)
+        if self.random_state is None:
+            rotation_matrix = np.eye(n_cols)
+        else:
+            rotation_matrix = sp.stats.ortho_group(n_cols, seed=self.random_state).rvs()
 
         # default alpha level
         alpha = 1
@@ -413,7 +422,10 @@ def _orthogonal(self, loadings, method):
 
         # initialize the rotation matrix
         _, n_cols = arr.shape
-        rotation_matrix = np.eye(n_cols)
+        if self.random_state is None:
+            rotation_matrix = np.eye(n_cols)
+        else:
+            rotation_matrix = sp.stats.ortho_group(n_cols, seed=self.random_state).rvs()
 
         # default alpha level
         alpha = 1
@@ -490,7 +502,10 @@ def _varimax(self, loadings):
 
         # initialize the rotation matrix
         # to N x N identity matrix
-        rotation_mtx = np.eye(n_cols)
+        if self.random_state is None:
+            rotation_mtx = np.eye(n_cols)
+        else:
+            rotation_mtx = sp.stats.ortho_group(n_cols, seed=self.random_state).rvs()
 
         d = 0
         for _ in range(self.max_iter):