Massive multiuser multiple-input multiple-output (massive MU-MIMO) technology is considered as a promising enabler to fulfill the rapid growth of traffic requirement for wireless mobile communications. The massive MU-MIMO system can achieve unlimited capacity when the base station (B-S) has accurate channel state information (CSI). In time-division-duplex (TDD) mode, the BS estimates CSI by receiving pilot signals sent from user terminals (UEs). However, because of using non-orthogonal pilots, pilot contamination happens to degrade the quality of the CSI estimation. To deal with pilot contamination problem, a low-complexity subspace minimum mean square error (MMSE) estimation method is proposed in this paper. Specifically, our approach operates the MMSE estimation in a low-dimensional subspace to avoid large matrix manipulation. Meanwhile, subspace projection helps to discriminate the desired signal and interfering signals in the power domain. Interference analysis shows the MMSE estimation can achieve interference-free estimation even in a low-dimensional subspace with a large number of BS antennas, and non-overlapping angles of arrival (AoAs) between desired and interfering UEs. Furthermore, thanks to the low-rank property of the channel covariance matrix in massive MU-MIMO systems, a two-step covariance matrix subspace projection method is proposed for further computational complexity reduction. The complexity analysis and simulation results indicate that our proposed approach has better channel estimation accuracy with lower complexity than the conventional MMSE estimation when the number of BS antennas is large. INDEX TERMS Massive MU-MIMO, channel estimation, pilot contamination, MMSE, subspace projection.