In this paper, we present a variational inference algorithm that decomposes a signal into multiple groups of related spectral lines. The spectral lines in each group are associated with a group parameter common to all spectral lines within the group. The proposed algorithm jointly estimates the group parameters, the number of spetral lines within a group, and the number of groups exploiting a Bernoulli-Gamma-Gaussian hierarchical prior model which promotes sparse solutions. Aiming to maximize the evidence lower bound (ELBO), variational inference provides analytic approximations of the posterior probability density functions (PDFs) and also gives estimates of the additional model parameters such as the measurement noise variance. While the activation variables of the groups and the associated group parameters (such as fundamental frequencies and the corresponding higher order harmonics) are estimated as point estimates, the remaining parameters such as the complex amplitudes of the spectral lines and their precision parameters are estimated as approximate posterior PDFs. We demonstrate the versatility and performance of the proposed algorithm on three different inference problems. In particular, the proposed algorithm is applied to the multi-pitch estimation problem, the radar signal-based extended object estimation problem, and variational mode decomposition (VMD) using synthetic measurements and to real multi-pitch estimation problem using the Bach-10 dataset. The results show that the proposed algorithm outperforms state-of-the-art model-based and pre-trained algorithms on all three inference problems.