{'fraction_non_alphanumeric': 0.03696077114761926, 'fraction_numerical': 0.011295641404419774, 'mean_word_length': 4.276176980003149, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On stochastic mirror descent with interacting particles: convergence properties and variance reduction

An open problem in optimization with noisy information is the computation of an exact minimizer that is independent of the amount of noise. A standard practice in stochastic approximation algorithms is to use a decreasing step-size. However, to converge the step-size must decrease exponentially slow, and therefore this approach is not useful in practice. A second alternative is to use a fixed step-size and run independent replicas of the algorithm and average these. A third option is to run replicas of the algorithm and allow them to interact. It is unclear which of these options works best. To address this question, we reduce the problem of the computation of an exact minimizer with noisy gradient information to the study of stochastic mirror descent with interacting particles. We study the convergence of stochastic mirror descent and make explicit the tradeoffs between communication and variance reduction. We provide theoretical and numerical evidence to suggest that interaction helps to improve convergence and reduce the variance of the estimate.


Introduction
Optimization models that arise in artificial intelligence and statistical learning applications often include noisy estimates of the function and its gradient. This is the case when for example the gradient is computed over a subset (or mini-batch) of the data. In such a situation it is known that the optimization algorithm will converge to a neighborhood of the minimizer [12]. The size of the neighborhood depends on the amount of noise. In addition, for constrained optimization problems noise can violate the constraints making the situation even more complex.
In various applications it can be beneficial to be able to control the fluctuations around the true minimum. The conventional way to control the error is to decrease the step size. Theoretical analysis suggests step sizes which are slow in practice. An alternative is to use heuristics such as to increase the batch size over time (see e.g. [37]); this however increases the computational costs and is difficult to tune. Another option is to run independent replicas of the algorithm. We will refer to each of these runs as a particle. The question we address in this paper is whether it is beneficial to allow these particles to interact with each other. We study this question using the general framework of Stochastic Mirror Descent (SMD) (see [26]). SMD can be used to solve constrained and unconstrained problems, and is known to be an optimal algorithm for certain classes of optimization problems ( [4]).
The aim of this paper is to demonstrate the advantages of using tthe interacting particle system, in comparison to independent copies of SMD. We will show that interaction reduces the variance of the estimation of x * and can improve the convergence properties including the convergence rate to the stationary distribution. In particular, we show that for a (strongly) convex objective the distance to the optimum is bounded by a term related to the standard optimization error in which the noise variance is reduced by a factor of N and an interaction term measuring the deviation of each particle from the system average. This latter term is bounded, so that under certain assumptions consensus and consequently convergence to the optimum is achieved. Using logarithmic Sobolev inequalities we can furthermore improve on this rate and show that the convergence to a stationary distribution can be achieved at an exponential rate. In addition, we believe that the ability to impose communication constraints through the interaction A will be beneficial also from a practical point of view as it will decrease communication costs in a parallel implementation. In this paper we will work with the continuous time formulation for SMD and its proposed interacting version. The main reason behind this is that the analysis for continuous time dynamics provides a clear and complete picture for the benefit of interaction. In our numerical experiments we used a simple forward Euler discretization scheme and we expect that the results can only improve with more sophisticated schemes. A detailed error analysis of the discretization scheme is beyond the scope of this paper.

Related Literature
It is a well-known problem that stochastic optimization algorithms converge to a neighborhood of the optimal solution. The size of this neighborhood is proportional to the noise variance (or the second moment of the sample gradient). Traditionally, common strategies for mitigating this include using a decreasing stepsize or attaining vanishing noise variance by increasing the batch size. Both of these come at additional computational costs. In the context of SMD, the authors of [23] study the convergence of SMD under the assumption that either the variance of the noise decreases over time or the step size is reduced exponentially slow. Various strategies have been also proposed to mitigate the effect of this noise variance to enable convergence closer to the optimum. The effect of sampling strategies such as importance sampling [25] can decrease the effect of the noise variance. Various variance reduction methods have been proposed, e.g. [15], [9], [11]. The work of [17] and [41] studied an accelerated version of SMD, however the distance to the optimum remains a function of the noise.
In this work we propose to use interacting particles to achieve a similar variance reduction. The work most closely related to ours is that of [31] in continuous time and [10] in discrete time. The authors in [31] study continuous and interacting SMD however only attained convergence rates for linear dynamics in the mirror domain and did not present any numerical results. The work of [10] covers the deterministic and noisy gradient setting in discrete time in a distributed setting; our work is in continuous time with all particles optimizing the same objective and uses different assumptions on the noise dynamics.
Our work also has parallels with the vast distributed optimization literature from which we list some indicative recent references: [19], [34], [16], [33]. We remark that our objective is different from that in distributed optimization. Our goal here is to study how interaction can improve the convergence properties of SMD and analyze the variance reduction effects, while in distributed optimization the task is to optimize efficiently an objective function that is distributed across different nodes. For example a distributed version of the interacting SDE above would require using a different drift −∇f i (x t t ) in the dynamics. There are also parallels with the work on consensus and synchronization, see e.g. [5], [36], [45], [43] and [44]. Our work generalizes some of these results to the SMD setting and combine it with results on convergence to the minimum. There is also parallel work using stochastic gradient descent (SGD) for constrained sampling problems considered in [14], [1], [42]. While there are many similarities, the goals of sampling and optimization are different and different dynamics are used in each case [20]; specifically, in optimization one wants to converge to the optimizer, while in sampling the objective is to converge to the correct invariant measure.

Contributions and Organization
Our main focus is on the convergence properties of stochastic mirror descent with interacting particles (ISMD). We derive regret bounds for the proximity of f (x i t ) to the optimal value using Lyapunov-based arguments. The application of Lyapunov techniques to optimization problems is an established approach. Our results have similarities with [31] and [23] bar a number of differences. We use different Lyapunov functions and analyze a general convex cost function. We show that in the case of ISMD there is a tradeoff between communication cost, i.e. how many particles interact, and variance reduction. In particular, with a fixed learning rate and non-vanishing noise variance, interaction between the particles can reduce the size of the neighborhood around the optimum to which the algorithm converges. Furthermore, we show that the particles converge to an area around the optimizer at an exponential rate using log-Sobolev inequalities and Bakry-Emery theory (see [2] and [3]). Although such results are standard for the analysis of interacting SDEs with convex potentianls (e.g. [21,39]), they have not been considered previously in the context of SMD. In addition, in our setting the particles do are not restricted to all interact simultaneously, so our results deviate from the mean field type analysis like in [21,39].
Our contributions can be summarized as follows: • We propose interaction between N particles as a way of controlling the convergence of a stochastic optimization algorithm -and in particular the distance to the minimum. In the presence of noise in the gradients interacting particles are an effective alternative to vanishing learning rates or noise variances.
• We establish that interaction leads to variance reduction. In particular under strong convexity and smoothness assumptions (to be specified later) we show that there exist positive constants K 1 and K 2 such that, where λ is the first nonzero eigenvalue of the graph Laplacian and κ is the strong convexity constant related to the mirrored objective. A similar average regret bound can be established for the convex case (and κ = 0). This result, together with the technical assumptions, proofs and extensions will be presented in Propositions 14-15 and 16-17. The distance to optimality is thus bounded by three terms: 1) a term which decays exponentially with time, 2) a noise term which decreases as the number of particles increases 3) a term which arises from the distance between the particle values which is smaller for a connected graph, a stronger interaction and a larger strong convexity constant.
• We provide explicit rates of convergence and concentration inequalities for the particle system using log Sobolev inequalities. This will be presented in Section 4.2. These results establish exponential rate of convergence in time for the law of the particle system, so that at equilibrium the samples z i t will oscillate around the optimum z * . The strong convexity constant plays a similar role as in the previous result, with a higher κ resulting in a smaller distance to the optimum.
• Finally, we show the benefits of using interaction in convergence speed and variance reduction through different numerical experiments. In a distributed optimization setup we manage to achieve comparable performance and convergence to the full batch gradient descent, when using a "mini-batch" approach that uses N particles with a mini-batch size that is 1/N times the total number of data-points or summands in f (see Figure 8 in Section 5 for more details).
The organization of this paper is as follows: In Section 2 present background material on MD and SMD together with some basic convergence results. Most of the material in this section is known, but we use this opportunity to set our framework and to relate the continuous time formulation with common discrete time implementations. We then propose ISMD in Section 3 and provide a detailed convergence analysis on how the particles approach the minimum as time increases. In Section 4 we apply the Bakry-Emery theory for the law of the corresponding particle system. In Section 5 we demonstrate the performance of ISMD in practice and demonstrate the effect of interaction on decreasing the variance in a variety of examples. Finally, we provide some concluding remarks in Section 6.

Notation
We will use the following notations: we denote by I d the d-dimensional identity matrix and 1 d the d-dimensional vector of ones. Let Diag(a) with a ∈ R d refer to a matrix with diagonal elements [a 1 , ..., a d ]. The Kronecker product is denoted by ⊗. Given an arbitrary norm || · || on R d , we will define B ||·|| := {v ∈ R d : ||v|| ≤ 1}. The dual norm || · || * is defined as ||z|| * := sup{z T v : v ∈ B ||·|| }. Denote with ||A|| 2 the spectral norm if A is a matrix. Similarly, let ||A|| F denote the Frobenius norm given by where λ i denote the eigenvalues of A. Assume that the dual norm is compatible with the spectral norm, i.e. ||Az|| * ≤ ||A|| 2 ||z|| * . Unless specified otherwise K, K etc. denote generic constants whose value may change according to context.
We will define a µ strongly convex function f (w.r.t. a norm · ) if there exists m > 0 such that for all x, y: The case m = 0 corresponds to the convex case. Note that equivalently strong convexity can be defined using Furthermore, for a Lipschitz function we let When looking at the Hessian of a convex function f we will denote with Hess(f ) λI d that Hess(f ) − λI d is positive semidefinite with λ ≥ 0.
For random variables X, Y we say X d = Y to denote equality in distribution. For a measure π and measurable function f we use π(f ) = E µ [f ] = f (x)π(dx). For any measurable space (Z, B(Z)) we use P(Z) to denote the space of all probability measures on Z and P 2 (Z) the one for finite second moments. A probability measure ν satisfies a Log-Sobolev inequality with constant C is for any smooth function f we have, where the entropy is defined as, For µ, ν ∈ P(Z) we will denote the Kullback-Leibler divergence or relative entropy as, otherwise.
The 2-Wasserstein distance is defined as where Π ∈ P 2 (Z ×Z) is the coupling of µ, ν with Π (Z, ·) = ν and Π (·, Z) = µ. We will assume that the optimization algorithm is started from a fixed deterministic point, and let E[·] refer to the expected value conditional on the initial value.

Preliminaries
We are interested in computing x * = arg min x∈X {f (x)} under the assumption of smoothness and convexity for f and X . Throughout the paper we will assume continuity and smoothness of f : We allow f and ∇f to be Lipschitz w.r.t arbitrary norms. As all norms are equivalent up to a proportionality constant, a natural choice without loss of generality is the Euclidean norm, which is also self dual, i.e. || · || * = || · || 2 (see [27,Lemma 2.1.2]). Despite this in the remainder we will distinguish between the chosen norm and its dual mainly in order to emphasize the difference when operating in the primal or mirror domain.
We will use a mirror map Φ : X → R d to convert the constrained optimization problem to an unconstrained one and adopt the following standard assumptions (e.g. [4, Assumption 9.3]): Assumption 2. Φ : X → R d is µ strongly convex and continuously differentiable.
The mirror map and its conjugate will be used to pass between the constrained and unconstrained space. Define to be the Legendre-Fenchel convex conjugate, so that Note that one has Note when Φ is strongly convex with constant µ, the conjugate correspondence theorem (see Section 5.3.1 in [4]) gives that the gradient of Φ * is Lipschitz continuous: We furthermore make the additional assumption that the conjugate of the mirror map maps R d to X : This simplifies the analysis as ∇Φ * maps directly to X and this avoids the need for projections. Extending our results without using Assumption 3 is possible by following a route similar to [23].
We will often use the notation z * for z * ∈ R n such that x * = ∇Φ * (z * ) ∈ arg min x∈X f (x).
The Bregman divergence is defined as, The Bregman divergence is meant to quantify how far a point x is from y and Φ can be thought as a distance generating function that adapts to the geometry or structure of X . Two well-known examples of the Bregman divergence are the Euclidean distance with Φ(x) = 1 2 ||x|| 2 2 , and the simplex constraint . An important property of using a mirror map with the Bregman divergence is that for some values z t and z * (x t and x * ) in unconstrained (constrained) space. In addition, the Φ-diameter of X is defined as We further let the anti-derivative be defined as, i.e. ∇V = ∇f • ∇Φ * . We note that knowledge of V here is mainly used to improve exposition, convey intuition and in practice we work only with ∇V. The main use of V itself will be to describe invariant distributions of the mirror descent SDEs as defined in Section 2.2. For the analysis, we note that Assumptions 1 and 2 imply ∇V is Lipschitz and we will use the following convexity assumptions in Section 3.3 to bound the fluctuations of the particle system and in Section 4 to establish rates of convergence. We will distinguish between the convex and strongly convex cases for V. Assumption 7. The function V is convex. Assumption 8. The function V is κstrongly convex with κ > 0.

Mirror Descent
We begin the presentation of mirror descent (MD) with some preliminaries on gradient descent. We provide a brief account of the main approaches in discrete time and later discuss the continuous time setting, which is of interest in this paper.

Projected Gradient Descent
One way of finding the solution of a constrained optimization problem is through projected gradient descent (GD), where k here denotes discrete time, η is the learning rate and we define the projection using the Euclidean norm, i.e. Π X (y) = arg min x∈X ||y − x|| 2 2 . The main drawback of this method is its slow convergence: when d, the dimension of x, increases the convergence can slow down. This is due to the fact that this scheme is tied to the Euclidean geometry of R d through the projection operator Π X . Consider a setting in which ||∇f (x)|| ∞ ≤ 1. This implies ||∇f (x)|| 2 < √ d, so projected gradient descent will converge at a rate d/t (see e.g. Example 9.17 from [4]).

Mirror Descent in Discrete and Continuous Time
A generalization of the projected gradient descent method is mirror descent, which is given by: The last step can be rewritten as, The current point x k is thus mapped into the dual space ∇Φ(x k ); this is then updated to ∇Φ(y k+1 ) by stepping in the direction of the negative gradient −η∇f (x k ) and then mapped back to x k+1 . Observe that for Φ(x) = 1 2 ||x|| 2 2 , we have ∇Φ(x) = x and the Bregman divergence is simply given by D Φ (x, y) = 1 2 ||x − y|| 2 2 , so that the algorithm is equivalent to projected gradient descent. Consider Φ(x) = n i=1 x(i) log x(i) and let X := ∆ d , the d-dimensional simplex. Observe that ∇ i Φ(x) = log x(i) + 1, and the projection of y with respect to this Bregman divergence onto the simplex ∆ d amounts to a scaling with the L 1 -norm. The benefit of MD is that it can adapt to the structure or geometry of X and hence result in a much better dependence of the rate of convergence on the problem dimension d (see Chapter 9 in [4]).
Using Assumption 3 the continuous time version of mirror descent is given by, where we have set η = 1 for simplicity.

Convergence of MD in continuous time
We will proceed with presenting two well-known convergence results for deterministic mirror descent and the dynamics of (5); see [8], [26]) for more details. These results are based on using Bregman divergence in a Lyapunov function V t and are included to ease exposition and to facilitate the comparison to the stochastic case. Note that D Φ,X is defined in (2).
Lemma 9 (Convergence of MD for a convex objective). Assume that Assumptions 1 -3 and 5 hold. We then have the following result, Proof. Consider, Observe that V (z * ) = 0 and V (z t ) ≥ 0 for z * = z t . Furthermore, where in the second equality we used x t = ∇Φ * (z t ) and in last inequality we have used the convexity of f . Integrating dV (z t )/dt from t = 0 to t = T , using the convexity of f and the relationship and the claim follows.
In the strongly convex case the convergence speed can be increased to an exponential rate. Lemma 10 (Convergence of MD for a strongly convex objective). Assume that Assumptions 1 -3 and 6 hold. We then have the following result, Proof. Since f is µ-strongly convex with respect to Φ, it holds, Observe that, where in the second equality we have used the definition of the Bregman divergence, in the second-to-last inequality we have used (6) and the fact that , and in the last inequality we have used the convexity of f , Then, where in the second-to-last inequality we have used (7). Then, since dV (t,zt) This implies,

Stochastic mirror descent
Consider adding noise to (4): In the limit of the time step ↓ 0 one recovers the following Itô stochastic integral equation [30]: Here B t is a d-dimensional Brownian motion and a local time construction on the boundary of X is used for u t to ensure that X t ∈ X almost surely; see [30,38] and the references therein for details. Understanding in detail the ergodicity properties or long time behaviour of (8) is largely unexplored in the literature. We note that at the level of the evolution of densities, reflecting boundary conditions are commonly used for the forward or Fokker-Planck equation, e.g. [29, Section 4.2.2], but a detailed study of (8) along these lines goes beyond the scope of this article. In the remainder we will consider the mapping to be on the mirror domain.
The stochastic mirror descent (SMD) is given by the SDE, where x t = ∇Φ * (z t ) and as before in (5) we use η = 1. Note that when σ = 0 we obtain a deterministic variant of continuous time mirror descent. The long time behavior of such dynamics are very well understood. The law of z t will converge to an invariant distribution with density proportional to exp(− 2 σ 2 V) where the antiderivative of ∇f • ∇Φ * is defined in (3). We will discuss this in more detail later in Section 4. The downside of this dynamics is that SMD cannot converge to the exact solution, i.e. the noise keeps the algorithm from fully converging exactly to x * (see [31], [23]). The latter would require the equilibrium distribution being a Dirac/atomic measure, which in turn requires using a decreasing noise in time; see [23] for more details.
Another interesting remark is that the dynamics in the dual space is then equivalent to that of SGD with objective function V. It is well-known that for a convex objective SGD methods converge linearly in time to a neighbourhood of the optimum the size of which is proportional to the noise variance [12]. A similar result can be shown to hold for SMD (see e.g. [31], [40], [23]).

Convergence of SMD
We proceed with presenting the analogs of Lemmata 9-10 for the SMD case. Proposition 11 (Convergence of SMD for a convex objective). Assume that Assumptions 1 -3 and 5 hold. We then have the following result, Proof. We define the following Lyapunov function, We then have the following results. By Itô's lemma we obtain, where for the last term we have σ( Integrating the expression we obtain, Rearranging we then obtain, Using the fact that the process Φ,X is smaller than D ||·|| 2 2 ,X then MD achieves a faster convergence rate than gradient descent. Proposition 12 (Convergence of SMD for a strongly convex objective). Assume that Assumptions 1 -3 and 6 hold. Then, We have through Itô's lemma, Integrating the expression we obtain, Observe again that for a µ-strongly convex function D f (x * , x t ) ≥ µD Φ (x * , x t ) and by the properties of the mir- Taking expected values we obtain, where we have used that the expected value of an Itô integral is zero and the relationship Again in the strongly convex case the convergence speed can be increased. However, as expected for both the convex and strongly convex results, as T → ∞, the gap to optimality is bounded from above by a quantity proportional to the noise variance σ 2 . The solution will not converge exactly to the minimum, but oscillate around it.

Interacting Stochastic Mirror Descent
It is clear that when using a fixed learning rate and a constant noise variance, SGD does not converge to the optimum. In fact, the distance to the optimum is controlled by the amount of noise. Various strategies such as using a vanishing noise variance by increasing the batch size or variance-reduced SGD (see e.g. [15], [9], [11]) have been proposed. In this work, we consider an alternative approach for controlling the distance from the optimum; namely by using interactions between the particles [31].
We consider as an alternative to (9) and we consider the following interacting particle dynamics, where B i t are independent Brownian motions and A = {A ij } N i,j=1 is an N × N doubly-stochastic matrix representing the interaction weights. As mentioned in Section 1 the matrix A represents an interaction graph which can also impose communication constraints on the particles. Each particle i will be influenced only from particles j for which A ij = 0.
The discretized version of (10) is then, where ζ t ∼ N (0, I d ) and is the discretization parameter. In the absense of interactions, i.e. when A ij = 0, we obtain a discretized version of SMD.
We remark here that the interacting mirror descent algorithm is equivalent to an algorithm in which the particle interaction is defined using a Bregman divergence instead of the L 2 distance. In particular for the time discretized version (11)- (12) we have the following result. Lemma 13 (Interacting mirror descent as Bregman distance interaction in primal space). Let σ = 0 and = 1, then interacting mirror descent in (11)-(12) can be rewritten as: Proof. We have, where in the last equality we have used ∇Φ(∇Φ * (x)) = x. Continuing, where in the first equality we have used the evolution of z i t+1 (note: the z i t terms cancel out due to the double stochasticity of A ij ). Then, where in the first equality we have used z j t = ∇Φ(x j t ), and -due to the minimum being taken over x -replaced Φ(x i t+1 ) with Φ(x j t ) and added the term ∇Φ(x j t ) T x j t and in the last equality we have used the definition of the Bregman divergence.
Unlike the standard L 2 consensus algorithm where here the particles interact in the Bregman distance. Lemma 13 can be combined with the analysis in [30] to extend (8) to an interacting reflected SDE like (10) x j t )dt as an interaction term. This is useful for conveying intuition, but the analysis in the mirror domain is much simpler. One does not have to deal with terms like u t in (8), and the interaction term in (10) is symmetric, which simplifies the analysis, see [21] for a more detailed discussion.

A useful reparameterization
It will often be useful to represent the dynamics of the particle system presented in (10) as a vector SDE with all particles stacked in a single vector variable We define the graph Laplacian as L := Diag(A1 N ) − A, and let L := L ⊗ I d , where ⊗ is the Kronecker product. Using the Laplacian, we can rewrite the evolution of the z i t -s in vector form as, Note that the interaction is linear and the distinctive properties of ISMD compared to N independent copies of SMD are contained L. Note that given A is doubly stochastic we have so the smallest eigenvalue of L is zero and L1 dN = 0. We will return to this point in Section 4. We assume throughout that the network graph corresponding to the Laplacian is connected, which in turn implies that L has eigenvalues ( [24]): We thus denote by λ the smallest non-zero eigenvalue of the Laplacian.

Convergence of ISMD: a general bound
As we saw earlier, the SMD dynamics does not converge to the minimum due to the noise in the optimization algorithm.
Here we show that introducing interaction between the particles can reduce the effect of noise.
We will decompose each particle as a sum of the particle average and a fluctuation term, z i t :=z N t +z i t , where we let and also define y N t := ∇Φ * (z N t ). We begin by deriving the evolution of the average of the particlesz N t . Observe, where we have used the fact that 1 Proposition 14 (Convergence of ISMD for a convex objective). Assume that Assumptions 1 -3 and 5 hold. We then have, Proof. Observe that, where we have used the L-Lipschitz continuity of f and the conjugate correspondence theorem for ∇Φ * . Furthermore, where we have again used the Lipschitz continuity of both f and ∇Φ * . We furthermore have, by convexity of f , where we used in the last inequality that from Assumption 5 it holds ||∇f (x)|| * ≤ L and the Lipschitz continuity of ∇Φ * . Using Lemma 23 (proved in the Appendix Section A), combining inequalities (14), (15) and (16), and taking expected values the result follows.
A similar result can be obtained for a strongly convex f .
Proposition 15 (Convergence of ISMD for a strongly convex objective). Assume that Assumptions 1 -4 and 6 hold. Then it holds, Proof. Observe that, using strong convexity, we have, Therefore, where we have used the Lipschitz continuity of ∇Φ * , the assumption that ||∇f (x)|| * ≤ L, the assumption that , and Lemma 24 (proved in Section A of the Appendix). Multiplying by e −µT , using the assumption that ||∇f (x)|| * ≤ L and the Lipschitz continuity of ∇Φ * , combining inequalities (14), (15) and (17), and taking expected values the result follows.
The deviation from the minimum is upper-bounded by four terms. The first two terms are the standard optimization errors, where we observe that the noise variance is reduced by a factor of N . The third and fourth terms are penalties incurred due to each of the particles having different values. These two terms measure the deviation of each individual particle from the particle average. Loosely speaking, if the fluctuation term ||z i t || * is bounded and small then the interaction will drive the particle system to a state near consensus and the gradient terms will direct the particles towards the minimum. There is furthermore a tradeoff between the interaction and the variance. If the interaction term is bounded and not increasing with N , the more particles, the smaller the distance to the optimum, as is witnessed by the term σ 2 2N ||∇Φ * || ∞ . This implies that the variance of the particles x i t is also smaller, and they lie closer around the optimal value. Comparing Proposition 14 in the convex case and Proposition 15 in the µ-strongly convex case we observe that the strong convexity of f can speed up convergence, mainly due to the additional factor 1 µ . We remark that the results presented in this section thus show that the expected value of the distance of the objective function evaluated at the time average to the objective function evaluated in the optimum decreases as N increases. In addition to this, as we will present in our numerical results later in Section 5 that both the variance of the cost function as well as the variance of the samples are smaller when using interaction.

A comparison with averaging
Consider the case of the convex objective. The interaction in the previous section is introduced in order to obtain a way to control the effects of noise. In particular we observed that if the term ||z i t || * was bounded and not an increasing function of N , the variance could be reduced by having interacting particles. We can compare the setting with interaction to one in which we simply average the trajectories of the particles.
If the following holds, then we have from Lemma 23 (in Appendix A) and the convexity of f that Notice that the Assumption in (18) holds for a function of the form f (x) = a T x + b + Φ(x) which is the particular loss function considered in [31]. Thus in the setup of this kind structure in the loss function and gradient, averaging the particles can decrease the effect of the noise.
In a general setting where the loss function is nonlinear, one can decrease the effects of noise by averaging the particles trajectories. Combining (15), (16) and Lemma 23 we have, In the nonlinear case the average of the gradients is not the gradient of the average and therefore the additional term representing the deviations of each individual particle from the particle average plays a role.

Bounding the fluctuation term
In this section we present a bound on the fluctuation term. When the fluctuation term is sufficiently small all particles have approximately the same value; we refer to this as the particles having achieved consensus. The dynamics of the fluctuation term in a vectorized form is given by, where we have used that Proof. Define the Lyapunov function, By Itô's lemma it then follows that, We can rewrite the drift term of dV t as, Note that,z By Assumption 5, Furthermore, by the properties of the graph Laplacian, it holds L1 dN = 0, so that, where we have additionally used thatz T dB t d = ||z|| 2 dB t . Rearranging, using the property of the Itô integral, and using V T ≥ 0 we find, We note that in the above result we used the convexity of V but we did not exploit the interaction between the particles. From the above statement we observe that in a deterministic setting, i.e. if σ = 0, consensus can be achieved even without interaction. Each particle is driven towards the optimum x * by the gradient term. Consensus is in this case achieved exactly at optimality. In the stochastic setting exact consensus can no longer be achieved and the distance to optimality remains bounded by a function of the noise. Proposition 17 (Bounding the fluctuation term: the strongly convex case). Let Assumptions 1 -3 and 8 hold. Then, Proof. Define again the Lyapunov function, Tz .

It is clear that
or in other words consensus has been achieved and the fluctuation term is zero. By Itô's lemma it then follows that, We can rewrite the drift term of dV t as, Note that,z By Assumption 8,z T (∇V(z N t ) − ∇V(z t )) ≤ −κz Tz . Furthermore, by the properties of the Laplacian, it holds L1 dN = 0, so that, Therefore, Integrating the expression and using the above bound we obtain, where we have additionally used thatz T dB t d = ||z t || 2 dB t . Then, taking expectations we have: The statement follows from norm equivalence, i.e. || · || 2 * ≤ K|| · || 2 2 .
Remark 18 (Constant in the norm equivalence). If the Euclidean norm is used above in Assumptions 1 -2, then K = 1. In certain cases however, K can be dimension-dependent.
We furthermore note, by the bound on E 1 ||z i t || * with the dependence on N remaining the same. This is formalized in the following lemma. Lemma 19. It holds for some constant K, Proof. Let || · || be an arbitrary norm. Observe that by Hölder's inequality where we have applied Young's inequality in the second inequality. We obtain, The statement follows using Proposition 17.
The above statements show that for a strongly convex objective approximate consensus can be achieved for a sufficiently high t and for a sufficiently large κ + λ. In the strongly convex case approximate consensus can be achieved even in the case of no interaction if κ is sufficiently large. For the convex case we have κ = 0. In this case approximate consensus can be achieved if λ is sufficiently high. For a convex objective without interaction we would be left with the noise term similar to the result in Proposition 16. Without interaction the strong convexity of the objective thus determines in how much the effect of noise can be reduced. With interaction the interaction type itself, i.e. the λ, also plays a role and for a sufficiently high λ the effects of noise on consensus can be minimized. We remark furthermore that we bound the dual norm by the L 2 norm using norm equivalence; the downside of this is that an additional dimensionality-dependence can be obtained. This could be mitigated by working with directly the dual norm as in [31] but in this case the obtained result is rather limited to a specific set of objective functions; alternatively one can work in the primal space directly. The latter idea will be addressed in future work.

Extensions to distributed optimization
A common setting is to consider minimizing f (x) = N i=1 f i (x) and each particle corresponds to a computing worker or processor having access only to f i and thus ∇f i • ∇Φ * (z i t ) would be used in (10) (instead of ∇f • ∇Φ * (z i t )). The proof of Proposition 16 and 17 is based the fact that each particle has access to the full objective function f . As a result terms that would be relevant in a distributed setup related to variability of each f i vanish. Clearly, this results in a better bound and thus a better convergence result, but the current results can be extended for the distributed case by increasing the interaction strength and multiplying A in (10) by a sufficiently high constant to ensure consensus arises. The details are left for future work.

On the tradeoff between noise reduction and interaction
From the results in Section 3.2 we observed that for interacting particles the convergence can be improved since the effect of noise in the term 1 N σ||∆Φ(z i t )|| ∞ is reduced by a factor 1/N . We however remark that this will only be achieved if the term ||z i t || * is bounded and non-increasing with N . As we showed in Section 3.3 we can present a bound on the interaction term. In particular, we showed that the fluctuation is bounded and non-increasing with N . The fluctuation is bounded in the L 2 norm, and a noise term does remain; as discussed, this noise term can be controlled by the strong convexity of the objective or by imposing an interaction strength, which means replacing A with ϑA for ϑ > 1. In other words, when λ + κ is sufficiently high the fluctuation term is sufficiently small. Specifically, in the strongly convex case for a large N, T the fluctuation term is bounded by dK (κ+λ) σ 2 for some constant K. We conclude this section by saying that as long as the decrease in the value of 1 N σ||∆Φ(z i t )|| ∞ is larger than the increase in the value of ||z i t || * , interaction with N particles achieves a closer convergence to the optimum. In other words, the amount of particles can be seen as an alternative to a decreasing learning rate or vanishing noise variance (e.g. the latter can be achieved by increasing the mini-batch size).

Understanding performance from convergence to stationarity
We proceed with exploiting tools from the analysis of SDEs and in particular the rate of convergence to their invariant distribution. Studying the convergence of the Law(z t ) or Law(z t ) is a cornerstone in the analysis of sampling schemes and can also provide valuable insights into the optimization problems; see [32,35] for recent works in this direction. In particular, we will discuss the convergence of (9) or (10) to the corresponding stationary distributions. Note that in the spirit of dual (or Nesterov) averaging the SDE in (10) is fairly standard as it evolves only in the mirror space. In the analysis so far we were investigating how close a time average 1 ) is to f (x * ) using a Lyapunov method based on Bregman divergence. We will complement these results with rates of convergence based on logarithmic Sobolev inequalities based on the celebrated Bakry-Emery approach, see [2], [3] for details.

Stationary distributions and particle correlations
We first consider the case when particles do not interact. One approach for finding x * is to run N independent copies of (9) for a long time and then select best particle i * as arg min i∈{1,...,N } f (x i t ). This intuitive approach is based on min f being equivalent to max exp(− 2 σ 2 V) and on the convergence of dLaw(zt) When stacking the particles together in z, then one can postpone choosing the best particle and considers min x instead, with each particle i here sharing the same properties and dynamics. To make this a bit clearer for our setup, denote with η t = Law(z t ) the law of N independent particles each following (9); or equivalently each z i t evolving as (10) with A = 0. It is well-known (see e.g. [29]) that under appropriate assumptions on the objective V, the stationary distribution of this SDE is given by, with Z being a finite normalizing constant. It is clear that finding the mode of η ∞ is equivalent to solving the following optimization problem: An important observation here is that independence of each i is not crucial and one can modify the cost V(z i ) by adding terms that do not affect the minimizer. This is true when using instead the following cost function for z: Note that the quadratic term due to interaction terms acts as a regularizer, which aims to impose consensus i.e. z i =z N . We can make the this claims more precise in the following lemma. Lemma 20 (Interaction preserves the minimum). Let z * := arg min z V(z), and let W(z) be as defined in (20). Then Proof. From the definition of z * and the optimality conditions we have ∇V(z * ) = 0 and that there exists β > 0 such that Hess(V)(z * ) βI d . We have, Clearly, when z i = z * for every i we then have ∇W(z) = 0. Furthermore: and This implies so that Hess(W)(z * ) βI dN and z * minimizes W since β > 0.
We will now extend the previous discussion on invariant distributions and z * for ISMD and (10) with A being nonzero. We will denote the distribution of the joint particle system and each marginal as follows: Recall z t evolves as in (13), whose the invariant distribution is given as where Z N is the normalization constant. Similar as in the non-interacting setting, finding the mode of η N ∞ is equivalent to solving the following optimization problem:z which we saw earlier in Lemma 20 is equivalent to solving (19).
Despite introducing correlations between the particles through L, the mode of the invariant distribution of the interacting particle system is exactly the minimizer of the objective V and is achieved when all particles are at consensus. The presence of the noise means that exact consensus or synchronization of all particles cannot be achieved. Similar to SMD stochasticity implies that at stationarity the z i t will behave as samples from the invariant distribution and thus oscillating together around the optimum z * . In this context it is convenient to use convergence results for η i,N t for deriving performance bounds for optimization or assessing the level of consensus based on differences between η i,N t and η j,N t . In addition to establishing positive curvature at z * , using diagonal dominance in the final steps in the proof of Lemma 20 implies positive definiteness for W. This is not surprising given the sum of quadratics in the interaction term acts as a convex regularizer. We summarize this in the corollary below: Corollary 21. Let Assumptions 1 -3 and 8 hold. Then W is strictly convex.
Proof. Assumption 8 implies that Hess(V)(z) κI d for some κ > 0. Note that (21)-(22) hold for any z and The result is somewhat negative in that to obtain strict convexity of W one requires strong convexity of V, that is adding z T Lz is not sufficient to obtain strong convexity when V is only convex. This is particularly relevant later when we will apply Bakry-Emery theory and derive log Sobolev inequalities.

Log-Sobolev inequalities for the particle system
We proceed to present the exponential convergence properties for η N t . Our first result in this section is to apply Bakry-Emery theory to derive the Log-Sobolev inequality for the particle system. Similar to the result in Lemma 3.5, Corollary 3.7 in [21] we have the following result: Proposition 22 (Convergence of ISMD using Bakry-Emery theory). Let Assumptions 1 -3 and 7 hold (w.r.t the Euclidean norm). Assume also η N 0 satisfies a log-Sobolev inequality with constant C 0 . Then η N t satisfies a log-Sobolev inequality for any smooth function f and with L N = L µ + ||L|| where ||L|| is any matrix norm on L.
Proof. For the proof we will rely on the Bakry-Emery theory. To state the Bakry-Emery criterion for (13), we need to define the following differential operators: The Bakry-Emery criterion requires to verify Γ 2 (f ) ρΓ(f, f ). We have: so need to verify that holds for ρ = σ 2 2 κ from Assumption 7 or Assumption 8, Corollary 21 and L1 = 0. For the second part we explicitly require ρ > 0, which follows straightforwardly from Assumption 8. This gives directly C t and we note that the second term in C t is due to the random initialization (e.g. see [21,Corollary 3.7]). When Assumption 8 holds and ρ > 0 we have For the mean-mode result, we first need to establish η N ∞ ∈ P 2 (R dN ). This can be shown using Theorem 25 in Appendix B.2 with the log Sobolev constant for any 1-Lipschitz function f , so consider f (z) = z and the moment generating function to get the following moment bounds for s = 0 so moments can be bound using exponential integrability and η N ∞ has finite second moment. Combining this with log det Σ N ≤ T r Σ N − I gives log det Σ N ≤ dN (2C ∞ − 1) .

We can then gather Corollary 28 and Proposition 29 in Appendix B.2 we get
The statement then follows by noting that, , and using the convergence rate from the result derived above in (23) for the first term and the mean-mode bound for the second.
This result shows that the rate of convergence to the invariant measure is exponential in time with factor ρ. Using the relationship between the mode of the invariant measure and the minimizer of the objective function, this result implies that with exponential convergence the samples z i t will oscillate around the optimum z * for a strongly convex objective. Note that the discrepancy between the mean and the mode at stationarity seen in (24) means this oscillation is not centred around z * . Equation (24) also shows that when the width of η N ∞ is small, i.e. σ 2 is low, then z t will lie closer to the minimizer. To identify z * in practice, one needs to look at the histogram of z i t (w.r.t i) or the occupation measure of W(z t ) after some burn-in time t 0 or a combination of the two. The amount of required burn-in can be lower when fast convergence to stationarity occurs. From the expression for ρ we can see that the rate of convergence ρ increases with the noise variance and the strong convexity coefficient of the objective function. Finally, the dependence on N arrives from dimensionality dependence of certain bounds involving Gaussian integrals; see Proposition 29 in Appendix B.2.
At the level of the marginal η i,N ∞ a similar bound to (24) can be obtained, where dN is replaced with d. Convergence of η N t to η N ∞ implies convergence of the marginals so one can deduce a similar result to Proposition 22 for the marginals and a similar concentration inequality for V(z i t ) showing that the dimension dependence of the optimization is with d and not dN . The details are omitted as the proof is very similar to Proposition 22.

On the tradeoff between noise reduction and interaction: a sampling perspective
One could attempt to compare Proposition 22 with the results obtained in Section 3.2, although they are quite different in nature. In Propositions 14 and 15 we showed that the expected time average converges to a neighborhood of the optimizer at rate e −µt for the strongly convex objective (with µ being the convexity constant for f ). Proposition 22 and specifically the result in (24) looks directly at W(z t ) − W(z * ) and uses instead κ, the convexity constant of V, which was also used to bound the fluctuation terms in Section 3.3. While this means direct comparisons are case specific and depend on the choice of the mirror map Φ, the convergence rates arein both settings exponential. In other words, for the strictly convex case both Proposition 15 and Proposition 22 lead to exponential rate albeit with different constants due to working with f and V respectively. Furthermore, in terms of the tradeoffs between interaction and convergence, both bounds give a similar result. Observe that from Proposition 15 combined with Proposition 17 we obtain, where we remark that we have simplified notation to shorten the presentation and we assume T 0 and denote with x i ∞ the corresponding exponential time average. The expectation is furthermore taken with respect to the initial condition. Furthermore, from (24), From the above bounds -while as mentioned different in nature -we observe that both consist of a term which decreases with N , and a term which is of O(1) in N . Both bounds additionally decreases with an increase in κ, or the strong convexity constant of the objective function.

Mapping back to the primal space
Convergence in the mirror domain and in z i t should imply convergence in the primal domain and x i t . Given Φ is invertible the law of X t follows from standard change of variable/coordinates of continuous random variables. This is we have for all t ≥ 0: Recall an implication of Proposition 22 is that W 2 (η N t , η N ∞ ) ≤ K e −ρt , so this is equivalent to In the context of Bakry-Emery theory and Markov semigroups we refer the the interested reader to Section 1.15.1 of [3] and to [14] for a rigorous treatment related to sampling methods including equivalence in total variation norm.

Numerical results
In this section we report numerical experiments using two standard benchmark problems. The first problem is based on linear regression and is an adapted version of the experimental setup that appeared in [4]. The second problem is based on a variation of the classic traffic assignment problem [23]. We consider compare between (Stochastic) Gradient Descent with Euclidean projections ((S)GD), (Stochastic) Mirror Descent with a Bregman divergence generated by the entropy function. In our numerical implementations of the optimization algorithm we use the Euler discretization of the continuous SMD dynamics: The parameters that we have to choose are the time step , the learning rate η, the number of training iterations T N and the noise level σ (set to zero in case of GD or MD).

Linear system with simplex constraints
We consider a first similar set-up as in Example 9.19 from [4] to compare mirror descent with projected gradient descent. Consider the problem, where X = ∆ n , the unit simplex, W ∈ R m×d and b ∈ R m . Unless otherwise mentioned, we set = 0.1 and T N = 2000. We generate A randomly with some given condition number κ, in particular a well-conditioned problem with κ(W ) = 1 and an ill-conditioned one with κ(W ) = 100. We let b i ∼ N (0, 1). We furthermore set A ij = 1 N , i.e. a mean-field interaction matrix. We set m = 100 and d = 100.
Single particle optimization We begin with a comparison of the convergence speed of GD, MD and SGD and SMD, to show that MD/SMD can outperform GD/SGD, as shown in the theoretical results in Lemmata 9-10 and Propositions 11-12. We observe from Figure 1 that mirror descent algorithms attain a much faster convergence than gradient descent in the system where the matrix W has a high condition number. For a low condition number the performance of mirror descent is similar to that of gradient descent. Mirror descent seems more robust to a higher noise level, and is able to converge even in a high noise setting. We set σ = 0.05.
Interacting particle optimization The theoretical results in Section 3.2 show that the expected value of the distance between the objective evaluated at the time average and the objective evaluated at the optimum is smaller for a larger number of particles. In this section we analyze this numerically and observe that ISMD is indeed able to converge closer to the minimum. Here we consider the linear system optimized with one, 10 and 100 particles. We set σ = 0.05 and η = 0.01 for κ(W ) = 10 and η = 0.001 for κ(W ) = 100. We remark that we thus use fixed learning rates. In Figure 2 we show the initial convergence speed. In the setting with a high condition number the convergence speed using interaction can be significantly faster than when considering just a single particle.
In Figure 3 we show the convergence speed as well as the distribution of the loss f (x i t ) − f (x * ) for the linear system with different condition numbers We observe from the plots and histograms that the noise in the interacting sampler is smaller so that it able to converge to closer to the global minimum than single-particle optimization. More specifically, the more particles used, then the smaller the variance of the samples is, so that the values are closer to the true minimizer.  Convergence: vanishing noise against interaction In this example we study how using a vanishing noise compares to using a fixed learning rate with and without interaction. Often in practice, a vanishing noise is achieved by using a larger batch size, so that the gradient noise is smaller. In Figure 4 the results for the distance to the optimum, f (x i t ) − f (x * ), are shown. As usual, the convergence speed using interaction is higher and using interaction can achieve the same effect as using a vanishing noise in terms of the distance to the optimum. We remark here that using an interaction strength, i.e. letting the interaction term be given by θ can further help in controlling the noise variance. The effects and necessity of such an interaction term will be studied in future work. An additional benefit of having noise is that noise allows to escape local minima in a non-convex problem. Being able to control this noise without decreasing the convergence speed is a significant benefit.
Variance reduction As implied by the theory in Section 3.2 and the numerical results in the last section, the expected value of the distance between the objective evaluated at the time average and the objective evaluated at the optimum is smaller for a larger number of particles. In addition to this we can study the variance of the particles as well as of the loss. In this section we show that the variance is decreased so that using more particles results in losses which lie closer together as well as particle values which are closer. We can study the dependence of the variance in the samples on the number of particles N in more detail. We keep the learning rate fixed at λ = 0.001 and set the condition number to 100. To make a fair comparison we sample N i.i.d. copies of the SMD algorithm, and use N particles in the ISMD algorithm. We first plot the variance of f (x i t ) − f (x * ) after convergence (t > 1000) for the i.i.d. an interacting setting in the left-hand side of Figure 5. Clearly, in the non-interacting setting the variance is not decreased when increasing the number of particles and due to the algorithm not having converged properly the variance is fluctuating. In the interacting setting more particles results in a lower loss variance. To understand the effects of interaction and the number of particles on the fluctuation term we compute ||x i T −x N t || 2 , whose results are shown in the right-hand side of Figure 5. As expected from the theoretical analysis the deviation from the mean is smaller in the interacting case Figure 3: A comparison between SMD and ISMD for the linear system with condition number 10 (L) and condition number 100 (R). The histogram shows the distribution of the losses f (x i t ) − f (x * ) after convergence (t > 1000). We observe from the loss and the histogram that more particles leads to a convergence closer to the true minimum due to the interaction reducing the effect of noise. Figure 4: The initial distance to the optimum (L) and the distance to the optimum after convergence (R) for the linear problem with condition number 100 for SMD with a fixed noise σ = 0.05, SMD with a vanishing noise σ = σ/(t + 1) 1/10 and ISMD with 10 particles. than in the i.i.d. case. Furthermore, as N increases the deviation remains constant, validating the theoretical results in Section 3.3 that the fluctuation is bounded and a non-increasing function of N .

Traffic Assignment Problem
The objective of the traffic assignment problem is to compute the optimal path between two nodes in a graph. To save space we refer the reader to [23] for a precise description of the problem. We only mention that the problem is a convex optimization problem with a simplex constraint and therefore fits our framework. Our experiments below are based on a random geometric graph G(n, r) i.e. we placed n points uniformly at random in [0, 1] 2 and connect any two points whose distance is at most some w. The weights of the edges in the graph also act as the weight of each edge. We randomly chose two nodes to act as the origin and destination nodes. To set up the optimization problem we also need to compute all the simple paths between the origin and destination node. We denote the length of the maximum path considered in this phase as r max . The size of r max decides the dimensionality of the problem. For the stochastic version of the problem we added Gaussian N (0, σ g ) noise on each edge. We used a constant step-size for all algorithms. For mirror descent (all variations) we used η = 0.2 and for gradient descent we manually selected the best step-size for the problem (see below for additional remarks regarding GD for this problem).
Comparison between MD and GD An interesting difference between the regression problem above and the traffic assignment problem is that the solution of the latter is sparse. As a result the solution is on the boundary with many paths having zero load. In this case the gradient descent algorithm performed considerably worse than all variations of mirror descent (i.e. with many or single particles) and in all variations of the problem we considered (small/large problems, with and without noise). In order to illustrate the difference in performance between mirror descent and gradient descent we plot the results for one example in Figure 6. The graph used to generate the problem above has 100 nodes, and the resulting optimization problem has 5831 possible paths i.e. d = 5831. We begin with an illustration of the performance of standard MD against GD in Figure 6; GD is significantly slower to converge than MD. Given the large difference in performance between MD and GD we will not report any more results from GD for the traffic assignment problem. Figure 6: A comparison between MD and GD. GD with Euclidean projections is considerably slower than MD for this problem because the solution is sparse. IMD performs the same as MD for deterministic problems.
Comparison between SMD and ISMD Similar to the experiments run for the linear regression we compare the iterations of the single particle mirror descent i.e. (SMD) and its interacting variant (ISMD) with a different number of particles (N = 10, 100, 1000). In the left of Figure 7 we show the iteration history of (SMD) and (ISMD) with On stochastic mirror descent with interacting particles: convergence properties and variance reduction N = 100 particles. To visualize the differences between (SMD) and (ISMD) and the impact of the number of particles we plot a histogram in the right plot in Figure 7 of the iterations after convergence (typically t > 100). Clearly there is a considerable reduction in the variance of the iterates between (MD) and (ISMD). We attempt to quantify the reduction in variance in the next set of experiments. Figure 7: Comparison between SMD and ISMD for a single path (L) and a histogram of the samples for the comparison between SMD and ISMD with 10, 100 and 1000 particles (R). With more particles the variance of the samples is lower.
In order to further quantify the variance reduction between ISMD and SMD in a fair way it is important to compare the reduction in variance between N i.i.d. copies of SMD and ISMD with N particles. In Table 1 we present the variance reduction results for different particles and observe that interaction can significantly reduce the variance.
Nodes r max d N  Table 1: A comparison of the variance of N i.i.d. copies of SMD and N particles in ISMD. We obtain a variance reduction using interaction.

Mini-batch optimization
In machine learning the optimization objective typically consists of a sum over data samples, i.e.
where m is the sample size. The gradient is typically computed over a subset of the data (a mini-batch), since for large m computing the full gradient is too costly. The downside of computing the gradient over a mini-batch is that the gradient contains noise. In our work we proposed interaction as a way of reducing the effects of this gradient noise. In this section we show how interaction helps in reducing the effect of noise in mini-batch gradient descent, and consider in more detail the computational costs related to interaction against an increased batch size. The ideas presented here can be generalized to popular machine learning techniques such as neural networks.
Consider a similar setup as in Section 5.1, In every iteration, the gradient is computed over a subset of the data, f S (x) = 1 |S| i∈S f i (X), where |S| refers to the size of the batch. The optimization algorithm is then given by, The noise is thus implicit in the gradient ∇f S . In Figure 8  with η = 0.1. Using N particles infers a computational cost of N times that of a single particle. However, as observed in the plot, using N particles allows to use a significantly smaller batch, i.e. due to the interaction the optimization can converge even with a higher noise. We conclude that interaction can be a convenient methodology for improving the convergence of the algorithm when using stochastic gradients.

Conclusion
In this work we analyzed stochastic mirror descent with interacting particles with the aim of understanding the tradeoffs between computation, communication and variance reduction. With a fixed learning rate or a constant noise variance convergence to the exact minimizer is not possible. Prior results showed that decreasing the learning rate [35] or using a vanishing noise variance [23] can result in good convergence properties of stochastic optimization. In our work we showed that using N interacting particles can in a similar fashion decrease the effect of noise and result in a closer convergence to the optimum. In addition using interaction decreases the variance of the samples further than that of N independent replicas of SMD. We therefore argue that interaction between N particles is a viable method and can control the convergence of the algorithm.
We analyzed the convergence from an optimization perspective exploiting Lyapunov functions. We showed that the deviation of an individual particle from the particle mean is bounded and showed that ISMD achieves a linear convergence rate in the convex case and an exponential one in the strongly convex setting. We furthermore presented the relationship between sampling from the invariant measure and converging to the optimizer of the objective function.
Using log-Sobolev inequalities we obtained an explicit convergence rate to the invariant measure, implying that the particles converge to a neighborhood of the optimizer at an exponential rate. Our numerical results show that interaction is beneficial in cases where the solution is sparse and therefore the behavior of the algorithm at the boundary is important. This is a case where the role of the mirror map becomes more critical, which is also reflected in our assumptions.
Interesting extensions include the study of SMD in the non-convex setting where the noise can help escape from local minima (see e.g. [47], [13]). In this case the interaction can play the role of a regularizer due to the convexification of the loss surface [46]. Additionally,even in a non-convex setting using interaction enforces the particles to converge to consensus in a local minimum. The interaction can also act as a pre-conditioner and its statistical properties may help speed up the convergence or minimize the variance . An open question is the study of the optimal interaction matrix beyond just mixing times [7]. Alternatively of interest is to study the choice of the interaction such that communication is minimal, but variance reduction and consensus is still guaranteed.
It is also interesting to study such algorithms in a truly distributed setting, i.e. when each particle has access to a subset of the data. This is the case in learning with privacy constraints or federated learning approaches. Some preliminary results on this have been presented in [6]. In this setting also the effect of an interaction strength will be of interest and we aim to address the necessity of considering an additional interaction strength θ, i.e. using θ N i=1 A ij (z j t − z i t ), in future work. More specifically, in a distributed setting the particles optimize different objectives, and -unlike in our setting of a same objective -convergence to consensus cannot be guaranteed without imposing a sufficiently high interaction strength.
Lastly, the effect of time discretization was not considered here. This is a well-understood topic. We refer the interested reader to [22], [39] for details related to the Euler discretization and [41] for more advanced schemes which also include acceleration.
Proof. Consider the Lyapunov function defined as Applying Itô's lemma and using the fact that σ(y N t − x * ) T 1 Integrating, and using the standard bounds we obtain, Similar to the convex case, in Proposition 15 we used the following lemma. Lemma 24. Under the continuous stochastic mirror descent dynamics (9) it holds,

Theorem 26. Assume
HessV ρI d , ρ ∈ R Then ν t satisfies (1) with constant If in addition ρ > 0, we get, Note that requiring V being strongly convex is crucial to obtain explicit rates of convergence to equilibrium. This is a particular instance of the celebrated (and more general) Bakry-Emery criterion on the curvature of Markov semigroups.
To state the more general version for any Markov semigroup we need to define the following differential operators defined on appropriate domains: The Bakry-Emery criterion consists of verifying, for any smooth f , which is equivalent to (27) for the dynamics in (26).
In addition, point 3. in Theorem 25 also implies a convergence in W 2 , i.e.

B.2 Empirical risk bounds for Proposition 22
The main purpose of this section is to provide so called empirical risk bounds for E η N ∞ (W) − W(z * ). We will base our derivation on Proposition 3.4 in [32] and specifically focus on certain aspects of the particle system with invariant distribution We will define the differential entropy for a probability measure η ∈ P 2 (R d ) as H(η) = − log dη dz (z)η(dz). From the definition of H we have for η N In what follows we will provide upper bounds for H(η N ∞ ) and lower bounds for log Z N to reach an empirical risk bound. Lemma 27. Let η ∈ P 2 (R d ) and let Σ denote the covariance matrix of η. Then H(η) ≤ d 2 (log 2π + 1) + 1 2 log det Σ.
Next we provide a lower bound for log Z N . Proposition 29. Let Assumptions 1, 2 hold w.r.t the Euclidean norm. We have that ∇W is L N -Lipschitz and as a result Proof. The first claim is trivial and one has L N = L µ + L where L is any matrix norm on L. Smoothness of ∇W implies the following quadratic upper bounds: Then we have, where in the last step we used the Gaussian integral exp − 1 2 z T Az + b T z dz = (2π) dN det A exp 1 2 b T A −1 b with A = 2L N σ 2 I dN and b = Az * . Remark 30. When non-Euclidean norms are used L N = L µ K + L , with K arising from norm equivalence.