Over the last decades, large systems of interacting particles (or agents) are widely used to investigate selforganization and collective behavior. They frequently appear in modeling phenomena such as biological swarms [14,18], crowd dynamics [2,7,8], self-assembly of nanoparticles [38] and opinion formation [3,35,47]. Similar particle models are also used in metaheuristics [1, 6, 10, 31], which provide empirically robust solutions to tackle hard optimization problems with fast algorithms. Metaheuristics are methods that orchestrate an interaction between local improvement procedures and global/high level strategies, and combine random and deterministic decisions, to create a process capable of escaping from local optima and performing a robust search of a solution space. Starting with the groundbreaking work Ref. [51] of Rastrigin on Random Search in 1963, numerous mechanisms for multi-agent global optimization have been considered, among the most prominent instances we recall the Simplex Heuristics [48], Evolutionary Programming [28], the Metropolis-Hastings sampling algorithm [34], Genetic Algorithms [36], Particle Swarm Optimization (PSO) [41, 50], Ant Colony Optimization (ACO) [23] and Simulated Annealing (SA) [37,42]. Despite the tremendous empirical success of these techniques, most metaheuristical methods

still lack proper mathematical proof of robust convergence to global minimizers. In fact, due to the random component of metaheuristics, their asymptotic analysis would require to discern the stochastic dependencies, which is a daunting task, especially for those methods that combine instantaneous decisions with memory mechanisms. Recent work Refs. [13,49] on Consensus-Based Optimization (CBO) focuses on instantaneous stochastic and deterministic decisions in order to establish a consensus among particles on the location of the global minimizers within a domain. In view of the instantaneous nature of the dynamics, the evolution of the system can be interpreted as a system of first order stochastic differential equations, whose large particle limit is approximated by a deterministic partial differential equation of mean-field type. The large time behavior of such a deterministic PDE can be analyzed by classical techniques of large deviation bounds and the global convergence of the mean-field model can be mathematically proven in a rigorous way for a large class of optimization problems. Certainly CBO is a significantly simpler mechanism with respect to more sophisticated metaheuristics, which can include different features including memory of past exploration. Nevertheless, it seems to be powerful and robust enough to tackle many interesting nonconvex optimizations [15, 29, 33], which would be harder to solve by gradient descent methods that have been dominating the field of optimization, especially in relevant applications such as machine learning. In fact in many of these problems the objective function is not differentiable and CBO do not use any gradient information for their exploration. Moreover, in certain models, such as training of deep neural networks, the gradient tends to explode or vanish [9]. Finally, although there exist situations where global optimization by gradient descent methods can be ensured under ad hoc conditions, see, e.g. Refs. [16,44,45], they do not offer in general guarantees of global convergence. Instead, CBO has potential of a rather general and rigorous global asymptotic/convergence analysis. Some theoretical gaps remain open in the analysis of CBO though, in particular the rigorous derivation of the mean-field limit, due to the difficulty in establishing bounds on the moments of the probability distribution of the particles [13]. Because of this lack of compactness, a direct proof of convergence of the stochastic particle method has been recently derived in Ref. [33], which does not require the mean-field limit, but at the same time it is not able to provide a quantitative error estimate with respect to the number of particles. Motivated by these theoretical gaps and several potential applications in machine learning, the main task of the present work is to develop a CBO approach to solve the following constrained optimization problem argmin(1.1)

where 0 is a given continuous cost function, which we wish to minimize over a hypersurface Γ. Here we assume that Γ is a connected smooth compact hypersurface embedded in , which is represented as the 0-level set of a signed distance function with = dist(v, Γ). This means that

If Γ = , we also assume for simplicity that 0 on the interior of Γ and 0 on the exterior. The gradient is then the outward unit normal on Γ wherever is defined. Moreover, we assume that there exists an open neighborhood Γ of Γ such that Γ). Such setting of hypersurface Γ has been used in Ref. [22].

Example 1.1. Examples of hypersurfaces Γ in this setting are

In particular, we consider a system of N interacting particles (()satisfying the following stochastic differential dynamics expressed in Itˆo’s form

= )())))2 ())))(1.2)

where 0 is a suitable drift parameter, 0 a diffusion parameter,

is the empirical measure of the particles (is the Dirac measure at ), and

For simplicity, we write to mean for any vector , and |v| is the standard Euclidean norm. This stochastic system is considered complemented with independent and identical distributed (i.i.d.) initial data Γ with i = 1, and the common law is denoted by (Γ). The trajectories ((denote N independent standard Brownian motions in . In (1.2) the projection operator ) is defined by

In the case of sphere we know ) = , so one has P(v) = , and (1.2) corresponds to a stochastic Kuramoto-Vicsek type model [12,43,53].

The choice of the weight function in (1.4) comes from the well-known Laplace’s principle [21,46,49], a classical asymptotic method for integrals, which states that for any probability measure (Γ), it holds

Let us discuss the mechanism of the dynamics. The right-hand-side of the equation (1.2) is made of three terms. The first deterministic term )())dt imposes a drift to the dynamics towards , which is the current consensus point at time t as approximation to the global minimizer. This drift term is the projection onto the tangent space to the hypersurfaces at of the vector . Hence gets drifted in the direction of , which is not normal to the hypersurface, and the term disappears when = . In fact, the consensus point is explicitly computed as in (1.4) and it may lay in general outside Γ. This choice of an embedded weighted barycenter is very simple, compatible with fast computations, and, for compact manifolds Γ, it is a good proxy for a minimizer . One could alternatively consider the computation of a weighted barycenter on the manifold

where is a (Riemannian) distance on Γ, 1, and is the particle distribution. However, the computation of is in general not explicit (because the distance function might not be simple or explicit) and one may have to solve at each time t an optimization problem over the manifold in order to compute . In view of the potentially variable curvature of Γ, the distance ) could locally behave in general as a geodesically convex or nonconvex function, rendering such a global optimization rather difficult in general. These are reasons for choosing the embedded alternative (1.4). The second stochastic term ))introduces a random decision to favor exploration, whose variance is a function of the distance of the particles to the current consensus points. In particular, as soon as consensus is reached, then the stochastic component vanishes. The last term ())))dt combined with ) is needed to ensure that the dynamics stays on the hypersurface despite the Brownian motion component. Intuitively, it simply compensates the normal direction of the noise produced by the isotropic Brownian motion and leaves the tangential component active. In fact, we will initially analyze the well-posedness of the system (1.2) as defined in the whole space instead of being immediately considered constrained on the hypersurface Γ, and we will check afterwards that, if the initial data ()is set on the hypersurface, then the particles (()remain there at all times. We further notice that the dynamics does not make use of any derivative of E, but only of its pointwise evaluations, which appear integrated in (1.4). Hence, the equation can be in principle numerically implemented at discrete times also for cost functions E which are just continuous and with no further smoothness. We require below more regularity to E exclusively to ensure formal well-posedness of the evolution.

The main results of this paper are about the well-posedness of (1.2) and its rigorous mean-field limit - which is an open issue for unconstrained CBO [13] - to the following nonlocal, nonlinear Fokker-Planck equation

with the initial data (Γ). Here (Γ) is a Borel probability measure on Γ and

The operators and ∆denote the divergence and Laplace-Beltrami operator on the hypersurface Γ respectively. The mean-field limit will be achieved through the coupling method [11,26,39,40,52] by introducing an auxiliary mono-particle process, satisfying the self-consistent nonlinear SDE

with the initial data distributed according to (Γ). Here we require = law(V), and we will show that as a measure concentrated on Γ solves the PDE (1.7). We call the SDE (1.8) the mean-field dynamic, and the PDE (1.7) is called the mean-field PDE. In our follow up paper Ref. [29], for more regular datum , we also prove existence and uniqueness of distributional solutions ([0(Γ)) at any finite time T > 0. This model and the analysis we provide in this paper are very much inspired by the homogeneous version of the kinetic Kolmogorov-Kuramoto-Vicsek model on the sphere Γ = , which was formally derived by Degond and Motsch [20] as a mean-field limit of the discrete Vicsek model [4, 17, 53]. Recently, the stochastic Vicsek model has received extensive attention in the mathematical community to establish its mean-field limit, hydrodynamic limit, and phase transitions. Bolley, Ca˜nizo and Carrillo [12] have rigorously justified the mean-field limit in case of smoothed force field, and Gamba and Kang [30] and Figalli, Kang, and Morales [27] extended the global existence and uniqueness of weak solutions of the kinetic Kolmogorov-Kuramoto-Vicsek model to the case of singular mean-field force field. Recently, Degond, Frouvelle and Liu [19] provided also a complete and rigorous description of phase

Figure 1: The Ackley function for d = 3 in the constrained case on two example hypersurfaces. The half sphere (left) and the torus with R = 1 and r = 0.5 (right). The global minimum on corresponds to the point = (0, 0, 1), over the torus to the point = (0.

transitions such as the number and nature of equilibria, stability, convergence rate, phase diagram and hysteresis. Our results build very much upon these developments.

Before starting our analysis of (1.2), (1.8), and (1.7), let us illustrate numerically the behavior of the dynamics in the case of the minimum solution of the Ackley function for d = 3 constrained over an hypersurface Γ

with A = 20, a = 0.2, b = 3, B = 20 and v = (with ) = 0.

The global minimum is at . In Figure 1 we report in the case d = 3 the Ackley function constrained over the half sphere 0 and over the torus with external radius R = 1 and inner radius r = 0.5. The colors over the surface represent the values attained by the function.

In the reported simulations we initialize the particles with a uniform distribution over the hypersurface (this typically requires the use of a suitable algorithm, see Ref. [29]) and employ a simple Euler-Maruyama scheme with projection and time step ∆t > 0. We report in Figure 2 the case of the half-sphere together with the particle trajectories for [0, 5] with N = 20, ∆t = 0.05, = 0.25 and = 50. On the left we consider the case with minimum at = (0, 0, 1), on the right the case with minimum at = (1. The time evolution of the particle distribution ) in the numerical mean field limit for N = 10is also reported in the upper part of the same figure.

In Figure 3 we report an analogous simulation in the constrained case on the torus. The numerical parameters are the same as in the previous case. On the left the minimum is located at = (0, whereas on the right is located at = (0.5, 0, 0).

As these figures show it, the dynamics consistently converges to the global minimizers totally regardless of the many local minimizers and the large particle system does converges numerically to its mean-field approximation. In the related paper Ref. [29] we rigorously analyze the large time behavior of the system for Γ = , we prove its convergence to global minimizers in a suitable sense, and we provide several applications in machine learning. A proof of convergence to global minimizers for more general hypersurfaces is the subject of current investigation.

Figure 2: Particles trajectories along the simulation for the Ackley function on the half-sphere in the case d = 3, N = 20 with minimum at = (0, 0, 1)(left) and = ((right). On the top corresponding time evolution of the particle distribution ) in angular coordinates at t = 1 and t = 2.5 for N = 10. The simulation parameters are ∆t = 0.05, = 0.25 and = 50. Although we run the simulation on the entire sphere, we display the result on the half-sphere because, eventually, all particles do converge to one point, hence, all particles do eventually belong to a half sphere (mind that the stochastic process is a time continuous function). Displaying half-sphere is simply aesthetically more pleasing and it allows a better reading of the figures. On the other hand, as soon as one generates initial particles on the half-sphere containing a minimizer, then they - empirically - tend generically to remain on that half sphere, hence there is no need of displaying the rest of the sphere and we do not use any confinement mechanism rather than the given dynamics. Also boundary conditions are actually not imposed.

Let us conclude this introduction by mentioning that the optimization on the hypersurface offers further numerous advantages, besides allowing a rigorous proof of mean-field limit thanks to the compactness of the hypersurface. First of all a vast class of optimization problems can be reduced to constrained optimizations over hypersurfaces: in a related paper Ref. [29], where we focus on the numerical implementation of the method on the sphere and its asymptotic/convergence analysis to global minimizers, we present also a few relevant and challenging applications in signal processing and machine learning. Due to compactness of the hypersurface, local smoothness and boundedness requirements on E are necessarily a uniform and global property. However, against these properties that greatly simplify the analysis of the well-posedness of the system and its mean-field limit, the specific topology of the hypersurface may make it harder to prove asymptotic convergence of the dynamics to global miminizers, requiring major technical variations with respect to the approach of unconstrained CBO [13]. We refer to Ref. [29] for the details of the large time analysis of solutions ) of (1.7) on the sphere and examples of applications in high-dimension. There we needed to develop ad hoc arguments to prove the asymptotic convergence. A similar consensus dynamics could be in principle defined also on more general compact manifolds [25], but

Figure 3: Particles trajectories along the simulation for the Ackley function constrained on the torus with R = 1 and r = 0.5 in the case d = 3, N = 20 with minimum at = (0(left) and = (0.5, 0, 0)(right). The simulation parameters are ∆t = 0.05, = 0.25 and = 50.

we still lack a general approach to deal with its (global) asymptotic analysis and this research direction needs certainly further investigations.

The rest of the paper is organized as follows: in Section 2 we show that the equations (1.2), (1.7) and (1.8) are well-posed, that is, there is a unique solution for any initial distribution , which depends continuously on the initial datum, and in Section 3 we address the mean-field limit.

This section focuses on proving the well-posedness for the particle system (1.2), the mean-field dynamic (1.8) and the mean-field PDE (1.7). To ensure this, we make the following smoothness assumption on the objective function E.

Assumption 2.1. The objective function 0 is locally Lipschitz continuous.

2.1 Well-posedness for the interacting particle system (1.2)

Recall that the particle system (1.2) is assumed to evolve in the whole space instead of on the hypersurface Γ directly. We choose this embedding because it provides an explicit and computable representation of the system and it allows for a global description. The difficulty in showing first the well-posedness of (1.2) in the ambient space is that the projection ), ∆) and ) may not be well defined when Γ. For example, in the case of the sphere , the projection ) = and ∆)) = (1) are not defined for = 0. In order to overcome this problem, we regularize the diffusion and drift coefficients, that is, we replace them with appropriate functions and respectively: let be a matrix valued map on with bounded derivatives such that ) = P(v) for all be a valued map on R with bounded derivatives such that ) = ∆) if Γ, and be a valued map on , again with bounded derivatives such that ) = ) if Γ. It is also useful to mention that

hold for any . For later use, let us denote [N] = {1, . . . , N}. Additionally, we regularize the locally Lipschitz function E: Let us introduce ˜E(v) satisfying the following assumptions

Assumption 2.2. The regularized extension function ˜is globally Lipschitz continuous and satisfies the properties

1. ˜E(v) =

2. ˜˜for all for a suitable global Lipschitz constant L > 0;

˜E := inf ˜˜sup ˜E =: ˜

Remark 2.1. Here ˜E is introduced as an auxiliary function for the proof of well-posedness and mean-field limit only, and it does not play any role in actual optimization problem, which is defined on Γ. Indeed, as we can see in Theorem 2.1 below, particles stay on the hypersurface Γ all the time, which means that certainly Γ, so one has ˜).

We note here that E and ˜E appearing the sequel always satisfy Assumptions 2.1 and 2.2. Given such and ˜E, we introduce the following regularized particle system

= )())))2 ())))(2.3)

for ], where

Next we can deduce that the coefficients in (2.3) are locally Lipschitz continuous and have linear growth. More precisely, we have the following result.

Lemma 2.1. Let 0 be arbitrary and ˜E satisfy Assumption 2.2. Then for any , and corresponding empirical measures = and it holds

and

where . Here we used the notations for norms of vectors = supand = .

Proof. First, one notes that

Then we have

where := . Next we split the error

For , the previous computation yields

Notice that it holds

Thus we conclude that

Collecting estimates (2.9) and (2.10) completes the proof.

Theorem 2.1. Under Assumptions 2.1 and 2.2 let be a probability measure on Γ and, for every , ()be N i.i.d. random variables with the common law . For every , there exists a pathwise unique strong solution (()to the particle system (1.2) with the initial data (). Moreover it holds that Γ for all ] and any t > 0.

Proof. Given and ˜E, the SDE (2.3) has locally Lipschitz coefficients, so it admits a pathwise unique local strong solution by standard SDE well-posedness result [24, Chap. 5, Theorem 3.1]. Moreover, it follows from Itˆo’s formula (5.1) (choosing ) = )), as long as Γ (the process is continuous and we can nevertheless consider a smooth extension of out of Γ), that

) = ) )()))) )12)))dt

where ) is the k-th component of ), is the unit vector along the kth dimension, and we have used the property of ) as in (2.2). Hence ) = ) = 0 for all t > 0, which ensures that the solution keeps bounded at finite time, hence we have a global solution. Since all Γ, the solution to the regularized system (2.3) is a solution to (1.2), which provides the global existence of solutions to (1.2).

To show pathwise uniqueness let us consider two solutions to (1.2) for the same initial distribution and Brownian motion. According to the above computation these two solutions stay on the hypersurface for any 0, hence they are solutions to the regularized system (2.3), whose solutions are pathwise unique due to the locally Lipschitz continuous coefficients. Hence we have uniqueness for solutions to (1.2).

2.2 Well-posedness for the mean-field dynamic (1.8)

For readers’ convenience, we give a brief introduction of the Wasserstein metric in the following definition, we refer to Ref. [5] for more details.

Definition 2.1 (Wasserstein Metric). Let 1 and ) be the space of Borel probability measures on with finite p-moment. We equip this space with the Wasserstein distance

where Π() denotes the collection of all Borel probability measures on with marginals and in the first and second component respectively. The Wasserstein distance can also be expressed as

where the infimum is taken over all joint distributions of the random variables Z, Z with marginals respectively.

Notice now that, for any )

A direct application of above leads to

p < C = C(C, α, L) > 0.

Proof. Let us compute the difference

where Π() is an arbitrary coupling of and . We can write the integrand as follows

where the last two terms on the right hand side can be rewritten as

Under Assumption 2.2, one can obtain the bounds

and

Similar to the estimate of , we estimate

Notice that ), so for any 1 one has

Collecting the above estimates, we obtain

where C depends only on and . Lastly, optimizing over all couplings yields (2.16).

The following theorem states the well-posedness for the mean-field dynamic (1.8).

Theorem 2.2. Let E and ˜E satisfy Assumptions 2.1 and 2.2. Then there exists a unique process C([0), T > 0, satisfying the nonlinear SDE (1.8)

in strong sense for any initial data Γ distributed according to (Γ), where

and = law(V) for all [0, T]. Moreover Γ for all [0, T].

Proof. The proof in the following is based on the Leray-Schauder fixed point theorem, see, e.g., [32, Chapter 10], and it will be carried out through 5 steps.

• Step 1: For some given ([0), a distribution on Γ and with law , we can uniquely solve the SDE

because the coefficients are locally Lipschitz continuous and is independent of V . In addition, following the same argument as in (2.11) we have ) = 0, so that Γ for all time. This introduces = law(V) and ([0)). Setting := ([0) we define the map

which we will show to be compact in the next step.

• Step 2: For any 0 , one has

Then applying Itˆo’s isometry yields

where the constants only depend on . Therefore, by defi-nition of the Wasserstein distance we have

By an application of Lemma 2.2 one obtains

This provides the H¨older continuity of ). Thus one has T (C([0([0C([0), which implies the compactness of the map T .

• Step 3: Let us define the set

For , there exists some satisfying (2.22) with its law ([0)) such that ). Due to (2.15), we have that for any [0, T]

since has a compact support, which verifies the boundedness of the set A. Applying the Leray-Schauder fixed point theorem, there exists a fixed point for the mapping T and thereby a solution of

• Step 4: In this step, we shall prove the uniqueness. Suppose we have two fixed points and , and their corresponding process satisfying (2.22) respectively, denote by and their laws. We compute the difference := and applying Itˆo’s isometry yields

where C depends on . According to Lemma 2.2, one has

which implies that

Therefore, applying Gronwall’s inequality with ] = 0, one obtain ] = 0 for all [0, T]. This leads to by (2.33). Hence, we obtain the uniqueness for solutions to (2.31).

• Step 5: Similar to the argument in Theorem 2.1, the unique solution to the regularized SDE (2.31) is also the unique solution to the nonlinear SDE (1.8) due to the fact that ) = 0 for all [0, T].

2.3 Well-posedness for the PDE (1.7)

Let Γ, and let be the solution to (1.8) obtained in the last section with the initial data distributed according to . For any (), it follows from Itˆo’s formula (5.1) that

where is the Hessian and A : B := Tr(). Taking expectation on both sides of (2.35), the law of as a measure on satisfies

As we have proved that Γ, almost surely, that is, the density is concentrated on Γ for any t, we have supp(Γ. Let us now define the restriction of on Γ by

is unique for . We know such exists since Γ), see for example [22, Section 2.1]. Let now Φ (Γ) and define a function () such that

Then defined above is 0-homogeneous in v in the strip Γ, so that ) = 0 for all v in the support Γ of , which leads to ) : = 0. Hence,

Let us now relate the Euclidean differential operators to corresponding operators on Γ, so that for Γ it holds ) = ) and ∆) = ∆). Therefore, we obtain

where

Thus by this construction we obtain a weak solution to the PDE (1.7).

Next we prove the uniqueness of solutions to (1.7). Assume that and are two solutions to (1.7) with the same initial data , and at each time t we treat them as measures on concentrated on the hypersurface Γ. Then we construct two linear process ((i = 1, 2) satisfying

with the common initial data distributed according to . Let us denote law(V) = ¯(i = 1, 2) as measures on , which are solutions to the following linear PDE

Since the uniqueness for the above linear PDE holds and is also a solution to the above PDE on (see (2.36)), it follows that ¯= (i = 1, 2). Consequently, the process (are solutions to the nonlinear SDE (1.8), for which the uniqueness has been obtained. Hence (and (are equal, which implies = ¯= ¯= . Thus the uniqueness is obtained.

The well-posedness of (1.2), (1.7) and (1.8) obtained in the last section provides all the ingredients we need for the mean-field limit. Let ((be N independent copies of solutions to (1.8). They are i.i.d. with the same distribution . Assume that (()is the solution to the particle system (1.2). Since Γ for all i and t, ((and (()are solutions to the corresponding regularized systems (2.31) and (2.3) respectively. We denote below = and = law(V).

Before stating our theorem on the mean-field limit, let us introduce the following lemma on a large deviation bound.

Lemma 3.1. Let E and ˜E satisfy Assumptions 2.1 and 2.2. Let ((be the solution to the mean-field dynamics (2.31), which are i.i.d. with common distribution ([0)). Then there exists a constant C depending only on Γ and such that

Proof. We start by estimating

Here we used explicitly the upper and lower bounds ˜E, ˜E. Denote

Here we have used the fact that Γ, so it holds for some 0 depending on Γ. Thus we conclude

By following the same argument it also holds that

Combining the above two estimates and (3.2), one has

Thus we have completed the proof.

Then we get the following mean-field limit result by the classical Sznitman’s theory.

holds for all 0 .

Remark 3.1. The estimate above guarantees the weak convergence of the empirical measure towards , in the following sense

Indeed, one has

We refer to Refs. [11,52] for more details.

Proof. Notice that ((and (()are also solutions to the corresponding regularized systems (2.31) and (2.3) respectively, and

By applying Itˆo’s formula one has

where () represents the k-th row of the matrix ). Using Lemma 2.1 it is easy to compute that )())

≤∥

C(α, , L, C, C)

Hence it yields

where we have used Cauchy’s inequality and C depends only on and , and

where we have used Lemma 2.1 and the fact that . Here C depends only on , and . This leads to

where C depends only on and . This implies that

where we have used the fact that

and Lemma 3.1 in the last inequality. Applying Gronwall’s inequality with ] = 0, one concludes

for all [0, T], which completes the proof.

Remark 3.2. The constant C > 0 appearing in the estimate (3.9) may depend on the dimension through the norm of or , see (3.14). Nevertheless we expect this dependency to scale at most linearly as 1. In fact, for the case of sphere Γ = it is ) = ∆) = . We conclude that in general there is no curse of dimensionality involved in the estimate (3.9).

We presented a new consensus-based model for global optimization on hypersurfaces, which is inspired by the kinetic Kolmogorov-Kuramoto-Vicsek equation on the sphere. The main results of this paper are the well-posedness of the stochastic particle system and its mean-field limit, which is obtained by the coupling method through introducing an auxiliary self-consistent nonlinear SDE. The presented mean-field limit is not affected by curse of dimension, i.e., the rate of convergence is of order in the particle number N. This favorable theoretical rate is confirmed by numerical experiments in very high dimension (3000) [29]. In the related paper Ref. [29] we analyze the large time behavior of the system for Γ = , we prove its convergence to global minimizers in a suitable sense, and we provide several applications in machine learning. We consider the results of this paper as a preliminary step towards the formulation of consensus-based optimization on compact manifolds [25], which requires a more general approach for the analysis of the global large time asymptotics.

Theorem 5.1 (Multidimensional Itˆo’s formula). Let

be an d-dimensional Itˆo process, where , and is -dimensional Brownian motion. Assume ) be a map from to R, then it holds that

with ) being the i-th row of the matrix v(s).

Massimo Fornasier and Hui Huang acknowledge the support of the DFG Project ”Identification of Energies from Observation of Evolutions” and the DFG SPP 1962 ”Non-smooth and Complementarity-based Distributed Parameter Systems: Simulation and Hierarchical Optimization”. The present project and Philippe S¨unnen are supported by the National Research Fund, Luxembourg (AFR PhD Project Idea “Mathematical Analysis of Training Neural Networks” 12434809). Lorenzo Pareschi acknowledges the support of the John Von Neumann guest Professorship program of the Technical University of Munich and the PRIN Project 2017KKJP4X ”Innovative numerical methods for evolutionary partial differential equations and applications”.

[1] Emile Aarts and Jan Korst. Simulated Annealing and Boltzmann Machines: A Stochastic Approach to Combinatorial Optimization and Neural Computing. John Wiley & Sons, Inc., New York, NY, USA, 1989.

[2] G. Albi, N. Bellomo, L. Fermo, S.-Y. Ha, J. Kim, L. Pareschi, D. Poyato, and J. Soler. Vehicular traffic, crowds, and swarms: From kinetic theory and multiscale methods to applications and research perspectives. Mathematical Models and Methods in Applied Sciences, 29(10):1901–2005, 2019.

[3] Giacomo Albi, Lorenzo Pareschi, Giuseppe Toscani, and Mattia Zanella. Recent advances in opinion modeling: control and social influence, pages 49–98. Springer, 2017.

[4] Maximino Aldana and Cristian Huepe. Phase transitions in self-driven many-particle systems and related non-equilibrium models: A network approach. Journal of Statistical Physics, 112:135–153, 07 2003.

[5] Luigi Ambrosio, Nicola Gigli, and Giuseppe Savar´e. Gradient Flows: In Metric Spaces and in the Space of Probability Measures. Springer Science & Business Media, 2008.

[6] Thomas Back, David B. Fogel, and Zbigniew Michalewicz, editors. Handbook of Evolutionary Computation. IOP Publishing Ltd., Bristol, UK, UK, 1st edition, 1997.

[7] Nicola Bellomo and Christian Dogbe. On the modeling of traffic and crowds: A survey of models, speculations, and perspectives. SIAM review, 53(3):409–463, 2011.

[8] Nicola Bellomo, Benedetto Piccoli, and Andrea Tosin. Modeling crowd dynamics from a complex system viewpoint. Mathematical Models and Methods in Applied Sciences, 22(supp02):1230004, 2012.

[9] Yoshua Bengio, Patrice Simard, Paolo Frasconi, et al. Learning long-term dependencies with gradient descent is difficult. IEEE Transactions on Neural Networks, 5(2):157–166, 1994.

[10] Christian Blum and Andrea Roli. Metaheuristics in combinatorial optimization: Overview and conceptual comparison. ACM computing surveys (CSUR), 35(3):268–308, 2003.

[11] Fran¸cois Bolley, Jos´e A Canizo, and Jos´e A Carrillo. Stochastic mean-field limit: non-Lipschitz forces and swarming. Mathematical Models and Methods in Applied Sciences, 21(11):2179–2210, 2011.

[12] Fran¸cois Bolley, Jos´e A Ca˜nizo, and Jos´e A Carrillo. Mean-field limit for the stochastic vicsek model. Applied Mathematics Letters, 25(3):339–343, 2012.

[13] Jos´e A Carrillo, Young-Pil Choi, Claudia Totzeck, and Oliver Tse. An analytical framework for consensus-based global optimization method. Mathematical Models and Methods in Applied Sciences, 28(06):1037–1066, 2018.

[14] Jos´e A Carrillo, Massimo Fornasier, Jes´us Rosado, and Giuseppe Toscani. Asymptotic flocking dynamics for the kinetic Cucker–Smale model. SIAM Journal on Mathematical Analysis, 42(1):218– 236, 2010.

[15] Jos´e A Carrillo, Shi Jin, Lei Li, and Yuhua Zhu. A consensus-based global optimization method for high dimensional machine learning problems. ESAIM: Control, Optimisation and Calculus of Variations, 2020.

[16] Yuxin Chen, Yuejie Chi, Jianqing Fan, and Cong Ma. Gradient descent with random initialization: fast global convergence for nonconvex phase retrieval. Mathematical Programming, 176(1):5–37, 2019.

[17] ID Couzin, J Krause, R James, GD Ruxton, and NR Franks. Collective memory and spatial sorting in animal groups. Journal of Theoretical Biology, 218:1 – 11, 2002.

[18] Felipe Cucker and Steve Smale. Emergent behavior in flocks. IEEE Transactions on automatic control, 52(5):852–862, 2007.

[19] Pierre Degond, Amic Frouvelle, and Jian-Guo Liu. Phase transitions, hysteresis, and hyperbolicity for self-organized alignment dynamics. Archive for Rational Mechanics and Analysis, 216(1):63–115, 2015.

[20] Pierre Degond and S´ebastien Motsch. Continuum limit of self-driven particles with orientation interaction. Mathematical Models and Methods in Applied Sciences, 18(supp01):1193–1215, 2008.

[21] Amir Dembo and Ofer Zeitouni. Large Deviations Techniques and Applications. Springer-Verlag Berlin Heidelberg, 2010.

[22] Alan Demlow and Gerhard Dziuk. An adaptive finite element method for the Laplace-Beltrami operator on implicitly defined surfaces. SIAM Journal on Numerical Analysis, 45(1):421–442, 2007.

[23] Marco Dorigo and Christian Blum. Ant colony optimization theory: A survey. Theoretical computer science, 344(2-3):243–278, 2005.

[24] Richard Durrett. Stochastic calculus: a practical introduction. CRC press, 2018.

[25] M. ´Emery and P.A. Meyer. Stochastic calculus in manifolds. Universitext (1979). Springer-Verlag, 1989.

[26] Razvan C Fetecau, Hui Huang, and Weiran Sun. Propagation of chaos for the Keller–Segel equation over bounded domains. Journal of Differential Equations, 266(4):2142–2174, 2019.

[27] Alessio Figalli, Moon-Jin Kang, and Javier Morales. Global well-posedness of the spatially homoge- neous Kolmogorov–Vicsek model as a gradient flow. Archive for Rational Mechanics and Analysis, 227(3):869–896, 2018.

[28] David B. Fogel. Evolutionary Computation: Toward a New Philosophy of Machine Intelligence (IEEE Press Series on Computational Intelligence). Wiley-IEEE Press, 2006.

[29] Massimo Fornasier, Hui Huang, Lorenzo Pareschi, and Philippe S¨unnen. Consensus-based optimization on the sphere: Convergence to global mininizers and machine learning. arXiv:2001.11988, 2020.

[30] Irene M Gamba and Moon-Jin Kang. Global weak solutions for Kolmogorov–Vicsek type equations with orientational interactions. Archive for Rational Mechanics and Analysis, 222(1):317–342, 2016.

[31] Michel Gendreau and Jean-Yves Potvin. Handbook of Metaheuristics. Springer Publishing Company, Incorporated, 2nd edition, 2010.

[32] David Gilbarg and Neil S Trudinger. Elliptic partial differential equations of second order. springer, 2015.

[33] Seung-Yeal Ha, Shi Jin, and Doheon Kim. Convergence of a first-order consensus-based global optimization algorithm. Mathematical Models and Methods in Applied Sciences, 2020.

[34] W. K. Hastings. Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1):97–109, 1970.

[35] Dirk Helbing. Quantitative Sociodynamics: Stochastic Methods and Models of Social Interaction Processes. Springer Science & Business Media, 2010.

[36] John H. Holland. Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control and Artificial Intelligence. MIT Press, Cambridge, MA, USA, 1992.

[37] Richard Holley and Daniel Stroock. Simulated annealing via Sobolev inequalities. Communications in Mathematical Physics, 115(4):553–569, 1988.

[38] Darryl D Holm and Vakhtang Putkaradze. Formation of clumps and patches in self-aggregation of finite-size particles. Physica D: Nonlinear Phenomena, 220(2):183–196, 2006.

[39] Hui Huang and Jian-Guo Liu. Error estimate of a random particle blob method for the Keller–Segel equation. Mathematics of Computation, 86(308):2719–2744, 2017.

[40] Hui Huang, Jian-Guo Liu, and Jianfeng Lu. Learning interacting particle systems: Diffusion param- eter estimation for aggregation equations. Mathematical Models and Methods in Applied Sciences, 29(01):1–29, 2019.

[41] James Kennedy. Particle swarm optimization. Encyclopedia of machine learning, pages 760–766, 2010.

[42] Scott Kirkpatrick, C Daniel Gelatt, and Mario P Vecchi. Optimization by simulated annealing. Science, 220(4598):671–680, 1983.

[43] Yoshiki Kuramoto. Self-entrainment of a population of coupled non-linear oscillators. In International symposium on mathematical problems in theoretical physics, pages 420–422. Springer, 1975.

[44] Jason D. Lee, Ioannis Panageas, Georgios Piliouras, Max Simchowitz, Michael I. Jordan, and Ben- jamin Recht. First-order methods almost always avoid strict saddle points. Mathematical Programming, 176(1):311–337, 2019.

[45] Shengchao Liu, Dimitris Papailiopoulos, and Dimitris Achlioptas. Bad global minima exist and sgd can reach them. arXiv:1906.02613, 2019.

[46] Peter David Miller. Applied Asymptotic Analysis, volume 75. American Mathematical Soc., 2006.

[47] Sebastien Motsch and Eitan Tadmor. Heterophilious dynamics enhances consensus. SIAM review, 56(4):577–621, 2014.

[48] John A. Nelder and Roger Mead. A simplex method for function minimization. Computer Journal, 7:308–313, 1965.

[49] Ren´e Pinnau, Claudia Totzeck, Oliver Tse, and Stephan Martin. A consensus-based model for global optimization and its mean-field limit. Mathematical Models and Methods in Applied Sciences, 27(01):183–204, 2017.

[50] Riccardo Poli, James Kennedy, and Tim Blackwell. Particle swarm optimization. Swarm intelligence, 1(1):33–57, 2007.

[51] L. A. Rastrigin. The convergence of the random search method in the external control of many- parameter system. Automation and Remote Control, 24:1337–1342, 1963.

[52] Alain-Sol Sznitman. Topics in propagation of chaos. In XIX—1989, pages 165–251. Springer, 1991.

[53] Tam´as Vicsek, Andr´as Czir´ok, Eshel Ben-Jacob, Inon Cohen, and Ofer Shochet. Novel type of phase transition in a system of self-driven particles. Physical review letters, 75(6):1226, 1995.