Matrix factor models have been growing popular dimension reduction tools for large-dimensional matrix time series. However, the heteroscedasticity of the idiosyncratic components has barely received any attention. Starting from the pseudo likelihood function, this paper introduces a Generalized Principal Component Analysis (GPCA) method for matrix factor model which takes the heteroscedasticity into account. Theoretically, we first derive the asymptotic distributions of the GPCA estimators by assuming the separable covariance matrices are known in advance. We then propose adaptive thresholding estimators for the separable covariance matrices and derive their convergence rates, which is of independent interest. We also show that this would not alter the asymptotic distributions of the GPCA estimators under certain regular sparsity conditions in the high-dimensional covariance matrix estimation literature. The GPCA estimators are shown to be more efficient than the state-of-the-art methods under certain heteroscedasticity conditions. Thorough numerical studies are conducted to demonstrate the superiority of our method over the existing approaches. Analysis of a financial portfolio dataset illustrates the empirical usefulness of the proposed method.