We present the first fully simulation-based hierarchical analysis of the light curves of a population of low-redshift type Ia supernovæ (SNæ Ia). Our hardware-accelerated forward model, released in the Python package slicsim, includes stochastic variations of each SN’s spectral flux distribution (based on the pre-trained BayeSN model), extinction from dust in the host and in the Milky Way, redshift, and realistic instrumental noise. By utilizing truncated marginal neural ratio estimation (TMNRE), a neural network-enabled simulation-based inference technique, we implicitly marginalize over 4000 latent variables (for a set of ≈100 SNæ Ia) to efficiently infer SN Ia absolute magnitudes and host-galaxy dust properties at the population level while also constraining the parameters of individual objects. Amortization of the inference procedure allows us to obtain coverage guarantees for our results through Bayesian validation and frequentist calibration. Furthermore, we show a detailed comparison to full likelihood-based inference, implemented through Hamiltonian Monte Carlo, on simulated data and then apply TMNRE to the light curves of 86 SNæ Ia from the Carnegie Supernova Project, deriving marginal posteriors in excellent agreement with previous work. Given its ability to accommodate arbitrarily complex extensions to the forward model, e.g. different populations based on host properties, redshift evolution, complicated photometric redshift estimates, selection effects, and non-Ia contamination, without significant modifications to the inference procedure, TMNRE has the potential to become the tool of choice for cosmological parameter inference from future, large SN Ia samples.