The paper presents an iterative scheme that incrementally minimizes the Kullback-Leibler divergence between joint posterior and nuisance parameter distribution, and directly minimizes the Kullback-Leibler divergence between joint posterior and parameter-of-interest distribution in alternating order. This so-called variational Bayes incremental expectation maximization (VB-IEM) algorithm has two virtues: first, it non-decreases an objective function at each iteration and secondly, it adopts the convergence properties of the exact EM algorithm if the parameter estimate is consistent and maximum a posteriori. When the VB-IEM algorithm is applied to joint multi-user decoding and multi-channel estimation in flat Rayleigh fading, Monte Carlo simulations show that the proposed scheme supports extremely high effective system loads.