The analysis of optical images of galaxy-galaxy strong gravitational lensing systems can provide important information about the distribution of dark matter at small scales. However, the modeling and statistical analysis of these images is extraordinarily complex, bringing together source image and main lens reconstruction, hyper-parameter optimization, and the marginalization over small-scale structure realizations. We present here a new analysis pipeline that tackles these diverse challenges by bringing together many recent machine learning developments in one coherent approach, including variational inference, Gaussian processes, differentiable probabilistic programming, and neural likelihood-to-evidence ratio estimation. Our pipeline enables: (a) fast reconstruction of the source image and lens mass distribution, (b) variational estimation of uncertainties, (c) efficient optimization of source regularization and other hyperparameters, and (d) marginalization over stochastic model components like the distribution of substructure. We present here preliminary results that demonstrate the validity of our approach.