Strong lensing images provide a wealth of information about both the magnified source and the mass distribution in the lens, allowing dark matter models to be constrained. However, due to the degeneracies inherent to lensing, making inferences about substructure requires very accurate and precise, yet flexible, reconstruction of the source. In anticipation of future high-resolution datasets, in this work we leverage a range of recent developments in machine learning to present a new end-to-end differentiable GPU-accelerated Bayesian strong lensing image analysis pipeline. We have also developed a new statistically principled source model based on an efficient approximation to Gaussian processes that also takes into account pixelisation effects. Using variational inference and stochastic gradient descent, we simultaneously derive approximate posteriors for tens of thousands of lens and source parameters, while also optimising hyperparameters. Besides efficient and accurate parameter estimation and uncertainty quantification, the main aim of the pipeline is the generation of training data for targeted simulation-based inference of dark matter substructure.