We study the problem of high-dimensional covariance matrix estimation from partial observations. We consider covariance matrices modeled as Kronecker products of matrix factors, and rely on observations with missing values. In the absence of missing data, observation vectors are assumed to be i.i.d multivariate Gaussian. In particular, we propose a new procedure computationally affordable in high dimension to extend an existing permuted rank-penalized least-squares method to the case of missing data. Our approach is applicable to a large variety of missing data mechanisms, whether the process generating missing values is random or not, and does not require imputation techniques. We introduce a novel unbiased estimator and characterize its convergence rate to the true covariance matrix measured by the spectral norm of a permutation operator. We establish a tight outer bound on the square error of our estimate, and elucidate consequences of missing values on the estimation performance. Different schemes are compared by numerical simulations in order to test our proposed estimator.