A Stochastic Second-Order Proximal Method
for Distributed Optimization

Stochastic optimization algorithms have been flourishing recently due to their appealing efficiency in machine learning [1, 2]. In the context of large-scale machine learning, parallel stochastic algorithms are often used to process large datasets[3, 4]. However, such methods need a central node to ensure the consistency of the variables from all the nodes, so that the communication burden of the central node becomes the bottleneck that restricts the algorithm performance.

On the other hand, a collection of distributed optimization algorithms have been proposed over the past decade in order to tackle various network control and resource allocation problems, where agents in a network only communicate with their neighbors and do not rely on any central coordination, eliminating potential communication bottlenecks in the computing infrastructure [5]. Typical methods include the distributed gradient descent (DGD) method [5], the decentralized exact first-order algorithm (EXTRA) [6], and the distributed gradient tracking algorithms [7].

Inheriting the merits of the above two algorithm types, distributed stochastic optimization algorithms have been attracting a lot of recent interest. For smooth, strongly convex optimization problems, [8] develops a distributed stochastic gradient descent (DSGD) method based on DGD, which is shown to achieve the optimal sublinear convergence rate (independent of the network) of a centralized stochastic gradient descent (SGD) method [9, 10]. The same convergence rate is attained by the exact diffusion method with adaptive step-sizes (EDAS) in [11]. In addition, the distributed stochastic gradient tracking (DSGT) method proposed in [12] is guaranteed to linearly converge to a neighborhood of the optimal solution in expectation. For non-convex problems, [13] designs a distributed primal-dual SGD algorithm with both fixed step-sizes (DPD-SGD-F) and adaptive step-sizes (DPD-SGD-T), where the former achieves linear convergence to suboptimality under the Polyak-Łojasiewics (PL) condition.

The aforementioned distributed stochastic algorithms are all constructed upon deterministic first-order methods by nature and only use stochastic gradients to evolve. As second-order information often leads to more accurate approximation and accelerates problem solving, we endeavor to develop a second-order distributed stochastic optimization algorithm.

To this end, we choose SoPro [14], a deterministic distributed second-order proximal algorithm, as the cornerstone. SoPro is developed by virtue of a decentralized second-order approximation of the augmented Lagrangian function in the classic method of multipliers [15], and its convergence performance outperforms that of various distributed first-order methods in the deterministic setting. In this paper, we adapt SoPro to the stochastic scenario. Specifically, instead of letting each agent compute exact local gradient and Hessian matrix determined by all its local data as SoPro does, we allow each agent to update using stochastic approximations of its local gradient and Hessian, which come from two randomly and uniformly chosen batches of samples from its local loss function. Such a stochastic variant of SoPro can significantly enhance the computational and memory efficiency of the agents. We call this algorithm a stochastic second-order proximal algorithm, referred to as St-SoPro.

Under the assumptions that the local loss functions are smooth and convex, and that their sum (i.e., the global objective function) is globally restricted strongly convex, we show that our proposed St-SoPro algorithm linearly converges to a neighborhood of the optimal solution in expectation over undirected networks. In particular, we provide an explicit upper bound on its ultimate suboptimality, and illustrate that this upper bound can be made arbitrarily small as long as the parameters are properly set. Finally, we validate the superior performance of St-SoPro in comparison with several recent distributed stochastic optimization methods over some real datasets for classification problems arising in machine learning with respect to convergence speed, communication load, computational efficiency, and classification accuracy.

The paper is organized as follows. Section II formulates the optimization problem to be solved, and Section III describes the proposed St-SoPro algorithm. Section IV provides the convergence analysis, Section V presents the numerical results, and Section VI concludes the paper.

Notation: For any differentiable function $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$, its gradient at $x\in\mathbb{R}^{d}$ is denoted by $\nabla f(x)$, and if $f$ is twice-differentiable, we use $\nabla^{2}f(x)$ to denote its Hessian matrix at $x$. For any set $S$, $|S|$ represents the cardinality of $S$. In addition, $\otimes$ is the Kronecker product, $\langle\cdot,\cdot\rangle$ is the Euclidean inner product, and $\|\cdot\|$ is the $\ell_{2}$ norm. We use $\mathbf{0}_{d},\mathbf{1}_{d},\mathbf{O}_{d},I_{d}$ to denote the $d$-dimensional all-zero vector, all-one vector, zero matrix, and identity matrix, respectively. Also, $\operatorname{diag}(A_{1},\ldots,A_{n})$ represents the block diagonal matrix whose diagonal blocks are sequentially $A_{1},\ldots,A_{n}$. Given a matrix $A\in\mathbb{R}^{d\times d}$, we write $A\succeq\mathbf{O}_{d}$ if it is positive semidefinite and $A\succ\mathbf{O}_{d}$ if it is positive definite. For any $A\succeq\mathbf{O}_{d}$ and $\mathbf{x}\in\mathbb{R}^{d},\|\mathbf{x}\|_{A}^{2}:=\mathbf{x}^{T}A\mathbf{x}$, $\lambda_{max}(A)$ and $\lambda_{min}(A)$ are the largest and smallest real eigenvalues of $A$, respectively, and $A^{\dagger}$ is $A$ ’s pseudoinverse.

Consider a set $\mathcal{V}=\{1,\ldots,N\}$ of agents, where the agents are connected through the link set $\mathcal{E}\subseteq$ $\{\{i,j\}\subseteq\mathcal{V}\times\mathcal{V}\mid i\neq j\}$. We model such a network as a connected undirected graph $(\mathcal{V},\mathcal{E})$, and denote the set of each agent $i$’s neighbors by $\mathcal{N}_{i}=\{j\in\mathcal{V}\mid\{i,j\}\in\mathcal{E}\}$. Suppose each agent $i$ observes a finite number of local samples $\xi_{i,j}\in\mathbb{R}^{m}$, $j=1,\ldots,\mathcal{C}_{i}$ that are independent random vectors, and attempts to solve the following optimization problem:

In problem (1), $f_{i}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is the local loss function of agent $i$, which is the average of every sample’s loss $l_{i,j}(x_{i},\xi_{i,j}):\mathbb{R}^{d}\rightarrow\mathbb{R}$ associated with agent $i$, i.e.,

Below, we impose the assumptions on problem (1).

Assumption 1 leads to the following inequalities. For any given $x,y\in\mathbb{R}^{d}$ and $\xi_{i,j}$, $m_{i}\|x-y\|^{2}\leq\langle\nabla l_{i,j}(x,\xi_{i,j})-\nabla l_{i,j}(y,\xi_{i,j}),x-y\rangle\leq M_{i}\|x-y\|^{2}$ and

for some $m_{i}\in[0,M_{i}]$. Also, the globally restricted strong convexity in Assumption 1a) guarantees that the optimal solution $x^{\star}$ is unique.

Problem (1) requires that the agents reach a consensus while minimizing all the sample losses throughout the network. Indeed, a wide range of real-world problems can be cast into the form of problem (1), such as distributed model predictive control [16], distributed spectrum sensing [17], and logistic regression [18]. Under many circumstances, these engineering problems involve huge datasets. Thus, we focus on solving problem (1) in a fully decentralized and stochastic fashion. Specifically, we only allow each agent to communicate with its neighbors and use a randomly chosen subset of its local samples to compute.

In this section, we develop a distributed stochastic algorithm for solving problem (1) over undirected networks.

To this end, we first provide a brief review of the distributed (deterministic) second-order proximal algorithm (SoPro) in [14]. Note that problem (1) is equivalent to

where $\mathbf{x}=(x_{1}^{T},\ldots,x_{N}^{T})^{T}$, $W=P\otimes I_{d}\succeq\mathbf{O}_{Nd}$, and

with $p_{ij}=p_{ji}>0$ $\forall\{i,j\}\in\mathcal{E}$, so that the null space of $P$ is $\operatorname{span}\{\mathbf{1}_{N}\}$. Also, the unique optimal solution of problem (3) is $\mathbf{x}^{\star}=((x^{\star})^{T},\ldots,(x^{\star})^{T})^{T}$.

The application of the method of multipliers [15] to solve (3) gives the following: Starting from any $\mathbf{v}^{0}\in\mathbb{R}^{Nd}$,

where $\mathbf{x}^{k}=((x_{1}^{k})^{T},\ldots,(x_{N}^{k})^{T})^{T}$ and $\mathbf{v}^{k}$ are the primal and dual variables, respectively, and $L_{\beta}(\mathbf{x},\mathbf{v}):\mathbb{R}^{Nd}\times\mathbb{R}^{Nd}\rightarrow\mathbb{R}$ is an augmented Lagrangian function given by $L_{\beta}(\mathbf{x},\mathbf{v})=f(\mathbf{x})+\mathbf{v}^{T}W^{\frac{1}{2}}\mathbf{x}+\frac{\beta}{2}\|W^{\frac{1}{2}}\mathbf{x}\|^{2}$, $\beta>0$.

Since (4)–(5) cannot be implemented in a distributed way, the SoPro algorithm in [14] introduces a decentralized second-order proximal approximation of $L_{\beta}\left(\mathbf{x},\mathbf{v}^{k}\right)$ in (4) and applies a variable change to (5). Specifically, it replaces $L_{\beta}\left(\mathbf{x},\mathbf{v}^{k}\right)$ with its second-order Taylor’s expansion at $\mathbf{x}^{k}$. Then, it replaces the remaining coupling term $\frac{1}{2}(\mathbf{x}-\mathbf{x}^{k})^{T}\rho W(\mathbf{x}-\mathbf{x}^{k})$ in the primal update with $\frac{1}{2}(\mathbf{x}-\mathbf{x}^{k})^{T}D(\mathbf{x}-\mathbf{x}^{k})$, where $D=\operatorname{diag}(D_{1},\ldots,D_{N})$ is a symmetric block diagonal matrix satisfying $\nabla^{2}f_{i}(x)+D_{i}\succ\mathbf{O}_{d}$ $\forall x\in\mathbb{R}^{d}$ $\forall i\in\mathcal{V}$. Furthermore, we define $\mathbf{q}^{k}=((q_{1}^{k})^{T},\ldots,(q_{N}^{k})^{T})^{T}=W^{\frac{1}{2}}\mathbf{v}^{k}$ as a substitute for $\mathbf{v}^{k}$, and $\mathbf{q}^{k}$ can be ensured to identically stay in the range of $W^{\frac{1}{2}}$ by letting $\sum_{i\in\mathcal{V}}q_{i}^{0}=0$. To summarize, SoPro takes the following form: Starting from $\mathbf{q}^{0}$ satisfying $\sum_{i\in\mathcal{V}}q_{i}^{0}=0$,

where $\nabla f(\mathbf{x}^{k})=(\nabla f_{1}(x_{1}^{k})^{T},\ldots,\nabla f_{N}(x_{N}^{k})^{T})^{T}$ and $\nabla^{2}f(\mathbf{x}^{k})=\text{diag}(\nabla^{2}f_{1}(x_{1}^{k}),\ldots,\nabla^{2}f_{N}(x_{N}^{k}))$ satisfying $\nabla^{2}f(\mathbf{x}^{k})+D\succ\mathbf{O}_{Nd}$.

The primal update of SoPro (6) requires that each agent uses up all its local samples. However, the agents may only be able to access or process a portion of their local samples at one time, especially in the big data scenario. Motivated by this, we consider approximating the gradient $\nabla f(\mathbf{x}^{k})$ and the Hessian $\nabla^{2}f(\mathbf{x}^{k})$ in (6) via a stochastic gradient $g(\mathbf{x}^{k})$ and a stochastic Hessian $h(\mathbf{x}^{k})$ given by

Here, for each agent $i\in\mathcal{V}$, $\mathcal{G}_{i}^{k}$ and $\mathcal{S}_{i}^{k}$ are two independent random sample sets uniformly chosen from $\{1,\ldots,\mathcal{C}_{i}$} without replacement, so that $g(\mathbf{x}^{k})$ and $h(\mathbf{x}^{k})$ are unbiased, i.e.,

for all $\mathbf{x}^{k}$. Due to (2), $m_{i}I_{d}\preceq h_{i}(x_{i}^{k})\preceq M_{i}I_{d}$ $\forall i\in\mathcal{V}$ and

where $\Lambda_{m}=\operatorname{diag}\left(m_{1},\ldots,m_{N}\right)\otimes I_{d}\succeq\mathbf{O}_{Nd}$ and $\Lambda_{M}=\operatorname{diag}\left(M_{1},\ldots,M_{N}\right)\otimes I_{d}\succ\mathbf{O}_{Nd}$.

Using the above randomly sampled gradient and Hessian, we obtain the following stochastic variant of SoPro: Starting from any $\mathbf{q}^{0}$ such that $\sum_{i\in\mathcal{V}}q_{i}^{0}=0$,

where each $D_{i}$ satisfies $h_{i}(x)+D_{i}\succ\mathbf{O}_{d}$ $\forall x\in\mathbb{R}^{d}$, i.e., $h(\mathbf{x})+D\succ\mathbf{O}_{Nd}$ $\forall\mathbf{x}\in\mathbb{R}^{Nd}$, so that (10) is well-posed. The above initialization and updates compose our proposed stochastic second-order proximal (St-SoPro) method. The distributed implementation of St-SoPro over the undirected network $(\mathcal{V},\mathcal{E})$ is described in Algorithm 1, in which $y_{i}^{k}$ $\forall i\in\mathcal{V}$ are auxiliary variables for better presentation.

This section provides the convergence analysis of St-SoPro.

We first impose an assumption on the expected deviation of each sample loss $l_{i,j}(x_{i},\xi_{i,j})$ from each entire local loss $f_{i}(x_{i})$ in terms of their gradients.

To simplify the notation, below we let $\mathcal{C}_{i}=\mathcal{C}$ $\forall i\in\mathcal{V}$ and $\mathcal{G}_{i}^{k}$ $\forall i\in\mathcal{V}$ $\forall k\geq 0$ be of the same size $\mathcal{G}$. We also abbreviate $\mathbf{E}_{\mathcal{G}_{i}^{k}}[\cdot]$ and $\mathbf{E}_{\mathcal{S}_{i}^{k}}[\cdot]$ to $\mathbf{E}[\cdot]$. According to [19, Chapter 2],

This is consistent with the fact that when computing the stochastic gradient $g(\mathbf{x}^{k})$, if we reduce the number $\mathcal{G}$ of randomly selected samples, then the discrepancy between $g(\mathbf{x}^{k})$ and $\nabla f(\mathbf{x}^{k})$ becomes larger in expectation.

Next, for the sake of presenting the convergence result, we introduce the following notation and definitions. According to [14], any $\mathbf{v}\in\mathbb{R}^{Nd}$ satisfying $\nabla f\left(\mathbf{x}^{\star}\right)=-W^{\frac{1}{2}}\mathbf{v}$ is a dual optimum of problem (3). Thus, we define

as a particular dual optimum of (3). Also, throughout this section, we let $\mathbf{v}^{k}=(W^{\dagger})^{\frac{1}{2}}\mathbf{q}^{k}$, $\mathbf{z}^{k}=((\mathbf{x}^{k})^{T},(\mathbf{v}^{k})^{T})^{T}$, and $\mathbf{z}^{\star}=((\mathbf{x}^{\star})^{T},(\mathbf{v}^{\star})^{T})^{T}$. Since such $\mathbf{v}^{k}$ satisfies (5),

In addition, we define $f_{\beta}(\mathbf{x})=f(\mathbf{x})+\frac{\beta}{4}\|\mathbf{x}\|_{W}^{2}$. Based on [14, Lemma 1], for any $\mathbf{x}\in\mathbb{R}^{Nd}$,

where $\zeta:(0,\infty)\rightarrow\mathbb{R}$ is given by $\zeta(\gamma)=\min\{\frac{m_{\bar{f}}}{N}-2M\gamma,\beta\frac{\lambda_{W}}{2(1+1/\gamma^{2})}\}$, $m_{\bar{f}}$ is given in Assumption 1a), $M=\max_{i\in\mathcal{V}}M_{i}>0$, and $\lambda_{W}$ is the smallest non-zero eigenvalue of $W$. It can be shown that $\zeta(\gamma)>0$ if and only if $\gamma\in(0,m_{\bar{f}}/(2MN))$, and its maximum value is attained at the unique positive root of the cubic equation $4MN\gamma^{3}+(\beta N\lambda_{W}-2m_{\bar{f}})\gamma^{2}+4MN\gamma-2m_{\bar{f}}=0$. We denote the maximum value of $\zeta(\gamma)$ by $m_{\beta}$, which is the convexity parameter of $f_{\beta}$ (and indeed can be taken as any positive $\zeta(\gamma)$).

Our convergence analysis relies on the following parameter condition. Suppose there exists $\eta_{s}\in(0,1)$ such that

Let $R=\frac{\Lambda_{m}+\Lambda_{M}}{2}+D$ and $Q=\operatorname{diag}(\beta R,I_{Nd})$, which are guaranteed to be positive definite due to (15). Furthermore, it follows from (9) and (15) that the condition $h(\mathbf{x})+D\succ\mathbf{O}_{Nd}$ required in Section III holds.

We provide our main result in the theorem below.

It can be shown that larger $\min_{i\in\mathcal{V}}m_{i}$, smaller $\max_{i\in\mathcal{V}}M_{i}$, or smaller $\lambda_{W}$ leads to faster convergence (i.e., larger $\delta_{s}$) of St-SoPro. Such analysis follows the idea in [14] and is thus omitted here due to space limitation. In addition, it can be seen from (17) that the error bound of $\operatorname{lim\;sup}_{k\rightarrow\infty}\mathbf{E}\left[\|\mathbf{z}^{k}-\mathbf{z}^{\star}\|_{Q}^{2}\right]$ drops with the decrease of $\tau$. Hence, essentially, larger random sample sets for computing the stochastic gradients lead to smaller optimality error.

In fact, the expected distance between $\mathbf{x}^{k}$ and $\mathbf{x}^{\star}$ can eventually reach an arbitrarily small value under proper parameter settings. To see this, for simplicity, let $m_{i}=m>0$ and $M_{i}=M\geq m$ for all $i\in\mathcal{V}$, and pick any $\eta_{s}\in(0,1)$, $c_{0}\in(0,2\eta_{s}m_{\beta})$, and $c_{1}>0$. We choose, for example, $D=\alpha I_{Nd}$ with $\alpha=(\frac{1}{2}+\lambda_{max}(W))\beta+\mu$, $\beta>0$, and $\mu>\tfrac{M-3m}{2}+\tfrac{{M}}{2(1-\eta_{s})}+\tfrac{({M}-{m})^{2}}{4c_{0}}$, so that (15) holds. This also results in $\kappa_{c_{0},\eta_{s}}=\mu-\frac{M-3m}{2}-\frac{{M}}{2(1-\eta_{s})}-\frac{({M}-{m})^{2}}{4c_{0}}>0$. From (17), we obtain $\operatorname{lim\;sup}_{k\rightarrow\infty}\mathbf{E}[\|\mathbf{x}^{k}-\mathbf{x}^{\star}\|^{2}]\leq\frac{\Gamma N\tau\sigma^{2}}{\delta_{s}\beta[(m+M)/2+\beta(1/2+\lambda_{max}(W))+\mu]}$. It can thus be shown that as $\beta\rightarrow\infty$, such an upper bound on $\operatorname{lim\;sup}_{k\rightarrow\infty}\mathbf{E}\left[\|\mathbf{x}^{k}-\mathbf{x}^{\star}\|^{2}\right]$ goes to zero. Since this bound is continuous at $\beta$, given any $\epsilon>0$, the above parameter setting with a sufficiently large $\beta$ guarantees $\operatorname{lim\;sup}_{k\rightarrow\infty}\mathbf{E}\left[\|\mathbf{x}^{k}-\mathbf{x}^{\star}\|^{2}\right]<\epsilon$.

This section compares the practical convergence performance of St-SoPro with several state-of-the-art distributed stochastic optimization algorithms.

In the numerical experiment, we intend to learn linear classifiers by solving $l_{2}$-regularized logistic regression of the following form over a randomly generated, undirected, and connected network:

where $\lambda>0$ is the regularization parameter and $\{a_{i,j},b_{i,j}\}$ are the data samples. Our experiment is conducted on two standard real datasets a4a and mushrooms from the LIBSVM library [20]. Table I lists the values of the problem and network parameters corresponding to these two datasets, including the problem dimension $d$, the number $N$ of agents, the network’s average degree $d_{a}=\sum_{i\in\mathcal{V}}\frac{|\mathcal{N}_{i}|}{N}$, the total number $\mathcal{C}_{i}$ of samples that we assign to each agent $i$, the sizes $|\mathcal{G}_{i}^{k}|$ and $|\mathcal{S}_{i}^{k}|$ of the random sample sets that each agent $i$ chooses per iteration, as well as the regularization parameter $\lambda$.

The simulations include DSGD [8], EDAS [11], DSGT [12], and DPD-SGD-T [13], which are all first-order methods, for comparison with our proposed St-SoPro. We fine-tune all the algorithm parameters so that the algorithms reach a given accuracy ($2\times 10^{-1}$ for a4a and $10^{-1}$ for mushrooms) within fewest possible iterations.

Figures 1(a)–(c) and 2(a)–(c) plot the evolutions of the optimality error $\frac{1}{N}\sum_{i\in\mathcal{V}}\|{x}_{i}^{k}-{x}^{\star}\|^{2}$ generated by the aforementioned algorithms over a4a and mushrooms with respect to the number of iterations, the number of communication bits (set as $32$ times the number of transmitted real scalars according to [21]), and computation time. Observe that St-SoPro converges faster than the other algorithms to reach the given accuracy, validating its computational and communication efficiency. It is worth mentioning that although St-SoPro is a second-order method, its computational cost can be comparable with the first-order methods when addressing such common machine learning problems. Figures 1(d) and 2(d) present the correct rates of classification on the test sets upon completing each iteration of these algorithms as training, whereby St-SoPro outperforms the others in the training effect.

We have developed St-SoPro, a distributed stochastic second-order proximal method, for addressing strongly convex and smooth optimization over undirected networks. Different from the existing first-order distributed stochastic algorithms, St-SoPro incorporates a second-order approximation of an augmented Lagrangian function and randomly samples each local gradient and Hessian. We show that St-SoPro linearly converges to a neighborhood of the optimal solution in expectation, and the neighborhood can be arbitrarily small. Simulations over two real datasets demonstrate that St-SoPro is both computationally and communication-wise efficient.