Robustly inferring the properties of a dark matter subhalo in a strong lensing image requires marginalizing over uncertainties in the main lens mass and source light distributions. This is an extremely difficult problem due to the high dimensionality of lensing observations. Here we present a new multi-stage method for performing subhalo inference that combines the strengths of parametric lensing models and simulation-based inference (SBI), a tool that leverages neural networks to directly obtain marginal posteriors from observations. In the first stage, we use our novel Gaussian process-inspired lensing model to closely fit an observation, obtaining approximate posteriors for all model parameters. We use the obtained posteriors to generate variations of the target image with different subhalo realizations. In the second stage we train a targeted inference network on these images to produce precision posteriors for the subhalo's parameters using the technique of neural likelihood-to-evidence estimation. The final inference is performed by applying the trained network to the original observation of interest. We present results of a mock analysis showing we can accurately reconstruct a subhalo's position and mass in a realistic, high-resolution observation, marginalizing over more than 100,000 lens and source parameters. The whole analysis can be performed rapidly using a single graphical processing unit and scales rather favorably (linearly) with models complexity.