The identification of relevant collective coordinates is crucial for the interpretation of coherent nonlinear spectroscopies of complex molecules and liquids. Using an h expansion of Liouville space generating functions, we show how to factorize multitime nonlinear response functions into products of lower-order correlation functions of collective coordinates, and derive closed expressions for linear, second- and third-order response functions. In addition to providing systematic quantum corrections, h offers a convenient bookkeeping device even for the purely classical response, since including quantum fluctuations allows to circumvent the expensive computation of stability matrices which is a major bottleneck in molecular dynamics simulations. The existing classical simulation strategies, including mode coupling in k space and in real space, Langevin equations, and instantaneous normal modes are compared from a unified viewpoint. (C) 2003 American Institute of Physics.