What do the first $k$ factors from factor analysis maximize? In principal components analysis, the first $k$ principal components are the $k$ orthogonal directions with the maximum variance. In other words, the first principal component is chosen to be the direction of maximum variance, the second principal component is chosen to be the direction orthogonal to the first with the maximum variance, and so on.
Is there a similar interpretation for Factor Analysis? For example, I'm thinking that the first $k$ factors are the factors that best explain the off-diagonal components of the original correlation matrix (in the sense of, say, squared error between the original correlation matrix and the correlation matrix defined by the factors). Is this true (or is there something similar that we can say)?
 A: PCA is primarily a data reduction technique where the objective is to obtain a projection of data onto a lower dimensional space. Two equivalent objectives are to either iteratively maximize variance or to minimize reconstruction error. This is actually worked out in some details in the answers to this previous question.
In contrast, factor analysis is primarily a generative model of a $p$-dimensional data vector $X$ saying that 
$$X = AS + \epsilon$$ 
where $S$ is the $q$ dimensional vector of latent factors, $A$ is $p \times k$ with $k < p$ and $\epsilon$ is a vector of uncorrelated errors. The $A$ matrix is the matrix of factor loadings. This yields a special parametrization of the covariance matrix as
$$\Sigma = AA^T + D$$
The problem with this model is that it is overparametrized. The same model is obtained if $A$ is replaced by $AR$ for any $k \times k$ orthogonal matrix $R$, which means that the factors themselves are not unique. Various suggestions exist for solving this problem, but there is not a single solution that gives you factors with the kind of interpretation you ask for. One popular choice is the varimax rotation. However, the criterion used only determines the rotation. The column space spanned by $A$ does not change, and since this is part of the parametrization, it is determined by whatever method is used to estimate $\Sigma$ - by maximum likelihood in a Gaussian model, say. 
Hence, to answer the question, the chosen factors are not given automatically from using a factor analysis model, so there is no single interpretation of the $k$ first factors. You have to specify the method used to estimate (the column space of) $A$ and the method used to choose the rotation. If $D = \sigma^2 I$ (all errors have same variance) the MLE solution for the column space of $A$ is the space spanned by the leading $q$ principal component vectors, which may be found by a singular value decomposition. It is, of course, possible to choose not to rotate and report these principal component vectors as the factors.  
Edit: To emphasize how I see it, the factor analysis model is a model of the covariance matrix as a rank $k$ matrix plus a diagonal matrix. Thus the objective with the model is to best explain the covariance with such a structure on the covariance matrix. The interpretation is that such a structure on the covariance matrix is compatible with an unobserved $k$ dimensional factor. Unfortunately, the factors can not be recovered uniquely, and how they may be chosen within the set of possible factors does not relate in any way to the explanation of data. As is the case with PCA, one can standardize the data upfront and thus fit a model that attempts to explain the correlation matrix as a rank $k$ plus a diagonal matrix.   
A: @RAEGTIN, I believe that you think right. After extraction and prior rotation, each successive factor does account for less and less of covariation/correlation, just like each successive component accounts for less and less of variance: in both cases, columns of a loading matrix A go in the order of fall of sum of squared elements (loadings) in them. Loading is correlation bw factor and variable; therefore one may say that the 1st factor explains the greatest portion of "overall" squared r in R matrix, the 2nd factor is second here, etc. The difference between FA and PCA, though, in predicting correlations by loadings is as follows: FA is "calibrated" to restore R quite finely with just m extracted factors (m factors < p variables), while PCA is rude in restoring it by m components, - it needs all p components to restore R without error.
P.S. Just to add. In FA, a loading value "consists" of clean communality (a part of variance responsible for correlating) while in PCA a loading is a mixture of communality and uniqness of the variable and therefore grabs variability.
