Analog network coding is a promising technique which can greatly improve the transmission efficiency of wireless communications. In two-way relay systems with multiple subchannels, multiple user pairs and multiple relays, however, the optimal joint relay selection and subchannel allocation problem has not been studied in a systematic way. In this paper, a random combinatorial optimization approach, referred to as the weighted random bipartite graph (WRBG) based minimum weighted matching (MWM) method, will be proposed to solve this problem. By analyzing the properties of the MWM on WRBG, we shall derive the outage probability and diversity-multiplexing tradeoff of each user after relay selection and channel allocation. Theoretical results will demonstrate that the outage probability, cooperative diversity, and frequency diversity of the proposed WRBG based MWM method for multi-user two-way relay systems is the same as that of two-way relay systems with only one user pair. The proposed algorithm for MWM also enjoys a low computation complexity of O(log 2 N) for parallel implementations, where N is the number of subchannels. Simulation results will illustrate the potential of the proposed method and also verify the theoretical derivations.