We extend the LCAO (Linear Combination of Atomic Orbitals) method to excited states by constructing a particularly simple basis in the space of orbital products. The residual error of our procedure vanishes exponentially with the number of products and our procedure avoids auxiliary sets of fitting functions and their intrinsic ambiguities. As an application of our technique, we compute the Kohn--Sham density response function $\chi_{0}$ for a molecule consisting of $N$ atoms in $O(N^{2}N_{\omega })$ operations, with $N_{\omega }$ the number of frequency points. Our construction of $\chi_{0}$ allows us to compute molecular spectra directly from the equations of Petersilka--Gossmann--Gross in $O(N^{2}N_{\omega })$ operations rather than from Casida's equations which takes $O(N^{3})$ operations. Ongoing work indicates that our method is well suited to a computation of the GW self-energy $\Sigma=\mathrm{i}GW$ and we expect a similar situation for the Bethe--Salpeter equation.