Is there any good reason to use PCA instead of EFA? Also, can PCA be a substitute for factor analysis? In some disciplines, PCA (principal component analysis) is systematically used without any justification, and PCA and EFA (exploratory factor analysis) are considered as synonyms.
I therefore recently used PCA to analyse the results of a scale validation study (21 items on 7-points Likert scale, assumed to compose 3 factors of 7 items each) and a reviewer asks me why I chose PCA instead of EFA. I read about the differences between both techniques, and it seems that EFA is favored against PCA in a majority of your answers here. 
Do you have any good reasons for why PCA would be a better choice? What benefits it could provide and why it could be a wise choice in my case?
 A: As you said, you are familiar with relevant answers; see also: So, as long as "Factor analysis..." + a couple of last paragraphs; and the bottom list here. In short, PCA is mostly a data reduction technique whereas FA is a modeling-of-latent-traits technique. Sometimes they happen to give similar results; but in your case - because you probably feel like constructing/validating latent traits as if real entities - using FA would be more honest and you shouldn't prefer PCA in hope that their results converge. On the other hand, whenever you aim to summarise/simplify the data - for subsequent analysis, for example - you would prefer PCA, as it doesn't impose any strong model (which might be irrelevant) on the data.
To reiterate other way, PCA gives you dimensions which may correspond to some subjectively meaningful constructs, if you wish, while EFA poses that those are even covert features that actually generated your data, and it aims to find those features. In FA, interpretation of the dimensions (factors) is pending - whether you can attach a meaning to a latent variable or not, it "exists" (FA is essentialistic), otherwise you should drop it from the model or get more data to support it. In PCA, the meaning of a dimension is optional.
And yet once again in other words: When you extract m factors (separate factors from errors), these few factors explain (almost) all correlation among variables, so that the variables are not left room to correlate via the errors anyhow. Therefore, so long as "factors" are defined as latent traits which generate/bind the correlated data, you have full clues to interpret that - what is responsible for the correlations. In PCA (extract components as if "factors"), errors (may) still correlate between the variables; so you can't claim that you've extracted something enough clean and exhaustive to be interpreted in that way.
You may want to read my other, longer answer in the current discussion, for some theoretical and simulation experiment details about whether PCA is a viable substitute of FA. Please pay attention also to outstanding answers by @amoeba given on this thread.

Upd: In their answer to this question @amoeba, who opposed there, introduced a (not well-known) technique PPCA as standing halfway between PCA and FA. This naturally launched the logic that PCA and FA are along one line rather than opposite. That valuable approach expands one's theoretical horizons. But it can mask the important practical difference about that FA reconstructs (explains) all the pairwise covariances with a few factors, while PCA cannot do it successfully (and when it occasionally does it - that is because it happened to mime FA).
A: (This is really a comment to @ttnphns's second answer)
As far as the different type of reproduction of covariance assuming error  by PC and by FA is concerned, I've simply printed out the loadings/components of variance which occur in the two procedures; just for the examples I took 2 variables.
We assume the construction of the two items as of one common factor and itemspecific factors. Here is that factor-loadingsmatrix:
  L_fa: 
          f1       f2      f3         
  X1:   0.894    0.447     .             
  X1:   0.894     .       0.447              

The correlation matrix by this is
  C:
         X1       X2 
  X1:   1.000   0.800
  X2:   0.800   1.000

If we look at the loadings-matrix L_fa and interpret it as usual in FA that f2 and f3 are error terms/itemspecific error, we reproduce C without that error, receiving
 C1_Fa 
        X1       X2 
 X1:  0.800   0.800
 X2:  0.800   0.800

So we have perfectly reproduced the off-diagonal element, which is the covariance (and the diagonal is reduced)
If we look at the pca-solution (can be done by simple rotations) we get the two factors from the same correlation-matrix:
 L_pca : 
         f1        f2
 X1:   0.949      -0.316
 X2:   0.949       0.316

Assuming the second factor as error we get the reproduced matrix of covariances
  C1_PC : 
        X1      X2
 X1:   0.900   0.900
 X2:   0.900   0.900

where we've overestimated the true correlation. This is because we ignored the correcting negative partial covariance in the second factor = error.
Note that the PPCA would be identical with the first example.
With more items this is no more so obvious but still an inherent effect. Therefore there is also the concept of MinRes-extraction (or -rotation?) and I've also seen something like maximum-determinant extraction and...

[update] As for the question of @amoeba:     
I understood the concept of "Minimal Residuals" ("MinRes")-rotation as a concurring method to the earlier methods of CFA-computation, to achieve the best reproduction of the off-diagonal elements of a correlation matrix. I learned this in the 80'ies/90'ies and didn't follow the development of factor-analysis (as indepth as before in the recent years), so possibly "MinRes" is out of fashion.
To compare it with the PCA-solution: one can think of finding the pc-solution by rotations of the factors when they are thought as axes in an euclidean space and the loadings are the coordinates of the items in that vectorspace.
Then for a pair of axes say x,y the sums-of-squares from the loadings of the x-axis and that of the y-axis are computed.
From this one can find a rotation angle, by which we should rotate, to get the sums-of-squares in the rotated axes maximal on the x° and minimal on the y°-axis (where the little circle indicates the rotated axes).
Doing this for all pairs of axes (where only always the x-axis is the left and the y-axis is the right (so for 4 factors we have only 6 pairs of rotation)) and then repeat the whole process to a stable result realizes the so-called "Jacobi-method" for the finding of the principal components solution: it will locate the first axis such that it collects the maximum possible sum of squares of loadings ("SSqL") (which means also "of the variance") on one axis in the current correlational configuration.
As far as I understood things, "MinRes" should look at the partial correlations instead of the SSqL; so it does not sum up the squares of the loadings (as done in the Jacobi-pc-rotation) but is sums up the crossproducts of the loadings in each factor - except of the "crossproducts" (=squares) of the loadings of each item with itself.
After the criteria for the x and for the y-axis are computed it proceeds the same way as described for the iterative Jacobi-rotation.
Since the rotation-criterion is numerically different from the maximum-SSqL-criterion the result/the rotational position shall be different from the PCA-solution. If it converges it should instead provide the maximum possible partial correlation on one axis in the first factor, the next maximal correlation on the next factor and so on. The idea seems to be, then to assume so many axes/factors such that the remaining/residual partial covariance becomes marginal.
(Note this is only how I interpreted things, I've not seen that procedure explicitly written out (or cannot remember at the moment); a description at mathworld seems to express it rather in terms of the formulae like in amoeba's answer) and is likely more authoritative. Just found another reference in the R-project documentation and a likely very good reference in the Gorsuch book on factoranalysis, page 116, available via google-books)
A: Disclaimer: @ttnphns is very knowledgeable about both PCA and FA, and I respect his opinion and have learned a lot from many of his great answers on the topic. However, I tend to disagree with his reply here, as well as with other (numerous) posts on this topic here on CV, not only his; or rather, I think they have limited applicability.

I think that the difference between PCA and FA is overrated.
Look at it like that: both methods attempt to provide a low-rank approximation of a given covariance (or correlation) matrix. "Low-rank" means that only a limited (low) number of latent factors or principal components is used. If the $n \times n$ covariance matrix of the data is $\mathbf C$, then the models are:
\begin{align}
\mathrm{PCA:} &\:\:\: \mathbf C \approx \mathbf W \mathbf W^\top \\
\mathrm{PPCA:}  &\:\:\: \mathbf C \approx \mathbf W \mathbf W^\top + \sigma^2 \mathbf I \\
\mathrm{FA:}  &\:\:\: \mathbf C \approx \mathbf W \mathbf W^\top + \boldsymbol \Psi
\end{align}
Here $\mathbf W$ is a matrix with $k$ columns (where $k$ is usually chosen to be a small number, $k<n$), representing $k$ principal components or factors, $\mathbf I$ is an identity matrix, and $\boldsymbol \Psi$ is a diagonal matrix. Each method can be formulated as finding $\mathbf W$ (and the rest) minimizing the [norm of the] difference between left-hand and right-hand sides.
PPCA stands for probabilistic PCA, and if you don't know what that is, it does not matter so much for now. I wanted to mention it, because it neatly fits between PCA and FA, having intermediate model complexity. It also puts the allegedly large difference between PCA and FA into perspective: even though it is a probabilistic model (exactly like FA), it actually turns out to be almost equivalent to PCA ($\mathbf W$ spans the same subspace).
Most importantly, note that the models only differ in how they treat the diagonal of $\mathbf C$. As the dimensionality $n$ increases, the diagonal becomes in a way less and less important (because there are only $n$ elements on the diagonal and $n(n-1)/2 = \mathcal O (n^2)$ elements off the diagonal). As a result, for the large $n$ there is usually not much of a difference between PCA and FA at all, an observation that is rarely appreciated. For small $n$ they can indeed differ a lot.
Now to answer your main question as to why people in some disciplines seem to prefer PCA. I guess it boils down to the fact that it is mathematically a lot easier than FA (this is not obvious from the above formulas, so you have to believe me here):


*

*PCA -- as well as PPCA, which is only slightly different, -- has an analytic solution, whereas FA does not. So FA needs to be numerically fit, there exist various algorithms of doing it, giving possibly different answers and operating under different assumptions, etc. etc. In some cases some algorithms can get stuck (see e.g. "heywood cases"). For PCA you perform an eigen-decomposition and you are done; FA is a lot more messy.
Technically, PCA simply rotates the variables, and that is why one can refer to it as a mere transformation, as @NickCox did in his comment above. 

*PCA solution does not depend on $k$: you can find first three PCs ($k=3$) and the first two of those are going to be identical to the ones you would find if you initially set $k=2$. That is not true for FA: solution for $k=2$ is not necessarily contained inside the solution for $k=3$. This is counter-intuitive and confusing.
Of course FA is more flexible model than PCA (after all, it has more parameters) and can often be more useful. I am not arguing against that. What I am arguing against, is the claim that they are conceptually very different with PCA being about "describing the data" and FA being about "finding latent variables". I just do not see this is as true [almost] at all. 
To comment on some specific points mentioned above and in the linked answers:


*

*"in PCA the number of dimensions to extract/retain is fundamentally subjective, while in EFA the number is fixed, and you usually have to check several solutions" -- well, the choice of the solution is still subjective, so I don't see any conceptual difference here. In both cases, $k$ is (subjectively or objectively) chosen to optimize the trade-off between model fit and model complexity.

*"FA is able to explain pairwise correlations (covariances). PCA generally cannot do it" -- not really, both of them explain correlations better and better as $k$ grows. 

*Sometimes extra confusion arises (but not in @ttnphns's answers!) due to the different practices in the disciplines using PCA and FA. For example, it is a common practice to rotate factors in FA to improve interpretability. This is rarely done after PCA, but in principle nothing is preventing it. So people often tend to think that FA gives you something "interpretable" and PCA does not, but this is often an illusion.
Finally, let me stress again that for very small $n$ the differences between PCA and FA can indeed be large, and maybe some of the claims in favour of FA are done with small $n$ in mind. As an extreme example, for $n=2$ a single factor can always perfectly explain the correlation, but one PC can fail to do it quite badly.

Update 1: generative models of the data
You can see from the number of comments that what I am saying is taken to be controversial. At the risk of flooding the comment section even further, here are some remarks regarding "models" (see comments by @ttnphns and @gung). @ttnphns does not like that I used the word "model" [of the covariance matrix] to refer to the approximations above; it is an issue of terminology, but what he calls "models" are probabilistic/generative models of the data:
\begin{align}
\mathrm{PPCA}: &\:\:\: \mathbf x = \mathbf W \mathbf z + \boldsymbol \mu + \boldsymbol \epsilon, \; \boldsymbol \epsilon \sim \mathcal N(0, \sigma^2 \mathbf I) \\
\mathrm{FA}: &\:\:\: \mathbf x = \mathbf W \mathbf z + \boldsymbol \mu + \boldsymbol \epsilon, \; \boldsymbol \epsilon \sim \mathcal N(0, \boldsymbol \Psi)
\end{align}
Note that PCA is not a probabilistic model, and cannot be formulated in this way.
The difference between PPCA and FA is in the noise term: PPCA assumes the same noise variance $\sigma^2$ for each variable, whereas FA assumes different variances $\Psi_{ii}$ ("uniquenesses"). This minor difference has important consequences. Both models can be fit with a general expectation-maximization algorithm. For FA no analytic solution is known, but for PPCA one can analytically derive the solution that EM will converge to (both $\sigma^2$ and $\mathbf W$). Turns out, $\mathbf W_\mathrm{PPCA}$ has columns in the same direction but with a smaller length than standard PCA loadings $\mathbf W_\mathrm{PCA}$ (I omit exact formulas). For that reason I think of PPCA as "almost" PCA: $\mathbf W$ in both cases span the same "principal subspace".
The proof (Tipping and Bishop 1999) is a bit technical; the intuitive reason for why homogeneous noise variance leads to a much simpler solution is that $\mathbf C - \sigma^2 \mathbf I$ has the same eigenvectors as $\mathbf C$ for any value of $\sigma^2$, but this is not true for $\mathbf C - \boldsymbol \Psi$.
So yes, @gung and @ttnphns are right in that FA is based on a generative model and PCA is not, but I think it is important to add that PPCA is also based on a generative model, but is "almost" equivalent to PCA. Then it ceases to seem such an important difference.

Update 2: how come PCA provides best approximation to the covariance matrix, when it is well-known to be looking for maximal variance?
PCA has two equivalent formulations: e.g. first PC is (a) the one maximizing the variance of the projection and (b) the one providing minimal reconstruction error. More abstractly, the equivalence between maximizing variance and minimizing reconstruction error can be seen using Eckart-Young theorem.
If $\mathbf X$ is the data matrix (with observations as rows, variables as columns, and columns are assumed to be centered) and its SVD decomposition is $\mathbf X=\mathbf U\mathbf S\mathbf V^\top$, then it is well known that columns of $\mathbf V$ are eigenvectors of the scatter matrix (or covariance matrix, if divided by the number of observations) $\mathbf C=\mathbf X^\top \mathbf X=\mathbf V\mathbf S^2\mathbf V^\top$ and so they are axes maximizing the variance (i.e. principal axes). But by the Eckart-Young theorem, first $k$ PCs provide the best rank-$k$ approximation to $\mathbf X$: $\mathbf X_k=\mathbf U_k\mathbf S_k \mathbf V^\top_k$ (this notation means taking only $k$ largest singular values/vectors) minimizes $\|\mathbf X-\mathbf X_k\|^2$.
The first $k$ PCs provide not only the best rank-$k$ approximation to $\mathbf X$, but also to the covariance matrix $\mathbf C$. Indeed, $\mathbf C=\mathbf X^\top \mathbf X=\mathbf V\mathbf S^2\mathbf V^\top$, and the last equation provides the SVD decomposition of $\mathbf C$ (because $\mathbf V$ is orthogonal and $\mathbf S^2$ is diagonal). So the Eckert-Young theorem tells us that the best rank-$k$ approximation to $\mathbf C$ is given by $\mathbf C_k = \mathbf V_k\mathbf S_k^2\mathbf V_k^\top$. This can be transformed by noticing that $\mathbf W = \mathbf V\mathbf S$ are PCA loadings, and so $$\mathbf C_k=\mathbf V_k\mathbf S_k^2\mathbf V^\top_k=(\mathbf V\mathbf S)_k(\mathbf V\mathbf S)_k^\top=\mathbf W_k\mathbf W^\top_k.$$
The bottom-line here is that 
$$ \mathrm{minimizing} \; 
\left\{\begin{array}{ll} 
\|\mathbf C-\mathbf W\mathbf W^\top\|^2 \\ \|\mathbf C-\mathbf W\mathbf W^\top-\sigma^2\mathbf I\|^2 \\ \|\mathbf C-\mathbf W\mathbf W^\top-\boldsymbol\Psi\|^2\end{array}\right\} \; 
\mathrm{leads \: to} \; 
\left\{\begin{array}{cc} \mathrm{PCA}\\ \mathrm{PPCA} \\ \mathrm{FA} \end{array}\right\} \;
\mathrm{loadings},$$
as stated in the beginning.

Update 3: numerical demonstration that PCA$\to$FA when $n \to \infty$
I was encouraged by @ttnphns to provide a numerical demonstration of my claim that as dimensionality grows, PCA solution approaches FA solution. Here it goes.
I generated a $200\times 200$ random correlation matrix with some strong off-diagonal correlations. I then took the upper-left $n \times n$ square block $\mathbf C$ of this matrix with $n=25, 50, \dots 200$ variables to investigate the effect of the dimensionality. For each $n$, I performed PCA and FA with number of components/factors $k=1\dots 5$, and for each $k$ I computed the off-diagonal reconstruction error $$\sum_{i\ne j}\left[\mathbf C - \mathbf W \mathbf W^\top\right]^2_{ij}$$ (note that on the diagonal, FA reconstructs $\mathbf C$ perfectly, due to the $\boldsymbol \Psi$ term, whereas PCA does not; but the diagonal is ignored here). Then for each $n$ and $k$, I computed the ratio of the PCA off-diagonal error to the FA off-diagonal error. This ratio has to be above $1$, because FA provides the best possible reconstruction.

On the right, different lines correspond to different values of $k$, and $n$ is shown on the horizontal axis. Note that as $n$ grows, ratios (for all $k$) approach $1$, meaning that PCA and FA yield approximately the same loadings, PCA$\approx$FA. With relatively small $n$, e.g. when $n=25$, PCA performs [expectedly] worse, but the difference is not that strong for small $k$, and even for $k=5$ the ratio is below $1.2$.
The ratio can become large when the number of factors $k$ becomes comparable with the number of variables $n$. In the example I gave above with $n=2$ and $k=1$, FA achieves $0$ reconstruction error, whereas PCA does not, i.e. the ratio would be infinite. But getting back to the original question, when $n=21$ and $k=3$, PCA will only moderately lose to FA in explaining the off-diagonal part of $\mathbf C$.
For an illustrated example of PCA and FA applied to a real dataset (wine dataset with $n=13$), see my answers here:


*

*What are the differences between Factor Analysis and Principal Component Analysis? 

*PCA and exploratory Factor Analysis on the same data set
A: In my view, the notions of "PCA" and "FA" are on a different dimension from that of notions of "exploratory", "confirmatory" or maybe "inferential".    So each of the two mathematical/statistical methods can be applied with one of the three approaches.    
For instance, why should it be unsensical to have a hypothese, that my data have a general factor and also the structure of a set of principal components (because my experiment with my electronical apparatus gave me nearly errorfree data) and I test my hypothese, that the eigenvalues of the subsequent factors occur with ratio of 75% ? This is then PCA in a confirmatory framework.      
On the other hand, it seems ridiculous that in our research team we create with much work an item battery for measuring violence between pupils and assuming 3 main behaves (physical agression, depression, search for help by authorities/parents) and putting the concerning questions in that battery ... and "exploratorily" work out how many factors we have... Instead to look, how well our scale contains three recognizable factors (besides neglectable itemspecific and possibly even spuriously correlated error). And after that, when I've confirmed, that indeed our item-battery serves the intention, we might test the hypothese, that in the classes of younger children the loadings on the factor indicating "searching-help-by-authorities" are higher than that of older pupils. Hmmm, again confirmatory...      
And exploratory? I have a set of measures taken from a research on microbiology from 1960 and they had not much theory but sampled everything they could manage because their field of research was just very young, and I re-explore the dominant factorstructure, assuming (for example), that all errors are of the same amount because of the optical precision of the microscope used (the ppca-ansatz as I have just learnt). Then I use the statistical (and subsequently the mathematical) model for the FA, but in this case in an explorative manner.       
This is it at least how I understand the terms.
Maybe I'm completely on the wrong track here, but I don't assume it.       

Ps. In the 90'ies I wrote a small interactive program to explore the method of PCA and factoranalysis down to the bottom. It was written in Turbo-Pascal, can still only be run in a Dos-Window ("Dos-box" under Win7) but has a really nice appeal: interactively switching factors to be included or not, then rotate, separate itemspecific error-variance (according to the SMC-criterion or the equal-variances-criterion (ppca?)), switch the Kaiser-option on and off, the use of the covariances on and off - just all while the factorloadingsmatrix is visible like in a spreadsheet and can be rotated for the basic different rotation-methods.
It is not highly sophisticated: no chisquare-test for instance, just intended for self-learning of the internal mathematical mechanics. It has also a "demo-mode", where the program runs itself, showing explanative comments on the screen and simulating the keyboard-inputs, which the user normally would do.
Whoever is interested to do selfstudy or teaching with it can download it from my small software-pages inside-(R).zip Just expand the files in the zip in a directory accessible by the Dos-Box and call "demoall.bat" In the third part of the "demoall" I've made a demonstration how to model itemspecific errors by rotations from an initially pca-solution... 
A: Just one additional remark for @amoebas's long (and really great) answer on the character of the $\Psi$-estimate.           
In your initial statements you have three $\Psi$: for PCA is $ \Psi = 0$, for PPCA is $ \Psi=\sigma ^2 I $ and for FA you left $\Psi$ indeterminate.       
But it should be mentioned, that there is an infinite number of various possible $\Psi$ (surely restricted) but exactly a single one which minimizes the rank of the factor matrix. Let's call this $\Psi_{opt}$ The standard (automatical) estimate for $\Psi_{std}$ is the diagonalmatrix based on the SMC's, so let's write this as $\Psi_{std}= \alpha^2 D_{smc}$ (and even some software (seem to) do not attempt to optimize $\alpha$ down from $1$ while $ \alpha \lt 1$ is (generally) required to prevent Heywood-cases/negative-definiteness). And moreover, even such optimized $\alpha^2$ would not guarantee minimal rank of the remaining covariances, thus usually we have this not equal: in general $\Psi_{std} \ne \Psi_{opt}$.
To really find $\Psi_{opt}$ is a very difficult game, and as far as I know (but that's no more so "far" as, say, 20 years ago when I was more involved and nearer to the books) this is still an unsolved problem.       

Well this reflects the ideal, mathematical side of the problem, and my distinction  between $\Psi_{std} $ and $\Psi_{opt}$ also might be actually small. A more general caveat is however, that it discusses the whole factorization machinery from the view that I study only my sample or have data of the whole population; in the model of inferential statistics, where I infer from an imperfect sample on the population, my empirical covariance- and thus also the factormatrix is only an estimate, it's only a shadow of the "true" covariance-/factormatrix. Thus in such a framework/model we should even consider that our "errors" are not ideal, and thus might be spuriously correlated. So in fact in such models we should/would leave the somehow idealistic assumption of uncorrelated error, and thus of a strictly diagonal form of $\Psi$, behind us.
