Distributed support-vector-machine over
dynamic balanced directed networks
Mohammadreza Doostmohammadian, Alireza Aghasi,  Themistoklis Charalambous, and  Usman A. Khan
MD is with the Faculty of Mechanical Engineering at Semnan University, Semnan, Iran, and with the School of Electrical Engineering at Aalto University, Espoo, Finland, . AA is with the Institute for Insight,  Robinson College of Business at Georgia State University, GA, USA, . TC is with the School of Electrical Engineering at Aalto University,
Espoo, Finland, . UAK is with the Electrical and Computer Engineering Department at Tufts University, MA, USA, . The work of UAK is partially supported by NSF under awards #1903972 and #1935555.

In this paper, we consider the binary classification problem via distributed Support-Vector-Machines (SVM), where the idea is to train a network of agents, with limited share of data, to cooperatively learn the SVM classifier for the global database. Agents only share processed information regarding the classifier parameters and the gradient of the local loss functions instead of their raw data. In contrast to the existing work, we propose a continuous-time algorithm that incorporates network topology changes in discrete jumps. This hybrid nature allows us to remove chattering that arises because of the discretization of the underlying CT process. We show that the proposed algorithm converges to the SVM classifier over time-varying weight balanced directed graphs by using  arguments from the matrix perturbation theory.

upport Vector Machine, constrained distributed convex optimization, matrix perturbation theory.

§ INTRODUCTION

Machine learning has been an area of significant research in recent signal processing and control literature [1, 2, 3, 4]. Among the topics of interest, Support Vector Machines (SVMs) are supervised-learning methods with several applications ranging from image/video processing to bioinformatics. Motivated by the recent progress in computing hardware and wireless communication, we are interested in developing distributed solutions for SVM classification. The basic idea is to process the raw data at each node in order to train a local classifier and then fuse these classifiers among the neighboring nodes. D-SVM (distributed SVM) finds applications where a subset of the data is acquired by different nodes/servers/agents possibly at different geographic locations, privacy is of concern, and communicating data to a fusion center (FC) is practically infeasible.

In binary classification, SVM defines the maximum-margin hyperplane (as the classifier) determined by the closest data samples known as the Support-Vectors (SVs). The preliminary work on D-SVM (referred as Distributed Parallel SVM (DP-SVM) [5] and Parallel SVM (P-SVM) [6]) is focused on local calculations and sharing of the SVs, as the representatives of discriminant information of
the dataset [7, 5, 6, 8, 9]. ${\mbox{In~[7, 6, 8]}}$, these local SVs are updated via a FC to improve the D-SVM performance. Ref. [5] implements D-SDV on multi-agent networks and requires a Hamiltonian cycle that visits every agent exactly once. A FC-free approach is considered in [9], where every agent locally solves a (coupled) convex optimization sub-problem via alternating direction method of multipliers (ADMM), which is not computationally efficient. A main drawback is that these approaches require sharing raw data over the communication network, raising data privacy and information security issues. More recently, consensus-based distributed optimization methods are proposed in [10, 11, 12, 13, 14, 15, 16, 17, 18], where instead of raw data, agents share processed information, which in case of leakage to unauthorized parties reveals little information about the original data. Among these, the solution in [16] requires distributed computation of the Hessian inverse, which is not practical since even for a sparse Hessian matrix, its inverse is not necessarily sparse. Penalty-based approaches are proposed in [17, 18], where the constrained convex cost function is reformulated by adding a new penalty term on consensus constraint violation. It is shown that there is a gap of $\mc{O}(\kappa)$ (with $\kappa$ as the penalty constant) between the optimal penalty-based solution and the original constrained one [19].
Other methods include finite/fixed-time algorithms [11, 12, 13, 14, 15] that are prone to
steady-state oscillations (known as chattering) due to non-Lipschitz dynamics.

In this paper, a D-SVM method is proposed that overcomes the challenges of (semi-centralized) FC-based solution and the chattering phenomena. Moreover, in contrast to Refs. [10, 11, 12, 13, 14, 15, 16, 17, 18], where either continuous-time (CT) or discrete-time (DT) protocols are considered, we propose a hybrid algorithm to address the topology switching of the multi-agent network in DT incorporated in a CT gradient-descent update [20]. Our hybrid approach enables more flexibility in considering mixed-dynamics[21, 20], which allows solving D-SVM via CT protocols over general dynamic digraphs in DT domain. To analyze the proposed hybrid model, we use matrix perturbation theory [22] to characterize the eigenspectrum of the proposed dynamics. The proposed D-SVM is fully distributed, as opposed to FC-based approaches, and does not require solving convex sub-problems unlike [7, 5, 6, 8, 9]. Due to Lipschitz-continuity of the proposed CT approach, it's DT approximation is free of the aforementioned chattering inherent to the non-Lipschitz dynamics [11, 12, 13, 14, 15]. Note that we directly solve the original constrained optimization  free of penalty-based approximation inaccuracies in [17, 18, 19].

We now describe the rest of the paper. Section <ref> recaps some preliminaries on algebraic graph theory while Section <ref> formulates the D-SVM problem. Section <ref> states our CT gradient descent method to address D-SVM whereas the convergence analysis over dynamic WB-digraphs is available in Section <ref>. Section <ref> provides an illustrative example, and finally, Section <ref> concludes the paper with some future research directions.

§ PRELIMINARIES ON ALGEBRAIC GRAPH THEORY

We represent the multi-agent network by a strongly-connected directed graph (SC digraph) $\mc{G}$.
Assuming a positive weight $w_{ij}$ for every link (from node $j$ to node $i$) and zero otherwise,  the irreducible adjacency matrix of $\mc G$ is $W=\{w_{ij}\}$, and the Laplacian matrix $\overline{W}=\{\overline{w}_{ij}\}$ is defined as,

\begin{eqnarray} \label{eq_laplac}
\overline{w}_{ij}=\left \{ \begin{array}{ll} w_{ij} & ,~ i\neq j, \\ -\sum_{j=1}^n w_{ij} & ,~ i=j. \end{array}\right.
\end{eqnarray}

The SC property of the graph is directly related to the rank of its Laplacian matrix as given in the next lemma.

[23] The given Laplacian $\overline{W}$ in (<ref>) for a SC digraph
has  eigenvalues whose real-parts are non-positive with one isolated eigenvalue at zero.

Next, we define a WB-digraph as an SC digraph with equal weight-sum of incoming and outgoing links at every node $i$, i.e., $\sum_{j=1}^n w_{ji} = \sum_{j=1}^n w_{ij}$,
implying the following lemma.

[23] For the Laplacian $\overline{W}$ of a WB-digraph, the vectors $\mb{1}_n^\top$ and $\mb{1}_n$ are respectively the left and right eigenvector associated with the zero eigenvalue, i.e., $\mb{1}_n^\top \overline{W}= \mb{0}_n$ and $\overline{W}\mb{1}_n=\mb{0}_n$, where $\mb{1}_n$ and $\mb{0}_n$ are the column vectors of $1$'s and $0$'s of size $n$, respectively.

In the rest of the paper, ${\lVert A\rVert_{\infty}}$ denotes the infinity norm of a matrix, i.e., $\lVert A\rVert_{\infty}= \max_{1\leq i\leq n} \sum_{j=1}^n |a_{ij}|.$

§ PROBLEM STATEMENT

Consider a binary classification problem for a given set of $N$ data points $\boldsymbol{\chi}_i \in \mathbb{R}^{m-1}$, $i=1,\ldots,N$, each belonging to one of two classes labeled by ${l_i \in \{-1,1\}}$. Using the entire training set, the SVM problem is to find a hyperplane ${\boldsymbol{\omega}^\top \boldsymbol{\chi} - \nu =0}$,  for ${\boldsymbol{\chi}\in\mbb R^{m-1}}$, based on the maximum margin linear classification to partition the data into two classes. Subsequently, a new test data point $\widehat{\boldsymbol{\chi}}$ belongs to the class labeled  as ${g(\widehat{\boldsymbol{\chi}})= \text{sgn}(\boldsymbol{\omega}^\top \widehat{\boldsymbol{\chi}} - \nu)}$.  In case the data points are not linearly separable, the input data is first projected into a high-dimensional space $\mc{F}$ via a nonlinear mapping $\phi(\cdot)$. This mapping is such that the inner products of two projected data points can be computed via a kernel function $K(\cdot)$, i.e., $K(\boldsymbol{\chi}_i,\boldsymbol{\chi}_j)=\phi(\boldsymbol{\chi}_i)^\top \phi(\boldsymbol{\chi}_j)$. By proper selection of $\phi(\cdot)$, a linear optimal  hyperplane defined by $\boldsymbol{\omega}$ and $\nu$ can be found in $\mc{F}$ such that ${g(\widehat{\boldsymbol{\chi}})= \text{sgn}(\boldsymbol{\omega}^\top \phi(\widehat{\boldsymbol{\chi}}) - \nu)}$ determines the class of $\widehat{\boldsymbol{\chi}}$. In centralized SVM, all the data points are sent to a central computation entity (the FC) that finds the optimal $\boldsymbol{\omega}~$ and $\nu$ by minimizing the following  convex loss function [24]:

\begin{equation} \label{eq_svm_cent}
\begin{aligned}
\displaystyle
& \min_{\boldsymbol{\omega},\nu}
~ &  \boldsymbol{\omega}^\top \boldsymbol{\omega} + C \sum_{j=1}^{N} \max\{1-l_j( \boldsymbol{\omega}^\top \phi(\boldsymbol{\chi}_j)+\nu),0\}^q
\end{aligned}
\end{equation}

where ${q = \{1,2,\ldots\}}$ defines smoothness of the loss function and its derivatives, and the positive constant $C$ determines the trade-off between increasing the margin size and ensuring that the projected data  $\phi(\boldsymbol{\chi}_j)$ lies on  correct side of the hyperplane. We note that the SVM loss function (<ref>) is not continuously twice-differentiable for $q=\{1,2\}$. Therefore, it is common to approximate $\max\{z,0\}^q$ for ${q=1}$ by the smooth function ${L(z,\mu)=\frac{1}{\mu}\log (1+\exp(\mu z))}$, which is the integral of the well-known sigmoid function [12]. It can be shown that by setting $\mu$ large enough ${L(z,\mu)}$ becomes arbitrarily close to $\max\{z,0\}$; see [25] for more smooth loss functions, e.g., for logistic regression with cross-entropy loss.

In distributed SVM (D-SVM), the data points are available over a network of $n$ agents and each agent $i$ possesses a local dataset with $N_i$ data points denoted  by $\boldsymbol{\chi}^i_j,{j=1,\ldots,N_i}$.  Since each agent has access to partial data, the locally found values $\boldsymbol{\omega}_i$ and $\nu_i$, obtained by solving (<ref>) over the local dataset $\boldsymbol{\chi}^i_j,{j=1,\ldots,N_i}$,

may differ for each agent $i$. The idea behind D-SVM is thus to develop a distributed mechanism to learn the global classifier parameters by making sure that no agent reveals its local data to any other agent. The corresponding distributed optimization  problem is given by:

\begin{equation} \label{eq_svm_dist}
\begin{aligned}
\displaystyle
\min_{\boldsymbol{\omega}_1,\nu_1,\ldots,\boldsymbol{\omega}_n,\nu_n}
\quad &  \sum_{i=1}^{n} f_i(\boldsymbol{\omega}_i,\nu_i) \\
\text{subject to} \quad&  \boldsymbol{\omega}_1 = \dots = \boldsymbol{\omega}_n,\qquad\nu_1 = \dots =\nu_n,
\end{aligned}
\end{equation}

where each local cost ${f_i:\mathbb{R}^m\rightarrow\mbb R}$ is approximated as (with ${z=1-l_j( \boldsymbol{\omega}_i^\top \phi(\boldsymbol{\chi}^i_j)+\nu_i)}$ and large enough $\mu>0$)

\begin{equation*}
f_i(\boldsymbol{\omega}_i,\nu_i)=\boldsymbol{\omega}_i^\top \boldsymbol{\omega}_i + C \sum_{j=1}^{N_i} \tfrac{1}{\mu}\log (1+\exp(\mu z)).
\end{equation*}

Let ${\mb{x}_i = [\boldsymbol{\omega}_i^\top;\nu_i]}\in\mathbb{R}^m$ and let $\mb{x} \in \mathbb{R}^{mn}$ be the global vector concatenating all $\mb{x}_i$'s, i.e.,
${\mb{x} = [\mb{x}_1;\mb{x}_2;\dots;\mb{x}_n]}$, where the symbol `$;$' denotes the column concatenation of the $\mb{x}_i$ vectors. Then, Problem (<ref>) takes the following form:

\begin{align}\nonumber
\min_{\mb{x} \in \mathbb{R}^{mn}}
F(\mb{x}),\qquad F(\mb x) = \sum_{i=1}^{n} f_i(\mb{x}_i)\\\label{eq_prob} \text{subject to} ~ \mb{x}_1 = \mb{x}_2 = \dots = \mb{x}_n.
\end{align}

We next provide the following lemma on the local costs.

[12] Each local cost $f_i$ is twice differentiable and strictly convex, i.e., the $m\times m$ Hessian matrix $\boldsymbol{ \nabla}^2 f_i(\mb{x}_i)$ is positive definite, for all non-zero $\mb x_i\in\mathbb R^m$.

Clearly, any solution $\mb x_i^*,{i=1,\ldots,n}$, of (<ref>) must satisfy ${\sum_{i=1}^{n} \boldsymbol{ \nabla} f_i(\mb{x}^*_i) = \mb{0}_m}$, such that ${\mb{x}^*_1=\ldots=\mb{x}^*_n=\overline{\mb{x}}^*}$, for some ${\overline{\mb{x}}^*\in\mathbb{R}^m}$. In other words, the optimality condition ${\boldsymbol{ \nabla} F(\mb{x}^*) = \mb{0}_{mn}}$ must hold for some ${\mathbf x^*\in\mathbb{R}^{mn}}$ such that ${\mb{x}^*=\mb{1}_n \otimes \overline{ \mb{x}}^*}$, where ${\boldsymbol{\nabla} F:\mathbb R^{mn}\rightarrow\mathbb R^{mn}}$ is the gradient of ${F:\mbb R^{mn}\rightarrow\mbb R}$.

§ PROPOSED ALGORITHM:
DYNAMICS AND AUXILIARY RESULTS

We now provide a distributed solution to  the D-SVM problem.
Let $\mb{x}_i(t)\in\mathbb R^m$ to be the state of agent $i$ at time $t$, where $t\geq 0$ is the continuous-time variable. To solve problem (<ref>), we consider the following continuous-time linear dynamics for all $\mb{x}_i(t) \in \mathbb{R}^{m},i\in\{1,\ldots,n\}$,

\begin{eqnarray} \label{eq_xdot}
\dot{\mb{x}}_i = -\sum_{j=1}^{n} w_{ij}(\mb{x}_i-\mb{x}_j)-\alpha \mb{y}_i,
\end{eqnarray}

where ${\dot{\mb{x}}_i=\frac{d{\mb{x}}_i}{dt}}$, ${W=\{w_{ij}\}}$ is the weighted adjacency matrix associated with  $\mc{G}$, and $\alpha>0$ is the stepsize.  We note that instead of the standard descend direction of $\nabla f_i(\mathbf{x}_i)$, the $\mathbf{x}_i$-update in (<ref>)
descends in the direction of an auxiliary variable $\mb y_i(t)\in\mathbb{R}^m$. The variable $\mb y_i(t)$ in fact tracks the sum of local gradients, asymptotically, and is updated via the following dynamics (see [1, 2, 3] for similar DT methods):

\begin{eqnarray} \label{eq_ydot}
\dot{\mb{y}}_i &=& -\sum_{j=1}^{n} a_{ij}(\mb{y}_i-\mb{y}_j) + \frac{d}{d t} \boldsymbol{ \nabla} f_i(\mb{x}_i),
\end{eqnarray}

where $\dot{\mb{y}}_i=\frac{d{\mb{y}}_i}{dt}$ and the matrix $A=\{a_{ij}\}$ is the weighted adjacency matrix with the same structure as $W$. In (<ref>),

\begin{eqnarray} \label{eq_dtdf}
\frac{d}{d t} \boldsymbol{ \nabla} f_i(\mb{x}_i)=  \boldsymbol{ \nabla}^2 f_i(\mb{x}_i)  \dot{\mb{x}}_i.
\end{eqnarray}

Note that the proposed algorithm, (<ref>) and (<ref>), is in continuous-time. However, the structure of the underlying graph $\mathcal G$ may change in time instances, that we consider as time steps in a discrete-time framework. This makes the proposed dynamics hybrid where the state variables, $\mathbf x$ and $\mathbf y$, evolve in CT over DT switching of the network topology.
For the ease of notation, we define an auxiliary global variable ${\mb{y}=[\mb{y}_1;\mb{y}_2;\dots;\mb{y}_n]\in\mbb{R}^{mn}}$ that concatenates the local $\mb{y}_i(t)$'s.
We make the following assumption on the weight matrices $W$ and $A$.

The weights ${W=\{w_{ij}\}}$ and ${A=\{a_{ij}\}}$ are associated to a WB-digraph with ${w_{ij}\geq0}$ and ${a_{ij}\geq0}$, respectively.

Further, $\sum_{j=1}^n w_{ij}<1$ and  ${\sum_{j=1}^n a_{ij}<1}$.

Following Assumption <ref>, we obtain from (<ref>) and (<ref>):

\begin{eqnarray}
\sum_{i=1}^n \dot{\mb{y}}_i
% 	&=& -\sum_{i,j=1}^{n} a_{ij}(\mb{y}_i-\mb{y}_j) + \sum_{i=1}^n\frac{d}{d t} \boldsymbol{ \nabla} f_i(\mb{x}_i)\\
&=& \sum_{i=1}^n\frac{d}{d t} \boldsymbol{ \nabla} f_i(\mb{x}_i),  \label{eq_sumydot} \\ \label{eq_sumxdot}
\sum_{i=1}^n \dot{\mb{x}}_i
% 	&=& -\sum_{i,j=1}^{n} w_{ij}(\mb{x}_i-\mb{x}_j)-\alpha \sum_{i=1}^n\mb{y}_i \\
&=& -\alpha \sum_{i=1}^n\mb{y}_i.
\end{eqnarray}

Integrating (<ref>) with respect to $t$ and initializing the auxiliary variable $\mb{y}(0)=\mb{0}_{nm}$, we have

\begin{eqnarray} \label{eq_sumxdot2}
\sum_{i=1}^n \dot{\mb{x}}_i = -\alpha \sum_{i=1}^n\mb{y}_i = -\alpha \sum_{i=1}^n \boldsymbol{ \nabla} f_i(\mb{x}_i),
\end{eqnarray}

which shows that the time-derivative of the sum of states $\mb{x}_i$'s is towards sum gradient, and therefore, the equilibrium ($\dot {\mathbf x_i}=0_m$) of the dynamics (<ref>)-(<ref>) is $\mb{x}^*$ satisfying ${(\mathbf 1_n^\top \otimes I_m) \boldsymbol{ \nabla} F(\mb{x}^*) = \mb{0}_m}$ ($I_m$ as the identity matrix of size $m$), which is the optimal state of  problem (<ref>) [10].

Initializing from any $\mb{x}(0) \neq \mb{1}_n \otimes  {\mb{x}}_0$, for some non-zero $\mathbf{x}_0\in\mathbb{R}^m$, and $\mb{y}(0)=\mb{0}_{nm}$, the  state $[\mb{x}^*;\mb{0}_{nm}]$ with $(\mathbf 1_n^\top \otimes I_m) \boldsymbol{ \nabla} F(\mb{x}^*) = \mb{0}_m$ is an invariant equilibrium point of the dynamics (<ref>)-(<ref>).

From (<ref>),

the following uniquely holds at ${\mb{x}=\mb{x}^*=\mb{1}_n \otimes \overline{ \mb{x}}^*}$,

\[\sum_{i=1}^n \dot{\mb{x}}_i = -\alpha (\mathbf 1_n^\top \otimes I_m) \boldsymbol{ \nabla} F(\mb{x}^*) =  \mb{0}_m.
\]

Further, from (<ref>) we have $\dot{\mb{x}}_i = \mb{0}_m$
and from (<ref>) and (<ref>),

\[\dot{\mb{y}}_i = \frac{d}{d t} \boldsymbol{ \nabla} f_i(\overline{ \mb{x}}^*)=  \boldsymbol{ \nabla}^2 f_i(\overline{ \mb{x}}^*) \dot{\mb{x}}_i = \mb{0}_m,
\]

which shows that $[\mb{x}^*;\mb{0}_{nm}]$ is an invariant equilibrium point of the dynamics (<ref>)-(<ref>).

The above lemma only shows that the state $[\mb{x}^*;\mb{0}_{nm}]$, with $\mb{x}^*$ as the optimal point of problem (<ref>), is the equilibrium of the proposed networked dynamics (<ref>)-(<ref>).

Note that the first term in Eq. (<ref>) drives the agents to reach consensus on $\mb{x}_i$'s while the second term along with the dynamics (<ref>) implements the gradient correction.

§ PROOF OF CONVERGENCE

In this section, we show that  dynamics (<ref>)-(<ref>) converge to the equilibrium state described in Lemma <ref>.
As it is the case in Section <ref>, in which the Laplacian matrix for $W$ is $\overline{W}$, the Laplacian matrix for $A$ is denoted by $\overline{A}$.
Define the $nm$-by-$nm$  Hessian matrix $H\coloneqq\text{blockdiag}[\boldsymbol{ \nabla}^2 f_i(\mb{x}_i)]$.
The dynamics (<ref>)-(<ref>) can be written in compact form as

\begin{eqnarray} \label{eq_xydot}
\left(\begin{array}{c} \dot{\mb{x}} \\ \dot{\mb{y}} \end{array} \right) = M (t,\alpha) \left(\begin{array}{c} {\mb{x}} \\ {\mb{y}} \end{array} \right),
\end{eqnarray}

\begin{eqnarray} \label{eq_M}
M(t,\alpha) = \left(\begin{array}{cc} \overline{W} \otimes I_m & -\alpha I_{mn} \\ H(\overline{W}\otimes I_m) & \overline{A} \otimes I_m - \alpha H
\end{array} \right).
\end{eqnarray}

Networked dynamics (<ref>)-(<ref>) represents a hybrid dynamical system because: (i) the matrix $H$ varies in CT; and (ii) the structure of $\overline{W}$ and $\overline{A}$ may change in DT in case of dynamic network topology. In this direction, towards convergence analysis, we first evaluate the stability properties of the matrix $M$ at every time-instant, and then generalize the convergence to the entire time horizon. In the rest of this paper for notation simplicity, we drop the dependence of $M$ on $(t,\alpha)$, unless where needed, despite the fact that it is a function of both time $t$ and stepsize $\alpha$.

[26, 27]
Let $P(\alpha)$ be an $n$-by-$n$ matrix depending smoothly on a real parameter $\alpha \geq 0$. Assume $P(0)$ has $l<n$ equal eigenvalues, denoted by $\lambda_1=\ldots=\lambda_l$, associated with right eigenvectors $\mb{v}_1,\ldots,\mb{v}_l$ and left eigenvectors $\mb{u}_1,\ldots,\mb{u}_l$ which are linearly independent.
Let $\lambda_i(\alpha)$ denote the eigenvalues of $P(\alpha)$, as a function of $\alpha$, corresponding to ${\lambda_i, i \in \{1,\ldots,l\}}$, and $P' = \frac{dP(\alpha)}{d\alpha}|_{\alpha=0}$. Then, $\frac{d\lambda_i}{d\alpha}|_{\alpha=0}$ are the eigenvalues of the following $l$-by-$l$ matrix,

\[\left(\begin{array}{ccc}
\mb{u}_1^\top P' \mb{v}_1 & \ldots & \mb{u}_1^\top P' \mb{v}_l \\
& \ddots & \\
\mb{u}_l^\top P' \mb{v}_1 & \ldots & \mb{u}_l^\top P' \mb{v}_l
\end{array} \right).
\]

Let Assumption <ref> hold. For a sufficiently small $\alpha$, all eigenvalues of $M$ have non-positive real-parts and the algebraic multiplicity of the zero eigenvalue is $m$.

Let $M=M_0+\alpha M_1$ with

\begin{eqnarray}\nonumber
M_0 &=&  \left(\begin{array}{cc} \overline{W} \otimes I_m & \mb{0}_{mn\times mn} \\ H(\overline{W} \otimes I_m) & \overline{A}\otimes I_m \end{array} \right),\\\nonumber
M_1 &=& \left(\begin{array}{cc} \mb{0}_{mn\times mn} & - {I_{mn}} \\ {\mb{0}_{mn\times mn}} & - H \end{array} \right),
\end{eqnarray}

where $\mb{0}_{mn\times mn}$ is the zero matrix of size $mn$.
Since matrix $M_0$ is block (lower) triangular we have,

\begin{eqnarray}
\sigma(M_0) = \sigma(\overline{W} \otimes I_m) \cup \sigma(\overline{A} \otimes I_m),
\end{eqnarray}

where $\sigma(\cdot)$ represents the eigenspectrum of the matrix. From Lemma <ref>, both matrices  $\overline{W}$ and $\overline{A}$ have $n-1$ eigenvalues in the LHP (left-half plane) and one isolated eigenvalue at zero. Therefore,  matrix $M_0$ has $m$ sets of eigenvalues associated with $m$ dimensions of vector states $\mb{x}_i$ i.e.,

$$\operatorname{Re}\{\lambda_{2n,j}\} \leq \ldots \leq \operatorname{Re}\{\lambda_{3,j}\} < \lambda_{2,j} = \lambda_{1,j} = 0,$$
where $j=\{1,\ldots,m\}.$
Using Lemma <ref>, we analyze  the spectrum of $M$ by considering it as the perturbed version of $M_0$ via the term $\alpha M_1$. We check the variation of the zero eigenvalues $\lambda_{1,j}$ and $\lambda_{2,j}$ by adding the (small) perturbation $\alpha M_1$. Denote these perturbed eigenvalues by $\lambda_{1,j}(\alpha)$ and $\lambda_{2,j}(\alpha)$. To apply Lemma <ref>, define the right eigenvectors corresponding to $\lambda_{1,j}$ and $\lambda_{2,j}$ as,

\begin{equation} \label{eq_V}
V = [V_1~ V_2] =\left(\begin{array}{cc}
\mb{1}_n& \mb{0}_n \\
\mb{0}_n & \mb{1}_n
\end{array} \right)\otimes I_m,
\end{equation}

Similarly, the left eigenvectors are $V^\top$. Note that these eigenvectors are defined using Lemma <ref> and satisfy $V^\top V=I_{2mn}$. Recall that,  $\frac{dM(\alpha)}{d\alpha}|_{\alpha=0}=M_1$ and following Lemma <ref>,

\begin{eqnarray} \label{eq_dmalpha}
V^\top M_1 V= \left(\begin{array}{cc}
\mb{0}_{m\times m}	& \mb{0}_{m\times m} \\
-nI_m	& -(\mb{1}_n \otimes I_m)^\top H (\mb{1}_n \otimes I_m)
\end{array} \right).
\end{eqnarray}

Following the definition of the Hessian matrix $H$,

\begin{equation} \label{eq_sum_df}
-(\mb{1}_n \otimes I_m)^\top H (\mb{1}_n \otimes I_m)= -\sum_{i=1}^n \boldsymbol{ \nabla}^2 f_i(\mb{x}_i) \prec 0,
\end{equation}

where the last inequality follows the strict convexity of the local functions $f_i(\mb{x}_i)$ (Lemma <ref>).

Recall that from Lemma <ref> the derivatives $\frac{d\lambda_{1,j}}{d\alpha}|_{\alpha=0}$ and $\frac{d\lambda_{2,j}}{d\alpha}|_{\alpha=0}$ depend on the eigenvalues of (<ref>) which clearly form a lower triangular matrix with $m$ zero eigenvalues and $m$ negative eigenvalues (following (<ref>)). Therefore, $\frac{d\lambda_{1,j}}{d\alpha}|_{\alpha=0} = 0$ and $\frac{d\lambda_{2,j}}{d\alpha}|_{\alpha=0}<0$,

which implies that considering $\alpha M_1$ as a perturbation, the $m$ zero eigenvalues $\lambda_{2,j}(\alpha)$ of $M$ move toward the LHP while $\lambda_{1,j}(\alpha)$'s remain zero. We recall that the eigenvalues are a continuous functions of the matrix elements [22], and therefore, for sufficiently small $\alpha$ we have,

\begin{equation}
\begin{aligned}
\operatorname{Re}\{\lambda_{2n,j}(\alpha)\} &\leq \ldots \leq \operatorname{Re}\{\lambda_{3,j}(\alpha)\} \\
&\leq \lambda_{2,j}(\alpha) < \lambda_{1,j}(\alpha) = 0.
\end{aligned}
\end{equation}

The proof is completed.

Theorem <ref> proves that for sufficiently small $\alpha$ the matrix $M$ defined by (<ref>) has $m$ zero eigenvalue each associated with a dimension of variable $\mb{x}_i$,  while all other eigenvalues remain in the LHP. Next, we determine an  upper-bound on $\alpha$ guaranteeing  the results of Theorem <ref>. In this direction,  we provide some relevant concepts regarding $\sigma(M_0)$ and $\sigma(M)$ as the spectrum of $M_0$ and $M$ matrices. Define the optimal matching distance as [28],
\[d(\sigma(M),\sigma(M_0)) = \min_{\pi} \max_{1\leq i\leq 2nm} (\lambda_i - \lambda_{\pi(i)}(\alpha)),
\]

with $\pi$ representing all possible permutations over the set $\{1,\ldots,2nm\}$. In fact, $d(\sigma(M),\sigma(M_0))$ is the smallest-radius circle among the circles centered at $\lambda_{1,j},\ldots,\lambda_{2n,j}$ (eigenvalues of $M_0$) which includes all the eigenvalues of $M$ denoted by $\lambda_{1,j}(\alpha),\ldots,\lambda_{2n,j}(\alpha)$. Loosely speaking, $d(\sigma(M),\sigma(M_0))$ represents the furthest distance between the eigenvalues of matrices $M$ and $M_0$. Recall that from Theorem <ref> we know that the first $2m$ eigenvalues of the perturbed matrix $M$ are $\lambda_{1,j}(\alpha)=0$ and $\lambda_{2,j}(\alpha)<0$. To show that all the other $(2n-2)m$ eigenvalues $\lambda_{3,j}(\alpha),\ldots,\lambda_{2n,j}(\alpha)$ remain in the LHP, it is sufficient  that,

\begin{eqnarray} \label{eq_lambda3}
\end{eqnarray}

where $\lambda_{\min}=\min_{1\leq j\leq m}|\operatorname{Re}\{\lambda_{3,j}\}|$. This guarantees that the distance between the $(2n-2)m$ eigenvalues of $M_0$ and $M$ is less than $\lambda_{\min}$ and therefore all the $(2n-2)m$ eigenvalues of $M$ remain in the LHP. In this direction, the following lemma provides a useful bound on $d(\sigma(M),\sigma(M_0))$.
[28]
The following holds for ${M = M_0 +\alpha M_1}$:

\begin{eqnarray*}
d(\sigma(M),\sigma(M_0))\leq 4(\lVert M_0\rVert_{\infty}+\lVert M\rVert_{\infty})^{1-\frac{1}{nm}} \lVert \alpha M_1\rVert_\infty^{\frac{1}{nm}}.
\end{eqnarray*}

From (<ref>) and Lemma <ref>, $\alpha$ can be bounded as follows.

Define ${\gamma = \max_{1\leq i\leq nm} \sum_{j=1}^{nm} |H_{ij}|}$ and ${\lambda_{\min}=\min_{1\leq j\leq m}|\operatorname{Re}\{\lambda_{3,j}\}|}$. The real-part of the eigenvalues ${\operatorname{Re}\{\lambda_{3,j}(\alpha)\},\ldots,\operatorname{Re}\{\lambda_{2n,j}(\alpha)\}<0}$ if $0<\alpha<\overline{\alpha}$ where for $\gamma<1$,

\begin{equation} \label{eq_alphabar1}
\overline{\alpha} = \argmin_{\alpha>0} |4(\max\{4+4\gamma+\alpha\gamma, 4+2\gamma+\alpha\})^{1-\frac{1}{nm}} \alpha^{\frac{1}{nm}}-\lambda_{\min}|,
\end{equation}

and for $\gamma\geq 1$,

\begin{equation} \label{eq_alphabar2}
\overline{\alpha} = \argmin_{\alpha>0} |4(4+4\gamma+\alpha\gamma)^{1-\frac{1}{nm}} (\alpha \gamma)^{\frac{1}{nm}}-\lambda_{\min}|.
\end{equation}

From Assumption <ref> and  Lemmas <ref>-<ref>,
$$||M_0||_\infty\leq 2(1+\gamma).$$

This is because, from Assumption <ref> the row sum of the absolute values of matrix $\overline{W}$ is at most $2$ and similarly, using Lemma <ref>, the row sum of the absolute values of matrix $H(\overline{W} \otimes I_m)$ is at most $2\gamma$. By similar reasoning,

\[||M||_\infty\leq \max\{2+\gamma(2+\alpha), 2+\alpha\},
\]

\[||\alpha M_1||_\infty\leq \max\{\alpha \gamma, \alpha\}.\]

Following (<ref>), for $\gamma<1$,

\[ 4(2(1+\gamma)+\max\{2+\gamma(2+\alpha), 2+\alpha\})^{1-\frac{1}{nm}} \alpha^{\frac{1}{nm}}<\lambda_{\min},
\]

and for  $\gamma\geq1$,

\[ 4(4+\gamma(4+\alpha))^{1-\frac{1}{nm}} (\alpha \gamma)^{\frac{1}{nm}}<\lambda_{\min}.
\]

Since the functions on the left-hand-side of the above inequalities are monotonically increasing for $\alpha>0$, the largest $\alpha$ satisfying the above inequalities is given  by (<ref>)-(<ref>).

Lemma <ref> gives a conservative upper-bound on $\alpha$ which guarantees the rest of the eigenvalues, other than $\lambda_{1,j}(\alpha)=0$ and $\lambda_{2,j}(\alpha)<0$, remain in the LHP and Theorem <ref> is valid.	However, the eigenvalues of $M$ may still remain in the LHP for a possible less-conservative choice of $\alpha>\overline{\alpha}$. In general, for a proper $\alpha$, matrix $M$ has $m$ zero eigenvalues associated with the eigenvectors $V_1$ given in (<ref>), and the null space of the time-varying matrix $M$ is,

\begin{equation}
\mc{N}(M) = \text{span}\{\left(\begin{array}{c}
\mb{1}_n \\
\mb{0}_n
\end{array} \right)\otimes I_m\},
\end{equation}

which is independent of time.

Let the conditions in Lemma <ref>, Lemma <ref>, and Theorem <ref> hold.  The proposed dynamics (<ref>)-(<ref>) converges to $[\mb{x}^*;\mb{0}_{nm}]$ with $\mb{x}^*$ as the optimal value of problem (<ref>).

Define the following proper positive-definite Lyapunov function proposed in [29],

\begin{eqnarray}
V(\boldsymbol{\delta}) = \frac{1}{2} \boldsymbol{\delta}^\top \boldsymbol{\delta} =  \frac{1}{2}\lVert\boldsymbol{\delta}\rVert_2^2
\end{eqnarray}

with $\delta \in \mathbb{R}^{2mn}$ defined as the difference of system state and the optimal state,

\begin{eqnarray}
\boldsymbol{\delta} = \left(\begin{array}{c} {\mb{x}} \\ {\mb{y}} \end{array} \right)-\left(\begin{array}{c} {\mb{x}^*} \\ {\mb{0}_{mn}} \end{array} \right).
\end{eqnarray}

Since from Lemma <ref> $[\mb{x}^*;\mb{0}_{nm}]$ is an invariant state of the networked dynamics (<ref>)-(<ref>), we have
$\dot{\boldsymbol{\delta}} = M \boldsymbol{\delta}$.
Then, the time-derivative of the proposed Lyapunov function is as follows,

\[\dot{V} = \boldsymbol{\delta}^\top \dot{\boldsymbol{\delta}}=  \boldsymbol{\delta}^\top M \boldsymbol{\delta},
%\boldsymbol{\delta}^\top M^\top \boldsymbol{\delta}=
%\boldsymbol{\delta}^\top \frac{M+M^\top}{2} \boldsymbol{\delta}
\]

Following Theorem <ref>, let ${\lambda}_{1,j}(\alpha)=0$, $Re\{{\lambda}_{i,j}(\alpha)\}<0$ for $2\leq i\leq2n,1\leq j\leq m$ represent the real-parts of the eigenspectrum of $M$.

It is known that [23],

\begin{eqnarray} \label{eq_Re2}
\boldsymbol{\delta}^\top M \boldsymbol{\delta} \leq \max_{1\leq j\leq m}\operatorname{Re}\{{\lambda}_{2,j}(\alpha)\} \boldsymbol{\delta}^\top  \boldsymbol{\delta}.
\end{eqnarray}

Since $M$ varies in time, $\max_{1\leq j\leq m}\operatorname{Re}\{{\lambda}_{2,j}(\alpha)\}$ also changes in time. However, from Theorem <ref>, it is always negative, implying that $\dot{V}< 0$ for $\delta \neq \mb{0}_{2mn}$.

We have,

\[\dot{V} =  0 \Leftrightarrow \boldsymbol{\delta} = \mb{0}_{2mn},
\]

and, from LaSalle’s invariance principle, convergence to the invariant set $\{\boldsymbol{\delta} = \mb{0}_{2mn}\}$ follows (see [29] Section 4.1).

Following (<ref>), the convergence rate of the dynamics (<ref>)-(<ref>) depends on $\operatorname{Re}\{{\lambda}_{2,j}(\alpha)\}$. Note that $\operatorname{Re}\{{\lambda}_{2,j}(\alpha)\}$ is tightly related to the parameter $\alpha$ and therefore, to improve the convergence rate parameter $\alpha$ needs not to be very small.

§ SIMULATION: NONLINEAR SVM EXAMPLE

For simulation consider the academic example given in [30] (page $747$). Consider $N=60$ uniformly distributed sample data points, shown in Fig. <ref>(Left), represented in two classes:  blue `*'s and  red `o's. Clearly, these points $\boldsymbol{\chi}_i=[{\chi}_i(1);{\chi}_i(2)]$ are not linearly separable in $\mathbb{R}^2$. The nonlinear mapping $\phi(\boldsymbol{\chi}_i) = [{\chi}_i(1)^2;{\chi}_i(2)^2;\sqrt{2}{\chi}_i(1){\chi}_i(2)]$, proposed by [30], properly maps the data points to $\mathbb{R}^3$ such that the projected data points are linearly separable via a hyperplane as shown in Fig. <ref>(Right). It can be shown that  the associated kernel function is
$K(\boldsymbol{\chi}_i,\boldsymbol{\chi}_j)=(\phi(\boldsymbol{\chi}_i)^\top \phi(\boldsymbol{\chi}_j))^2$.

(Left) Training data and the optimal nonlinear classifier (the ellipse) in 2D. (Right) The same points mapped into 3D space via a nonlinear mapping. Linear SVM optimally classifies the data points via the gray hyperplane which represents the  ellipse in the left figure by inverse mapping.

We evaluate the proposed dynamics (<ref>)-(<ref>) (with $\alpha=10$) for D-SVM over a network of $n=5$ agents each having access to $50\%$ random selection of the data points.
The loss function $f_i(\mb{x}_i)$ follows the smooth approximation discussed in Section <ref> with $\mu=3$ and $C=1.5$. Every agent finds the optimal separating hyperplane defined by $\mb{x}_i = [\boldsymbol{\omega}_i^\top;\nu_i]$ ($\boldsymbol{\omega}_i \in \mathbb{R}^3$) and shares this value along with the auxiliary variable $\mb{y}_i$ with its direct neighbors in $\mc{G}$. The agents' network $\mc{G}$ is considered as the union of a directed cycle and a $2$-hop digraph
(see examples in [2]). To satisfy the weight-balanced condition the link weights in each network are equal and randomly chosen in the range $(0,0.5)$. Using MATLAB's  function, we randomly change the permutation of the nodes in the network and the link weights every $0.05$ seconds to simulate  a dynamic network in DT domain. The time-evolution of $\mb{x}_i  = [\boldsymbol{\omega}_i^\top;\nu_i] \in \mathbb{R}^4$, loss function $F(\mb{x})$, and sum of the gradients $\sum_{i=1}^{5} \boldsymbol{ \nabla} f_i(\mb{x}_i) \in \mathbb{R}^4$ are shown in  Fig. <ref>.

The time-evolution of the SVM classifier parameters  ${\omega}_i$ and $\nu_i$ (at all $5$ agents) under dynamics (<ref>)-(<ref>)  along with overall loss function $F(\mb{x})$ and sum of the gradients $\sum_{i=1}^{5} \boldsymbol{ \nabla} f_i(\mb{x}_i)$. The optimal values based on the centralized SVM are also shown for comparison.

The agents reach consensus on the optimal value ${\overline{\mb{x}}^*=[\overline{\omega}(1),\overline{\omega}(2),\overline{\omega}(3),\overline{\nu}]^\top}$ as the parameters of the separating hyperplane in $\mathbb{R}^3$. Via the inverse mapping, the hyperplane in $\mathbb{R}^3$ represents an ellipse formulated as ${\overline{\omega}(1)z_1^2 + \overline{\omega}(2)z_2^2 -\overline{\nu}=0}$ ($z_1$ and $z_2$ as the Cartesian coordinates in $\mathbb{R}^2$), which separates the original data points $\boldsymbol{\chi}_i$'s in $\mathbb{R}^2$. The calculated separating ellipses by  all $5$ agents are shown in Fig. <ref> at two different time-instants.

Optimal elliptical SVM classifiers computed by all the agents under the distributed optimization  dynamics (<ref>)-(<ref>)  at time $t=0.1$ (Left) and ${t=2}$ (Right).

§ CONCLUSION AND FUTURE RESEARCH

In this work, a Lipschitz dynamics is proposed to solve D-SVM over a dynamic WB-digraph in a hybrid setting. We adopt matrix perturbation analysis to prove convergence of the CT dynamics (<ref>)-(<ref>) whose parameters vary due to  switching network topology in DT domain. In particular, our proposed distributed optimization in D-SVM setup enables the agents to cooperatively learn the classifier over a dynamic network via local information, improving classical D-SVM methods in terms of data privacy [7, 6, 8] and computational complexity [9].

As future research direction, one can extend the results to the DT counterpart.

For example, for Euler-Forward method, the DT version of matrix $M$ in (<ref>) is $M_d = I+TM$ with $T$ as the sampling period. Then, the explicit upper bound on $T$ such that a stable CT dynamics from Theorem <ref>-<ref> remains stable after discretization can be defined. Additionally, extensions to time-delayed networks, online D-SVM, and sparse digraphs are directions of interest.

[1]
R. Xin, S. Kar, and U. A. Khan,
“Decentralized stochastic optimization and machine learning: A
unified variance-reduction framework for robust performance and fast
IEEE Signal Processing Magazine, vol. 37, no. 3, pp. 102–113,

[2]
R. Xin, S. Pu, A. Nedić, and U. A. Khan,
“A general framework for decentralized optimization with first-order
Proceedings of the IEEE, vol. 108, no. 11, pp. 1869–1889,

[3]
A. Simonetto, E. Dall'Anese, S. Paternain, G. Leus, and G. B. Giannakis,
“Time-varying convex optimization: Time-structured algorithms and
Proceedings of the IEEE, vol. 108, no. 11, pp. 2032–2048,

[4]
M. Doostmohammadian, A. Aghasi, and T. Charalambous,
“Fast-convergent dynamics for distributed resource allocation over
sparse time-varying networks,”
arXiv preprint arXiv:2012.08181, 2020.

[5]
Y. Lu, V. Roychowdhury, and L. Vandenberghe,
“Distributed parallel support vector machines in strongly connected
IEEE Transactions on Neural Networks, vol. 19, no. 7, pp.
1167–1178, 2008.

[6]
E. Y. Chang, K. Zhu, H. Wang, and H. Bai,
“Psvm: Parallelizing support vector machines on distributed
in Foundations of Large-Scale Multimedia Information Management
and Retrieval, pp. 213–230. Springer, 2011.

[7]
A. Navia-Vázquez, D. Gutierrez-Gonzalez, E. Parrado-Hernández, and
J. J. Navarro-Abellan,
“Distributed support vector machines,”
IEEE Transactions on Neural Networks, vol. 17, no. 4, pp. 1091,

[8]
A. Bordes, S. Ertekin, J. Weston, L. Botton, and N. Cristianini,
“Fast kernel classifiers with online and active learning,”
Journal of Machine Learning Research, vol. 6, no. 9, 2005.

[9]
P. A. Forero, A. Cano, and G. B. Giannakis,
“Consensus-based distributed support vector machines.,”
Journal of Machine Learning Research, vol. 11, no. 5, 2010.

[10]
B. Gharesifard and J. Cortés,
“Distributed continuous-time convex optimization on weight-balanced
IEEE Transactions on Automatic Control, vol. 59, no. 3, pp.
781–786, 2014.

[11]
B. Ning, Q. Han, and Z. Zuo,
“Distributed optimization for multiagent systems: An edge-based
fixed-time consensus approach,”
IEEE Transactions on Cybernetics, vol. 49, no. 1, pp. 122–132,

[12]
K. Garg, M. Baranwal, A. O. Hero, and D. Panagou,
“Fixed-time distributed optimization under time-varying
communication topology,”
arXiv preprint arXiv:1905.10472, 2019.

[13]
S. Rahili and W. Ren,
“Distributed continuous-time convex optimization with time-varying
cost functions,”
IEEE Transactions on Automatic Control, vol. 62, no. 4, pp.
1590–1605, 2017.

[14]
M. Doostmohammadian,
“Single-bit consensus with finite-time convergence: Theory and
IEEE Transactions on Aerospace and Electronic Systems, vol. 56,
no. 4, pp. 3332–3338, 2020.

[15]
Z. Li and Z. Ding,
“Time-varying multi-objective optimisation over switching graphs via
fixed-time consensus algorithms,”
International Journal of Systems Science, vol. 51, no. 15, pp.
2793–2806, 2020.

[16]
P. Armand and R. Omheni,
“A globally and quadratically convergent primal–dual augmented
lagrangian algorithm for equality constrained optimization,”
Optimization Methods and Software, vol. 32, no. 1, pp. 1–21,

[17]
P. Srivastava and J. Cortés,
“Distributed algorithm via continuously differentiable exact penalty
method for network optimization,”
in IEEE Conference on Decision and Control (CDC). IEEE, 2018,
pp. 975–980.

[18]
F. Mansoori and E. Wei,
“A fast distributed asynchronous newton-based optimization
IEEE Transactions on Automatic Control, vol. 65, no. 7, pp.
2769–2784, 2019.

[19]
K. Yuan, Q. Ling, and W. Yin,
“On the convergence of decentralized gradient descent,”
SIAM Journal on Optimization, vol. 26, no. 3, pp. 1835–1854,

[20]
R. Goebel, R. G. Sanfelice, and A. R. Teel,
“Hybrid dynamical systems,”
IEEE control systems magazine, vol. 29, no. 2, pp. 28–93,

[21]
D. L. Ly and H. Lipson,
“Learning symbolic representations of hybrid dynamical systems,”
The Journal of Machine Learning Research, vol. 13, no. 1, pp.
3585–3618, 2012.

[22]
G. W. Stewart and J. Sun,
“Matrix perturbation theory,”

[23]
R. Olfati-Saber and R. M. Murray,
“Consensus problems in networks of agents with switching topology
and time-delays,”
IEEE Transactions on Automatic Control, vol. 49, no. 9, pp.
1520–1533, Sept. 2004.

[24]
O. Chapelle,
“Training a support vector machine in the primal,”
Neural computation, vol. 19, no. 5, pp. 1155–1178, 2007.

[25]
D. Jurafsky and J. H. Martin,
Speech and Language Processing,
Prentice Hall, 2020.

[26]
A. P. Seyranian and A. A. Mailybaev,
Multiparameter stability theory with mechanical applications,
vol. 13,
World Scientific, 2003.

[27]
K. Cai and H. Ishii,
“Average consensus on general strongly connected digraphs,”
Automatica, vol. 48, no. 11, pp. 2750–2761, 2012.

[28]
R. Bhatia,
Matrix analysis,
Springer Science & Business Media, 2013.

[29]
M. Mesbahi and M. Egerstedt,
Graph theoretic methods in multiagent networks, vol. 33,
Princeton University Press, 2010.

[30]
S. Russell and P. Norvig,
“Artificial intelligence: a modern approach,”