Multi-scale fuzzy entropy can well measure the complexity of the vibration signal, but it lacks the effective use of other channel information. To make full use of the vibration information of other channels, the multivariate sample entropy theory that characterizes the multivariate complexity of synchronized multi-channel data is applied to the bearing fault diagnosis. To accurately extract fault features of bearing signals, a bearing multi-fault diagnosis method based on complete ensemble empirical mode decomposition with adaptive noise（CEEMDAN） and refined composite generalized multivariate multiscale fuzzy entropy（RCGmvMFE） is proposed. First, CEEMDAN is used to decompose multi-channel raw signals to obtain IMF without mode mixing. Then the correlation analysis method is used to screen the IMF components, and the IMF sensitive to the fault characteristics is selected as the multi-channel data to constitute the multivariate variable, and the RCGmvMFE is calculated to constitute the fault feature. Then, t-distributed stochastic neighbor embedding（t-SNE） is used to reduce the dimensionality of high-dimensional features. Finally, the whale optimization algorithm（WOA）is used to optimize the kernel extreme learning machine（WOA-KELM） so as to classify the low-dimensional fault features. Experimental results show that this method can effectively diagnose different fault severity of bearings, and provides a supplementary method for fault diagnosis of rolling bearings.