Neural network approach to data-driven estimation of chemotactic sensitivity in the Keller-Segel model

1.   Introduction

Interacting particle systems play an important role in many applications in science - on the one hand as a modeling framework for social and biological systems, on the other as a tool for computational algorithms used in data science. In the latter case the collective behavior of interacting particle systems is used to solve high-dimensional problems, often resulting from non-convex optimization tasks in data science. Well known algorithms include particle swarm optimization (PSO) ^[1], ant colony optimization ^[2] or evolutionary ^[3] and genetic algorithms ^[4].

PSO was first introduced in ^[1] and has been successfully used in engineering applications ^[5]. Each particle in a PSO algorithm adjusts its position due to information of the global best, personal best and a noise term that allows for exploration of its neighborhood. Consensus-based optimization (CBO) ^[6,7] combines the idea of swarm intelligence with consensus formation techniques ^[8,9,10] to obtain a global optimization algorithm for non-convex high-dimensional problems. On the one hand particles explore the state space via an amplitude modulated random walk. On the other a drift term convects them towards the weighted global best. The method was first introduced in ^[6] and analyzed at the mean-field level in ^[7]. Recent developments of CBO include component-wise diffusion and utilize random mini-batch ideas to reduce the computational cost of calculating the weighted average ^[11]. Other contributions investigate a CBO dynamic that is restricted to the sphere ^[12,13]. Also, convergence and error estimates for time-discrete consensus-based optimization algorithms have been discussed ^[14]. CBO-type systems are related to large interacting particle systems, in which the dynamics are driven by weighted average quantities, see ^[15,16,17].

The model proposed in the work is based on the component-wise diffusion variant introduced in ^[11] and combines it with personal best information. This adjustment is motivated by the original work on PSO by Eberhart and Kennedy ^[1], where the particles move towards a (stochastic) linear combination of their personal best or the common global best position. The new information leads to an additional drift term in the dynamics. We investigate two types of memory effects - either using a weighted personal best over time or the personal best value in the past. The latter corresponds to record processes, see ^[18] for an overview. The former is used in the presented analysis and approximates the personal best of each particle. We expect by arguments similar to the Laplace principle that the weighted mean converges towards the personal best.

The proposed stochastic dynamics with weighted personal best fall into the class of stochastic functional differential equations. These equations are in general non-Markovian and their mean-field limit has been investigated in special cases only. For example, Gadat and Panloup ^[19] investigated a non-Markovian process with memory, which corresponds to the weighted average of the drift all along the particle's trajectory. This memory term is of a special form allowing them to rewrite the system as a 2-dimensional non-homogeneous Markovian dynamical system. Moreover, they exploit this special structure to analyze the existence and long time behavior of solutions as well as the mean-field limit. The strategy of increasing the dimension to get around the non-Markovian nature goes back to the Mori-Zwanzig formalism, see ^[20]. This strategy was recently adapted for non-Markovian interacting dynamics by Duong and Pavliotis in ^[21]. Another interesting work by Kuntzmann, see ^[22], investigates the ergodic behavior of self-interacting diffusions depending on the empirical mean of the process. The proposed generalization of CBO with weighted personal best does not fall into this category, hence the derivation and analysis of the respective mean-field dynamics, which often give useful insights into the dynamics, is to the best of the authors' knowledge open. This applies as well for personal best, where the update of the best function value corresponds to a record process. Hence, we focus on the well-posedness of the stochastic system as well as a detailed computational investigation of the dynamics.

This paper is organized as follows: we introduce the particle dynamics with (weighted) personal best in Section 2 and illustrate its dynamics with first toy examples. Section 3 discusses well-posedness and existence of solutions to the SDE model with weighted personal best. Section 4 presents extensive computational experiments of various benchmark optimization problems.

2.   Consensus based optimization with personal best

In this section we discuss how personal best information can be included in consensus based optimization algorithms as proposed by Carrillo and co-workers in ^[6,7,11]. We start by introducing the notation before continuing with the modeling.

2.1.   Notation

We refer the euclidean norm by $|x| = (x_1^2 + \dots + x_d^2)^{1/2}$ for $x \in \mathbb{R}^d$ and $|Y| = (\sum_{i, j = 1}^{dN} Y_{ij}^2)^{1/2}$ for matrices $Y \in \mathbb R^{dN \times dN}$. The set of natural numbers without $0$ is denoted by $\mathbb N^* = {1, 2, 3, \dots}$ and the half-line $[0, \infty)$ by $\mathbb R^+$. A vector valued function or vector $x \in \mathbb R^{dN}$ is assumed to be of the form $x = (x^1, \dots, x^N)$ with $x^i \in \mathbb R^d.$ When discussing the stochastic systems we follow the notation of ^[23]: $(\Omega, \mathcal F, \mathbb P, \{ \mathcal F_t \}_{t\ge 0})$ corresponds to the stochastic basis with sample space $\Omega,$ filtration $\mathcal F$ and probability function $\mathbb P.$ Moreover, $S_d^p[0, T]$ is the space of (equivalence classes of) $\mathcal P$ -measurable continuous stochastic processes $X \colon \Omega \times [0, T] \rightarrow \mathbb R^{dN}$ such that

Two processes $X, Y$ are called equivalent if $(X_t = Y_t \forall t \in [0, T]) \; \mathbb P$-almost surely ($\mathbb P$-a.s.). Furthermore, $S_d^p$ is the space of (equivalence classes of) $\mathcal P$-measurable continuous stochastic processes $X \colon \Omega \times \mathbb R_+ \rightarrow \mathbb R^{d}$ such that for all $T > 0$ the restriction $X_{|[0, T]}$ of $X$ to $[0, T]$ belongs to $S_d^p[0, T]$. Analogously, we define $\Lambda_d^p(0, T)$ as the space of (equivalent classes) of $\mathcal P$-measurable processes $X \colon [0, T] \rightarrow \mathbb R^d$ such that

We refer to $\Lambda_d^p$ as the space of (equivalence classes of) $\mathcal P$-measurable continuous stochastic processes $X \colon \Omega \times (0, +\infty) \rightarrow \mathbb R^{d}$ for which for all $T > 0$ the restriction $X_{|[0, T]}$ of $X$ to $[0, T]$ belongs to $\Lambda_d^p(0, T).$ Moreover, for any $\phi \in C(\mathbb R_+, \mathbb R^{dN})$ we define

2.2.   The model

We wish to approximate the global minimum

of a given non-negative, continuous objective function $f: \mathbb{R}^d \rightarrow \mathbb{R}$. In doing so we consider $N \in \mathbb{N}$ particles and denote the position of the $i$-th particle at time $t$ by $X_t^i : = X^i(t) \in \mathbb{R}^d$, $i = 1, \ldots N$. Note that we use $X_t = X(t) = (X^1(t), \ldots X^N(t)) \in \mathbb{R}^{dN}$, when referring to the positions of all particles at time $t$. In CBO particles compare their current function value with a weighted mean value based on the current information of the whole system. A particle moves towards the position of the weighted mean, if the function value of the weighted mean is lower. Following the ideas of ^[6,7], we use the weighted average

with $\alpha > 0$, to approximate the global best, that is, the particle with the lowest function value. Note that even though the weighted average uses only information of the current time step, it is assumed to approximate the global best over time as well, since a particle that is close to the weighed average experiences only small drift and diffusion. The parameter $\alpha$ scales the influence of local and global minima in the weighted mean. In fact, for $\alpha = 0$ the weights are independent of the function values, and all particles are weighted equally. For $\alpha > 0$ the particle with the best function value has the largest weight. Moreover, the Laplace principle from large deviations theory ^[24] assures that $v_f(t)$ converges to the global best, as $\alpha \rightarrow \infty.$ For more details on the Laplace principle in the CBO context, we refer to ^[6,7].

In the original version of PSO, see ^[25], particles compare their current position with the global best as well as their personal best value up to that time. We propose two different approaches how to include the personal best $p_i$ of the i-th particle. First, we consider the true personal best by setting

Moreover, the personal best can be approximated similarly to the global best, $v_f(t)$, defined in (2.2). Hereby, we use the entire trajectory in the past and refer to this trajectory by $X = (X^1, \dots, X^N)$ with $X ^i \in C(\mathbb R_+, \mathbb R^d)$ for all $i = 1, \dots, N$. Let $X^i_0$ denote the initial position of the $i$-th particle at time $t = 0$, the weighted mean over time of the $i$-th particle is defined by

with $\beta > 0$. Note that the well-posedness result presented in Section 3 holds for the weighted personal best (2.4) only. Again, by the Laplace principle, we expect that $p_f^i(t) \rightarrow P_f^i$ as $\beta \rightarrow \infty$.

We recall that particles either move towards the global or personal best state. The respective CBO dynamics for the $i$-th particle, $i = 1, \ldots N$ are then given by the following SDE:

where

The function $\text{H}$ corresponds to the Heaviside function and $\sigma > 0$ denotes the standard deviation. System (2.5) is supplemented with the initial condition $X_0^i = \xi_i$, $i = 1, \dots, N.$ The drift and diffusion are motivated by the following considerations:

1. If the global best $v_f$ is better than the current position $X^i_t$ and the personal best $p_f^i$, the particle moves towards the current global best $v_f$.

2. If the personal best $p_f^i$ is better than the current position $X^i_t$ and the global best $v_f$, the particle moves towards the personal best $p_f^i$.

3. If none of the above holds, the particle still explores the function landscape via Brownian motion until it reaches the global best $v_f$.

Note that the drift coefficients depend on the past of each particle, hence system (2.5) is non-Markovian. The form of the memory does not allow us to use existing results, such as ^[23] to rewrite the system. Hence the existence and form of the respective mean-field model is, up to the authors' knowledge, not known.

Remark 1. (1) The CBO version proposed in ^[11] can be recoverd by setting

(2) Note that the personal best (2.3) and weighted personal best (2.4) can be computed very efficiently due to their accumulative structure; this does not significantly increase the computational cost.

Throughout this manuscript we will refer to the dynamics defined by (2.5) with (2.6) as CBO, and to (2.5) with (2.3) or (2.4) as personal best (PB) or weighted personal best (wPB), respectively.

2.3.   Toy example: CBO vs. PB dynamics

In the following we will illustrate the differences between CBO and (w)PB using a 1D toy objective function $f$ and $3$ particles. We consider a double well-type $f$ of the form:

For this function, shown in Figure 1 the global and local minimum, located at $x = -1.00125$ and $x = 0.998748$ respectively, are very close. In the following we perform $1000$ Monte Carlo (MC) simulations with deterministic initial conditions $\xi$. We count a run as run successful, if the final position of the particles satisfy $|v_f(T) - X^i(T)| < 0.4$ for all $i = 1, \dots, N.$ The final time is set to $T = 100,$ the time step size $dt = 10^{-3}$ and $\beta$ in (2.4) to $\beta = 30$. We study the dynamics for the following two initial conditions:

(IC1) Initialize 2 particles near the local minimizer and 1 particle near the global minimizer.

(IC2) Initialize 1 particle near the local minimizer and 2 particles near the global minimizer.

The initial positions (IC1) and (IC2) of the particles correspond to the gray dots in Figure 1 and in Figure 3, respectively. We discuss (IC1) first. In this situation the weighted average, $v_f(0),$ is located near $x = 0.9,$ thus, the Heaviside functions are zero and the system would be in a stationary state for $\sigma = 0$.

For $\sigma > 0,$ the diffusion term drives the dynamic and the particles are exploring their neighborhood. Due to the multiplicative factor, the particle on the left is exposed to more diffusion than the particles on the right. In case of the CBO scheme, the particle on the left has a high probability of jumping out of the basin of the global minimum. Then, all particles concentrate near the local minimum. For one run, this behavior is illustrated by the positions of $v_f(T)$ shown in Figure 1 (left). This alone does not reflect the concentration which becomes apparent in Figure 2. In fact, the orange lines show fluctuations for small times but stabilize quickly indicating that no diffusion is present and thus that all particles are concentrated. This behavior changes when personal best information is included. Here, particles still explore their neighborhood, however at some point their current positions are worse than their personal best, and hence the drift starts pulling them back towards their personal best. This behavior is also illustrated by the success rates stated in Table 1. We see that (w)PB outperform CBO for large and small values of $\alpha$. The 'pull-back' effect slows down the convergence of (w)PB - we observe that the respective energies decrease slower than for CBO in Figure 2. Nevertheless, they find the global minimum.

Next we consider initial condition (IC2). Again, in the deterministic case $\sigma = 0$ the initial configuration is stationary. For $\sigma > 0$ the particles on the left are less diffusive than the particle on the right. Therefore, it is more likely that the particle on the right jumps into the basin of the global minimum. This is illustrated in Figure 3 and confirmed by the success rates in Table 2. Again, CBO converges faster than (w)PB, see Figure 4. Nevertheless, the function values at the point of concentration are smaller for (w)PB which means that the slower algorithms find better approximations. Note that the scale of the time step-axis is much smaller than in Figure 2.

3.   Well-posedness results

In the following we discuss well-posedness of the wPB model. We begin by considering CBO with component-wise diffusion, which was proposed in ^[11].

3.1.   Well-posedness of CBO with component-wise diffusion

Theorem 1. Let $f$ be locally Lipschitz and $N \in \mathbb{N}$. Then system (2.5) with $\lambda(t, X) \equiv \lambda, \; \mu(t, X) \equiv 0$ admits a unique strong solution for any initial condition $\xi = (\xi_1, \dots, \xi_N)$ satisfying $\mathbb E\, |\xi|^2 < \infty$.

A detailed proof can be found in the Appendix. Let us just emphasize that the estimates in the proof of Theorem 1 are independent of the dimension, $d,$ as was already highlighted in ^[11] for the mean-field setting. This is in contrast to ^[6,7], where the estimates depend on $d.$

3.2.   Well-posedness in case of weighted personal best

Next, we present an existence and uniqueness result for the proposed SDE model with weighted personal best and smoothed Heaviside functions. Note that the structure of the weighted personal best suggests the idea of introducing new variables for the numerator and the denominator. This re-formulation converts the non-Markovian process into a Markovian system for times $t > 0$, but violates the initial condition. We therefore use different proofs for properties of SDEs with local Lipschitz conditions as well as path-dependent SDEs that can be found in the literature ^[23]. To the authors' knowledge none of them covers the case of path-dependent SDEs with local Lipschitz conditions. In the following we present a proof which combines the two techniques to obtain a well-posedness result.

We assume that the regularized Heaviside function $H^\epsilon$ satisfies the following conditions:

(A1) Let $0 \le H^\epsilon(x) \le 1$ for all $x \in \mathbb{R}$.

(A2) There exists a constant $C > 0$ such that

This corresponds to the following regularized problem

with

Moreover, we assume that the objective function $f$ satisfies the following properties:

(A3) Positivity: it holds $0 \le f(x) \quad \text{ for all } x \in \mathbb{R}^d,$

(A4) Quasi-local Lipschitz condition: for any $n < \infty$ and $|x|, |y| \le n$ it holds

with a constant $L_f > 0$ depending on $n$ only.

Remark 2. The well-known regularization of the Heaviside function

satisfies the assumptions (A1) and (A2). Note that in the context of optimization problems, the positivity assumption on $f$ is not too restrictive. Since $f$ corresponds to a minimization functional it is naturally bounded from below and can be shifted to satisfy the positivity constraint.

The following proof is based on a combination of arguments of Theorem 3.17 and Theorem 3.27 in ^[23] - this yields well-posedness of (3.2). We begin with two lemmata providing necessary estimates. The first lemma is concerned with properties of the weighed averages. Note that the global best $v_f$ depends on the current state of the process and the personal best $p_f$ on the whole trajectory. We therefore write $v_f[\varphi(t)] = v_f(t)$ and $p_f[\varphi] = p_f$.

Lemma 1. Let $f$ satisfy (A3) and (A4), $N \in \mathbb N$ and $\varphi = (\varphi^1, \dots \varphi^N) \in \mathcal C(\mathbb R_+, \mathbb{R}^{dN}).$ Then

Moreover, the averages satisfy the local Lipschitz conditions:

for all $t \in [0, \infty)$ with $|\varphi|_t, |\hat \varphi|_t \le n$ with constants

Here $L_f$ is the Lipschitz constant of $f$ in $B_n = \{ x \colon |x| \le n\}$ and $\underline f, \overline f$ correspond, respectively, to the minimal and maximal values of $f$ on $B_n.$

The proof of Lemma 1 can be found in the Appendix. Using Lemma 1 we show that the drift and diffusion terms satisfy local Lipschitz and linear growth conditions. These properties allow us to apply the existence and uniqueness result later on.

Lemma 2. Let (A1)–(A4) hold. Then

given by

and

with $\varphi^i = (\varphi_{(i-1)d+1}, \dots, \varphi_{(i-1)d+d})$ satisfy the following conditions for all $\varphi, \psi \in C(\mathbb R_+, \mathbb R^{dN})$ and all $R > 0:$

(i) $|b(t, \varphi) - b(t, \psi)| \le L_R \|\varphi - \psi\|_t,$       (iii) $|\Sigma(t, \varphi(t)) - \Sigma(t, \psi(t))| \le \ell_R |\varphi(t)- \psi(t)|,$

(ii) $|b(t, \varphi)| \le a_R \|\varphi\|_t$,             (iv) $|\Sigma(t, \varphi(t))| \le b_R |\varphi(t)|$,

where $L_R, \ell_R, a_R, b_R \in \mathbb R.$

Proof. To show (ⅰ) we calculate

where

From Lemma 1 we know that $|v_f[\psi(t)] - v_f[\varphi(t)] |^2 \le C_1 |\varphi(t) - \psi(t)|^2,$ and $|p_f^i[\varphi](t) - p_f^i[\psi](t)|^2 \le C_2 \|\varphi^i - \psi^i\|_t^2$ with constants $C_1$ and $C_2$ given by (3.6) with $n = R$. This yields (ⅰ) since

The Lipschitz bound (ⅲ) follows from similar arguments using the diagonal structure of $\Sigma$:

The last two inequalities hold due to

Equipped with this lemma, we have everything at hand to prove the main theorem.

Theorem 2. Let (A1)–(A4) be satisfied and $\xi \in L^0(\Omega, \mathcal F_0, \mathbb P, \mathbb R^{dN})$ with $\mathbb E|\xi|^p < \infty$ for each $p > 0.$ Then, there exists a unique strong global solution to (3.2). Moreover, there exists a constant $C_{p, T, L_r, \ell_R}$ such that

The proof combines arguments of Theorem 3.17 (path-dependent SDE) and Theorem 3.27 (SDE with local Lipschitz coefficients) in ^[23].

Proof. We start by proving uniqueness. Let $X, \hat X \in S_{dN}^0$ be two solutions to (3.2) corresponding to initial data $\xi, \hat \xi \in L^0(\Omega, \mathcal F_0, \mathbb P, \mathbb R^{dN}),$ respectively. Then it holds $\mathbb E\, \xi = \mathbb E\, \hat \xi$. Define the stopping time $\tau_n(\omega) = \inf \{ t \ge 0 \colon |X_t(\omega)| + |\hat X_t(\omega)| \ge n\}.$ Then the two solutions satisfy

and $\tau_n \rightarrow \infty$ for $n \rightarrow \infty.$ Note that we have global Lipschitz constants for $b$ and $\sigma$ for all times $t \in [0, \tau_n]$ with $n$ arbitrary but fixed. This allows us to use Theorem 3.8 in ^[23], see proof of Theorem 3.27 in ^[23] for more details, to obtain

for some $V^R$ depending on the local Lipschitz constants $L_R, \ell_R, a_R, b_R$ and $p \ge 2$ arbitrary. Hence, the uniqueness of the solution on $[0, \tau_n].$ As $\tau_n \rightarrow \infty$ for $n \rightarrow \infty,$ this allows us to conclude the global uniqueness of $X \in S_{dN}^0.$

Next, we show the existence of solutions. Let $M \in \mathbb N^*$ and $0 = T_0 < T_1 < \dots < T_M = \tau_n$ with $T_i = \frac{i \tau_n}{M}$. It holds

We employ a fixed point argument for the mapping $\Gamma \colon S_{dN}^p [0, T_1] \rightarrow S_{dN}^p[0, T_1]$ given by

where we need no stopping times due to $t < \tau_n$ on $[0, T_1].$ Indeed, the mapping $\Gamma$ is well-defined since for all $\varphi \in C(\mathbb R_+, \mathbb R^{dN})$

is satisfied. Because of the Lipschitz continuity, both stochastic processes $b(\cdot, U)$ and $\sigma(\cdot, U_s)$ are progressively measurable for all $U \in S_{dN}^p[0, \tau_n]$ and $b(\cdot, U) \in L^p(\Omega, L^1(0, \tau_n))$ and $\sigma(\cdot, U) \in \Lambda^p_{dN \times dN}(0, \tau_n).$ Therefore,

We will show that the operator $\Gamma$ is a strict contraction on the complete metric space $S_{dN}^p[0, T_1]$ for sufficiently large $M$ (where $S_{dN}^p[0, T_1]$ is equipped with the usual distance $d_{p, M}(U, V) = (\mathbb E\| U -V\|_{T_1}^p)^{1/{p \vee 1}}$). Let $U, V \in S_{dN}^p[0, T_1].$ By the Burkholder-Davis-Gundy inequality we have

Let $M_0 \in \mathbb N^*$ such that $(1 \vee 2^{p-1}) \alpha(\frac{\tau_n}{M_0}) \le \left(\frac{1}{2} \right)^{1 \vee p}$. Then $\Gamma$ is a strict contraction in $S_{dN}^p[0, T_1]$ and thus (3.2) has a unique solution $X \in S_{dN}^p[0, T_1].$ We extend the solution to the interval $[0, T_2]$ by defining a mapping, again, called $\Gamma \colon S_{dN}^p[0, T_2] \rightarrow S_{dN}^p[0, T_2]$:

We repeat the argument $M_0$ times to be valid the whole interval $[0, \tau_n]$. Since $\tau_n \rightarrow \infty$ almost surely, the uniqueness of the solution implies that

Hence, the process $X \in S_d^0$ is defined by $X_t(\omega) = X_t^n(\omega)$ if $0 \le t \le \tau_n(\omega)$ and $\tau_n(\omega) > 0$ is the unique solution to the regularized problem (3.2).

4.   Numerical results

The numerical simulations are based on the direct simulation of system (2.5) using the Euler-Maruyama scheme, in which we do not approximate the Heaviside function $H$. Note that the smoothing of $H$ was only needed for analytic considerations. In practise, we want only one of the drift terms to effect the particles dynamics, which is why we have not pursued this option any further. The final time is set to $T = 15$, discretized into $3 \times 10^4$ time steps. All presented results are averaged over $M = 5000$ realizations. The standard deviation is set to $\sigma = 0.5$, while the number of agents depends on the dimension of the function space. In particular, we set the number of agents to $3, 5$ or $10$ times the space dimension. The initial positions of particles are drawn from a uniform distribution within a specific domain for each function. The parameters $\alpha$ and $\beta$ to compute the global and personal best are set to

A realization is successful if the average mean is close to the function minimum $f_{\min}$, in particular

We compare the performance of the CBO scheme with $\mu = 0$, PB and wPB for the following benchmark problems:

1. Alpine ^[26]: This non-convex differentiable function has a global minimum at $x_{\min} = (0, \ldots, 0)$

2. Ackley ^[27]: This function is continuous, non-differentiable and non-convex and has its global minimum at $x_{\min} = (0, \ldots, 0)$.

3. Rastrigin ^[28]: The Rastrigin function is continuous, differentiable and convex, has lots of local minima and a global minimum at $x_{\min} = (0, \ldots 0)$.

4. Xinsheyang2: ^[26] This function is continuous but not differentiable and non-convex with a global minimum at $x_{\min} = (0, \ldots 0)$.

The choice of these functions is based on the different characteristics they have, see Figure 5 for plots in 2D. Table 3 shows the results for the Alpine function (4.1) and the Ackley function (4.2). We observe that the success rate increases with the number of particles, and decreases for higher space dimension. Weighted personal best and personal best give comparable results, which is not surprising since wPB approximates PB for large values of $\beta$. A similar behavior can be seen in the case of the Rastrigin function (4.3) and the Xinsheyang function (4.4) in Table 4.

We conclude by investigating a 2D version of the toy problem considered in Section 2.3:

This function has a global minimum at $x_{\min} = (-1.00125, 0)$ and a local minimum at $(0.998748, 0)$. We wish to explore the dynamics of this 2D version for different number of particles. In doing so we consider $4$, $8$ or $16$ particles and start the particle schemes with $2$, $4$ and $8$ placed in each of the two wells with a random perturbation. Furthermore we set $\alpha = 10$ and $\beta = 20$. Table 5 shows the results as the number of particles increases. We observe that CBO and (w)PB perform equally well for large numbers of particles, and that (w)PB outperform CBO for few particles. This could be explained by the fact that the probability of all particles deviating from the global minimum decreases as their number increases.

5.   Conclusion

In this paper we introduced a consensus based global optimization scheme, which includes the personal best information of each particle. The proposed generalization is motivated by the original works on particle swarm algorithms, in which particles adjust their position as a linear combination of moving towards the current global best and their personal best value.

We discussed how information about the personal best can be included in consensus based optimization schemes, leading to a system of functional stochastic differential equations. A well-posedness result for the respective regularized non-Markovian SDEs was presented. New features of the algorithm with personal best were illustrated and compared in computational experiments. The numerical results indicate that information about the personal best leads to higher success rates in the case of few particles and that the corresponding weighted means are better approximations of the global function minima.

Acknowledgments

The authors would like to thank Christoph Belak (TU Berlin), Stefan Grosskinsky (University of Warwick), Greg Pavliotis (Imperial College London) and Oliver Tse (TU Eindhoven) for the helpful discussions and constructive input. CT was partly supported by the European Social Fund and by the Ministry Of Science, Research and the Arts Baden-Württemberg. MTW was partly supported by the New Frontier's grant NST-0001 of the Austrian Academy of Sciences ÖAW.

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendix

Proof. [Theorem 1] In the following we denote in abuse of notation by $v_f[X]$ the vector $(v_f, \dots, v_f) \in \mathbb R^{dN}.$ We rewrite (2.5) with $\lambda(X^i(t), v_f, p_f^i) \equiv \lambda, \mu(X^i(t), v_f, p_f^i) \equiv 0$ as

with $\text{diag}(X(t) - v_f[X(t)]) \in \mathbb R^{dN \times dN}$ being the diagonal matrix with $d_k = diag(X^k(t) - v_f) \in \mathbb R^{d \times d}$ for $k = 1, \dots, N$ and $dB_t$ a $dN$-dimensional Brownian motion. The argument follows the lines of the well-posedness in ^[7]. In fact, let $M[X(t)] = \text{diag}(X(t) - v_f[X(t)])$ and $n\in \mathbb N$ arbitrary. We have to check that there exists a constant $C_n$ such that

for every $|X(t)| \le n.$ Note that $f(X(t))$ is bounded for $|X(t)| \le n$ due to its local Lipschitz continuity. Hence, the estimate for the first term on the left-hand side is identical to the one in ^[7]. Indeed, we have

For the component-wise drift we obtain

Combining the two preceding estimates we obtain Eq (5.1). Now, employing [29,Ch 5.3,Thm 3.2] yields the desired result.

[Lemma 1] We start by showing continuity. For $p_f[\varphi]$ we need to check continuity as $t \rightarrow 0$. In fact, l'Hospital's rule yields

This directly implies the continuity of the vector $p_f[\varphi].$ Moreover, it is easy to see that

We are left to show the local Lipschitz continuity. Therefore, we estimate

to obtain

with $L_f$ being the global Lipschitz constant of $f$ on $B_n = \{ x \in R^d \colon |x| \le n\}$ and $\underline f, \overline f$ are the minimal and maximal value of $f$ on $B_n,$ respectively. Altogether, this yields the estimate

For the vectors $p_f(t) = (p_f^i[\varphi](t))_{i = 1, \dots, N}$ and $\hat p_f(t) = (\hat p_f^i[\varphi](t))_{i = 1, \dots, N}$ this implies

Similarly, we have

which satisfy

Thus, we get

Taking squares leads to the estimate