Currently, PCA decomposes the covariance matrix using the Jacobi eigenvalue algorithm. A singular value decomposition (SVD) would probably be preferable. It would likely be more numerically stable as it avoids computing the covariance matrix. It would also likely be in two phases: bidiagonalization and full diagonalization. This may help with the recursion depth issue.
Currently, PCA decomposes the covariance matrix using the Jacobi eigenvalue algorithm. A singular value decomposition (SVD) would probably be preferable. It would likely be more numerically stable as it avoids computing the covariance matrix. It would also likely be in two phases: bidiagonalization and full diagonalization. This may help with the recursion depth issue.