The temperature dependence of thermal expansion of hexagonal close-packed Mg was analyzed based on the quasi-harmonic approximation (QHA) using first-principles density functional theory calculations. To consider the structural anisotropy of a Mg single crystal, we introduced two individual structural parameters into the QHA scheme so that static total energy and lattice vibration frequencies of the system were numerically described by the approximate polynomials as a function of the lattice parameters a and c. We found that our approach can successfully reproduce the thermal expansion behavior of Mg over the wide temperature range by adopting the second- and higher-order polynomial to describe the lattice vibration, in a manner consistent with the experimental measurements. The nonlinear dependence of the lattice vibration frequencies on the axial strain was suggested to play an important role in understanding the thermal expansion anisotropy due to the anharmonicity in the interatomic interactions.