  Computing the matter power spectrum, $P(k)$, as a function of cosmological
parameters can be prohibitively slow in cosmological analyses, hence emulating
this calculation is desirable. Previous analytic approximations are
insufficiently accurate for modern applications, so black-box, uninterpretable
emulators are often used. We utilise an efficient genetic programming based
symbolic regression framework to explore the space of potential mathematical
expressions which can approximate the power spectrum and $\sigma_8$. We learn
the ratio between an existing low-accuracy fitting function for $P(k)$ and that
obtained by solving the Boltzmann equations and thus still incorporate the
physics which motivated this earlier approximation. We obtain an analytic
approximation to the linear power spectrum with a root mean squared fractional
error of 0.2% between $k = 9\times10^{-3} - 9 \, h{\rm \, Mpc^{-1}}$ and across
a wide range of cosmological parameters, and we provide physical
interpretations for various terms in the expression. We also provide a simple
analytic approximation for $\sigma_8$ with a similar accuracy, with a root mean
squared fractional error of just 0.4% when evaluated across the same range of
cosmologies. This function is easily invertible to obtain $A_{\rm s}$ as a
function of $\sigma_8$ and the other cosmological parameters, if preferred. It
is possible to obtain symbolic approximations to a seemingly complex function
at a precision required for current and future cosmological analyses without
resorting to deep-learning techniques, thus avoiding their black-box nature and
large number of parameters. Our emulator will be usable long after the codes on
which numerical approximations are built become outdated.
