We propose a new, likelihood-free approach to inferring the primordial matter
power spectrum and cosmological parameters from arbitrarily complex forward
models of galaxy surveys where all relevant statistics can be determined from
numerical simulations, i.e. black-boxes. Our approach builds upon approximate
Bayesian computation using a novel effective likelihood, and upon the
linearisation of black-box models around an expansion point. Consequently, we
obtain simple "filter equations" for an effective posterior of the primordial
power spectrum, and a straightforward scheme for cosmological parameter
inference. We demonstrate that the workload is computationally tractable, fixed
a priori, and perfectly parallel. As a proof of concept, we apply our framework
to a realistic synthetic galaxy survey, with a data model accounting for
physical structure formation and incomplete and noisy galaxy observations. In
doing so, we show that the use of non-linear numerical models allows the galaxy
power spectrum to be safely fitted up to at least $k_\mathrm{max} = 0.5$
$h$/Mpc, outperforming state-of-the-art backward-modelling techniques by a
factor of $\sim 5$ in the number of modes used. The result is an unbiased
inference of the primordial matter power spectrum across the entire range of
scales considered, including a high-fidelity reconstruction of baryon acoustic
oscillations. It translates into an unbiased and robust inference of
cosmological parameters. Our results pave the path towards easy applications of
likelihood-free simulation-based inference in cosmology.