Multidimensional item response theory (MIRT) models have generated increasing interest in the psychometrics literature. Efficient approaches for estimating MIRT models with dichotomous responses have been developed, but constructing an equally efficient and robust algorithm for polytomous models has received limited attention. To address this gap, this paper presents a novel Gaussian variational estimation algorithm for the multidimensional generalized partial credit model (MGPCM). The proposed algorithm demonstrates both fast and accurate performance, as illustrated through a series of simulation studies and two real data analyses. [This is the online version of an article published in "Psychometrika."]