  The mean field (MF) theory of multilayer neural networks centers around a
particular infinite-width scaling, where the learning dynamics is closely
tracked by the MF limit. A random fluctuation around this infinite-width limit
is expected from a large-width expansion to the next order. This fluctuation
has been studied only in shallow networks, where previous works employ heavily
technical notions or additional formulation ideas amenable only to that case.
Treatment of the multilayer case has been missing, with the chief difficulty in
finding a formulation that captures the stochastic dependency across not only
time but also depth.
  In this work, we initiate the study of the fluctuation in the case of
multilayer networks, at any network depth. Leveraging on the neuronal embedding
framework recently introduced by Nguyen and Pham, we systematically derive a
system of dynamical equations, called the second-order MF limit, that captures
the limiting fluctuation distribution. We demonstrate through the framework the
complex interaction among neurons in this second-order MF limit, the
stochasticity with cross-layer dependency and the nonlinear time evolution
inherent in the limiting fluctuation. A limit theorem is proven to relate
quantitatively this limit to the fluctuation of large-width networks.
  We apply the result to show a stability property of gradient descent MF
training: in the large-width regime, along the training trajectory, it
progressively biases towards a solution with "minimal fluctuation" (in fact,
vanishing fluctuation) in the learned output function, even after the network
has been initialized at or has converged (sufficiently fast) to a global
optimum. This extends a similar phenomenon previously shown only for shallow
networks with a squared loss in the ERM setting, to multilayer networks with a
loss function that is not necessarily convex in a more general setting.
