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

Question

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()

Answer 1

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$?

My answer would be yes.

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.