Let's define the correlation matrix by $C$. Since it is positive semi-definite, but not positive definite, its spectral decomposition looks something like
$$C = Q D Q^{-1},$$
where the columns of $Q$ consist of orthonormal eigenvectors of $C$ and
$$D = \begin{pmatrix}\lambda_1 & 0 & \cdots & \cdots &\cdots & \cdots& 0\\ 0 & \lambda_2 & \ddots & && &\vdots \\ \vdots & \ddots &\ddots & \ddots && &\vdots \\ \vdots & &\ddots &\lambda_n &\ddots &&\vdots \\ \vdots & & & \ddots &0 & \ddots& \vdots \\ \vdots & & & &\ddots & \ddots& 0\\ 0 & \cdots &\cdots & \cdots &\cdots & 0& 0\end{pmatrix}$$
is a diagonal matrix containing the eigenvalues corresponding to the eigenvectors in $Q$. Some of those are $0$.
Moreover, $n$ is the rank of $C$.
A simple way to restore positive definiteness is setting the $0$-eigenvalues to some value that is numerically non-zero, e.g. $$\lambda_{n+1}, \lambda_{n+2},... = 10^{-15}.$$ Hence, set
$$\tilde{C} = Q \tilde{D} Q^{-1},$$
where
$$\tilde{D} = \begin{pmatrix}\lambda_1 & 0 & \cdots & \cdots &\cdots & \cdots& 0\\ 0 & \lambda_2 & \ddots & && &\vdots \\ \vdots & \ddots &\ddots & \ddots && &\vdots \\ \vdots & &\ddots &\lambda_n &\ddots &&\vdots \\ \vdots & & & \ddots &10^{-15} & \ddots& \vdots \\ \vdots & & & &\ddots & \ddots& 0\\ 0 & \cdots &\cdots & \cdots &\cdots & 0& 10^{-15}\end{pmatrix}$$
Then, perform the factor analysis for $\tilde{C}.
In Matlab, one can obtain $Q,D$ using the command:
[Q,D] = eig(C)
Constructing $\tilde{C}$ is then just simple Matrix manipulations.
Remark: It would be hard to tell how this influences the factor analysis though; hence, one should probably be careful with this method. Moreover, even though this is a $C$ is a correlation matrix, $\tilde{C}$ may well be not. Hence, another normalisation of the entries might be necessary.
