Summary: We show that the formation of large-scale structures through gravitational instability in the expanding universe can be fully described through a path-integral formalism. We derive the action \(S[f]\) which gives the statistical weight associated with any phase-space distribution function \(f(\vec x,\vec p,t)\). This action \(S\) describes both the average over the Gaussian initial conditions and the Vlasov-Poisson dynamics. Next, applying a standard method borrowed from field theory we generalize our problem to an \(N\)-field system and we look for an expansion over powers of \(1/N\). We describe three such methods and we derive the corresponding equations of motion at the lowest nontrivial order for the case of gravitational clustering. This yields a set of nonlinear equations for the mean \(\overline{f}\) and the two-point correlation \(G\) of the phase-space distribution \(f\), as well as for the response function \(R\). These systematic schemes match the usual perturbative expansion on quasi-linear scales but should also be able to treat the nonlinear regime. Our approach can also be extended to non-Gaussian initial conditions and may serve as a basis for other tools borrowed from field theory.