  Searching for gravitational waves in pulsar timing array data is
computationally intensive. The data is unevenly sampled, and the noise is
heteroscedastic, necessitating the use of a time-domain likelihood function
with attendant expensive matrix operations. The computational cost is
exacerbated when searching for individual supermassive black hole binaries,
which have a large parameter space due to the additional pulsar distance, phase
offset and noise model parameters needed for each pulsar. We introduce a new
formulation of the likelihood function which can be used to make the Bayesian
analysis significantly faster. We divide the parameters into projection and
shape parameters. We then accelerate the exploration of the projection
parameters by more than four orders of magnitude by precomputing the expensive
inner products for each set of shape parameters. The projection parameters
include nuisance parameters such as the gravitational wave phase offset at each
pulsar. In the new scheme, these troublesome nuisance parameters are
efficiently marginalized over using multiple-try Markov chain Monte Carlo
sampling as part of a Metropolis-within-Gibbs scheme. The acceleration provided
by our method will become increasingly important as pulsar timing datasets
rapidly grow. Our method also makes sophisticated analyses more tractable, such
as searches for multiple binaries, or binaries with non-negligible
eccentricities.
