  We present a new method to compute short-time expectation values in large
collective spin systems with generic Markovian decoherence. Our method is based
on a Taylor expansion of a formal solution to the equations of motion for
Heisenberg operators. This expansion can be truncated at finite order to obtain
virtually exact results at short times that are relevant for metrological
applications such as spin squeezing. In order to evaluate the expansion for
Heisenberg operators, we compute the relevant structure constants of a
collective spin operator algebra. We demonstrate the utility of our method by
computing spin squeezing, two-time correlation functions, and
out-of-time-ordered correlators for $10^4$ spins in strong-decoherence regimes
that are otherwise inaccessible via existing numerical methods. Our method can
be straightforwardly generalized to the case of a collective spin coupled to
bosonic modes, relevant for trapped ion and cavity QED experiments, and may be
used to investigate short-time signatures of quantum chaos and information
scrambling.
