We perform a lattice QCD study of the local form factors governing the $B_{s}\to \mu^{+}\mu^{-}\gamma$ decay. To determine the $B_{s}$ meson form factors, we perform lattice simulations for several values of the heavy-strange meson masses $m_{H_{s}}$, within the range $m_{H_{s}} \in [ m_{D_{s}}, 2 m_{D_{s}} ]$, and extrapolate to the physical $B_{s}$ meson mass, $m_{B_{s}} \simeq 5.367~{\rm GeV}$, using heavy quark effective theory (HQET) scaling laws. For this calculation we employ the gauge configurations generated by the ETM Collaboration with $N_{f} = 2+1+1$ flavours of Wilson-Clover twisted-mass fermions at maximal twist. We explore the region of large di-muon invariant masses, $\sqrt{q^{2}} > 4.16\,{\rm GeV}$, and use our results to estimate the branching fraction for $B_{s}\to \mu^{+}\mu^{-}\gamma$, recently measured by LHCb in the region $\sqrt{q^{2}} > 4.9~{\rm GeV}$.