This paper provides a comprehensive estimation framework via nuclear norm plus $l_1$ norm penalization for high-dimensional approximate factor models with a sparse residual covariance. The underlying assumptions allow for non-pervasive latent eigenvalues and a prominent residual covariance pattern. In that context, existing approaches based on principal components may lead to misestimate the latent rank, due to the numerical instability of sample eigenvalues. On the contrary, the proposed optimization problem retrieves the latent covariance structure and exactly recovers the latent rank and the residual sparsity pattern. Conditioning on them, the asymptotic rates of the subsequent ordinary least squares estimates of loadings and factor scores are provided, the recovered latent eigenvalues are shown to be maximally concentrated and the estimates of factor scores via Bartlett's and Thompson's methods are proved to be the most precise given the data. The validity of outlined results is highlighted in an exhaustive simulation study and in a real financial data example.