This paper considers a multi-user multiple-input multiple-output (MU-MIMO) system with independent phase noise (PN) at user terminals (UTs) due to non-synchronous noisy local oscillators. The conventional factor graph of MU-MIMO involves an observation factor with multiple high-dimensional variables in the form of multiplication and summation. In this work, the MU-MIMO factor graph is represented by introducing auxiliary variables, which enables the decomposition of the observation factor into several simpler sub-factors and then use different message passing techniques to obtain approximate marginals. Specifically, belief propagation and expectation propagation (BP-EP) are used for the linear factors for random walk process, modulation and coding, while mean field (MF) is used for non-linear factor in exponential form. To reduce the computational complexity, the non-Gaussian MF messages are approximated to be Gaussian by using the second order Taylor expansion at its belief obtained i