# Factor Analysis does not give a better covariance estimate than the Empirical Covariance matrix?

## I do not see that Factor Analysis gives a better covariance estimate than the empirical covariance estimate, from the toy data simulation with explanation and code below. Am I doing something wrong?

### Generative model setup

From the CS229 lecture notes on Factor Analysis, I built a (random) generative model for Factor Analysis with observation dimension $$N$$, latent dimension $$K$$ and number of samples $$M$$, and a ground truth mean of $$\mu=0$$:

$$x=Wz+\epsilon$$

with

$$W=\Lambda\sim U(0,1)^{N\times K}\in\mathbb{R}^{N\times K}\\ z\sim\mathcal{N}(0,I_K)\in\mathbb{R}^{K}\\ \epsilon \sim\mathcal{N}(0,\Psi)\in\mathbb{R}^{N}\\ \Psi=diag(\psi)\\ \psi\sim U(0,1)^{N}$$

From this we can construct a (derived) ground truth covariance matrix $$\Sigma\in\mathbb{R}^{100\times100}$$. I use an observation dimension of $$N=100$$ and latent dimension of $$K=10$$.

### Sampling

I sample data from this generative model for a different number of total observed samples $$m$$, yielding a data matrix $$X\in\mathbb{R}^{N\times M}$$ and for every specific number of samples I average the covariance estimation error across 20 different seeds. I compare the covariance matrix recovered by Factor Analysis $$\hat{\Sigma}_{FA}$$ to the empirical covariance matrix $$\hat{\Sigma}_{emp}$$. This is done in terms of the difference with respect to the true covariance matrix under a Frobenius norm: $$||\Sigma - \hat{\Sigma}||_F$$

### Results

Strangely enough, the Factor Analysis covariance estimate seems to not perform significantly better than the empirical covariance estimate, even though it does get the true latent dimension size ($$K=10$$) given:

To put the relative error (%) w.r.t. the true covariance in perspective (where |·| is taken elementwise):

Results with higher $$N$$, seem to show little (relative) change. Here for $$N=1000, K=10$$:

To put the relative error (%) w.r.t. the true covariance in perspective (where |·| is taken elementwise):

### Question

I would expect the Factor Analysis covariance estimate to have a significantly lower error w.r.t. the true covariance than the empirical covariance estimate, especially in the small sample setting. Am I doing something wrong?

### Code

The following code generates these results (for Frobenius norm):

import numpy as np
from sklearn.decomposition import PCA, FactorAnalysis
import matplotlib.pyplot as plt
from tqdm.notebook import tqdm
import scipy
from matplotlib.pyplot import figure

obs_dim = 100 # observation dimension
lat_dim = 10 # latent dimension
n_random_draws = 20 # number of trials to average over

fa_errors = []
fa_errors_error = []
emp_errors = []
emp_errors_error = []

n_samples_list = np.rint(np.logspace(1, 2, 10)).astype(int)

x = []
for n_samples in tqdm(n_samples_list):
fa_errors_seed = []
emp_errors_seed = []
for seed in range(n_random_draws):
np.random.seed(seed)
W = np.random.rand(obs_dim, lat_dim)/10
Z = np.random.multivariate_normal(np.zeros(lat_dim), np.eye(lat_dim), size=n_samples).T

gt_noise_variance = np.random.rand(obs_dim)/1000.
psi = np.diag(gt_noise_variance)
gt_errors = np.random.multivariate_normal(np.zeros(obs_dim), psi, size=n_samples).T

true_corrcov = W.dot(W.T)
gt_cov = true_corrcov + psi

X = W.dot(Z) + gt_errors

emp_cov = np.cov(X)

model_fa = FactorAnalysis(n_components=lat_dim)
model_fa.fit(X.T)
est_fa_cov = model_fa.get_covariance()

fa_error = np.linalg.norm(gt_cov - est_fa_cov, ord="fro")
fa_errors_seed.append(fa_error)

emp_error = np.linalg.norm(gt_cov - emp_cov, ord="fro")
emp_errors_seed.append(emp_error)

fa_errors.append(np.mean(np.array(fa_errors_seed)))
fa_errors_error.append(np.std(np.array(fa_errors_seed)))

emp_errors.append(np.mean(np.array(emp_errors_seed)))
emp_errors_error.append(np.std(np.array(emp_errors_seed)))

x.append(n_samples)

plt.rcParams.update({'font.size': 20})
plt.rcParams['text.usetex'] = True
from_end = 10

fig, ax = plt.subplots(figsize=(16, 8), tight_layout=True)
ax.errorbar(x[:from_end], fa_errors[:from_end], yerr = fa_errors_error[:from_end], label=r"Factor Analysis Covariance estimate error $$||\Sigma - \hat{{{\Sigma}}}_{{{FA}}}||_F$$", capsize=10)
ax.errorbar(x[:from_end], emp_errors[:from_end], yerr = emp_errors_error[:from_end], label=r"Empirical Covariance estimate error $$||\Sigma - \hat{{{\Sigma}}}_{{{emp}}}||_F$$", linestyle="dashed", capsize=10)
ax.set_ylabel(r"$$||\Sigma - \hat{{{\Sigma}}}||_F$$")
ax.set_xlabel("number of samples")
plt.grid()
plt.legend()
plt.show()

• I don't do this sort of thing in python, so I'm not 100% sure, but I think you need to make N larger to see a difference. Mar 3, 2022 at 17:38
• @JeremyMiles thanks for the suggestion! I added the result with 10x higher N, but there seems to be little change in the relative performance of FA and its empirical counterpart. Also, suggestions outside of Python are welcome to maybe solve the question differently? R? Mar 3, 2022 at 22:38
• I can't quite read the uncommented python code. Are you sure you are doing FA extraction and not PCA extraction? Mar 5, 2022 at 19:45
• Thanks for the comment @ttnphns! On top I import sklearn.decomposition. FactorAnalysis and then (1) initialize it as model_fa = FactorAnalysis(n_components=lat_dim) and (2) fit it on the data with model_fa.fit(X.T). Then I extract the covariance matrix with est_fa_cov = model_fa.get_covariance().The official documentation provides more info and examples: scikit-learn.org/stable/modules/generated/… Mar 6, 2022 at 20:53
• Are you sure your est_fa_cov = model_fa.get_covariance returns you the full reproduced covariance matrix (i.e., with original variances on the diagonal) and not the reduced reproduced one (i.e., with communalities on the diagonal)? Mar 7, 2022 at 1:23

I am assuming (without going to deal with python functions) that your post is in fact asking: Is covariance matrix $$\bf\hat R$$ reproduced by Factor analysis of a sample covariance matrix $$\bf R$$ a better (closer) estimate of the population ("true") covariance matrix $$\bf R_t$$ than $$\bf R$$ is? That is, is $$\bf ||R_t-\hat R||$$ smaller than $$\bf ||R_t-R||$$, where $$\bf R$$ is a random realization of $$\bf R_t$$ for some sample size $$n$$?

I've conducted (in SPSS) a series of simulations, each one randomly creating some population factor loading p x m matrix $$\bf A_t$$ (it is what you call W); $$p$$ variables = 50 and $$m$$ factors = 5. According to the factor theorem (footnote 1), this yields us the corresponding population covariance matrix $$\bf R_t = A_tA_t^{'}+\text{diag}(u_t{^2})$$, $$\bf u_t{^2}$$ being the p-length vector of some uniquenesses, also randomly generated. So, this $$\bf R_t$$ is the m common factor covariance matrix.

For each such population $$\bf R_t$$, I generated 1000 sample covariance matrices $$\bf R$$ under normal population assumption and the sample size, say, $$n$$=200 observations. I did not actually generated the data; rather, I used Wishart distribution to get $$\bf R$$s (you did equivalently by actually generating the data sample according to the factor model and then computing covariance matrix of your variables).

On each $$\bf R$$, I performed Factor analysis (like you did it on your sample data), by Principal Axis factoring method with initial communalities estimated as the images (this is the usual way), extracting precisely $$m$$ factors. The obtained loadings $$\bf A$$ then yielded the reproduced covariance matrix $$\bf \hat R = AA^{'}$$ with the final communalities on the diagonal. But we need the full, nonreduced matrix, so the diagonal was then replaced by the diagonal of $$\bf R$$. (Thus, diagonals of $$\bf R$$ and $$\bf \hat R$$ are equal.)

Departures $$\bf ||R_t-\hat R||$$ as well as $$\bf ||R_t-R||$$ were recorded for each sample, and their difference $$\bf ||R_t-\hat R|| - ||R_t-R||$$ computed.

This difference averaged over the 1000 samples was negative, - and so was observed in every simulation of population factor structure $$\bf R_t$$. Moreover, the difference was negative in all individual samples. We may conclude that Factor analysis of a sample covariance matrix (or sample data) yields an estimate of the population (true) covariances which is (always) a more accurate estimate than the sample covariance matrix itself. Of course, provided you are extracting in your FA the true number of factors $$m$$ (but who will tell you the true $$m$$ for certain, in practice, when all what you have is just one sample of data sized $$n$$, and n is not very large?)

[Note: the difference $$\bf ||R_t-\hat R|| - ||R_t-R||$$ is non zero only for off-diagonal covariances; since diagonals of $$\bf R$$ and $$\bf \hat R$$ are equal - see above.]

I've conducted the alike simulation series also with Principal Component analysis in place of FA, to extract factors and get the $$\bf \hat R$$. The difference $$\bf ||R_t-\hat R|| - ||R_t-R||$$ now was positive, not in every sample though, but in the great majority of samples; and positive overall (averaged). We may say that PCA of a sample covariance matrix (or sample data) yields an estimate of the population (true) covariances which is (most of the time) a less accurate estimate than the sample covariance matrix itself.

But PCA can approach FA as a latent structure discoverer under certain conditions. More comparison of FA and PCA with an extensive simultion study can be found in this big thread, particularly in my answer.

• First, thanks a lot for the detailed explanation and simulation! Three remarks/questions: (1) could you share numerical results/graphs? (2) I don't see how replacing the diagonal is motivated in FA; loadings^T*loadings + psi should be equivalent, no? (it is not in simulation :( ) (3) I am still underwhelmed by FA's covariance estimate in my simulations: only sometimes a few % better than empirical (in Frobenius norm), and on average 10-20% away from the true cov in a small sample setting (note that some (co)variances are estimated up to 100% wrong). I would expect like 1% relative error. Mar 7, 2022 at 17:47
• (1) Sorry, I will not give pics; the primary reason is that pics would be expressive if I did this time simulations on different numbers of n, p and m, but I had done on them only 200, 50, an 5, respectively. (2) In FA, the diagonal of AA' is called communalities and their complements to the diagonal of the analyzed cov. matrix - the variances - are called uniquenesses. So if you are comparing some full cov. matrices, i.e. having variances on their diagonals, you must have variances there in them both (and not to have variances in one matrix and communalities in the other, as the diagonal). Mar 7, 2022 at 18:15
• I might recommend you to do exactly like and what I did, if you wish. Will you arrive at the same conclusions or not? Mar 7, 2022 at 18:19
• Will you arrive at the same conclusions or not?. Yes. But, "positive" or "negative" only give directions, not sizes or percentages; covariance matrix items which are 10-20% off on average (compared to the true covariance) are still quite unsatisfactory, despite FA doing consistently better than the empirical covariance. That is why I am interested in the % w.r.t. true covariance. With these results it seems that FA is useless for covariance estimation, with its consistent but marginal improvement over the empirical covariance estimate. Mar 7, 2022 at 18:28
• P.S. To (2) point. FA does not estimate variances; it estimates only co-variances. The variances - it splits them each into two summands: the communality part and the uniqueness part. The communality part is the direct function of the loadings, so can be said is "estimated", like loadings are, but the uniqueness is simply the counterpart left. The sum of the two summonds always equals the input variance. Mar 7, 2022 at 18:29