  In Einstein's general relativity, gravity is mediated by a massless metric
field. The extension of general relativity to consistently include a mass for
the graviton has profound implications for gravitation and cosmology. Salient
features of various massive gravity theories can be captured by Galileon
models, the simplest of which is the cubic Galileon. The presence of the
Galileon field leads to additional gravitational radiation in binary pulsars
where the Vainshtein mechanism is less suppressed than its fifth-force
counterpart, which deserves a detailed confrontation with observations. We
prudently choose fourteen well-timed binary pulsars, and from their intrinsic
orbital decay rates we put a new bound on the graviton mass, $m_g \lesssim 2
\times 10^{-28}\,{\rm eV}/c^2$ at the 95% confidence level, assuming a flat
prior on $\ln m_g$. It is equivalent to a bound on the graviton Compton
wavelength $\lambda_g \gtrsim 7 \times 10^{21}\,{\rm m}$. Furthermore, we
extensively simulate times of arrival for pulsars in orbit around stellar-mass
black holes and the supermassive black hole at the Galactic center, and
investigate their prospects in probing the cubic Galileon theory in the near
future.
