In this study, we introduces a novel approach, called  the quasi maximum likelihood estimation (Q-MLE), for estimating large-dimensional matrix factor models. In contrast to the principal components based approach, Q-MLE takes into account heteroskedasticities of the diosyncratic error term, which are simultaneously estimated with other parameters. Theoretical analysis shows that the Q-MLE estimator of the factor loading matrices achieves faster convergence rates than  most existing estimators under similar conditions. We also present the asymptotic distributions of the Q-MLE estimators. Extensive numerical experiments demonstrate that the Q-MLE method performs better empirically, especially when heteroscedasticity exists. Furthermore, two real examples in finance and macroeconomics reveal factor patterns across rows and columns, which coincide with financial, economic, or geographical interpretations.