Vibration signals of rolling element bearings (REBs) contain substantial bearing motion state information. However, the property of nonlinear and nonstationary vibration signals decreases the diagnostic accuracy of REBs. To improve the accuracy of fault diagnosis for REBs, an ensemble approach based on ensemble empirical mode decomposition (EEMD), multi-scale permutation entropy (MPE), and backpropagation (BP) neural network optimized by genetic algorithm (GA) is proposed. Firstly, the REBs are decomposed into a set of intrinsic mode functions (IMFs) that contain various fault features by EEMD. The fault features of the first four IMFs are extracted by MPE, and the feature vectors are formed. Then, the BP neural network optimized by GA is utilized as a classifier for fault diagnosis to train and test the feature vector set, and the fault diagnosis of the REBs is realized in the form of probability output. Experimental results show that the proposed method can identify the fault pattern of the vibration signals of REBs precisely. Compared with the existing fault diagnosis methods, the proposed method can realize the fault diagnosis of REBs with 16 fault patterns, and demonstrates an excellent accuracy rate.