The standard diagram technique for impurity scattering is extended to the case in which the impurity potential linearly depends on the intensity of the local field. A Bethe-Salpeter equation is derived which describes diffusion of light in a nonlinear random medium. In order to cope diagrammatically with the nonlinear behavior of the impurities, one has to introduce vertices at which four propagators meet. As for the case of linear impurities, configurational averages can be taken at the diagrammatic level. However, the resulting dressed diagrams are of a much more complicated structure than usual. Furthermore, two diagrams that are topologically different from each other may represent the same expression. In decomposing the dressed diagrams, certain cuts can be performed only if other cuts are carried out first. This implies that the Bethe-Salpeter equation cannot be derived by simply cutting all dressed diagrams into irreducible pieces. Instead, a more sophisticated procedure has to be followed. It is developed by performing a topological classification of the dressed diagrams. In calculating the diagrams of the Bethe-Salpeter equation those with more than two vertices are discarded. This leads to three coupled nonlinear integral equations which describe radiative transport in a dilute suspension of weakly nonlinear Kerr particles. The dimension of the particles is small as compared to the wavelength of the light. The solution of the integral equations must satisfy a constraint that follows from energy conservation. Analytical work and predictions on experiments will be presented in a forthcoming paper.