  We develop and analyze a fault-tolerant quantum algorithm for computing
$n$-th order response properties necessary for analysis of non-linear
spectroscopies of molecular and condensed phase systems. We use a
semi-classical description in which the electronic degrees of freedom are
treated quantum mechanically and the light is treated as a classical field. The
algorithm we present can be viewed as an implementation of standard
perturbation theory techniques, focused on {\it ab initio} calculation of
$n$-th order response functions. We provide cost estimates in terms of the
number of queries to the block encoding of the unperturbed Hamiltonian, as well
as the block encodings of the perturbing dipole operators. Using the technique
of eigenstate filtering, we provide an algorithm to extract excitation energies
to resolution $\gamma$, and the corresponding linear response amplitude to
accuracy $\epsilon$ using
${O}\left(N^{6}\eta^2{{\gamma^{-1}}\epsilon^{-1}}\log(1/\epsilon)\right)$
queries to the block encoding of the unperturbed Hamiltonian $H_0$, in double
factorized representation. Thus, our approach saturates the Heisenberg
$O(\gamma^{-1})$ limit for energy estimation and allows for the approximation
of relevant transition dipole moments. These quantities, combined with
sum-over-states formulation of polarizabilities, can be used to compute the
$n$-th order susceptibilities and response functions for non-linear
spectroscopies under limited assumptions using
$\widetilde{O}\left({N^{5n+1}\eta^{n+1}}/{\gamma^n\epsilon}\right)$ queries to
the block encoding of $H_0$.
