Since its introduction in the late 19th century, symmetry breaking has been found to play a crucial role in physics. In particular, it appears as one key phenomenon controlling hydrodynamic and acoustic instabilities in problems with rotational symmetries. A previous paper investigated its desired potential application to the control of circumferential thermoacoustic modes in one annular cavity coupled with multiple flames (Bauerheim et al., J. Fluid Mech., vol. 760, 2014, pp. 431–465). The present paper focuses on a similar problem when symmetry breaking appears unintentionally, for example when uncertainties due to tolerances are taken into account. It yields a large uncertainty quantification (UQ) problem containing numerous uncertain parameters. To tackle this well-known ‘curse of dimensionality’, a novel UQ methodology is used. It relies on the active subspace approach to construct a reduced set of input variables. This strategy is applied on two annular cavities coupled by 19 flames to determine its modal risk factor, i.e. the probability of an azimuthal acoustic mode being unstable. Since each flame is modelled by two uncertain parameters, it leads to a large UQ problem involving 38 parameters. An acoustic network model is then derived, which yields a nonlinear dispersion relation for azimuthal modes. This nonlinear problem, subject to bifurcations, is solved quasi-analytically. Results show that the dimension of the probabilistic problem can be drastically reduced, from 38 uncertain parameters to only 3. Moreover, it is found that the three active variables are related to physical quantities, which unveils underlying phenomena controlling the stability of the two coupled cavities. The first active variable is associated with a coupling strength controlling the bifurcation of the system, while the two others correspond to a symmetry-breaking effect induced by the uncertainties. Thus, an additional destabilization effect appear caused by the non-uniform pattern of the uncertainty distribution, which breaks the initial rotating symmetry of the annular cavities. Finally, the active subspace is exploited by fitting the response surface with polynomials (linear, quadratic and cubic). By comparing accuracy and cost, results prove that 5 % error can be achieved with only 30 simulations on the reduced space, whereas 2000 are required on the complete initial space. It exemplifies that this novel UQ technique can accurately predict the risk factor of an annular configuration at low cost as well as unveil key parameters controlling the stability.