(1) As we want to have clusters of Variables which are correlated with
each other in one and the same factor, where exactly do we optimize
for this cohesity?
Factor analysis is not cluster analysis, it is latent variable modeling. It could be seen as a "clustering" method for correlated variables, however, the different "clusters" may show up only after rotation (such as varimax) to a simple structure where different subsets of items appear loaded by different factors. Yet, reconstruction of correlations holds even without any rotation. (See: A point in "Contrast list PCA vs FA" here; the starting section here.)
In other words, how do the loadings of the variables ... reflect the
proximity (i.e. their correlations) to each other...?
Reconstruction (= explaining) of correlations/covariances is the primary aim of FA. The reconstruction is formalized in fundamental factor theorem (see it in footnote 2 here). In short, correlation (or covariance) between two variables X and Y is restored, with (ideally small) residual, by those variables' m factor loadings: $r_{xy} \approx a_{x1}a_{y1}+a_{x2}a_{y2}+...+a_{xm}a_{ym}$.
(2) What are the Factors, numerically, when we try to optimize the
model (i.e. the equations above for minimizing)? I assume that if we
want to optimize for the loadings, we must already know the values for
the Fs – or not? I am thinking for instance about regression, where we
have the dependent and the independent variables of several
datasamples, and we try to optimize the regression line.
FA model is like a regressional model, still it is not regression analysis. In regression, the aim is to minimize the error. In FA, the main task is to restore correlations. Maximize variance explained by common factors (= hence minimize error or "unique variance") is second, not as principal, aim. Different extraction methods may put different accents in such two-fold (even 2+ fold) task.
We do the FA optimization on the level of loadings, the parameters a, not factor values. We do everything basic by analyzing the correlation/covariance matrix, not casewise dataset. Optimization tries to minimize error of the factor theorem (i.e. restore the matrix), not the "error" of the factor model. So no, we don't need to know F's values to optimize what we want and we don't estimate those values at the FA factor extraction phase optimization.
Btw, in regression we technically estimate its parameters b from correlations, too, without assistance from casewise data.
In regression (its simple forms), the parameters relate to the data directly because no latent feature modeling is involed. In FA, we may estimate "true" factor loadings very adequately sometimes, but the values of those factors remain black box for us; it is known as "factor values indeterminacy problem": computed factor scores are never exact true values (read, read) unlike in PCA (read answer and comments).
(3) When we optimize the model above, are there other criteria then
just reducing the distance between the estimated correlation matrix
and the original correlation matrix? Especially I always read that the
factors should each capture as much overall variance as possible – is
this just a byproduct or an actual criteria of the optimization
process?
This point was answered in my previous section at its beginning. It is a by-product. But often very desirable. In optimizations, we often can settle more than one objective function at once. (In this answer, I've tried to show visually on the 2nd picture how you can perform a factor extraction, and spoke of two objectives of "pulling the arrow by finger": one was to orthogonalize common factor to unique factors - and by that virtue to approximate correlation, and the other was to maximize the factor's variance. And here's a numerical example showing FA and PCA extractions step-by-step.)
(4) Determination of the Number of Factors: I know that after we do
have [have extracted] some number of Factors, we can use ... criteria
to get rid of those who do not capture enough variance. But what
determines how many factors we have before [the extraction]? Is this
initial number of factors determined by the dimensions of the
correlation matrix we start with? So we can say how many factors are
initially there even without having yet determined the loadings?
FA is fitting by modelling. This means we first select the number of factors m to extract - then extract them (and see the fit quality). You are not correct: we do not throw away some factors after extraction. If we don't like this or that factor (poor loadings, poor interpretation, etc) we usually reconsider m or revise variables to include in the analysis. Or simply ignore the factor, saying "this factor appears important for the overall fit but we have difficulty to describe/appreciate it". We don't "retain" factors in FA - unlike we do components in PCA; we might, however, in rare cases decide to rotate/interpret only a subset of the extracted factors.
To cite myself: After rotation, you may not say "this factor is more important than that" because they were rotated vis-a-vis each other (to be honest, in FA, unlike PCA, you may hardly say it even after the extraction because factors are modelled as already "important").
How do we determine m to extract? For this, please read some text on FA (including mayby this site's Q/A) - because the question is traditional/sacramental. Please find to read on 1) Kaiser's Eigenvalue >threshold criterion, 2) Cattell's Scree-plot elbow rule, 3) Very simple structure (VSS) and Minimum average partial (MAP) criterions, 4) Humphrey-Ilgen Parallel analysis, 5) post-extraction fit measures, such as size/pattern of correlation residuals, Chi-square goodness-of-fit, loadings st. errors, etc, 6) old good Percentage of variance explained, 7) Interpretability (sure!). And a number of other approaches.
