ABSTRACT We present a machine-learning-based model for the total density profiles of subhaloes with masses $M \gtrsim 7\times 10^8\, h^{-1}{\rm M}_\odot$ in the IllustrisTNG100 simulation. The model is based on an interpretable variational encoder (IVE) which returns the independent factors of variation in the density profiles within a low-dimensional representation, as well as the predictions for the density profiles themselves. The IVE returns accurate and unbiased predictions on all radial ranges, including the outer region profile where the subhaloes experience tidal stripping; here its fit accuracy exceeds that of the commonly used Einasto profile. The IVE discovers three independent degrees of freedom in the profiles, which can be interpreted in terms of the formation history of the subhaloes. In addition to the two parameters controlling the normalization and inner shape of the profile, the IVE discovers a third parameter that accounts for the impact of tidal stripping on to the subhalo outer profile; this parameter is sensitive to the mass loss experienced by the subhalo after its infall on to its parent halo. Baryonic physics in the IllustrisTNG galaxy formation model does not impact the number of degrees of freedom identified in the profile compared to the pure dark matter expectations, nor their physical interpretation. Our newly proposed profile fit can be used in strong lensing analyses or other observational studies which aim to constrain cosmology from small-scale structures.