Orthogonality of the basis in factor analysis I have been studying principal component analysis (PCA) and then I have gone up to factor analysis (FA). I understood that PCA seeks orthonormal basis, but I am not so sure if this is the case for factor analysis. It seems that if my code is correct, the basis it finds is not orthonormal.
Is it true that the basis is not orthonormal in factor analysis? What is the theory behind it?
 A: You seem to be familiar with the probabilistic PCA, so I will use it in my explanation.
In probabilistic PCA, the model of the data is $$\mathbf x|\mathbf z \sim \mathcal N(\mathbf W \mathbf z + \boldsymbol \mu, \sigma^2 \mathbf I),\hspace{5em}\mathbf z\sim \mathcal N(\mathbf 0, \mathbf I),$$ where $\mathbf z$ is lower-dimensional than $\mathbf x$.
Maximum likelihood solution is not unique, but one of these solutions is "special" and has an analytical expression in terms of standard PCA: columns of this $\mathbf W_\mathrm{PPCA}$ are proportional to the principal directions of the data matrix $\mathbf X$. In fact, they are principal directions scaled by the corresponding eigenvalues (these are PCA loadings), and then scaled  a bit down. We can call them PPCA loadings.
If you compute them like that, then they are of course orthogonal to each other (but do not have unit length, so not orthonormal). But note that $\mathbf W$ can be multiplied by any rotation matrix and will be an equally good (equally likely) solution, and if do that then its columns will stop being orthogonal. So if, instead of taking PCA loadings and converting them into PPCA loadings by the analytical formula, you were to use expectation-minimization (EM) algorithm to find optimal $\mathbf W$ then the algorithm would converge to some random solution and it will not necessarily have orthogonal columns.
Factor analysis has a very much related model $$\mathbf x|\mathbf z \sim \mathcal N(\mathbf W \mathbf z + \boldsymbol \mu, \boldsymbol \Psi),\hspace{5em}\mathbf z\sim \mathcal N(\mathbf 0, \mathbf I),$$ where $\boldsymbol \Psi$ is a diagonal matrix. In the comments you ask if PPCA is a special case of FA; no, I would rather say that it's a "restricted" FA. In any case, there is no analytical solution for FA. Using EM algorithm, you can find  maximum likelihood $\mathbf W_\mathrm{FA}$ (which is again not unique). Will it have orthogonal columns? No it will not, see above. It will also not have columns ordered by variance.
However, we can rotate this solution by doing SVD on $\mathbf W_\mathrm{FA} = \mathbf {USV}^\top$ and then taking $\mathbf {US}$ as our new $\mathbf W_\mathrm{FA}$. This new matrix will have columns ordered by factor variance and can be called the matrix of FA loadings. It is easy to see that its columns will also be orthogonal.
So there is no difference between FA and PPCA in this respect.
