CERN-LHC, ATLAS. Searches for the exclusive decays of Higgs and $Z$ bosons into a vector quarkonium state and a photon are performed in the $\mu^+\mu^- \gamma$ final state with a proton--proton collision data sample corresponding to an integrated luminosity of $139$ fb$^{-1}$ collected at $\sqrt{s}=13$ TeV with the ATLAS detector at the CERN Large Hadron Collider. The observed data are compatible with the expected backgrounds. The 95% confidence-level upper limits on the branching fractions of the Higgs boson decays into $J/\psi \gamma$, $\psi(2S) \gamma$, and $\Upsilon(1S,2S,3S) \gamma$ are found to be $2.0\times10^{-4}$, $10.5\times10^{-4}$, and $(2.5,4.2,3.4)\times10^{-4}$, respectively, assuming Standard Model production of the Higgs boson. The corresponding 95% CL upper limits on the branching fractions of the $Z$ boson decays are $1.2\times10^{-6}$, $2.4\times10^{-6}$, and $(1.1,1.3,2.4)\times10^{-6}$. An observed 95% CL interval of $(-133,175)$ is obtained for the $\kappa_c/\kappa_\gamma$ ratio of Higgs boson coupling modifiers, and a 95% CL interval of $(-37,40)$ is obtained for $\kappa_b/\kappa_\gamma$.