Graph-based linear scaling electronic structure theory for quantum-mechanical molecular dynamics simulations [J. Chem. Phys. 144, 234101 (2016)] is adapted to the most recent shadow potential formulations of extended Lagrangian Born-Oppenheimer molecular dynamics, including fractional molecular-orbital occupation numbers [J. Chem. Phys. 152, 104103 (2020); Eur. Phys. J. B 94, 164 (2021)], which enables stable simulations of sensitive complex chemical systems that may have unsteady charge solutions. The proposed formulation includes a preconditioned Krylov subspace approximation for the integration of the extended electronic degrees of freedom, which requires quantum response calculations for electronic states with fractional occupation numbers. For the response calculations we introduce a graph-based canonical quantum perturbation theory that can be performed with the same natural parallelism and linear scaling complexity as the graph-based electronic structure calculations for the unperturbed ground state. The proposed techniques are particularly well-suited for semi-empirical electronic structure theory and the methods are demonstrated using self-consistent charge density-functional tight-binding (SCC-DFTB) theory, both for the acceleration of self-consistent field calculations and for quantum molecular dynamics simulations. The graph-based techniques combined with the semi-empirical theory enable stable simulations of large, complex chemical systems, including tens-of-thousands of atoms.