I am trying to match the factor coefficient matrix from psych package for factor analysis with manual calculation. My book says that the principal axis method takes the first few principal components of $\textbf S - \hat \Psi_1$ where $\hat \Psi _1 = \text {diag}(\hat \psi^2_1,\dots,\hat \psi_d^2)$ and $\hat\psi_j^2=1/s^{jj}$ where $\textbf S^{-1}=[(s^{jk})]$. For the model $\textbf x = \mu+\Gamma\textbf f+\epsilon$ then $\Gamma = \textbf M_1 \Phi_1^{1/2}$ is a solution where $\textbf M_1$ are the first few positive standardized eigenvectors of $\textbf S-\hat\Psi_1$ and $\Phi_1$ is diagonal of the first few positive eigenvalues.
My code is
library(psych)
covariances <- ability.cov$cov
correlations <- cov2cor(covariances)
fact <- fa(correlations, nfactors=2, rotate="none", fm="pa")
fact
Psi1 <- diag(c(1/solve(covariances)[1,1], 1/solve(covariances)[2,2], 1/solve(covariances)[3,3],
1/solve(covariances)[4,4], 1/solve(covariances)[5,5], 1/solve(covariances)[6,6]))
Phi1 <- diag(eigen(covariances - Psi1)$values[1:2])
M1 <-eigen(covariances - Psi1)$vectors[, 1:2]
library(expm)
M1 %*% sqrtm(Phi1)
The results of fact and M1 %*% sqrtm(Phi1) don't match. I'm wondering what change can be made to make them the same.
