extfieldsofstudy:
Computer Science
Mathematics

s2fieldsofstudy:
Computer Science
Physics

text:
Momentum Doesn't Change the Implicit Bias

The momentum acceleration technique is widely adopted in many optimization algorithms. However, the theoretical understanding on how the momentum affects the generalization performance of the optimization algorithms is still unknown. In this paper, we answer this question through analyzing the implicit bias of momentum-based optimization. We prove that both SGD with momentum and Adam converge to the L2 max-margin solution for exponential-tailed loss, which is the same as vanilla gradient descent. That means, these optimizers with momentum acceleration still converge to a model with low complexity, which provides guarantees on their generalization. Technically, to overcome the difficulty brought by the error accumulation in analyzing the momentum, we construct new Lyapunov functions as a tool to analyze the gap between the model parameter and the max-margin solution.

INTRODUCTION
It is widely believed that the optimizers have implicit bias in terms of selecting output parameters among all the local minima on the landscape. The pioneer work (Soudry et al., 2018) shows that gradient descent would converge to the L 2 max-margin solution for the linear classification task with exponetial-tailed loss, which mirrors its good generalization property in practice. Since then, many efforts have been taken on analyzing the implicit bias of various local-search optimizers, including stochastic gradient descent (Nacson et al., 2019), AdaGrad (Qian & Qian, 2019) and optimizers for homogeneous neural networks (Lyu & Li, 2019;Wang et al., 2021;Qian & Qian, 2019).
However, though the momentum acceleration technique is widely adopted in the optimization algorithms in both convex and non-convex learning tasks (Sutskever et al., 2013;Vaswani et al., 2017;Tan & Le, 2019), the understanding on how the momentum would affect generalization performance of the optimization algorithms is still unclear. A natural question is:

Can we theoretically analyze the implicit bias of momentum-based optimizers?
In this paper, we take the first step to analyze the convergence of momentum based optimizers and unveil their implicit bias. Specifically, we study the classification problem with linear model and exponential-tailed loss using Stochastic Gradient Descent with Momentum (SGDM) and Adam optimizers. We consider the optimizers with constant learning rate and constant momentum hyperparameters, which are widely adopted in practice, e.g., the default setting in popular machine learning frameworks (Paszke et al., 2019) and in experiments (Xie et al., 2017). We note that Gradient Descent with Momentum (GDM) can be viewed as a special case of SGDM, and naturally share the properties for SGDM. Our main results are summarized in Theorem 1.
Theorem 1 (informal). With linear separable dataset S, for SGDM and Adam, the loss converges to 0 with rate O( 1 t ) where t is the number of iterations, the parameter norm diverges to infinity, and direction of parameters converges to the direction of the L 2 max-margin solution.
Theorem 1 states SGDM converges to the L 2 max-margin solution, which is the same as SGD, indicating that momentum does not affect the convergent direction. The good generalization behavior of the output parameters of SGDM is well validated as the margin of a classifier is positively correlated with its generalization error (Jiang et al., 2019). This is supported by existing experimental observations (c.f. Figure 1, (Soudry et al., 2018) and Figure 2, (Nacson et al., 2019)). Similar claims hold for Adam, which is also well supported by empirical results in Wang et al. (2021).
Our contributions are significant in terms of the following aspects: • We establish the implicit bias of the momentum based optimizers, an open problem since the initial work Soudry et al. (2018). The momentum based optimizers are widely used in practice and our theoretical characterization deepens the understanding on their generalization property, which is important by its own.
• Technically, we propose a new Lyapunov function to analyze the convergence of SGDM, which helps to bound the sum of squared gradients along the training trajectory. Compared to the usual one, the new Lyapunov function depends on a middle variable of an alternative update rule of SGDM, which helps to capture the historical dependence in the momentum update. To our knowledge, such a technique has not been exploited ever before and can be of independent interest for convergence analysis of momentum-based optimizers. We then construct a new Lyapunov function to bound the difference of learned parameters and the scaled max-margin solution, which finally leads to the direction convergence. This Lyapunov function provides a direct way to establish the convergence to the desired direction.
Organization of This Paper. Section 2 collects further related works on the implicit bias of first order optimizers and convergence of momentum-based optimizers. Section 3 shows basic settings and assumptions which will be used throughout this paper. Section 4 studies the implicit bias of GDM as a warm up, while Section 5 and Section 6 explore respectively the implicit bias of SGDM and Adam. Discussions of these results are put in Section 7.

FURTHER RELATED WORKS
Implicit Bias of First-order Optimization Methods. Soudry et al. (2018) prove that gradient descent on linear classification problem with exponential-tailed loss converges to the direction of the max L 2 margin solution of the corresponding hard-margin Support Vector Machine. Nacson et al.
(2019) extend the results in (Soudry et al., 2018) to the stochastic case, proving that the convergent direction of SGD is the same as GD almost surely. Qian & Qian (2019) go beyond the vanilla gradient descent methods and consider the AdaGrad optimizer instead. They prove that the convergent direction of AdaGrad has a dependency on the optimizing trajectory, which varies according to the initialization.  propose a primal-dual analysis framework for the linear classification models, and prove a faster convergent rate of the margin by increasing the learning rate according to the loss. Based on ,  design another algorithm with an even faster convergent rate of margin by applying the Nesterov's Acceleration Method on the dual space. However, the corresponding form of the algorithm on the primal space is no longer a Nesterov's Acceleration Method nor GDM, which is significantly different from our settings.
On the other hand, there is another line of work trying to extend the result in the linear case to deep neural networks. Ji & Telgarsky (2018); Gunasekar et al. (2018) study the deep linear network and Soudry et al. (2018) study the two-layer neural network with ReLU activation. Lyu & Li (2019) propose a framework to analyze the asymptotic direction of GD on homogeneous neural networks, proving that given there exists a time the network achieves 100% training accuracy, GD will converge to some KKT point of the L 2 max-margin problem. Wang et al. (2021) extend the framework of Lyu & Li (2019) to adaptive optimizers, and prove RMSProp and Adam without momentum have the same convergent direction as GD, while AdaGrad doesn't. The results (Lyu & Li, 2019;Wang et al., 2021) indicate that results in the linear model can be extended to deep homogeneous neural networks, and suggest that the linear model is a proper start point to study the implicit bias.
Convergence of Momentum-Based Optimization Methods. For convex optimization problems, the convergence rate of Nesterov's Acceleration Method has been proved in Nesterov (1983). In contrast, although GDM (Polyak's Heavy-Ball Method) was proposed in (Polyak, 1964) prior to the Nesterov's Acceleration Method, the convergence of GDM on convex loss with Lipschitz gradient was not solved until Ghadimi et al. (2015) provides an ergodic convergent result for GDM, i.e., the convergent result for the running average of the iterates. However, the ergodic result is undesired under many learning scenarios, e.g., in classification tasks, the optimization algorithms usually output the parameters of the last step. To the best of our knowledge, there are only two works on the non-ergodic analysis of GDM: Sun et al. (2019) and Tao et al. (2021). Sun et al. (2019) prove that if the training loss is coercive (the training loss goes to infinity whenever parameter norm goes to infinity), convex, and globally smooth, then, with constant momentum hyper-parameter, the training loss converges to minima with rate O(t −1 ). Tao et al. (2021) analyze a case when momentum coefficient increases to 1 and the gradient is bounded all over the parameter space, showing SGDM can improve the convergence rate of SGD by a factor log(t).
However, there is no existing work on the implicit bias of momentum-based optimizers for the classification problem, which is first analyzed by this paper.

PRELIMINARIES
In this paper, we focus on the linear model with exponential-tailed loss. We first derive the results for binary classification, then we show that the methodology can be easily extended to the multi-class classification problem.
Problem setting. The dataset used for training is defined as , where x i ∈ R d is the i-th feature, and y i ∈ R is the i-th label (i = 1, 2, · · · , N ). We will use the linear model to fit the label: for any feature x ∈ R d and parameter w ∈ R d , the prediction is given by w, x .
For binary classification, given any data z i = (x i , y i ) ∈ S, the individual loss for parameter w is given as ℓ(y i w, x i ). A common setting is to ensemble the feature and label together asx i = y i x i (there is a mapping T : (x, y) → yx, andx i is T ((x i , y i ))). The individual loss can be rewritten asl The optimization target is defined as the averaged loss: Without loss of generality, we consider the case with normalized data 1 , that is, Optimizer. Here we will introduce the update rules of SGDM and Adam. SGDM can be viewed as a stochastic version of GDM by randomly choosing a subset of the dataset to update. Specifically, the update rule of SGDM is given as where B(t) is a subset of S with size b which is sampled independently and uniformly with replace- We also define F t as the sub-sigma field such that {w(t)} ∞ t=1 is adapted with respect to Filtration {F t } ∞ t=1 . The Adam can be viewed as a variant of GDM in which the preconditioner is adopted, whose form is characterized as follows: where 1 √ν (t−1)+ε1 d ) is called the preconditioner at step t.
1 The proof can be naturally applied to unnormalized data by letting ℓ(x) = ℓ(maxx ∈T (S) x · x) and Assumptions: The analysis of this paper are based on three assumptions, respectively on the separability of the dataset, the individual loss behaviour at the tail , and the smoothness of the individual loss. We list them as follows: Assumption 1 (Linearly Separable Dataset). There exists one parameter w ∈ R d , such that Assumption 2 (Exponential-tailed Loss). The individual loss ℓ is exponential-tailed, i.e., • Differentiable and monotonically decreasing to zero, with its derivative also converging to zero, i.e., lim x→∞ ℓ(x) = lim x→∞ ℓ ′ (x) = 0, and ℓ ′ (x) > 0 ∀x; • Close to exponential loss when x is large enough, i.e., there exist positive constants c, a, K, µ + , µ − , x + , x − and x 0 , such that, We provide explanations of these three assumptions respectively. Based on Assumption 1, we can formally define the margin and the maximum margin solution of an optimization problem, which is introduced in Definition 1 Definition 1. Let the margin γ(w) of parameter w defined as the lowest score of the prediction of w over the dataset S, i.e., γ(w) = minx i∈T (S) w,x i . By Assumption 1 and the positive homogeneous of γ, γ(ŵ γ(ŵ) ) = 1, and thus we define the maximum margin solutionŵ and the margin γ of the dataset S as follows: Since · 2 is strongly convex and set {w : γ(w) ≥ 1} is convex,ŵ is uniquely defined.
Assumption 2 constraints the loss to be exponential-tailed, which is satisfied by many popular choices of ℓ, including the (ℓ exp (x) = e −x ) and the logistic loss (ℓ log (x) = log(1 + e −x )). Also, as c and a can be respectively absorbed by resetting the learning rate and data as η = cη and x i = ax i , without loss of generality, in this paper we only analyze the case that c = a = 1.
The globally smooth assumption (Assumption 3. (S)) is strictly stronger than the locally smooth assumption (Assumption 3. (D)). One can easily verify that both the exponential loss and the logistic loss meet Assumption 3. (D), and the logistic loss also meets Assumption 3. (S).

WARM UP: IMPLICIT BIAS OF GDM
In this section, we study the implicit bias of GDM as a warm up. The update rule of GDM is which is a determined case of SGDM, with batch size b = N and thus without randomness. Although the analysis of GDM is essentially easier than that of SGDM, it helps us gain a better understanding of this problem, and helps us to demonstrate several techniques which will be applied latter. The formal theorem of the implicit bias of GDM is as follows: Theorem 2. Let Assumptions 1, 2, and 3. (D) hold. Let β ∈ [0, 1) and η < 2 . Then, for almost every data set S, with arbitrary initialization point w(1), GDM (Eq. (5)) satisfies that w(t) − log(t)ŵ is bounded as t → ∞, and lim t→∞ Before we move on to the proof, we would like to provide some explanations of the results in Theorem 2. To begin with, Theorem 2 adopts an constant momentum hyper-parameter, which agrees with the practice use (e.g., β is set to be 0.9 (Xie et al., 2017)). Also, Theorem 2 puts no restriction on the range of β, which allows wider choices of hyper-parameter tuning. Furthermore, Theorem 2 shows the implicit bias of GDM agrees with GD in linear classification with exponential-tailed loss (c.f. Soudry et al. (2018) for results on GD), and this consistency can be verified by existing results (c.f. Section 7 for detailed discussions).
We then present a proof sketch of Theorem 2, which is divided into two parts: we first prove that the sum of squared gradients is bounded, which indicates both the loss and the norm of gradient converge to 0; we then show the difference between w(t) and log(t)ŵ is bounded, and therefore, the direction ofŵ dominates as t → ∞.
Stage I: Bound the sum of squared gradients. The core of Stage I is to select a proper Lyapunov function ξ(t), which is required to correlated with the training loss L and be non-increasing along the optimization trajectory. For GD, since L itself is non-increasing with properly chosen learning rate, we can just pick ξ(t) = L(t). However, as the update of GDM doesn't align with the direction of negative gradient, training loss L(t) in GDM is no longer monotonously decreasing, and the Lyapunov function requires special construction. A choice of ξ(t) is proposed as the follows: Remark 1. Although this Lyapunov function is obtained by (Sun et al., 2019) by directly examining the Taylor's expansion at w(t), the proof here is non-trivial as we only requires the loss to be locally smooth instead of globally smooth in (Sun et al., 2019), and the Taylor's expansion can only be applied to w(t + 1) if it's ensured that all parameters on the line αw(t) + (1 − α)w(t + 1) (α ∈ (0, 1)) have training loss no larger than L(w(0)).
To prove Lemma 1, we define w(t + α) = αw(t) + (1 − α)w(t + 1) for any t ∈ Z + and α ∈ (0, 1), and prove a generalized version of Eq. (6): The proof idea is that as long as Eq. (7) holds for all time in [1, t + α), the training loss across [1, t + α] will be smaller than L(w(1)), and the strict inequality in Eq. (7) holds. Consequently, there exists some small enough positive real ε, such that Eq. (7) holds for all time in [1, t+α+ε), and we are able to extend the feasible set where Eq. (7) holds. The formal description of the generalized lemma and its corresponding proof is deferred to Appendix B.1.1.
However, it no longer holds for GDM, as w(t + 1) − w(t) is a exponentially-decayed sum of gradients at w(s), s ≤ t. To handle this problem, we propose a novel Lyapunov's function for the direction convergence, concluded as the following Lemma.

Lemma 2. Let all conditions in Theorem 2 hold. Then, r(t) is bounded if and only if the function
) is upper bounded.
The first claim in Lemma 2 provides an alternative approach to verify that r(t) is bounded, while the second claim shows this approach can be fulfilled. The motivation of this lemma is to construct a easy-to-verify criterion (i.e., g(t) is upper bounded) of r(t) is upper bounded. To demonstrate how g(t) is constructed intuitively, we consider the continuous dynamics approximation of GDM, i.e., We then calculate the 1 2 r(t) 2 − 1 2 r(1) 2 by applying the integral of the dynamics as where the last equation is due to Integration by part. We ds , dw(s) ds ds) + 1 2 r(t) 2 , and it can be verified the derivative ofg(t) is r(t), − η 1−β ∇L(w(t)) − 1 tŵ . By replacing dw(t) dt as w(t) − w(t − 1), and dr(t) dt as r(t) − r(t − 1) ing(t), we obtain g(t), and g(t + 1) − g(t) has the form r(t), − η 1−β ∇L(w(t)) − log t+1 tŵ , which can be analyzed following the similar routine as (Soudry et al., 2018). The proof of the second claim is completed.

TACKLE THE DIFFICULTY BROUGHT BY RANDOM SAMPLING
In this section, we analyze the implicit bias of SGDM. The randomness introduced by random sampling makes the analysis of GDM no longer work for SGDM. As an example, if we want to follow the same routine of the GDM case to show L(w(t)) + β 2η E w(t) − w(t − 1) 2 is a Lyapunov function of SGDM, we can only have As the squared first momentum is smaller than the second momentum, EL(w(t)) + β 2η E w(t) − w(t − 1) 2 may no longer be monotonously decreasing for every β ∈ [0, 1). Furthermore, if the margin γ is small, β will be required to be upper bounded by a small number, which contradicts to the relatively large choice of β in practice, e.g. 0.9 (We defer a detailed discussion to Appendix B.2.3). To tackle this problem, we propose a new Lyapunov function, with which we derive the following theorem for the implicit bias of SGDM: Theorem 3. Let Assumption 1, 2, and 3. (S) hold. Let β ∈ [0, 1) and η < min{ 1−β 2H 2 N 3 β 2 }. Then, for almost every data set S, with arbitrary initialization point w(1), SGDM (Eq. (1)) satisfies w(t) − log(t)ŵ is bounded as t → ∞ and lim t→∞ w(t) w(t) =ŵ ŵ , almost surely (a.s.).
To the best of our knowledge, this is the first convergence analysis of SGDM with no restriction on the gradient norm and with constant learning rates. Also, as both 1−β are monotonously increasing with respect to batch size b, Theorem 3 also sheds light on the learning rate tuning, i.e., the larger the batch size is, the larger the learning rate is. Furthermore, similar to the GDM case, Theorem 3 shows the implicit bias of SGDM under this setting is consistent with SGD (c.f. (Nacson et al., 2019) for the implicit bias of SGD). This matches the observations in practice (c.f. Section 7 for details). Remark 2. In practice, mini-batch SGDM are more widely adopted, whose update rule differs from SGDM only by the way obtaining B(t). Specifically, within the T -th epoch E(T ) = {KT + 1, · · · , KT +T } (K is the length of one epoch), {B(t)} t∈E(T ) is a randomly and uniformly partition of S. One may wonder whether the same result hold for mini-batch SGDM. The answer is "Yes", with the a.s. condition removed in this case. We defer the detailed description of the corresponding theorem together with the proof to Appendix D.1.
The proof sketch follows the same framework as that for GDM by dividing the proof into two stages. However, the implementations in both stages differ, among which the proof of Stage I is significantly distinctive, as a new Lyapunov function is proposed.
By a simple rearrangement of the second equation, we have u(t+1) = w(t)+ 1 1−β (w(t+1)−w(t)), which differs from w(t + 1) only by a larger step size updated from w(t). The following lemma then indicates that EL(u(t)) is a proper selection of Lyapunov function. Lemma 4. Let all conditions in Theorem 3 hold. Then, we have EL(u(t)) may be not monotonously decreasing. However, the upper bound of EL(u(t)), i.e., L(u(1)) − t s=1η 4 'E ∇L(w(s)) 2 , is monotonously decreasing byη 4 E ∇L(w(t)) 2 at step t, and leads to the sum of expected squared gradients being finite along the trajectory (which further indicates the sum of squared gradients is finite, a.s.).
Stage II. Bound the difference between w(t) and log(t)ŵ. The Lyapunov function used by Stage II of SGDM is the same as the g(t) in the GDM case. As we have proved the sum of squared gradient along the trajectory is bounded, one can easily obtain that ∞ t=1 w(t + 1) − w(t) 2 is bounded, and r(t) being bounded is still equivalent to that g(t) is upper bounded. However, when it comes to the detailed calculation of analyzing g(t + 1) − g(t), it differs from the GDM case due to the random subset. We defer the proof to Appendix B.2.2.

ANALYZE THE EFFECT OF PRECONDITIONERS
When it comes to the analysis of Adam, the effect of the preconditioner should be taken into consideration when designing the corresponding Lyapunov functions. To tackle this problem, we incorporate the preconditioner into the Lyapunov function, and obtain the implicit bias of Adam as follows: Theorem 4. Let Assumption 1, 2, and 3. (D) hold. Let 1 > β 2 > β 4 1 ≥ 0, and the learning rate η is a small enough constant (The upper bound of learning rate is complex, and we defer it to Appendix C.1). Then, for almost every data set S, with arbitrary initialization point w(1), Adam (Eq. (2)) satisfies that w(t) − log(t)ŵ is bounded as t → ∞, and lim t→∞ Remark 3. Existing literature usually assume a time-decaying hyperparameter choice of β 1 or β 2 (c.f., (Kingma & Ba, 2014;Chen et al., 2018)). To the best of our knowledge, the only work analyzing Adam with constant β 1 and β 2 is (Reddi et al., 2019), which assumes β 2 > β 2 1 (stronger than our assumption). Our result indicates that the ranges of β 1 and β 2 can be broader in hyperparameter selection.
We still start the proof sketch by proving the sum of squared gradients is finite in Stage I. Compared to GDM, the difference is that in Stage I, we will also prove the loss converges to 0 with rate Θ(t −1 ). This property will be used in Stage II to analyze the effect of preconditioner on implicit bias.
Stage I: Bound the sum of squared gradients. The next lemma can be viewed as an extension of Lemma 1 by taking the preconditioner into consideration, and characterizes the one-step loss update in Adam. Lemma 5. Let all conditions in Theorem 4 hold. Then, for any t ≥ 1, The difference between the proof of Lemma 5 and Lemma 1 is that we need to handle the gap between the preconditioners at step t and step t + 1, this leads to a amplifying factor 4 , which, however, is smaller than the shrinking factor ) introduced by the gap between the coefficients of momentumm(t) and m(t − 1). As t → ∞, both β t 1 → 0 and (cβ 1 ) t → 0, and we can obtain that ξ(t) is a Lyapunov function. Based on Lemma 5, we can prove the following asymptotic rate of L. Lemma 6. Let all conditions in Theorem 4 hold. Then, The proof idea of Lemma 6 is that by the exponential-tailed property, when time is large enough, the training loss can be bounded by the gradient, and further bounded by Stage II. Bound the difference between w(t) and log(t)ŵ. The Lyapunov function for the direction difference is also modified according to the preconditioner, as introduced in Lemma 7:

Lemma 7. Let all conditions in Theorem 4 hold. Then, r(t) is bounded if and only if g(t) is upper bounded, where g(t) is defined as follows.
The first claim of Lemma 7 is similar to that in GDM and SGDM case. However, when we come to the second claim, simple calculation of g(t + 1) − g(t) leads to This is where we need Lemma 6, which bounds the gap between ε1 d +ν(t) and can be tackled by following the same routine as the GDM/SGDM case. Consequently, the proof of Lemma 7 is completed.

DISCUSSIONS
Extension to the Multi-Class Classification Problem. As mentioned in the introduction, despite all the previous analyses are aimed at the binary classification problem, they can be naturally extended to the analyses multi-class classification problem. Specifically, in the linear multi-class classification problem, for any (x, y) ∈ R dX × {1, · · · , C} in the sample set S, the (individual) logistic loss with parameter W ∈ R C×dX is denoted as Correspondingly, dataset S is separable if there exists a parameter W , such that ∀(x, y) ∈ S, we have W y, x > W i, x, ∀i = y. The multi-class L 2 max-margin problem is then defined as DenoteŴ as the L 2 max-margin solution, we have (mini-batch) SGDM and Adam still converges to the direction ofŴ , the proof of which is deferred to the Appendix D.2. Theorem 5. For linear multi-class classification problem using logistic loss and almost every separable data, with a small enough learning rate, and 1 > β 2 > β 4 1 ≥ 0 (for Adam), (mini-batch) SGDM and Adam converge to the multi-class L 2 max-margin solution (a.s. for SGDM).  (Lyu & Li, 2019) can be generalized, and successfully extend the Lyapunov function to AdaGrad, RMSProp and Adam (without momentum) on deep homogeneous neural networks. It will be interesting to see whether such a Lyapunov function can be constructed for GDM and Adam in homogeneous neural networks.

CONCLUSION
In this paper, we study the implicit bias of momentum-based optimizers in linear classification with exponential-tailed loss. Our results indicates that for SGD and Adam, adding momentum will not influence the implicit bias, and the direction of parameter converges to the L 2 max-margin solution.
Our theoretical results stands with existing experimental observations.

Supplementary materials for "Momentum Doesn't Change the Implicit Bias"
A PREPARATIONS This section collect definitions and lemmas which will be used throughout the proofs. To start with, we define support vectors and support set, which are two common terms in margin analysis.
Definition 2 (Support vectors and support set).
The following lemma deliversŵ as an linear combination of support vectors.
Lemma 8 (Lemma 12, Soudry et al. (2018)). For almost every datasets S, there exists a unique vector v = (v 1 , · · · , v N ), such thatŵ can be represented aŝ where v satisfies v i = 0 ifx i / ∈ T (S s ), and v i > 0 ifx i ∈ T (S s ). Furthermore, the size ofS s is at most d.
By Lemma 8, we further have the following corollary: Corollary 2. For almost every datasets S, the unique v given by Lemma 8 further satisfies that for any positive constant c, there exists a non-zero vectorŵ, such that,x i ∈ T (S s ), we have Proof. For almost every datasets S, any subsets with size d of S is linearly independent. SinceS s has size no larger than d (by Lemma 8), and Eq. (12) is equivalent to linear equations, the proof is completed.
For the stochastic case, we will also need the following lemma when we calculate the form of parameter at time t.

A.2 PREPARATIONS OF THE OPTIMIZATION ANALYSIS
This section collects technical lemmas which will be used in latter proofs. We begin with a lemma bounding the smooth constants if the loss is bounded.

Lemma 10. If loss ℓ satisfies (D) in Assumption 3, then for any
for any two parameters w 1 and w 2 close enough to w, and thus, which completes the proof that L is locally H s0 smooth at w.
Following the same routine as the locally smooth proof, we complete the second argument.
The proof is completed.
Based on Assumption 2, we also have the following lemma characterizing the relationship between loss ℓ and its derivative ℓ ′ when x is large enough.
Lemma 11. Let loss ℓ satisfy Assumption 2. Then, there exists an large enough x 0 and a positive real K, such that, ∀x > x 0 , we have Proof. By Assumption 2, there exists a large enough x 0 , such that ∀x > x 0 , we have On the other hand, as lim t→∞ ℓ(x) = 0, we have which by Eq. (13) leads to Therefore, setting K = 4 completes the proof.
The following lemma bridges the second moment of ∇L B(t) with its squared first moment.
Lemma 12. Let the dataset S satisfies the separable assumption 1. Let B be a random subset of S with size b sampled independently and uniformly without replacement. Then, at any point w, we have Proof. To start with, notice that Therefore, the first inequality can be directly obtained by Cauchy-Schwartz's inequality. To prove the second inequality, we first calculate the explicit form of ∇L B (w). Therefore, On the other hand, where Eq. (⋆) is due to ∀z ∈ x, −γŵ ≥ γ and ℓ ′ < 0. Therefore, The proof is completed by putting Eqs. (15) and (16) together.
In the following lemma, we show the updates of GDM, Adam, and SGDM are all non-zero. Lemma 13. Regardless of GDM, Adam, or SGDM, the updates of all steps are non-zero, i.e., Proof. We start with the alternative forms of the update rule of GDM, Adam, and SGDM using the gradients along the trajectory respectively. For GDM, by Eq. (5), the update rule can be written as Similarly, the update rule of SGDM can be written as while the update rule of Adam can be given as On the other hand, by the definition of empirical risk L, the gradient of L at point w can be given as By Eq. (20) and Eq. (17), we further have for GDM, By Assumption 1, there exists a non-zero parameterŵ, such that, ŵ,x i > 0, ∀i. Therefore, by executing inner product between Eq. (21) andŵ, we have where Eq. ( * ) is due to ℓ ′ < 0. This complete the proof for GDM.
Similarly, for SGDM, we have which completes the proof of SGDM.
For Adam, we have which completes the proof of Adam.
The proof is completed.

B IMPLICIT BIAS OF GD/SGD WITH MOMENTUM
This section collects the proof of the implicit bias of gradient descent with momentum and stochastic gradient descent with momentum.

B.1 IMPLICIT BIAS OF GD WITH MOMENTUM
This section collects the proof of Theorem 2.

B.1.1 PROOF OF THE SUM OF SQUARED GRADIENTS CONVERGES
To begin with, we will prove the sum of squared norm of gradients along the trajectory is finite for gradient descent with momentum. To see this, we first define the continuous-time update rule as We then prove a generalized case of Lemma 1 for any w(t + α). Lemma 14 (Lemma 1, extended). Let all conditions in Theorem 2 hold. We then have where c is a positive real such that η = 2 1−β Hs 0 c.
We further divide the proof into two cases depending on the value of α * .
By Lemma 10, we have L is H s0 smooth at w(t * ). Therefore, by Taylor's expansion for L at point w(t * ), we have for small enough α > 0 where Eq. ( * ) is due to a simple rearrangement of the update role of gradient descent with momentum (Eq. (5)), i.e., Inequality ( * * ) is due to Cauchy Schwarz's inequality and arithmetic-geometric average inequality, and Inequality (⋄) is due to Eq. (22) indicates Eq. (7) holds at (t * , α) for α > 0 is small enough, which contradicts to α * = 0.
The proof is completed.
By Lemma 1, one can easily obtain the sum of the squared norms of the updates across the trajectory converges.
By triangle inequality, we further have where Eq. (⋆) is due to Cauchy-Schwartz's inequality.
The proof is completed.
By the negative derivative of the loss and the separable data, we can finally prove Corollary 1.
By the exponential-tailed assumption of the loss (Assumption 2), we further have the following corollary.
Consequently, there exists an large enough time t 0 , such that, ∀t > t 0 , ∀i, we have w(t),x i > 0, and

B.1.2 BOUNDING THE ORTHOGONAL PART
To prove Theorem 2, we only need to show w(t) − log(t)ŵ (t ≥ 1) has bounded norm for any iteration t > 0. Letting c = η (1−β)N in Corollary 2, we obtain an constant vectorw satisfying Eq. (12). Define Asw is a constant vector, that w(t) − log(t)ŵ (t ≥ 1) has bounded norm is equivalent to r(t) has bounded norm. Lemma 2 then propose an equivalent proposition of r(t) is bounded, and further prove this proposition is fulfilled. As the proof is rather complex, we separate it into two sub-lemmas. We first prove r(t) is bounded if and only if function g(t) is upper bounded.

Lemma 15 (First argument in Lemma 2). Let all conditions in Theorem 2 hold. Then, r(t) is bounded if and only if function g(t) is upper bounded.
Proof. We start the proof by showing that A 1 (t) By the definition of r(t) , we have Therefore, the absolute value of A 1 (t) can be bounded as where Inequality (⋆) is due to the Inequality of arithmetic and geometric means, and Inequality (•) is due to Corollary 3 and log τ τ −1 = O( 1 τ ).

Therefore, g(t) is upper bounded is then equivalent to
is upper bounded, we will prove r(t) is bounded by reduction to absurdity. By Corollary 3, we have lim t→∞ w(t)− w(t− 1) = 0, and there exists a large enough time T , such that w(t) − w(t − 1) < 1 for any t ≥ T . On the other hand, since r(t) is unbounded from above, there exists an increasing time sequence k i > T , i ∈ Z + , such that lim i→∞ r(k i ) = ∞.
Therefore, we have which leads to contradictory, and completes the proof of necessity.
On the other hand, if r(t) is upper bounded, since w(t) − w(t − 1) is also upper bounded, we have 1 2 r(t) 2 + β 1−β r(t), w(t) − w(t − 1) is upper bounded, which completes the proof of sufficiency.
The proof is completed. Therefore, the last piece of this puzzle is to prove g(t) is upper bounded ∀t > 0.

Lemma 16 (Second argument in Lemma 2). Let all conditions in Theorem 2 hold. Then, we have g(t) is upper bounded.
Proof. We start the proof by calculating g(t + 1) − g(t). For any t ≥ 2, we have On the other hand, by simply rearranging the update rule Eq. (5), we have which further indicates Denote A 2 (t) = r(t + 1) − r(t) 2 , and A 3 (t) = r(t), − log t+1 t ŵ − η 1−β ∇L(w(t)) . We then prove respectively ∞ t=1 A 2 (t) and where Eq. (•) is due to Lemma 3 and log t+1 t = O( 1 t ). Then we only need to prove ∞ t=1 A 3 (t) < ∞. To begin with, by adding one additional term 1 tŵ into A 3 , we have On the one hand, by Corollary 3, w(t) = O( √ t), which further leads to On the other hand, by direct calculation of the gradient, we have We then analysis these two terms respectively. As for A 4 (t), due to ℓ ′ < 0, we have By Corollary 4, we further have ∀t > t 0 For each term r(t), − η we divide the analysis into two parts depending on the sign of r(t), x i .
We are now ready to prove Theorem 2.
Proof of Theorem 2. By Lemma 16, we have g(t) is upper bounded. Therefore, by Lemma 2, we have r(t) is bounded, which further indicates w(t) − log(t)ŵ is bounded. Therefore, the direction of w(t) can be calculated as The proof is completed.

B.2 IMPLICIT BIAS OF SGD WITH MOMENTUM
This section collects the proof of Theorem 3. Following the same framework as Appendix B.1, we will first prove that the sum of the squared gradient norms along the trajectory is finite. One may expect L(w(t)) + β 2η w(t) − w(t − 1) 2 is a Lyapunov function of SGDM. However, due to the randomness of the update rule of SGDM, L(w(t)) + β 2η w(t) − w(t − 1) 2 may no longer decrease (we will show this in the end of Appendix B.2, please see Appendix B.2.3 for explanation).

B.2.1 PROOF OF THE SUM OF GRADIENTS ALONG THE TRAJECTORY IS FINITE
We first provide a proof of Lemma 3.

Proof of Lemma 3.
We denote the parameter of the t-th step in Eq. (9) asw, while that in Eq. (1) as w(t). With w(1) =w(1), we will provew(t) = w(t), ∀t > 1 iteratively. Specifically, suppose for any 1 ≤ k ≤ t w(k) =w(k). Then, we have Therefore, which by iteration further indicates On the other hand, by Eq. (1), we have The proof is completed.
Using the alternative form given by Lemma 3, we then prove Lemma 4, which indicates L(u(t)) is a proper choice of Lyapunov function.
As L(u(1)) is upper bounded, we have the following corollary given by Lemma 4.
Proof. By Lemma 4, we have for any T > 1, which completes the proof of Eq. (38). The rest of claims follows immediately by Fubini's Theorem and Assumption 2.
The proof is completed.

B.2.2 BOUNDING THE ORTHOGONAL PART
Similar to the case of GDM, we definew as the solution of Eq. (12) with c = η (1−β)N . We also let n(t) be given by Lemma 9, and define r(t) in this case as Asw is a constant vector, and n(t) → 0 as t → ∞, we have w(t) − log(t)ŵ has bounded norm if and only if r(t) is upper bounded. Similar to the GDM case, we have the following equivalent condition of that r(t) is bounded.

Lemma 17. Let all conditions in Theorem 3 hold. Then, r(t) is bounded almost surely if and only if function g(t)
is upper bounded almost surely, where g : Z + → R is defined as Proof. To begin with, we prove that almost surely | t τ =2 r(τ ) − r(τ − 1), w(τ ) − w(τ − 1) | is upper bounded for any t. By Corollary 5, we have almost surly ∞ t=1 ∇L(w(t)) 2 < ∞.
On the other hand, for any w, we have Therefore, we have almost surely, which further leads to almost surely By the definition of r(t) (Eq. (39)), we further have where Inequality (⋆) is due to n(τ ) − n(τ − 1) = O( 1 τ ) and log τ +1 τ = O( 1 τ ). Therefore, g(t) is upper bounded almost surely is equivalent to 1 2 r(t) 2 + β 1−β r(t), w(t) − w(t − 1) is upper bounded, which can be shown to be equivalent with r(t) is bounded following the same routine as Lemma 2.
The proof is completed.
As the case of GDM, we only need to prove g(t) is upper bounded to complete the proof of Theorem 3.

Lemma 18. Let all conditions in Theorem 3 hold. We have g(t) is upper bounded.
Proof. Following the same routine as Lemma 2, we have where ∞ t=1 r(t + 1) − r(t) 2 is upper bounded.
On the other hand, by the definition of r(t) (Eq. (39)), we have while by Lemma 9, Combining the above two equations, we further have Therefore, we only need to prove By directly applying the form of ∇L(w(t)), we have We will investigate these two terms respectively. As w(t),x → ∞, ∀x ∈ T (S), a.s., we have a.s., there exists a large enough time t 0 , s.t., ∀t ≥ t 0 , ∀x ∈ T (S), x, w(t) > 0. Therefore, where Eq. (⋆) is due the definition of θ (Eq. (31)), and Eq. (•) is due to lim t→∞ e − w+n(t), On the other hand, A 6 (t) can be rewritten as If r(t),x i ≥ 0, we have for ε < 0.5, where Eq. (•) is due to n(t) = O( 1 t 0.5+ε ), and Eq. (⋄) is due to e − r(t),xi r(t),x i ≤ 1 e .
On the other hand, if r(t),x i < 0, we have If −2 > r(t),x i , then for large enough time t, e − r(t)+n(t),xi ≥ e 3 2 , 1 − e −µ− w(t),xi ≥ e − 1 2 , and 1 t e − w, Conclusively, if r(t),x i < 0, for large enough t, we have which further indicates, for large enough t, we have <∞.
The proof is completed.

B.2.3 EXPLANATION FOR PROPER LYAPUNOV FUNCTION
Based on the success of applying Lyapunov function L(w(t)) + β 2η w(t) − w(t − 1) 2 to analyze gradient descent with momentum, it is natural to try to extend this routine to analyze stochastic gradient descent with momentum. However, in this section, we will show such Lyapunov function is not proper to analyze SGDM as this will put constraints on the range of the momentum rate β. Specifically, at any step t, since the loss L is H smooth at w(t), we can expand the loss L in the same way as the GDM case: L(w(t + 1)) ≤L(w(t)) + w(t + 1) − w(t), ∇L(w(t)) + H 2 w(t + 1) − w(t) 2 .

C IMPLICIT BIAS OF ADAM
This section collects the proof of the convergent direction of Adam, i.e., Theorem 4. The methodology of this section bears great similarity with GDM, although the preconditioner of Adam requires specific treatment for analysis. The proof is still divided into two stages: (1). we first prove the sum of squared gradients along the trajectory is finite. Additionally, we prove the convergent rate of loss is O( 1 t ); (2). we prove w(t) − log(t)ŵ has bounded norm. Before we present these two stages of proof, we will first give the required range of η for which Theorem 3 holds.
The "sufficiently small learning rate" in Theorem 3 means .
To ensure η is well-defined, we need to prove and we introduce the following technical lemma: Proof. First of all, by definition, .
The left side of the above inequality is no smaller than 1, while the right side is no larger than 1, which completes the proof.
We are now ready to prove η is well-defined. First of all, for every t, we have where Eq. (⋆) is by Lemma 19 and cβ 1 = 4 √ β 2 < 1.
On the other hand, we have By Eq. (45) and Eq. (46), we obtain 1−(cβ1) t 1−(cβ1) t−1 is lower bounded by some positive constant across t, and η is well defined.

C.2 SUM OF GRADIENTS ALONG THE TRAJECTORY IS BOUNDED
We start with the following lemma, which indicates L(w(t))+ 4 ε1 d +ν(t)⊙(w(t)−w(t−1)) 2 is a proper Lyapunov function for Adam.

Preprint and under review
Eq. (47) contradicts the fact that the equality in Eq. (10) holds for 1 + α * , which completes the proof of t = 1.
This contradicts to that the equality in Eq. (10) holds for t * + α * .
The proof is completed.
and we only need to prove r(t) is bounded over time. We then prove Lemma 7, providing an equivalent condition of r(t) being bounded. As the GDM and SGDM case, we separate the proof into two sub-lemmas.
The proof is completed.

D.2 EXTENSION TO THE MULTI-CLASS CLASSIFICATION PROBLEM
Here we use several notations and lemmas from (Soudry et al., 2018). We define w = vec(W ), w = vec(Ŵ ), e i ∈ R C (i ∈ {1, · · · , C}) satisfying (e i ) j = δ ij , and A i = e i ⊗ I dX , where I dX is the identity matrix with dimension d X . We still consider the normalized data, i.e., x ≤ 1, ∀(x, y) ∈ S. Then, the individual loss of sample (x, y) can be then represented as ℓ(y, W x) = log e w,Ayx one can then easily verify all absolute value of the eigenvalues of HL(w) is no larger than 2, which indicates L is 2-globally smooth.
On the other hand, the separable assumption leads to ŵ, (A y − A i )x > 0, ∀y = i, which further indicates ∇L(w),ŵ > 0.
Let γ = 1 ŵ , following the similar routine as the binary case, we have for a random subset of S sampled uniformly without replacement with size b, we have Similarly, we have for any positive real series {a t } t2 t=t1 , γ t2 t=t1 a(t) ∇L(w(t)) ≤ t2 t=t1 a(t)∇L(w(t)) ≤ t2 t=t1 a(t) ∇L(w(t)) .
The proofs of Stage I can then be obtained with Lyapunov functions unchanged and by replacing the corresponding lemmas using Eq. (66) and Eq. (67).
The proof of Theorem 5 is then completed.


Achieving safe minimum circle circumnavigation around multiple targets: a dynamic compensation approach

Minimum circle circumnavigation is proposed in this paper, which is of special value in target monitoring, capturing and/or attacking. In this paper, a safe minimum circle circumnavigation of multiple targets based on bearing measurements is studied. In contrast with the traditional circumnavigation problem, with the new pattern, one agent is able to enclose multiple targets along a minimum circle with the desired enclosing distance and tangential speed. To achieve the minimum circle circumnavigation, an algorithm including dynamic compensators and a control protocol is proposed, by which collision is avoided between the agent and the multiple targets during the whole moving process. Moreover, the control protocol developed for a single agent is further extended to the scenarios of multiple agents by adding a coordination mechanism into the tangential velocity term, which drives the agents to distribute evenly on each expected circular orbit with the same radius or different radius. Simulations results illustrate the effectiveness of the proposed methods.

Introduction
Due to the wide applications in applications in satellite clusters, monitoring, rescue, etc, considerable effort has been made to study the coordination control of the multi-agent system including consensus, formation, containment control and so on [1][2][3][4]. In the area of cooperative control of multi-agent systems, circumnavigation is an important problem, which is to drive an agent/agents to enclose some targets or an area of interest. A classic method dealing with such problem is that the agent first localize the targets, and then encloses the targets with the desired radius and tangential speed. Previous literature taking this classic approach includes, but is not limited to, [5][6][7][8].
In the early studies on circumnavigations with only bearing information, [17] is a pioneering work, where an algorithm is proposed enabling an agent to enclose a target with the desired radius and tangential velocity. The target's position is directly estimated by constructing the errors of the target position through orthogonal projection and the errors converge to zero with the movement of the agent. In [18], this algorithm is further extended for a multi-target scene, where multiple agents enclose multiple targets with the desired radius.
In recent studies, other key issues are considered while solving the localization and circumnavigation problem. To closely detect or monitor the targets, one important problem is concerned with how to realize a circumnavigation pattern that enables an agent to enclose the targets closely. In addition, how to design an algorithm to avoid collision with the targets is also worth investigating when enclosing unknown targets with only bearing measurements.
There have been some research efforts to investigate such localization and circumnavigation problems. Different from the previous circular circumnavigation, the authors of [19] propose an elliptical circumnavigation pattern. Since elliptical circumnavigation is closer than the traditional circle circumnavigation, it shows a better performance in harsh environments. But collision avoidance with the targets is not considered and the distribution range of targets need to be known in advance. Although [20] proposes a safe circle circumnavigation pattern, where the distance between the agent and the targets is larger than a safe distance, the pattern is not close enough and the distribution of initial positions of the estimates of the targets is strictly limited as well. Subsequently, a new scalar estimator is designed in [21], which provides a collision avoidance mechanism with a loose initial condition. However, [21] does not consider multitarget scenarios. Inspired by [19,21], a safe convex-circumnavigation pattern is proposed in [22] for multiple targets. The pattern realizes the closer circumnavigation than elliptical circumnavigation under loose initial conditions. Moreover, there is no collision between the agent and the convex hull of multiple targets.
However, the convex-circumnavigation pattern is not good enough in engineering practice to design control law. The reason is that the curvature of the orbit in the convex hull vertex region is large, thus it is difficult to control the agent to turn in this region. It should be noted that the control law is different from the control protocol. The speed of the agent is regulated with the control protocol and the force acting on the agent is regulated with the control law. With this consideration, a smooth circumnavigation pattern is expected because smooth orbitals require simpler accelerations than complex orbitals. Inspired by the above analysis, in this paper, we consider the localization and safe minimum circle circumnavigation problem for multiple targets using bearing-only measurements, where the agent can safely circumnavigate a min-imum radius circle covering all targets with the desired enclosing distance and tangential velocity. To deal with the incomplete measurement information, dynamic compensators are first designed. Based on the dynamic compensation approach proposed in this paper, a control protocol is designed to make the agent achieve the minimum circle circumnavigation. To make the agent circumnavigate the targets closely enough, different from the control protocol proposed in [18,19], where the velocity of the agent is decomposed into the direction of the geometric center of the targets, we develop a new control protocol with the agent velocity pointing to the estimated minimum circle center and its orthogonal direction. In addition, by the developed method, there is no collision between the agent and the multiple targets during the whole moving process. Furthermore, the proposed method can be easily extended for the multiple agents scenarios, where an evenly spaced configuration along the desired minimum circle orbit is achieved.
The main contributions of this work are threefold. First, a minimum circle circumnavigation pattern is proposed. With this pattern, the agent can surround the targets more closely than the traditional circular circumnavigation, which can help the agent to detect the targets more thoroughly. Second, a dynamic compensation based anti-collision mechanism is established without knowing the distribution range of the targets, by which the minimum circle circumnavigation is achieved, and no collision happens under the loose initial constrained conditions. Third, compared with the convex-circumnavigation, under the minimum circle circumnavigation pattern, it is easier to design the control law and is more convenient to adjust the proposed algorithm for the multi-agent case.
The rest of this paper is structured as follows. In Section 2, we introduce some mathematical preliminaries and formally propose the problem of localization and safe minimum circle circumnavigation of multiple targets. In Section 3, we propose an algorithm to solve this problem and prove the convergence of the algorithm. In Section 4, we extend the proposed algorithm to multiple targets scenarios. Simulation results are shown in Section 5 and we demonstrate the conclusions and future works in Section 6.

Preliminaries and problem statement
In this section, we first define some notations used throughout the paper, and then briefly review some closely related theories.

Mathematical preliminaries and notations
Suppose there are n stationary targets with unknown positions and an agent moving in a 2D space. Both the targets and the agent are assumed to be modeled as points. Among them, the targets are represented by µ 1 , µ 2 ,..., µ n .
In the beginning, some definitions of mathematical symbols are introduced. ρ i (t) ∈ R ≥0 the distance between the agent and the target µ i , ρ i (t) the set of unit vectors from the agent to the targets at time t, the distance between the agent and the estimated position of the target µ i along the direction of ϕ i (t) ρ i (t) ∈ R the estimation error of the position of the target µ i ,ρ i (t) ϕ(t) ∈ U ϕ(t) rotates π/2 clockwise α ∈ R + the desired tangential enclosing speed r s ∈ R + the safety distancê r t (t) ∈ R ≥0 the radius of minimum circle of the estimated positions of n targets r t ∈ R + the radius of minimum circle of the real positions of n targets d ∈ R + the desired distance around the minimum circle r(t) ∈ R + the estimated radius of the desired orbit,r(t) r t (t) + d r ∈ R + the real radius of the desired orbit, r r t + d Here, R is the set of real numbers. R ≥0 is the set of non-negative real numbers.
R + is the set of positive real numbers. R 2 is the set of 2D column vectors. · is the Euclidian norm. U is the set of unit vectors, U {u ∈ R 2 | u = 1}.
Remark 1. Some points need to be explained.
(1)ρ i (t) is a variable that can be positive or negative.
(2) The minimum circle of a plane point set is unique, thusr t (t),ĉ(t) are unique at a certain time.
In this paper, the concept of the convex hull is used to analyze the security of the agent motion, for which the related definitions are introduced.
Definition 2 (Convex Hull of Finite Points Set). [22] Suppose is a subset of R 2 . We define conv(P ) as the convex hull of P , such that Definition 3 (Relationship Between a Point and a Convex Set). [22] Suppose Q ∈ R 2 is a set and a ∈ R 2 is a point. Define dist(a, Q) as the distance from a to Q, that is, If Q is compact and convex, then we can define prj(a, Q) as the minimum distance point in Q to a, which satisfies that If dist(a, Q) > 0, define uv(a, Q) as the unit vector pointing from a to Q, i.e., uv(a, Q) prj(a, Q) − a dist(a, Q) .

Remark 2.
In Definition 3, it should be noted that prj(a, Q) exists and unique.
With the concepts defined above, we give some definitions of the convex hull related to safe minimum circle circumnavigation.
Definition 4. Define: D(t) ∈ R ≥0 the distance from the agent to the convex hull of all of the targets, the distance from the agent to the convex hull of all of the targets, ψ(t) ∈ R 2 the unit vector pointing from the agent to the convex hull of all the unit vector pointing from the agent to the convex hull of all of the targets, ψ(t) uv(y(t),X c ) (ifD(t) > 0) The relationships among some notations defined above are illustrated in Figure 1 and Figure 2.
Then, some angle descriptions and comparison criteria of angle size are defined in advance for the convenience of the following description. First, u, v ∈ R ≥0 is defined as the angle between v and u.
is a subset of U. v r is defined as the rightmost unit vector of V denoted by ruv(V), and v l is defined as the leftmost unit vector of V denoted by luv(V).
Finally, we define the standard of angle size comparison. Note that the angle of v i is defined as the angle when v r is taken as the starting vector and rotated counterclockwise to v i .
Furthermore, we summarize some important results that will be used to analyze the convergence of dynamic compensators and control protocol.
Then we can assert that L P ≤ L Q . Proposition 3. [26] Suppose C ⊂ R 2 is a convex hull of a finite point set.
Define a function Let l uv(a, C). Then we can assert that f is differentiable and ∇f (a) = −l.  Proof. With Proposition 4, we know the center of the minimum circle is either the incenter of the triangle or the midpoint of the longest side of the triangle.
Obviously, the center of the circle can be expressed linearly by the three vertices of the triangle, which are all points on the convex hull. According to the definition of the convex set, we know that the center of the circle is also a point on the convex hull. Hence, the conclusion of this proposition follows readily.

Problem statements
Consider a multi-agent system containing multiple unknown static targets and an agent, the agent can only measure its own position as well as the bearing information from itself to each target, while the distance information from the agent to each target is unavailable. Assume that the agent is with single integrator dynamics, i.e., the velocity of the agent can be directly controlled.
The purpose of this paper is to design an algorithm for an agent to locate the unknown targets and enclose them with a safe minimum circle. Concretely, the objectives of the paper include: 1) An algorithm should be designed to locate the targets. That is, for all the targets, lim t→+∞ρi (t) = 0, i = 1, 2, ..., n.
2) A control protocol is designed to make the agent stably encircle the minimum circle of all the targets with the desired distance and tangential speed. That is, lim t→+∞ ρ(t) = r and lim t→+∞φ 3) Safety should be guaranteed. That is, the algorithm shall ensure that collision never happens.
Remark 3. The common definition of collision is that the agent does not collide with the targets. But the agent should avoid entering the area occupied by the targets to achieve real security in some application scenarios. Therefore, in this paper, we say collision happens at time t if D(t ) < r s .

Proposed algorithm
In this section, dynamic compensators are designed to compensate for the lack of measurement values of some state variables in the system, and then a control protocol is developed to achieve the safe minimum circle circumnavigation pattern through the action of the dynamic compensators.
First, the dynamics of the agent is described by the following nonlinear Here, x, u and z are the state, control input and control output of the agent respectively. f (·) and h(·) are unknown vector-valued functions. In fact, the speed of the agent is the control input, and the position of the agent is the system output. That is, u =ẏ and z = y.
Then, in the multi-agent circumnavigation control system, we expect to take the relative position of the agent and the targets as the state of the system.
Compared with the cartesian coordinates, in order to make full use of the known bearing information, we use the polar coordinates to describe the system state.
Therefore, we choose ϕ i and ρ i as the state variables of the system in this paper.
The state equation of this system is as follows In this paper, the agent does not know the position of the targets, i.e., the information of state variable ρ i are unknown, which are required for the design of the control input. Therefore, we need first to design dynamic compensators to compensate for the lack of the measured values of ρ i . The dynamic compensators are designed as follows Here,ρ i andx i are the dynamic compensator state variables and output respectively. In addition, sgn(·) is the sign function, u(·) is the step function and h is a positive constant.
Remark 4. In the design of dynamic compensators, some points need to be explained.
( (2) Since ϕ i and y can be measured, the term sgn(φ T iẏ )φ T iφ i is known to us until time t if we view the derivativeφ as left derivative. Therefore,ρ i can be solved.
Under the action of dynamic compensators, the control input (control protocol) is designed as follows Remark 5. As long asx i is determined, the estimated minimum circle is unique. So only know y andx i , one can getr t ,r,ρ,φ andφ.
Finally, to ensure that the agent can safely enclose the minimum circle of the targets, we need the following assumptions.
These assumptions are quite reasonable and weak. Assumption 1 is to ensure that collision does not happen in the beginning. Assumption 2 means that the desired enclosing distance is larger than the safety distance and the safety distance is larger than a small positive constant. Assumption 3 is a constraint on the initial value of each target's estimate and is loose.

Stability Analysis
In this section, we prove the stability and convergence property of the proposed algorithms. We first verify the safety can be guaranteed, and then prove that the agent eventually performs safe minimum circle circumnavigation around all of the targets successfully.

Safety
In this subsection, our goal is to confirm that collision does not happen in the whole time domain by using the proposed algorithm. To better explain that safety performance is guaranteed, the concept of the convex hull is applied.
However, it should be noted that in real implementation, it is not necessary to construct the convex hull. Denote Omitting "(t)" for simplicity, it follows thaṫ To ensure that the collision will not occur, we first give the important angle relations.
Then, based on the above lemma and under some restricted conditions, the following lemma is given to ensure the safety. Proof. Considering Proposition 2 and 0 <ρ i (t) ≤ ρ i (t), we knowD(t) ≤ D(t).
According to Lemma 1, from Assumption 2, we know that ψ T (t)φ(t) > 0 and ρ(t) <r(t) since d > r s , which implies thatḊ(t) > 0 if D(t) < r s . With Assumption 1, D(0) ≥ r s , then it follows that D(t) ≥ r s . This completes the proof.
Finally, the following two lemmas prove that the given constraint conditions are satisfied. Proof. It follows from (7)  Proof. First, we prove thatρ i (t) ≤ ρ i (t), that is,ρ i (t) ≤ 0, which is performed by three steps.
Step 1. Suppose there exists t 1 > 0 such thatρ i (t 1 ) = 0. Define the zero set Step 2. We prove that there isρ i (t) > h when t > t 2 . Apply the reduction to absurdity. Suppose there exists t 3 > t 2 such thatρ i (t 3 ) = h. Define the zero set B {t |ρ i (t) = h, t ≥ t 2 }. B is also a closed set. Define t 4 min{B}.

Circumnavigation
In this subsection, we prove that the agent finally performs minimum circle circumnavigation around all of the targets. The proof is divided into two steps.
First, we prove that the agent will circumnavigate the minimum circle of some points with the desired enclosing distance d and the tangential speed α. Then, we prove that these points are the real targets.
Subsequently, we prove that the agent will circumnavigate the minimum circle of some points with the desired enclosing distance and tangential velocity.
Here, c z is the center of the minimum circle of Z, and r z is the radius of the minimum circle of Z.
Proof. Since lim t→+∞xi (t) = z i , for any ε > 0, there exists t 1 > 0 such that Thus when t > t 1 , with Proposition 4, we know that the center of the minimum circle can be divided into the following two cases.
Case 1. The center of the circle is the midpoint of the longest side. We know lim t→+∞ĉ (t) = lim t→+∞ (x h (t) +x j (t))/2, wherex h (t),x j (t) are two points used to determine the estimated minimum circle. Similarly, Therefore, ĉ(t) − c z < ε for t > t 1 , which implies that lim t→+∞ĉ (t) = c z . It is obviously that lim t→+∞rt (t) = r z .
Case 2. The center of the circle is the incenter of a triangle. According to the definition of the incenter of a triangle, lim t→+∞ ĉ(t)−x k (t) = lim t→+∞ ĉ(t)− x l (t) = lim t→+∞ ĉ(t) −x p (t) , wherex k (t),x l (t),x p (t) are three points used to determine the estimated minimum circle. Then, Thusĉ(t) is also the incenter of the triangle by z k , z l and z p , which implies that c(t) is also the center of the minimum circle of Z. Therefore, lim t→+∞ĉ (t) = c z .
Then, it is obviously that lim t→+∞rt (t) = r z .
Therefore, the conclusion of this lemma follows readily.
Finally, we prove that the outputs of the dynamic compensator shall converge to the real targets and the agent performs the safe minimum circle circumnavigation around the real targets. Proof. The proof is divided into two steps.
Step 1. We prove that the localization tasks are achieved.
Suppose there exists i 0 , such that lim t→+∞ρi0 (t) = 0. Then with Lemma 6, we know ϕ i0 (t) must converge. With Theorem 2, it is obvious that the agent must circumnavigate the minimum circle of Z with the desired enclosing distance and tangential speed. Thus ϕ i0 (t) is not converging. Contradiction appears. Therefore, we have lim t→+∞ρi (t) = 0, for i = 1, 2, ..., n.
Step 2. We prove that the agent circumnavigates the minimum circle of the real targets with the desired distance d and tangential speed α.
With step 1, we have lim t→+∞xi (t) = x i , thus X = Z. With Theorem 2, we know that the conclusion of this theorem obviously holds.

Extension to multiple agents
In this section, we will show the aforementioned algorithm can be extended for the multiple agents case. Inspired by [24,25], for the multiple agents' case, it is required that the agents not only converge to the desired shaped orbit but also maintain an equiangular spaced formation around the targets.
However, different from [24,25], this paper studies a multi-target circumnavigation problem. We do not know the direction of the agent pointing to the desired orbital center. Therefore, we apply the direction pointing to the center of the minimum circle of the estimated targets instead. We define θ i as the angle ofφĉ i i , whereφĉ i i represents the unit vector pointing from agent i to the minimum circle center estimated by agent i. The reference axis of θ i is globally fixed in this paper.
Convert the variables to local versions by adding the subscript i. To be more practical, it is assumed that the agents cannot detect each other unless their distance is less than a threshold distance M . We denote the neighbor set of agent i as N i = {j | y j − y i ≤ M, j = i}. The bearing angle of agent j from agent i as θ j i being the angle of (y j − y i ), j ∈ N i and wrap It is worth pointing out that i + and i − can be determined only using bearings without communication. Suppose that positions ofĉ i + andĉ i − are available to agent i from its neighboring agent i + and i − via communication respectively. In addition, as long as the neighbor agents of agent i are within its detection range, the agent i can obtain the location of its neighbor agents. Then we define To achieve circumnavigation of multiple agents for the targets at equiangular spaced formation on the same circular orbit, in fact, only the tangential velocity term in (4) needs to be modified. Thus, we define the new control protocol aṡ whereρĉ i i (t) represents the distance between the agent i and the center of the minimum circle estimated by it,r i (t) represents the desired orbit radius estimated by agent i,r i,t (t) is the minimum circle radius estimated by agent i, and f : (0, +∞] → (0, 1] is any strictly increasing continuous function. For instance, we can let f (x) = 1 − exp(−x).
Define the angle of −φĉ i i as φ i = θ i + π mod 2π. For convenience, we label agents at time t according to their positions in a counterclockwise radial order around the targets, i.e., j > k ⇔ φ j > φ k . We also define δφ i = φ i+1 − φ i for i ∈ {1, ..., n − 1} (n ≥ 3) and define δφ n = φ 1 + 2π − φ n . Under this labeling, Proof. Based on the analysis in Section 3, because f (σ + i /σ − i )α ∈ (0, α], then for all i,ρĉ i i → r andĉ i = c as t → ∞. The security of the agents is also guaranteed. The condition M ≥ 2r sin(π/n) guarantees that neither N + i nor N − i is empty around the equilibrium point. Because the minimum circle is strictly convex, it can be easily proved that i + = i+1 and i − = i−1 around the equilibrium point.
Therefore, whenρĉ i i = r andĉ i = c for all i, we have(see Figure 3) σ + i = δφ i and σ − i = δφ i−1 . With the new control protocol (11), we know that the tangential velocity of each agent is  Remark 8. For the multiple agents' case, it is required that the agents not only converge to the different desired circle orbit but also maintain an equiangular spaced formation around the targets. Similarly, we just need to multiply the tangential velocity term in (11) by their respective desired radius.

Simulation
In this section, a series of simulation results for the safe minimum circle circumnavigation algorithm of four targets are presented. distance. Figure 5 describes the distance between the agent and the real minimum circle, which shows that the distance eventually converges to the desired enclosing distance d. Figure 6 describes the image ofφ T (t)ẏ(t), which means that the tangential velocity of the agent enclosing the minimum circle will eventually converge to the desired velocity α. Figure 7 depicts the estimation error ρ i (t) and shows that all the errors converge to 0, which means that the dynamic compensator outputs will converge to the real target's positions.
To sum up, it can be concluded that there is no collision in the whole circumnavigation process, and the dynamic compensators and control protocol eventually converge. Therefore, the safe minimum circle circumnavigation is successfully performed by the proposed algorithm.    Figure 8(a), which shows that multiple agents safely circumnavigate the minimum circle of multiple targets with the desired enclosing distance. The angular configuration is demonstrated by evolution of δφ i for i = 1, 2, 3, 4. It can be seen from Figure 8(b) that each δφ i converges to π/2. One can observe that all the four agents converge to the desired minimum circle of the real targets and achieve even distribution over the minimum circle as expected. Meanwhile, the security of the agents is also guaranteed.  Case 3. We further consider the case that multiple agents distributed in different radius with equiangular spaced formation. It can be seen in Figure 9(a) and Figure 9(b) that the four agents still converge to the desired minimum circle of the real targets and achieve even distribution over the different minimum circles as expected. Meanwhile, the security of the agents is also guaranteed.


SGI: automatic clinical subgroup identification in omics datasets

Abstract Summary The ‘Subgroup Identification’ (SGI) toolbox provides an algorithm to automatically detect clinical subgroups of samples in large-scale omics datasets. It is based on hierarchical clustering trees in combination with a specifically designed association testing and visualization framework that can process an arbitrary number of clinical parameters and outcomes in a systematic fashion. A multi-block extension allows for the simultaneous use of multiple omics datasets on the same samples. In this article, we first describe the functionality of the toolbox and then demonstrate its capabilities through application examples on a type 2 diabetes metabolomics study as well as two copy number variation datasets from The Cancer Genome Atlas. Availability and implementation SGI is an open-source package implemented in R. Package source codes and hands-on tutorials are available at https://github.com/krumsieklab/sgi. The QMdiab metabolomics data is included in the package and can be downloaded from https://doi.org/10.6084/m9.figshare.5904022. Supplementary information Supplementary data are available at Bioinformatics online.

Introduction
The identification of patient subgroups from high-dimensional molecular profiles has become a central approach in biomedical research, driven by the wide availability of modern 'multi-omics' datasets (Collisson et al., 2011;Eddy et al., 2020;Rouzier et al., 2005). The central idea is that genomics, transcriptomics, proteomics, metabolomics and other deep molecular phenotypes will inherently define groups of patients that are similar with respect to disease-relevant clinical outcomes. Note that, 'outcome' is here defined in a statistical sense, and includes parameters such as sex, prevalent disease and current BMI. Recent examples of molecularly defined subgroups include the identification of subtypes of lymphoma that severely impact survival (Nowakowski and Czuczman, 2015), subtypes of various cancers identified in the 'The Cancer Genome Atlas' (Sanchez-Vega et al., 2018) and patient stratification in allergic diseases (Agache and Akdis, 2019).
Various bioinformatics methods for the computational extraction of subgroups from omics data have been published in recent years. Most of those methods implement novel metrics for the pairwise similarity of samples in a multi-omics setting (Loh et al., 2019;Rappoport and Shamir, 2019;Speicher and Pfeifer, 2015;Wang et al., 2014) followed by standard clustering, or improved clustering approaches addressing the robustness of resulting subgroups (Nguyen et al., 2017).
Here, we present the 'SGI' (subgroup identification) package, which implements a new method for the automatic detection of omics-based subgroups. It provides the following novel features compared to previously published methods: (i) the algorithm works hierarchically and can thus identify subgroups of any granularity in complex patient cohorts. (ii) Any sample-wise distance metric can be used for hierarchical clustering, including classical Euclidian distance, or more advanced measures such as similarity network fusion (Wang et al., 2014) and multiple kernel learning (Speicher and Pfeifer, 2015). (iii) It can handle an arbitrary number of clinical variables of interest simultaneously. Combined with the hierarchical approach, this allows the user to explore the complex relationships between various outcomes. (iv) The algorithm is assumption-free and does not require any model fitting, since it only operates on a distance matrix across the omics samples. Moreover, the toolbox implements a comprehensive set of methods to visualize the associations for further interpretation.

SGI method
The SGI algorithm generates a hierarchical clustering of samples and runs a two-group association test against the analyzed clinical outcomes at each branching point in the tree. An example output

573
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Applications Note plot is shown in Figure 1, which will be further discussed in the application section below. The algorithm works as follows: (i) Clustering. A dendrogram of the samples is generated using standard hierarchical clustering on the input data matrix. The function accepts any hclust object, giving the user full control over the choice of distance and linkage functions. (ii) Generate valid cluster pairs. To avoid low-powered calculations in small clusters, the method enumerates all branching points where both left and right subclusters are above a user-defined size threshold. This results in a list of 'valid' cluster pairs. SGI runs with a default setting of 5% of the sample size. (iii) Association analysis. Run statistical tests with all clinical outcomes for all valid cluster pairs, i.e. compare left versus right subcluster at the respective branching points. SGI has builtin implementations for categorical outcomes (Fisher's exact tests), continuous outcomes (two-sample t-tests) and survival outcomes (log-rank tests). The appropriate test is automatically determined by the toolbox based on the data type of the clinical variable. Furthermore, the user can define arbitrary association functions for more complex data types. (iv) Multiple testing correction. Since a dendrogram clusters the samples into cascaded, non-overlapping groups, all statistical tests are strictly independent. Thus, SGI performs Bonferroni multiple testing correction by adjusting each Pvalue by a factor of the number of valid cluster pairs. In addition to its core functionality, SGI provides support functions to extract, plot, print and summarize all clustering and association results and test statistics, which allows the user to access all intermediate results to further analyze the subgroup results.

Visualization
Association results with multiple outcomes on a hierarchical tree are inherently complex to visualize. The SGI toolbox provides a variety of dedicated functions to visually inspect the statistical associations. This includes tree visualizations of all statistically significant outcome associations that are displayed at the respective branching points and heatmaps of the actual data (Fig. 1). Moreover, the user can generate plots to inspect specific associations, e.g. boxplots of a quantitative clinical outcome between two clusters.
This allows the user to obtain a quick overview of the correlation structure between the input omics dataset, the resulting patient groups and the clinical features that are analyzed. The simultaneous visualization of all data also allows to dissect the relationship Fig. 1. Application example. Blood metabolomics-based clustering of n ¼ 356 participants of the QMdiab study. White circles in the tree indicate the left/right splitting points of the samples in the data (note that, these are not centered if the subclusters are of unequal size). Markings on the tree indicate statistically significant associations of the parameter with the respective left and right subgroups at that split. Heatmap track below the tree shows individual values for selected parameters. Red circles between gaps indicate significant results for left versus right at that split and are horizontally aligned with their respective white circles on the tree. Bottom panel shows the metabolomics data matrix behind the clustering between confounding variables, avoiding the predefined choice of a list of confounder variables to correct for.

Multi-omics datasets
The SGI package provides clustering functionality for the analysis of multi-omics datasets, i.e. datasets where more than one omics layer has been measured for the same samples. To this end, SGI generates a joint samples X samples distance matrix from the individual distance matrices of each omics layer. Since different omics layers will have varying numbers of variables, the respective distance values are not at comparable scales. The toolbox thus normalizes each individual distance matrix by its maximum and defines D ¼ D 1 =maxðD 1 Þ þ D 2 =maxðD 2 Þ þ Á Á Á þ D l =maxðD l Þ, with D representing the final distance matrix, D 1 etc. the original distance matrices and l the number of omics datasets. The approach was adapted from Chavent et al. (2018), where it was originally introduced to generate a Ward-like clustering. Notably, the method also works for the normalization of multi-omics contributions for other types of linkages, such as average linkage and complete linkage, and works with any distance metric.
A detailed example of the multi-omics capabilities of SGI on combined plasma, urine and saliva metabolomics data can be found in the example R codes of the github repository.

Application examples
We will demonstrate the functionality of the SGI package on plasma metabolomics dataset from the 'QMdiab' diabetes case/control study with 356 participants (Do et al., 2018;Mook-Kanamori et al., 2014). We chose Euclidean distance with Ward linkage for hierarchical clustering. Outcome parameters were type 2 diabetes diagnosis and nine anthropometric and clinical biochemistry parameters: age, sex, BMI, HbA1c, albumin, hemoglobin, LDL cholesterol, total cholesterol and skin auto-fluorescence (AF score) (Mook-Kanamori et al., 2013). The goal was to determine clusters of study participants defined by their profiles of circulating metabolites, and how these metabolomic clusters correlate with the different clinical parameters. The resulting visualization is shown in Figure 1. The following lines of code directly generate that plot using the SGI package: # hierarchical clustering hc ¼ hclust (dist(sgi::qmdiab_plasma), method ¼ "ward.D2") # initialize SGI structure sg ¼ sgi_init(hc, outcomes ¼ sgi::qmdiab_clin) # run SGI as ¼ sgi_run(sg) # generate tree plot, show results for adjusted p-values <¼0.05 (gg_tree ¼ plot(as, p_th ¼ 0.05)) # plot overview, including clinical data and metabolomics data matrix plot_overview (gg_tree ¼ gg_tree, as ¼ as, outcomes ¼ sgi::qmdiab_clin, xdata ¼ sgi::qmdiab_plasma) In this R code, sgi::qmdiab_plasma and sgi::qmdiab_clin are data frames holding the metabolomics and clinical variables, respectively. These data frames are contained in the package. The tree shows how metabolomic profiles separate study participants into two major groups at the top level of the tree (clusters 2 versus 3), one with a higher proportion of males, higher BMI, higher HDL-C, as well as higher AF score, and the other group with the reversed effect directions. Inside those two groups, further subgroups were identified; for example, two clusters with different proportions of diabetes, which also associate with Hba1c, age and further diabetesrelated risk factors (clusters 4 versus 5).
In addition to the diabetes showcase, we ran SGI on two datasets from The Cancer Genome Atlas, in order to assess whether the algorithm can reproduce known molecular subtypes of cancer. In the first example, we recovered IDH-related mutational subgroups (Ceccarelli et al., 2015;Sanchez-Vega et al., 2018) using SGI on copy number variation data from low-grade glioma samples (Supplementary Material S1). The second example demonstrates how SGI identifies previously described subgroups in uterine corpus endometrial carcinoma samples (Levine et al., 2013;Sanchez-Vega et al., 2018), where both our groups and the originally reported groups were based on copy number variation measurements (Supplementary Material S2).

Conclusion
SGI provides a flexible, unbiased and data-driven way to automatically identify sample subgroups in omics profiles. It identifies and visualizes complex, hierarchical relationships for an arbitrary number of clinical outcomes in a visually intuitive way. The toolbox is easy to use, open source and comes with a series of examples in the online repository.