The Laser Interferometer Space Antenna (LISA) will observe (among other things) the gravitational-wave emission of every mHz compact binary system in the Milky Way, collectively known as Galactic binaries (GBs). This astrophysical population will be realized within LISA data in two ways: ~10,000 individually-observable deterministic signals ("resolved" GBs) and a prominent stochastic confusion noise from the remaining tens of millions of systems ("unresolved" GBs, also known as the Galactic foreground). Together, the resolved and unresolved systems comprise a complete sample of mHz GBs -- a rare opportunity in astronomy, full of exciting potential to unveil new aspects of Galactic and stellar astrophysics.
However, current population inference approaches employed to study compact binary merger populations in LIGO-Virgo-KAGRA data will fail in the signal-dominated LISA regime, confounded by the global nature of LISA data analysis and the circular dependence of the resolved and unresolved GBs on each other and their shared underlying population. Hidden within this apparent challenge is a promising opportunity to maximize the astrophysical information gained by characterizing the entire population holistically --- resolved GBs, the unresolved Galactic foreground, and the line between them all providing crucial levers for astrophysical inference. We present a novel method to perform GB population inference in this complex data analysis landscape and discuss applications of this approach to similar problems in near-future PTA and 3G terrestrial datasets.