This paper deals with mmWave overloaded multiuser multi-input multi-output (MU-MIMO) detection, where the number of receive antennas is less than that of transmitted streams. Belief propagation (BP) is well known strategy for achieving large-scale MU detection (MUD) with low-complexity and high-accuracy. However, in mmWave massive MUD, the BPbased signal detector is subject to ill convergence behavior of iterative detection due to under-determined problem induced by spatial overloading and strong correlation among user channels induced by narrow angular spread of receive signal and lineof-sight (LOS) environments. To alleviate these impairments, we propose a novel iterative MUD approach based on beam-domain subspace marginalized BP (SMBP). Exploiting the approximate sparsity of beam-domain channels, the maximum likelihood (ML) principle is used to combine the strongly correlated signal subspace with reduced dimension while the BP-based detection is used for the remaining complementary subspace. The space partitioning criterion is adaptively determined based on channel state information (CSI) so that the two subspaces are as orthogonal as possible. Numerical results show that the proposed method is able to serve a massive number of wireless connections with low computational complexity even in the LOS environment, while providing excellent BER performance.