This paper proposes a novel profile likelihood method for estimating the covariance parameters in exploratory factor analysis (EFA) with high-dimensional Gaussian data. By implementing a Lanczos algorithm and a limited-memory quasi-Newton method, we develop a matrix free algorithm (HDFA) which does partial singular value decomposition (partial SVD) for data matrix where number of observations $n$ is typically less than the dimension $p$ and it only requires limited amount of memory during likelihood maximization. We perform simulation study with both the randomly generated models and the data-driven models. Results indicate that HDFA substantially outperforms the EM algorithm in all cases. Furthermore, Our algorithm is applied to fit factor models for a fMRI dataset with suicidal attempters, suicidal nonattempters and a control group.