Global convergence of neuron birth-death dynamics

Introduction

As a consequence of the universal approximation theorems, sufficiently wide single layer neural networks are expressive enough to accurately represent a broad class of functions [Cyb89, Bar93, PS91]. The existence of a neural network function arbitrarily close to a given target function, however, is not a guarantee that any particular optimization procedure can identify the optimal parameters. Recently, using mathematical tools from optimal transport theory and interacting particle systems, it was shown that gradient descent [RVE18, MMN18, SS18, CB18b] and stochastic gradient descent converge asymptotically to the target function in the large data limit.

This analysis relies on taking a “mean-field” limit in which the number of parameters $n$ tends to infinity. In this setting, gradient descent optimization dynamics is described by a partial differential equation (PDE), corresponding to a Wasserstein gradient flow on a convex energy functional. While this PDE provides a powerful conceptual framework for analyzing the properties of neural networks evolving under gradient descent dynamics, the formula confers few immediate practical advantages. Nevertheless, analysis of this Wasserstein gradient flow motivates the interesting possibility of altering the dynamics to accelerate convergence.

In this work, we propose a dynamical scheme involving a parameter birth/death process. It can be defined on systems of interacting (e.g., neural network optimization) or non-interacting particles. We prove that the resulting modified transport equation converges to the global minimum of the loss in both interacting and non-interacting regimes (under appropriate assumptions), and we provide an explicit rate of convergence in the latter case for the mean-field limit. Interestingly—and unlike the gradient flow—the only fixed point of the dynamics is the global minimum of the loss function. We study the fluctuations of finite particle dynamics around this mean-field convergent solution, showing that they are of the same order throughout the dynamics and therefore providing algorithmic guarantees directly applicable to finite single-layer neural network optimization. Finally, we derive algorithms that converge to the birth-death PDEs and verify numerically that these schemes accelerate convergence even for finite numbers of parameters.

Global convergence and monotonicity of the energy with birth-death dynamics — We propose in Section 3 two distinct modifications of the original gradient flow that can be interpreted as birth-death processes. In this sense, the processes we describe amount to non-local mass transport in the equation governing the parameter distribution. We prove that the schemes we introduce guarantee global convergence and increase the rate of contraction of the energy compared to gradient descent and stochastic gradient descent for fixed $\mu$. We also derive asymptotic rates of convergence (Section 4).

Analysis of fluctuations and self-quenching — The birth-death dynamics introduces additional fluctuations that are not present in gradient descent dynamics. In Section 5 we calculate these fluctuations using tools from the theory of measure-valued Markov processes. We show that these fluctuations, for $n$ sufficiently large, are of order $O(n^{-1/2})$ and “self-quenching” in the sense that they diminish in magnitude as the quality as the optimization dynamics approaches the optimum.

Algorithms for realizing the birth-death schemes — In Section 6 we detail numerical schemes (and provide implementations in PyTorch) of the birth-death schemes described below. In the particular case of neural networks, the computational cost of implementing our procedure is minimal because no additional gradient computations are required. We demonstrate the efficacy of these algorithms on simple, illustrative examples in Section 7.

Related Works

Non-local update rules appear in various areas of machine learning and optimization. Derivative-free optimization [RS13] offers a general framework for optimizing complex non-convex functions using non-local search heuristics. Some notable examples include Particle Swarm Optimization [Ken11] and Evolutionary Strategies, such as the Covariance Matrix Adaptation method [Han06]. These approaches have found some renewed interest in the optimization of neural networks in the context of Reinforcement Learning [SHC+17, SMC+17] and hyperparameter optimization [JDO+17].

Our setup of non-interacting potentials is closely related to the so-called Estimation of Distribution Algorithms [BC95, LL01], which define update rules for a probability distribution over a search space by querying the values of a given function to be optimized. In particular, Information Geometric Optimization Algorithms [OAAH17] study the dynamics of parametric densities using ordinary differential equations, focusing on invariance properties. In contrast, our focus in on the combination of transport (gradient-based) and birth/death dynamics.

Dropout [SHK+14] is a regularization technique popularized by the AlexNet CNN [KSH12] reminiscent of a birth/death process, but we note that its mechanism is very different: rather than killing a neuron and replacing it by a new one with some rate, Dropout momentarily masks neurons, which become active again at the same position; in other words, Dropout implements a purely local transport scheme, as opposed to our non-local dynamics.

Finally, closest to our motivation is [WLLM18], who, building on the recent body of works that leverage optimal transport techniques to study optimization in the large parameter limit [RVE18, CB18b, MMN18, SS18], proposed a modification of the dynamics that replaced traditional stochastic noise by a resampling of a fraction of neurons from a base, fixed measure. Our model has significant differences to this scheme, namely we show that the dynamics preserves the same global minimizers and accelerates the rate of convergence. Finally, our interpretation of the modified dynamics in terms of a generalized gradient flow is related to the unbalanced optional transport setups of [KMV16, LMS18, CPSV18].

Mean-field PDE and Birth-death Dynamics

we see that, up to an irrelevant constant depending only on the data distribution, we arrive at (1) with

We also consider non-interacting objective functions in which $K=0$ in (1). Optimization problems that fit this framework include resource allocation tasks in which, e.g., weak performers are eliminated, Evolution Strategies, and Information Geometric Optimization [OAAH17].

In the case of gradient descent dynamics, the evolution of the particles $\boldsymbol{\theta}_{i}$ is governed for $i=1,\ldots,n$ by

To analyze the dynamics of this particle system, we consider the “mean-field” limit $n\to\infty$. As the number of particles becomes large, the empirical distribution of particles

leads to a deterministic partial differential equation at first order [RVE18, MMN18, CB18b, SS18],

and (8) should be interpreted in the weak sense in general:

where $C^{\infty}_{c}(D)$ denotes the space of smooth functions with compact support on $D$.

Interestingly, $V$ is the gradient with respect to $\mu$ of an energy functional $\mathcal{E}[\mu]$,

As a result, the nonlinear Liouville equation (8) is the Wasserstein gradient flow with respect to the energy functional $\mathcal{E}[\mu]$. Local minima of $V$ (where $\nabla V=0$) are clearly fixed points of this gradient flow, but these fixed points may not always be minimizers of the energy when $\operatorname{supp}\mu\subset D$. When the initial distribution of parameters has full support, neural networks evolving with gradient descent avoid these spurious fixed points under appropriate assumptions about their nonlinearity [CB18b, RVE18, MMN18].

2 Birth-Death augmented Dynamics

Here we consider a more general dynamical scheme that involves nonlocal transport of particle mass. As we shall see in Section 4, this dynamics avoids spurious fixed points and local minima, and converges asymptotically to the global minimum. Consider the following modification of the Wasserstein gradient flow above:

The additional term $-\alpha V\mu_{t}$ is a birth/death term that modifies the mass of $\mu$. If $V$ is positive, this mass will decrease, corresponding to the removal or “death” of parameters. If $V$ is negative, this mass will increase, which can be implemented as duplication or “cloning” of parameters. For a finite number of parameters, this dynamics could lead to changes in the architecture of the network. In many applications it is preferable to fix the total population, achieved by simply adding a conservation term to the dynamics,

where $\bar{V}\equiv\int_{D}Vd\mu_{t}$. This equation (like (12)) should in general be interpreted in the weak sense. Here we will focus on solutions of (13) for the initial condition $\mu_{0}\in\mathcal{M}(D)$, the space of probability measures on $D$, that satisfy

and $\boldsymbol{\Theta}(t,\boldsymbol{\theta})$ satisfies

Formula (14) can be formally established by solving (13) by the method of characteristics. In the non-interacting case, since $V(\boldsymbol{\theta},[\mu_{t}])=F(\boldsymbol{\theta})$, (14) is explicit and well-posed under appropriate assumptions on $F$ (see Assumption 4.1 below). In the interacting case, (14) is implicit since the right hand side depends on $\mu_{t}$. Following Chizat & Bach [CB18b], we know that under appropriate assumptions on $F$ and $K$ (see Assumption 4.4 below), solutions to (14) exist for all $t>0$ for appropriate initial $\mu_{0}$ that are compactly supported in $D$. Here we will assume global existence of solutions to this equation for $\mu_{0}$ such that $\operatorname{supp}\mu_{0}=D$ with $D$ open: if $\mu_{0}$ decays sufficiently fast at infinity, this assumption is supported by the alternative derivation of (12) based on a proximal gradient formulation given in Sec. 3.3.

Note that solutions of (12) that satisfy (14) are probability measures since they are positive by definition and we can set $\phi=1$ in (14) to deduce that $\mu_{t}(D)=1$. We can also show that the birth-death terms improve the rate of energy decay, as stated in the following proposition:

Let $\mu_{t}$ be a solution of (13) for the initial condition $\mu_{0}\in\mathcal{M}(D)$ that satisfies (14) for all $t\geq 0$. Then, $\mu_{t}(D)=1$ for all $t\geq 0$, and $E(t)=\mathcal{E}[\mu(t)]$ satisfies

Proof: (17) can be formally obtained by testing (13) against $V(\boldsymbol{\theta},[\mu_{t}])$ and using the chain rule to deduce that $d\mathcal{E}[\mu_{t}]/dt=\int_{D}V(\boldsymbol{\theta},[\mu_{t}])\partial_{t}\mu_{t}(d\boldsymbol{\theta})$. To complete the proof, we need to show that this testing is legitimate and the terms at the right hand side of (17) are well-defined; this is done in Appendix D by differentiating $C(t)$. $\square$

The birth-death term thus contributes to increase the rate of decay of the energy at all times. A natural question is whether such improved energy decay can lead to global convergence of the dynamics to the global minimum of the energy. As it turns out, the answer is yes: the fixed points of the birth-death PDEs (12) and (13) are the global minimizers of the energy $\mathcal{E}[\mu]$, as we prove in Section 4. How to implement a particle dynamics consistent with (13) is discussed in Sections 5 and 6.

We also note that there are several ways in which we can modify (13) to certain advantages: this is discussed in Appendix A.

3 Proximal formulation of birth-death dynamics

where $W_{2}(\mu,\mu_{k})$ denotes the $2$-Wasserstein distance between the probability measures $\mu$ and $\mu_{k}$. Interestingly, the birth-death PDE relies on a different measure of “distance”: the PDE

can be obtained as the time-continuous limit of the proximal optimization scheme: given an initial $\mu_{0}$ such that $\mathcal{E}[\mu_{0}]<\infty$, set

where the minimum is taken over all probability measures $\mu\in\mathcal{M}(D)$ and $D_{\text{KL}}(\mu||\mu_{k})$ is the Kullback-Leibler divergence

We verify this claim formally; notice that the Euler-Lagrange equation for the minimizer $\mu_{k+1}$, obtained by zeroing the first variation of the objective function in (20), reads

where $\lambda$ is a Lagrange multiplier added to enforce $\int_{D}d\mu_{k+1}=1$. (22) can be reorganized into

where $C$ is adjusted so that $\int_{D}d\mu_{k+1}=1$. (23) is the discrete equivalent of (14) If $\tau$ is small, we can expand the exponential to arrive at

Setting $\mu_{k+1}=\mu_{k}+O(\tau)$ in $V$ and expanding again gives

where we have also expanded $C$ and solved for it explicitly at leading order in $\tau$. Subtracting $\mu_{k}$ for both sides, dividing by $\tau$, and letting $\tau\to 0$ gives (19). The full PDE (13) can be obtained by alternating (18) and (20).

Convergence of Transport Dynamics with Birth-death

Here, we compare the solutions of the original PDE (8) with those of the PDE (13) with birth-death. We restrict ourselves to situations where $F$ and $K$ in (11) are such that $\mathcal{E}[\mu]$ is bounded from below. Our main technical contributions are results about convergence towards global energy minimizer as well as convergence rates as the dynamics approaches these minimizers. We consider separately the non-interacting and the interacting cases.

Under gradient descent dynamics, global convergence can be established with appropriate assumptions on the initialization and architecture of the neural network. [MMN18] establishes global convergence and provides a rate for neural networks with bounded activation functions evolving under stochastic gradient descent. Similar results were obtained in [CB18b, RVE18], in which it is proven that gradient descent converges to the globally optimal solution for neural networks with particular homogeneity conditions on the activation functions and regularizers. Closely related to the present work, [WLLM18] provides a convergence rate for a “perturbed” gradient flow in which uniform noise is added to the PDE (8). It should be emphasized that, unlike our formulation, the addition of uniform noise changes the fixed point of the PDE and convergence to only an approximate global solution can be obtained in that setting.

With no loss of generality we set $F(\boldsymbol{\theta}^{*})=0$ since adding an offset to $F$ in (13) does not affect the dynamics. We also denote by $H^{*}=\nabla\nabla F(\boldsymbol{\theta}^{*})$ the Hessian of $F$ at $\boldsymbol{\theta}^{*}$: recall that a Morse function is such that its Hessian is nondegenerate at all its critical points (where $\nabla F=0$) and it is coercive if $\lim_{\boldsymbol{\theta}\to\infty}F(\boldsymbol{\theta})=\infty$. Our main result is

Furthermore the rate of convergence becomes exponential in time asymptotically: for all $\delta>0$, $\exists t_{\delta}$ such that

The theorem is proven in Appendix B This proof shows that the additional birth-death terms in the PDE (12) allow the measure to concentrate rapidly in the vicinity of $\boldsymbol{\theta}^{*}$; subsequently, the transport term takes over and leads to the exponential rate of energy decay in (28). The proof also shows that, if we remove the transportation term $\nabla\cdot\left(\mu_{t}\nabla V\right)$ in the PDE (12), the energy only decreases linearly in time asymptotically. This means that the combination of the transportation and the birth-death terms accelerates convergence. A similar theorem can be proven for the PDE (55).

2 Interacting Case

Let us now consider the interacting case, when $V$ is given by (9) with $K\not=0$. We make

The set $D$ is a $k$-dimensional differentiable manifold which is either closed (i.e. compact, with no boundaries), or open (i.e. with no closed subset), or the Cartesian product of a closed and an open manifold.

This technical assumption typically holds for neural networks. Assumption 4.4 guarantees that the quadratic energy $\mathcal{E}[\mu]$ in (11) has a (unique) minimum value. While we cannot guarantee in general that this minimum is reached only by minimizers, below we will work under the assumption that minimizers exist. These are solutions in $\mathcal{M}(D)$ of following Euler-Lagrange equations:

where $\bar{V}[\mu]\equiv\int_{D}V(\boldsymbol{\theta},[\mu])\mu(d\boldsymbol{\theta})$. These equations are well-known [Ser15]: for the reader’s convenience we recall their derivation in Appendix C.

Minimizers of the energy should not be confused with fixed points of the dynamics. In particular, a well-known issue with the PDE (8) is that it potentially has many more fixed points than $\mathcal{E}[\mu]$ has minimizers: Indeed, rather than (30), these fixed points only need to satisfy

It is therefore remarkable that, if we pick an initial condition $\mu_{0}$ for the birth-death PDE (13) that has full support, the solution to this equation converges to a global minimizer of $\mathcal{E}[\mu]$:

Let $\mu_{t}$ denote the solution of (13) that satisfies (14) for the initial condition $\mu_{0}$ with $\operatorname{supp}\mu_{0}=D$. If $\mu_{t}\rightharpoonup\mu_{*}$ as $t\to\infty$ for some probability measure $\mu_{*}\in\mathcal{M}(D)$, then under Assumptions 4.3 and 4.4 $\mu_{*}$ is a global minimizer of $\mathcal{E}[\mu]$.

This theorem is proven in Appendix D. Note that the theorem holds under the assumption that $\mu_{t}$ converges to a fixed point $\mu_{*}$, which we cannot guarantee a priori but should be true for a wide class of $F$ and $K$ and initial conditions $\mu_{0}$ satisfying properties like $\mathcal{E}[\mu_{0}],\infty$—for more details on these conditions see the proof in Appendix D. One aspect of this proof is based on the evolution equation (17) for $\mathcal{E}[\mu_{t}]$. Since $d\mathcal{E}[\mu_{t}]/dt\leq 0$ and since $\mathcal{E}[\mu_{t}]$ is bounded from below by Assumption 4.4, by the bounded convergence theorem, the evolution must stop eventually. By assumption, this involves $\mu_{t}$ converging weakly towards some $\mu_{*}$. This happens when both integrals in (17) are zero, i.e. $\mu_{*}$ must satisfy the first equation in (30) as well as (31). What remains to be shown is that $\mu_{*}$ must also satisfy the second equation in (30), which we check in Appendix D.

Regarding the rate of convergence, we have the following result:

Under the same conditions as in Theorem 4.5, $\exists C>0$ and $t_{C}>0$ such that $E(t)=\mathcal{E}[\mu_{t}]-\mathcal{E}[\mu_{*}]\geq 0$ satisfies

The proof of this theorem is given in Appendix E where we show that

From Mean-field to Particle Dynamics with Birth-Death

In practice the number of parameters $n$ is finite, so we must verify that we can implement dynamics at finite particle numbers that is consistent with the PDEs with birth-death terms introduced in Sec. 3 in the mean-field limit $n\to\infty$. We must also ensure that the fluctuations arising from the discrete particles do not pose a problem for the optimization dynamics. In this section, we carry out this program in the context of the PDE (13). Analogous calculations can be performed in the case of (55). These results rely on the theory of measure-valued Markov processes [Daw06], and are detailed in Appendix F.

The dynamics of the particles $\{\boldsymbol{\theta}_{i}(t)\}_{i=1}^{n}$ is specified by a Markov process defined as follows: the birth-death part of the evolution is realized by equipping each particle $\boldsymbol{\theta}_{i}$ with an independent exponential clock with (signed) rate

Between these birth events the particles evolve by the GD flow (6).

Due to the interchangeability of the particles, the evolution of their empirical distribution $\mu^{(n)}_{t}$ defined in (7) is also Markovian: it is referred to in the probability literature as a measured-valued Markov process [Daw06]. We can write down the generator of this process, which specifies the evolution of the expectation of functionals of $\mu^{(n)}_{t}$, and analyze its behavior as $n\to\infty$. These calculations are performed in Appendix F, and they lead to:

Let the empirical distribution of the initial position of the particles be $\mu_{0}^{(n)}=n^{-1}\sum_{i=1}^{n}\delta_{\boldsymbol{\theta}_{i}(0)}$ and assume that $\mu^{(n)}_{0}\rightharpoonup\mu_{0}$ as $n\to\infty$. Then, for all for $t\in[0,\infty)$, $\mu_{t}^{(n)}=n^{-1}\sum_{i=1}^{n}\delta_{\boldsymbol{\theta}_{i}(t)}\rightharpoonup\mu_{t}$ in law as $n\to\infty$, where $\mu_{t}$ satisfies (13) with the initial condition $\mu_{t=0}=\mu_{0}$.

This statement verifies that, to leading order, the large particle limit recovers the mean-field PDE (13).

While the limit gives rise to the birth-death term of the PDE as expected, we can also quantify the scale and asymptotic behavior of the higher order fluctuations at finite $n$. This computation ensures that finite $n$ fluctuations do not overcome the convergence expected from the mean-field analysis. To do so, we we introduce the discrepancy distribution defined by the difference, scaled by $\sqrt{n}$, between the empirical distribution and its mean-field limit

where $\mu^{(n)}_{t}$ is the empirical distribution defined in (7) and $\mu_{t}$ is limit satisfying (54). We can then analyze the generator of the joint process $(\mu_{t},\omega^{(n)}_{t})$ and deduce the following proposition:

We should emphasize that these conclusions rely on $n$ being large enough that both the LLN and the CLT apply. In practical situations, it may be difficult to determine the threshold value of $n$ to reach this regime—it may grow with the dimension of $D$. At finite $n$, we also cannot rule out the possibility of some distinct dynamical regime in which the fluctuations grow with time—our results simply indicate that, in the regime where the LLN and CLT apply, the timescale for such a phenomenon would be diverging with $n$. These concerns are partially placated by the fact that our experiments show no signs of any such distinct dynamical regime and clearly indicate that birth-death helps accelerating convergence at moderate values of $n$.

Finally we want to stress that, while the calculations above indicate convergence with the birth-death dynamics alone when $n$ is large enough, the gradient flow probably plays a crucial part in accelerating the underlying optimization procedure, especially at moderate values of $n$. Without the transport term, the birth-death dynamics can only adjust the weight of existing neurons, which is clearly inefficient in some cases. That is, we do not advocate the use of birth-death dynamics alone, but rather to combine it with GD.

Algorithms

Numerical schemes that converge to the PDEs presented in Sec. 3 are both straightforward to design and easy to implement. In absence of the GD part of the dynamics, we could use Kinetic Monte Carlo (also called the Gillespie algorithm) to simulate birth-death without time-discretization error. However, in the large parameter regime, this would be computationally expensive: every particle has its own exponential clock, and the time between successive birth-death events scales like $1/n$. Because we must time-discretize the GD flow, we carry out the birth-death dynamics using the same time-discretization.

While this type of update is standard in machine learning, more accurate integration schemes could be used.

The corresponding particle system is a discretized version, both in particle number and time, of the PDE (13) and it converges to this equation as $n\to\infty$ and $\Delta t\to 0$. The error we make at finite $n$ is analyzed in Sec. 5; the error we make at finite $\Delta t$ can be deduced from standard results about time discretization of differential equations: with the Euler scheme used above, this error scales as $O(\Delta t)$.

In the case of neural network parameter optimization, the birth-death algorithm does not incur any significant computational cost beyond regular stochastic gradient descent. Denoting the parameters $\boldsymbol{\theta}_{i}=(c_{i},\boldsymbol{y}_{i}).$ and writing the neural network function as

the potential $V(\boldsymbol{\theta}_{i})=F(\boldsymbol{\theta}_{i})+n^{-1}\sum_{j=1}^{n}K(\boldsymbol{\theta}_{i},\boldsymbol{\theta}_{j})$ is given by

Since this quantity is computed in the SGD update, the only additional computation is the sum of $V_{P}$ over the $n$ particles. The cost of the algorithm is $O(nP)$ at every iteration.

For neural networks of the form given in Eq. (40) a particularly simple modification of Algorithm 1 enables particle creation from a prior distribution. The algorithm proceeds through the initial birth-death loop as in Algorithm 1. At the end of the initial loop, if the total population has decreased, then additional particle are sampled with configurations $(c,\boldsymbol{y})$ distributed according to the prior distribution

so that a reinjected particle has zero contribution to the total energy.

Finally, let us note that it is possible to design algorithms for the particles that mimic the proximal optimization scheme introduced in (20). For concreteness we focus on the cases of neural networks—the ideas below can be easily adapted to the others situations treated in this paper. Assume that the neural representation at iterate $k$ is

where $\boldsymbol{\theta}_{i}^{k}$ denotes the parameter in the network and $w_{i}^{k}\geq 0$ are extra weights satisfying $n^{-1}\sum_{i=1}^{n}w^{k}_{i}=1$—we will define a dynamics for these weights in a moment. Notice that (44) can be written as

1. Gradient step. Evolve the parameters $\boldsymbol{\theta}_{i}^{k}$ by GD (or SGD if we need to use the empirical loss) with the weights $w_{i}^{k}$ kept fixed. Do this for $m$ steps of size $\Delta t$ to obtain a new set of $\{\boldsymbol{\theta}_{i}^{k+1}\}_{i=1}^{n}$.

2. Proximal step. Evolve the weights $w_{i}^{k}$ with the parameter $\boldsymbol{\theta}_{i}^{k+1}$ fixed using a proximal step based on the particle equivalent of (20), i.e.

where the minimization is done under the constraint that $n^{-1}\sum_{i=1}^{n}w_{i}=1$. The equation for the minimizer $w_{i}^{k+1}$ is the discrete equivalent of (24)

where $C$ is a constant to be adjusted so that $n^{-1}\sum_{i=1}^{n}w^{k+1}_{i}=1$ and

(48) is implicit in $w_{i}^{k+1}$ and should be solved by iteration. Note that this proximal step is guaranteed to decrease the loss. In practice, this step could eventually lead to big variations of the weights. Should this happen, we add the additional step:

3. Resampling step. Resample the weights $\{w_{i}^{k+1}\}_{i=1}^{n}$ so as to keep them roughly equal to $1$ each, that is: eliminate the ones that are too small and transfer their weights to the others: split the remaining (large) weights into bits of size roughly $1$. There are standard ways to do this resampling step that are unbiased and preserve the population size exactly. This resampling step may increase the loss, though not to leading order. This step is the actual birth-death step in the scheme (and it is also the only random component of it if the exact loss is used).

If we set $\tau=\alpha m\Delta t$ and set $\Delta t\to 0$ and $n\to\infty$, the scheme above is formally consistent with the PDE

However, it is obviously not necessary to take either of these limits explicitly in practice, and, as explained above, the proximal step is guaranteed to decrease the loss. With a strict version of the the resampling step performed at every iteration, in which the weights are taken to be in $\{0,1\}$ the scheme above recovers the one described in Algorithm 1. The main difference is that in Algorithm 1 the proximal step (48) is solved in one iteration, by substituting $w^{k+1}_{i}$ by $w_{i}^{k}$ at the right hand side of (48).

Finally notice that if we were to implement the proximal step only and skip both the gradient and the resampling steps, the scheme above is a naive implementation of the lazy training scheme discussed in [CB18a]. This highlights again why using birth-death alone is not an efficient way to perform network optimization, and it should be combined with standard GD.

Numerical Experiments

We take as an illustrative example a mixture of Gaussians in dimension $d$,

which we approximate as a neural network with Gaussian nonlinearities with fixed standard deviation $\sigma<\min_{i}\sigma_{i}$,

denoting the parameters $\boldsymbol{\theta}_{i}=(c_{i},\boldsymbol{y}_{i}).$ This is a useful test of our results because we can do exact gradient descent dynamics on the mean-squared loss function:

In Fig. 1, we show convergence to the energy minimizer for a mixture of three Gaussians (details and source code are provided in the SM). The non-local mass transport dynamics dramatically accelerates convergence towards the minimizer. While gradient descent eventually converges in this setting—there is no metastability—the dynamics are particularly slow as the mass concentrates near the minimum and maxima of the target function. However, with the birth-death dynamics, this mass readily appears at those locations. The advantage of the birth-death dynamics with a reinjection distribution $\mu_{\textrm{b}}$ is highlighted by choosing an unfavorable initialization in which the particle mass is concentrated around $y=-2.$ In this case, both GD and GD with birth-death (12) do not converge on the timescale of the dynamics. With the reinjection distribution, new mass is created near $y=2$ and convergence is achieved.

2 Student-Teacher ReLU Network

As shown in Fig. 2, we find that the birth-death dynamics accelerates convergence to the teacher network. We emphasize that because the birth-death dynamics is stochastic at finite particle numbers, the fluctuations associated with the process could be unfavorable in some cases. In such situations, it is useful to reduce $\alpha$ as a function of time. On the other hand, in some cases we have observed much more dramatic accelerations from the birth-death dynamics.

Conclusions

The success of an optimization algorithm based on gradient descent requires good coverage of the parameter space so that local updates can reach the minima of the loss function quickly. Our approach liberates the parameters from a purely local dynamics and allows rapid reallocation to values at which they can best reduce the approximation error. Importantly, we have constructed the non-local birth-death dynamics so that it converges to the minimizers of the loss function. For a very general class of minimization problems—both interacting and non-interacting potentials—we have established convergence to energy minimizers under the dynamics described by the mean-field PDE with birth-death. Remarkably, for interacting systems with we can guarantee global convergence for sufficiently regular initial conditions. We have also computed the asymptotic rate of convergence with birth-death dynamics.

These theoretical results translate to dramatic reductions in convergence time for our illustrative examples. It is worth emphasizing that the schemes we have described are straightforward to implement and come with little computational overhead. Extending this type of dynamics to deep neural network architectures could accelerate the slow dynamics at the initial layers often observed in practice. Hyperparameter selection strategies based on evolutionary algorithms [SMC+17] provide another interesting potential application of our approach.

While we have characterized the basic behavior of optimization under the birth-death dynamics, many theoretical questions remain. First, we did not address generalization; understanding the role of the extra birth/death term in controlling the generalization gap is an important future question, in particular relating it to the lazy-training regime of [CB18a]. Next, we need to assume the existence of weak solutions through (14) with an initial measure $\mu_{0}$ that has full support, yet it may be possible to certify that the dynamics exist for all times if $\mu_{0}$ decays sufficiently fast. Besides, more explicit calculations of global convergence rates for the interacting case and tighter rates for the non-interacting case would be exciting additions. The proper choice of $\mu_{\textrm{b}}$ is another question worth exploring because, as highlighted in our simple example, favorable reinjection distributions can rapidly overcome slow dynamics. Finally, a mean-field perspective on deep neural networks would enable us to translate some of the guarantees here to deep architectures.

Acknowledgments

We would like to acknowledge the useful and detailed comments by Sylvia Serfaty and Yann Ollivier on previous versions of this manuscript.

References

Appendix A Generalizations of (13)

Here we mention two ways in which we can modify (13) to certain advantages. For example, we can replace this equation with

While the birth-death dynamics described above ensures convergence in the mean-field limit, when $n$ is finite, particles can only be created in proportion to the empirical distribution $\mu^{(n)}.$ In particular, such a birth process corresponds to “cloning” or creating identical replicas of existing particles. In practice, there may be an advantage to exploring parameter space with a distribution distinct from the instantaneous empirical particle distribution (7). To enable this exploration we introduce a birth term proportional to a distribution $\mu_{\textrm{b}}$ which we will assume has full support on $D$. In this case, the time evolution of the distribution is described by

where $\alpha,\alpha^{\prime}>0$, $(V-\bar{V})_{+}=\max(V-\bar{V},0)\geq 0$, $(V-\bar{V})_{-}=\max(\bar{V}-V,0)\geq 0$. That is, we kill particles in proportion to $\mu_{t}$ in region where $V>\bar{V}$ but create new particles from $\mu_{\textrm{b}}$ in regions where $V\leq\bar{V}$. We could also combine (54) with (55) to obtain other variants.

These alternative birth-death dynamical schemes also satisfy the consistency conditions of Proposition 3.1:

Proof: By considering again $1$ and $V(\cdot,[\mu_{t}])$ as a test function in (54) or (55), we verify that $\partial_{t}\mu_{t}(D)=0$. In addition, (54) implies that

which proves (56) for (55) since all the terms at the right hand side of this equation are negative. $\square$

Appendix B Convergence and Rates in the Non-interacting Case

Let us look first at the PDE satisfied by the measure $\mu$ in the non-interacting case, i.e. with $V=F$ satisfying Assumption 4.1, and without the transportation term:

Therefore, by plugging this last expression in equation (58), we obtain the explicit expression

where $G(\alpha t)$ is the function defined as:

At late times, the factor $e^{-\alpha tF(\boldsymbol{\theta})}$ focuses all the mass in the vicinity of the global minimum of $F$. Therefore, we can neglect the influence of the density $\rho_{0}$ in this integral. More precisely a calculation using the Laplace method indicates that

where $H^{*}=\nabla\nabla F(\boldsymbol{\theta}^{*})$ is the Hessian at the global minimum located at $\boldsymbol{\theta}^{*}$, and $\sim$ indicates that the ratio of both sides of the equation tend to 1 as $\alpha t\to\infty$. This shows that

B.2 Non-interacting Case with Transportation and Birth-death

We first prove the following intermediate result

Proof: By slightly abusing notation, we define

We consider the following Lyapunov function:

Observe that $0\leq f_{\delta}(t)<1$ because otherwise $F$ would be flat (in which case the energy is ). Also, we can assume wlog that $E(t)-\delta>0$, since otherwise the statement of the lemma is trivially verified. By plugging (67) and (B.2.1) into (66) we have

Finally, since $f^{-1}_{\delta}(t)\geq 0$, we have

which concludes the proof of the Lemma. $\square$

Proof of Theorem 4.2: In order to prove (27), we apply the previous lemma for $\delta\to 0$. Let $\theta^{*}=\arg\min V(\theta)$, We have $F(\boldsymbol{\theta}^{*})=0$, and $\|\nabla\nabla F(\boldsymbol{\theta})\|\leq\beta$ for some $\beta>0$. Then, for $\delta$ sufficiently small, the indicator function $\phi_{\delta}(\boldsymbol{\theta})$ is localized in the set

where $H^{*}=\nabla\nabla F(\boldsymbol{\theta}^{*})$. It follows that for sufficiently small $\delta$,

which implies that in order to reach an error $\epsilon$, we need

For large $t$, we can again use Laplace method to confirm that $\rho(t,\boldsymbol{\theta})$ concentrates near the absolute minimum of $F(\boldsymbol{\theta})$ located at $\boldsymbol{\theta}^{*}$. To see why notice that $\boldsymbol{\Theta}(t,\boldsymbol{\theta})$ converge, as $t\to\infty$, near local minima of $F$. Suppose that these minima are located at $\boldsymbol{\theta}_{1}^{*}=\boldsymbol{\theta}^{*}$, $\boldsymbol{\theta}_{2}^{*}$, etc. At these minima we have $\nabla F(\boldsymbol{\theta}^{*}_{j})=0$, and if in (79) we replace $F(\boldsymbol{\theta})$ by its quadratic approximation around any $\boldsymbol{\theta}_{j}^{*}$, $\tfrac{1}{2}\langle\boldsymbol{\theta}-\boldsymbol{\theta}_{j}^{*},H_{j}^{*}(\boldsymbol{\theta}-\boldsymbol{\theta}_{j}^{*})\rangle$ with $H_{j}^{*}=\nabla\nabla H(\boldsymbol{\theta}_{j}^{*})$ positive definite, the solution to this equation reads

This quantifies the late stages of the global convergence to the minimum and confirms the asymptotic decay rate in (28), thereby concluding the proof of Theorem 4.2. $\square$

Denote by $\boldsymbol{\Theta}(t,\boldsymbol{\theta})$ the solution of the ODE

Then under the conditions of Theorem 4.2, the solution $\mu_{t}$ of the PDE (12) has a density $\rho_{t}$ given by

where $G(\boldsymbol{\theta})=\Delta F(\boldsymbol{\theta})-\alpha F(\boldsymbol{\theta})$.

Proof: Since the initial $\mu_{0}$ has a density $\rho_{0}>0$, so does $\mu_{t}$ for all $t>0$ (but not in the limit as $t\to\infty$) and its density satisfies

If $\boldsymbol{\Theta}(t,\boldsymbol{\theta})$ satisfies

By using $\boldsymbol{\Theta}(t,\boldsymbol{\Theta}(s,\boldsymbol{\theta}))=\boldsymbol{\Theta}(t+s,\boldsymbol{\theta})$ and the normalization condition, this implies

This is (77) and terminates the proof of the lemma. $\square$

Appendix C Derivation of (30)

Let $\mu_{*}$ be a minimizer and compare its energy to that of any other probability measure $\mu$. Since the energy minimum is unique by convexity, we must have $\mathcal{E}[\mu]\geq\mathcal{E}[\mu_{*}]$. A direct calculation shows that

The last term at the right hand side is always non-negative. Focusing on the second term, if we denote $\operatorname{supp}\mu_{*}=D_{*}$, we can write it as

where we used $V(\boldsymbol{\theta},[\mu_{*}])=\bar{V}[\mu_{*}]$ on $D_{*}$ and $\mu_{*}=0$ on $D_{*}^{c}$. The only possibility to make this term nonnegative for all $\mu$ is to have $V(\boldsymbol{\theta},[\mu_{*}])\geq\bar{V}[\mu_{*}]$ on $D_{*}^{c}$.

Appendix D Proof of Theorem 4.5

We begin by noting that, if (14) holds for al $t>0$, then $\bar{V}[\mu_{t}]=-\alpha^{-1}d\log C(t)/dt$ must be well-defined at all times. From (15), this derivative is given by

Using (16) to replace $\dot{\boldsymbol{\Theta}}(t,\boldsymbol{\theta})$ by $-\nabla V(\boldsymbol{\Theta}(t,\boldsymbol{\theta}),[\mu_{t}])$ and (14) to express these integral as expectations against $\mu_{t}$ gives

Therefore the terms at right hand side of (17) must be well-defined and we must also have

Since $\mu_{t}\rightharpoonup\mu_{*}\in\mathcal{M}(D)$ by assumption, we can take the limit as $t\to\infty$ to deduce that

We will use these properties below, along with

which is require in order that both $\bar{V}[\mu_{t}]$ and $\mathcal{E}[\mu_{t}]$ be well-defined at all $t>0$ and in the limit as $t\to\infty$.

With these preliminaries, we now recall that the argument given after Theorem 4.5 implies that any fixed point $\mu_{*}$ of the PDE (13) must satisfy the first equation in (30). That is, we must have

Therefore, to prove Theorem 4.5, it remains to show that the second equation in (30) must be satisfied as well. We will argue by contradiction: Let $D_{*}=\operatorname{supp}\mu_{*}$, assume $D_{*}^{c}\not=\emptyset$, and suppose that there exists a region $N\subseteq D_{*}^{c}$ where $V(\boldsymbol{\theta},[\mu_{*}])<\bar{V}[\mu_{*}]$. If it exists, this region must have nonzero Hausdorff measure in $D$ since, by Assumption 4.4, $V(\boldsymbol{\theta},[\mu_{t}])\in C^{2}(D)$ for all $t\geq 0$ and $V(\boldsymbol{\theta},[\mu_{*}])\in C^{2}(D)$. $V(\boldsymbol{\theta},[\mu_{*}])-\bar{V}[\mu_{*}]$ must also reach a minimum value inside $D$ even if $D$ is open, for otherwise (16) would eventually carry mass towards infinity, which contradicts $\mu_{t}\rightharpoonup\mu_{*}$. This implies that, if we pick $\delta\in(0,\bar{V}[\mu_{*}]-\min_{\boldsymbol{\theta}}V(\boldsymbol{\theta},[\mu_{*}]))$ and let

then $N_{\delta}$ is not empty. Since $V(\boldsymbol{\theta},[\mu_{*}])$ is twice differentiable in $\boldsymbol{\theta}$, for $\delta$ close enough to $\bar{V}[\mu_{*}]-\min_{\boldsymbol{\theta}}V(\boldsymbol{\theta},[\mu_{*}])$, $N_{\delta}$ is also compact and such that

Given any solution $\mu_{t}$ of the PDE (13) that is supposed to converge to $\mu_{*}$ as $t\to\infty$, consider

Since $\mu_{t}$ is positive everywhere at any finite time, we must have $f_{\delta}(t)>0$ for $t\in(0,\infty)$ However, since $\mu_{t}\to\mu_{*}$, we must also have

where $\hat{n}(\boldsymbol{\theta})$ is the inward pointing unit normal to $\partial N_{\delta}$ at $\boldsymbol{\theta}$ and $\sigma_{t}$ is the probability measure on $\partial N_{\delta}$ obtained by restricting $\mu_{t}$ on this boundary: If $\phi_{\epsilon}\in C^{\infty}_{c}(D)$ is a sequence of test functions with $\operatorname{supp}\phi_{\epsilon}=N_{\delta}$ and converging towards the indicator set of $N_{\delta}$ as $\epsilon\to 0$, $\sigma_{t}$ is defined as

Restricting ourselves to $t>t_{+}$, we therefore have

where we used the definition of $N_{\delta}$. Looking at the last term, we can assess its magnitude using

where (using the compactness of $N_{\delta}$)

Since we work under the assumption that $\mu_{t}\rightharpoonup\mu_{*}$, $M(t)$ must tend to as $t\to\infty$. As a result, $\exists t_{\delta}>0$ such $\forall t>t_{\delta}$ we have $N(t)<\delta$, which, from (104), implies that $\forall t>\max(t_{+},t_{\delta})$ we have $\dot{f}_{\delta}(t)>0$, a contradiction with (95). Therefore the only fixed points accessible by the PDE (13) are those for which both equations in (30) hold, which proves the theorem.

Appendix E Proof of Theorem 4.6

Let $\mu_{*}=\lim_{t\to\infty}\mu_{t}$ be the stationary point reached by the solution of (13) and denote $E(t)=\mathcal{E}[\mu_{t}]-\mathcal{E}[\mu_{*}]\geq 0$. Then

where we used $\int_{D}V^{2}d\mu_{t}-\bar{V}^{2}=\int_{D}|V-\bar{V}|^{2}d\mu_{t}$. By convexity

In Lemma E.1 below we show that $\exists t_{+}>0$ such that

As a result, $dE^{-1}/dt\geq\alpha$ for $t>t_{+}$. Integrating this relation in time on $[t_{0},t]$ with $t_{+}<t_{0}\leq t$ gives

Note that the proof only takes into account the effects of birth-death terms; adding transport may accelerate the rate.

There exist $t_{+}>0$ such that (111) holds.

Proof: Let $\nu_{t}=\mu_{t}-\mu_{*}$ and for future reference note that $\nu_{t}$ is a signed measure on $D_{*}=\operatorname{supp}\mu_{*}$ but $\nu_{t}\geq 0$ on $D^{*}_{c}$. Denote

Recall that $V_{*}=\bar{V}_{*}$ on $\operatorname{supp}\mu_{*}$. As a result

We can combine these two equations to obtain

where we used $\int_{D}\bar{V}_{*}d\nu_{t}=\bar{V}_{*}\int_{D}(d\mu_{t}-d\mu_{*})=0$ to get the penultimate equality and $V_{\star}-\bar{V}_{*}=0$ on $D_{*}$ to get the last.

Proceeding similarly using again $V_{*}=\bar{V}_{*}$ on $\operatorname{supp}\mu_{*}$ as well as $\int_{D}d\nu_{t}=\int_{D}(d\mu_{t}-d\mu_{*})=0$, we can also obtain

Let us now compare the square of (118) to (120). Since $V_{*}-\bar{V}_{*}\geq 0$ and $\nu_{t}\geq 0$ on $D_{*}^{c}$, we have

Case 1: $\int_{D^{c}_{*}}(V_{*}-\bar{V}_{*})d\nu_{t}>0$ (which requires $D_{*}^{c}\not=\emptyset$). Since $\nu_{t}\rightharpoonup 0$ as $t\to\infty$ the last term in (116) is higher order. As a result, for any $\delta>0$, $\exists t_{1}>0$ such that

which also implies that (using again $\nu_{t}\geq 0$ on $D^{c}_{*}$)

Similarly, the first term at the right hand side of (120) dominates all the other ones as $t\to\infty$ in the sense that, for any $\delta>0$, $\exists t_{2}>0$ such that

Taken together, (124) and (125) imply the statement of the lemma with any $C>0$ (since $\nu_{t}(D^{c}_{*})\to 0$ as $t\to\infty$). As a result $\lim_{t\to\infty}tE(t)=0$ in this case since $\int_{D}|V-\bar{V}|^{2}d\mu_{t}/|\int_{D}(V-\bar{V})d\mu_{*}|^{2}\to\infty$.

Case 2: $\int_{D^{c}_{*}}(V_{*}-\bar{V}_{*})d\nu_{t}=0$ (i.e. $D_{*}^{c}=\emptyset$ or $V_{*}=\bar{V}_{*}$ on $D_{*}^{c}$ as well as $D_{*}$). In this case it is easier to use (119) via the inequality

where we use the fact that $R$ reduces to (using $V_{*}=\bar{V}_{*}$ and $\int_{D}V_{*}d\nu_{t}=\bar{V}_{*}\int_{D}(d\mu_{t}-d\mu_{*})=0$)

Since $\int_{D}K(\boldsymbol{\theta},\boldsymbol{\theta}^{\prime})\nu_{t}(d\boldsymbol{\theta}^{\prime})\not=0$ on $D_{*}$, the leading order terms in $\int_{D}|V-\bar{V}|^{2}d\mu_{*}$ and $\int_{D}|V-\bar{V}|^{2}d\mu_{t}$ are the same and given by

That is, for any $\delta>0$, $\exists t_{3}>0$ such that

Together with (126), this implies the statement of the lemma with $C=1$. $\square$

Appendix F Proof of Propositions 5.1 and 5.2

Here we give formal proofs Propositions 5.1 and 5.2 using tools from the theory of measure-valued Markov processes [Daw06].

where $\mu^{(n)}_{t-}=\lim_{\epsilon\to 0+}\mu^{(n)}_{t-\epsilon}$ Similarly if particle $\boldsymbol{\theta}_{i}(t)$ gets duplicated at time $t$ and particle $\boldsymbol{\theta}_{j}(t)$ gets killed, the change this induces on $\mu^{(n)}_{t}$ is

We can use the properties of the Dirac distribution to rewrite the generator in (135) as

The operator in (137) is now defined for any $\mu\in\mathcal{M}(D)$, and we will use it in this form in our developments below.

The generator (137) can be used to write an evolution equation for the expectation of functionals evaluated on $\mu^{(n)}_{t}$. That is, if we define

then this time-dependent functional satisfies the backward Kolmogorov equation (BKE)

The proof of Proposition 5.1 is based on analyzing the properties of this equation in the limit as $n\to\infty$, which we expand upon in Appendix F.1. The proof of Proposition 5.2 is based on writing a similar equation for an extended process in which we magnify the dynamics of $\mu^{(n)}_{t}$ around its limit, as shown in Appendix F.2.

If we take the limit of $({\mathcal{L}}_{n}{\Phi})[\mu^{(n)}]$ as $n\to\infty$ on a sequence such that $\mu^{(n)}\rightharpoonup\mu$, we deduce that $({\mathcal{L}}_{n}{\Phi})[\mu^{(n)}]\to({\mathcal{L}}{\Phi})[\mu]$ with

Correspondingly, in this limit the BKE (140) becomes

Since (141) is precisely the generator of process defined by the PDE (13), this shows that, if $\mu^{(n)}_{t=0}=\mu^{(n)}\rightharpoonup\mu$ as $n\to\infty$, then

where $\mu_{t}$ solves the PDE (13) for the initial condition $\mu_{t=0}=\mu$. This proves the weak version of the LLN stated in Proposition 5.1.

F.2 Proof of Proposition 5.2

To quantify the fluctuations around the LLN, let $\mu_{t}$ be the limit of $\mu^{(n)}_{t}$ (i.e. the solution to the PDE (13)) and define

and similarly for $D^{2}_{\omega}\hat{\Phi}$. The operator in $\mu$ in (145) is the same as in (141), confirming the LLN; the operator in $\omega$ is a second order operator, i.e. it is the generator of a stochastic differential equation. That is, we have established that, as $n\to\infty$,

where $\omega_{t}(d\boldsymbol{\theta})$ is Gaussian random distribution whose equation can be obtained from the generator in (146) Formally

where $\eta(t)$ is a white-noise term with covariance consistent with (146):

This equation should also be interpreted in the weak sense by testing it against some $\phi\in C_{c}^{\infty}(D\times D)$, and it can be seen that it conserves mass in the sense that $\Sigma_{t}(d\boldsymbol{\theta},D)=\Sigma_{t}(D,d\boldsymbol{\theta}^{\prime})=0$ for all $t>0$ since this is true initially and $\partial_{t}\Sigma_{t}(d\boldsymbol{\theta},D)=\partial_{t}\Sigma_{t}(D,d\boldsymbol{\theta}^{\prime})=0$.