In this paper, a generic extension of variational mode decomposition (VMD) algorithm for multivariate or multichannel data sets is presented. We first define a model for multivariate modulated oscillations that is based on the presence of a joint or common frequency component among all channels of input data. Using that model for multivariate oscillations, we construct a variational optimization problem that aims to extract an ensemble of band-limited modes containing inherent multivariate modulated oscillations present in multivariate input signal. The cost function to be minimized is the sum of bandwidths of all signal modes across all input data channels, which is a generic extension of the cost function used in standard VMD to multivariate data. Minimization of the resulting variational model is achieved through the alternate direction method of multipliers (ADMM) approach. That yields an optimal set of multivariate modes in terms of narrow bandwidth and corresponding center frequencies that are assumed to be commonly present among all channels of a multivariate modulated oscillation. We demonstrate the effectiveness of the proposed method through results obtained from extensive simulations involving test (synthetic) and real world multivariate data sets. Specifically, we focus on the ability of the proposed method to yield joint oscillatory modes in multivariate data which is a prerequisite in many real world applications involving nonstationary multivariate data. We also highlight the utility of the proposed method in two real world applications which include the separation of alpha rhythms in multivariate electroencephalogram (EEG) data and the decomposition of bivariate cardiotocographic signals that consist of fetal heart rate and maternal uterine contraction (FHR-UC) as its two channels.