Low-frequency intermolecular dynamics in liquids is studied by ultrafast four- and six-wave mixing. The theory of these nonlinear optical processes is given for electronically nonresonant optical interactions up to fifth order in the electric field. The Born-Oppenheimer approximation is used to separate the motional part of the response functions from coordinate independent electronic hyperpolarizabilities. A large variety of experiments, involving far-infrared absorption, ordinary Rayleigh-Raman or hyper Rayleigh-Raman scattering is covered by this theory. The response in nonresonant six-wave mixing comprises four dynamically different processes. It is shown that one of the terms contains information on the time scale(s) of intermolecular dynamics, that is not available from lower-order nonresonant experiments. For instance, homogenous and inhomogeneous contributions to line broadening can be distinguished. The optical response of harmonic nuclear motion is calculated for nonlinear coordinate dependence of the polarizabilities. Results for level-dependent and level-independent damping of the motion are compared. It is shown that level-dependent damping destroys the interference between different quantum mechanical pathways, yielding an extra contribution to the fifth-order response that has not been discussed before. When two or more nuclear modes determine the optical response, their relative contributions to the four- and six-wave mixing signals are in general different. These contributions are determined by the coordinate dependence of the electronic polarizability, which is usually not fully known, Model calculations are presented for the dynamic parameters of liquid CS2. The theory of this paper will be employed in Part II, to analyze experimental results on femtosecond four- and six-wave mixing.