Covariance matrix estimation plays an important role in risk management, asset pricing, and portfolio allocation. This task becomes challenging when the dimensionality is comparable or much larger than the sample size. A widely used approach for reducing dimensionality is based on multi-factor models. Although it has been well studied and quite successful in many applications, the quality of the estimated covariance matrix is often degraded due to a nontrivial amount of missing data in the factor matrix for both technical and cost reasons. In this paper, we consider a new matrix completion paradigm using the factor models directly and apply the proximal gradient method for the recovery. When the scale of problem is huge, the computational time and complexity is so large that high-dimensional covariance matrix estimation is challenging. Therefore, the random projection method is used to accelerate the speed. Numerical experiments show that our proposed randomized algorithm is promising.