We present a detailed theoretical analysis of very rare, exclusive hadronic decays of the electroweak gauge bosons V = W, Z from first principles of QCD. Our main focus is on the radiative decays V → Mγ , in which M is a pseudoscalar or vector meson. At leading order in an expansion in powers of Λ QCD /m V the decay amplitudes can be factorized into convolutions of calculable hard-scattering coefficients with the leading-twist light-cone distribution amplitude of the meson M . Power corrections to the decay rates arise first at order (Λ QCD /m V ) 2 . They can be estimated in terms of higher-twist distribution amplitudes and are predicted to be tiny. We include one-loop O α s $$ \mathcal{O}\left({\alpha}_s\right) $$ radiative corrections to the hard-scattering coefficients and perform the resummation of large logarithms ( α s  ln( m v 2 / μ 0 2 )) n (with μ 0 ∼ 1 GeV a typical hadronic scale) to all orders in perturbation theory. Evolution effects have an important impact both numerically and conceptually, since they reduce the sensitivity to poorly determined hadronic parameters. We present detailed numerical predictions and error estimates, which can serve as benchmarks for future precision measurements. We also present an exploratory study of the weak radiative decays Z → MW . Some of the decay modes studied here have branching ratios large enough to be accessible in the high-luminosity run of the LHC. Many of them can be measured with high accuracy at a future lepton collider. This will provide stringent tests of the QCD factorization formalism and enable novel searches for new physics.