The kernel polynomial method allows to sample overall spectral properties of a quantum system, while sparse diagonalization provides accurate information about a few important states. We present a method combining these two approaches without loss of performance or accuracy. We apply this hybrid kernel polynomial method to improve the computation of thermodynamic quantities and the construction of perturbative effective models, in a regime where neither of the methods is sufficient on its own. To achieve this we develop a perturbative kernel polynomial method to compute arbitrary order series expansions of expectation values. We demonstrate the efficiency of our approach on three examples: the calculation of supercurrent and inductance in a Josephson junction, the interaction of spin qubits defined in a two dimensional electron gas, and the calculation of the effective band structure in a realistic model of a semiconductor nanowire.