We present a practical nonperturbative method for exact treatment of intensity-dependent generalized nonlinear optical susceptibilities χ(ω) in intense polychromatic fields, valid for arbitrary laser intensities, detunings, and relaxation. By means of the many-mode Floquet theory, the time-dependent Liouville equation can be transformed into an equivalent time-independent infinite-dimensional Floquet–Liouville supermatrix (FLSM) eigenvalue problem. It is then shown that the nonlinear optical susceptibilities χ(ω) can be completely determined simply from the supereigenvalues and eigenfunctions of the Floquet–Liouvillian L̂F. In addition to this exact FLSM approach, we have also presented higher-order perturbative results, based on the extension of the Salwen’s nearly degenerate perturbation theory, appropriate for somewhat weaker fields and near-resonant multiphoton processes, but beyond the conventional perturbative or rotating wave approximation (RWA). In the case of two-level systems, for example, the implementation of Salwen’s method in the time-independent L̂F allows the reduction of the infinite-dimensional FLSM into a 4×4 dimensional effective Hamiltonian, from which essential analytical formulas for intensity-dependent χ(ω) can be obtained. These methods are applied to a detailed study of intensity-dependent spectral line shapes (such as hole burning and extra resonance peaks at the line center, and the effects of saturation, detuning, and radiative and collisional damping, etc.) and subharmonic structures in nonlinear multiple wave mixings χ[(m+1)ω1−mω2] for two-level systems in intense linearly polarized bichromatic fields.