2022MNRAS.517.3275E

The distribution of dark and luminous matter can be mapped around galaxies that gravitationally lens background objects into arcs or Einstein rings. New surveys will soon observe hundreds of thousands of galaxy lenses and current labour-intensive analysis methods will not scale up to this challenge. We develop an automatic Bayesian method, which we use to fit a sample of 59 lenses imaged by the Hubble Space Telescope. We set out to leave no lens behind and focus on ways in which automated fits fail in a small handful of lenses, describing adjustments to the pipeline that ultimately allows us to infer accurate lens models for all 59 lenses. A high-success rate is key to avoid catastrophic outliers that would bias large samples with small statistical errors. We establish the two most difficult steps to be subtracting foreground lens light and initializing a first approximate lens model. After that, increasing model complexity is straightforward. We put forward a likelihood cap method to avoid the underestimation of errors due to pixel discretization noise inherent to pixel-based methods. With this new approach to error estimation, we find a mean ∼1 per cent fractional uncertainty on the Einstein radius measurement, which does not degrade with redshift up to at least z = 0.7. This is in stark contrast to measurables from other techniques, like stellar dynamics and demonstrates the power of lensing for studies of galaxy evolution. Our PyAutoLens software is open source, and is installed in the Science Data Centres of the ESA Euclid mission.