A new theoretical approach to the problem of non-linear polarisation of matter is suggested. The density matrix equation for a low-density gas in a strong light field has been solved without using the phenomenological relaxation constants. Use is made of the general Hamiltonian that explicitly involves, besides the interaction of atoms with the external field, also their interactions with each other and with the vacuum electromagnetic field. The equation is solved by the canonical transformation method (a variant of the Bogolyubov-Mitropolski method of averages). It is generalised and used in the more complicated case of interacting multilevel atoms (molecules) in a strong light field. An expression for the cubic susceptibility is obtained, enabling analysis of the resonant saturation in both stationary and non-stationary cases. The proposed method provides correct expressions for the dynamic Stark shifts, amplitudes of multiphoton transitions and of radiational collisions, as well as for effective interatomic interactions caused by photon exchange. The proposed theoretical approach can also be used in studying non-linear optical processes in dense gases, liquids and solids.