# Learning Action Duration and Synergy in Task Planning for Human-Robot
Collaboration

Samuele Sandrini1, Marco Faroni1, Nicola Pedrocchi1 This work was partially
supported by ShareWork project (H2020, European Commission – G.A. 820807).1
STIIMA-CNR - Institute of Intelligent Industrial Technologies and Systems,
National Research Council of Italy<EMAIL_ADDRESS>
###### Abstract

A good estimation of the actions’ cost is key in task planning for human-robot
collaboration. The duration of an action depends on agents’ capabilities and
the correlation between actions performed simultaneously by the human and the
robot. This paper proposes an approach to learning actions’ costs and coupling
between actions executed concurrently by humans and robots. We leverage the
information from past executions to learn the average duration of each action
and a synergy coefficient representing the effect of an action performed by
the human on the duration of the action performed by the robot (and vice
versa). We implement the proposed method in a simulated scenario where both
agents can access the same area simultaneously. Safety measures require the
robot to slow down when the human is close, denoting a bad synergy of tasks
operating in the same area. We show that our approach can learn such bad
couplings so that a task planner can leverage this information to find better
plans.

###### Index Terms:

Human-Robot Interaction; Task And Motion Planning; Task planning; Learning for
task planning.

##  I Introduction

Human-Robot Collaboration (HRC) often requires the system to make decisions
based on human users’ observed or predicted behavior. In such a context,
planning and allocating tasks to the human and robot agents are two complex
problems, even for tasks composed of a few activities. Two modelling issues
are essential. First, human behavior is intrinsically unpredictable and
partially uncontrollable [1]. Second, the coupling between the human and the
robotic agents is characterized by a large variability. For example, the
feasibility and the duration of an action may vary because of the interference
between the human and the robot (_e.g._ , safety stops of the robot or even
due to path-switch in motion re-planning [2]).

In recent years, task planning and allocation problems for HRC have been
investigated. Existing works tried to model human preferences explicitly and
ergonomics into planning [3]. For example, [4, 5] proposed a hierarchical
agent-based task planner, where complex tasks can be decomposed into simpler
actions. This approach can improve the collaborative experience by considering
human preferences as social costs [6, 7]. Manufacturing-oriented works focus
on process throughput by minimizing the expected duration [8] or planning
contingent plans to reduce process errors [9, 10]. A promising approach to
deal with uncertainty on the task duration is timeline-based planning [11].
Timeline-based planning explicitly models the duration variability of tasks
and finds plans that are robust with respect to it. Timeline-based planning
was used in HRC with pre-computed robot motions [12], and online motion
planning [13], demonstrating robust plans allow for less frequent re-planning
and shorter average task duration.

Most planning-based methods assume a guess of the action cost (e.g., the
duration) is available from a domain expert. If this information is
unreliable, the task planner reasons on wrong assumptions so that the
resulting plans become sub-optimal or even infeasible (leading to frequent re-
planning when using a motion re-planner such as [2]). Moreover, the duration
guess does not consider possible couplings between the tasks executed
concurrently by the human and the robot. Therefore, the task planner neglects
each agent’s positive or negative effects on the others. For example, if
concurrent tasks require the human and the robot to work in the same area, the
robot would probably slow down because of safety. If the task planner knew
this effect, it could favour pairs of tasks that avoid robot safety slowdowns.

In this work, we propose a method to learn the expected duration of a joint
plan. We leverage a minimum-time formulation of the task planning/allocation
problem for HRC. Then, we estimate the expected duration of each task from
previous executions and learn a synergy coefficient that represents the effect
of one task over the other. A synergy coefficient greater than one means that
the human task causes a slow down of the robot task. We show that it is
possible to cast the estimation of the synergy coefficients into a set of
linear regression problems (one for each task). The resulting synergy
coefficients can be exploited in task planning and allocation method to select
advantageous task couplings. The proposed approach can be used offline – to
obtain a guess of the task duration and synergies – or iteratively as the
number of executions grows to refine the initial guess over time. We
demonstrate the proposed approach in an HRC scenario where a human and a
collaborative robot shall perform a sequence of pick-and-place operations.
Experiments show that pairs of tasks that drive the nearby agents lead to high
synergy coefficients, recognizing favorable and unfavourable pairs.

The paper is organized as follows. Section II defines the task planning
problem for HRC systems. Section III builds on that model to formulate the
task expected duration in terms of average duration and synergy coefficients.
Then, it casts the estimation of such coefficients into a set of linear
regression problems. Section IV describes the software architecture to collect
and process the data from task executions. Section V applies the proposed
methodology to the manufacturing example and shows that the resulting
coefficients reflect the good and bad coupling effects of robot safety stops.
Finally, conclusions and future works are discussed in Section VI.

##  II Preliminaries

###  II-A Task Planning Problem

A task planning problem can be formalized as an optimization problem. Given a
set $\\{H,R\\}$ of agents, i.e., human and robot, and a set of Tasks
$\mathcal{T}=\\{\tau_{i}\\}$, the objective of the problem is to obtain a task
plan and assignment $\pi$ that minimizes the duration of the process.

We denote a task by a tuple $\tau=(l,d,t)$, in which:

  * •

$l\in\\{H,R\\}$ is an assignment variable that specifies which agent can
perform the specific task, i.e., human, or robot;

  * •

$d\in\mathbb{R}^{+}$ is a guess of the task duration;

  * •

$t=\bigl{[}t^{\text{start}}\;;\;t^{\text{end}}\bigl{]}$ is an interval with
endpoints corresponding to the start and end time.

We refer to $\mathcal{T}^{H}$ and $\mathcal{T}^{R}$ as the subsets of
$\mathcal{T}$ such that $l=H$ and $l=R$, respectively.

In this context, it is possible to introduce a binary assignment variable that
defines the allocation of $\tau_{i}\in\mathcal{T}$ to the robot ($a^{R}_{i}$)
or to the human ($a^{H}_{i}$):

$\displaystyle a^{R}_{i}$ $\displaystyle=\left\\{\begin{aligned} &1,&&\text{if
task $\tau_{i}$ is assigned to the robot}\\\
&0,&&\text{otherwise}\end{aligned}\right.$ (1) $\displaystyle a^{H}_{i}$
$\displaystyle=\left\\{\begin{aligned} &1,&&\text{if task $\tau_{i}$ is
assigned to the human }\\\ &0,&&\text{otherwise.}\end{aligned}\right.$

The duration of a plan $\pi$ is the maximum between the duration of the
robot’s and the human’s plan, denoted by $d^{H}_{\pi}$ and $d^{R}_{\pi}$
respectively. The duration of each agent’s plan can be calculated as:

$\displaystyle d^{H}_{\pi}$ $\displaystyle=\sum_{i}d^{H}_{i}a^{H}_{i}$ (2)
$\displaystyle d^{R}_{\pi}$ $\displaystyle=\sum_{i}d^{R}_{i}a^{R}_{i}$

By defining a cost function $J$ that represents the duration of a plan $\pi$
as:

$J=\max\\{d_{\pi}^{H},d_{\pi}^{R}\\}$ (3)

the optimal plan $\pi^{*}$ is:

$\pi^{*}=\operatorname*{argmin}_{\pi}{J}$ (4)

subject to the constraint that $\tau_{i}\in\mathcal{T}$ can only be assigned
during $\pi$:

$a^{H}_{i}+a^{R}_{i}=1\quad\forall\
i\;\;\mbox{s.t.}\;\;\tau_{i}\in\mathcal{T}$ (5)

and assignment, temporal, and causal constraints owed to the process
requirements111If (4) only involves assignment constraints, the problem is a
task allocation; if causal and/or temporal constraints are considered, the
problem becomes a task planning&scheduling problem..

##  III Methodology

The task planning problem assumes that a guess of the task duration exists. A
domain expert usually provides this guess. However, this information is often
unreliable and comes with a sizeable unmodeled uncertainty. Moreover, it does
not account for the effect of the task performed by other agents. Indeed, the
task duration is a function of the tasks executed simultaneously by the other
agent, i.e.:

$\displaystyle d^{R}_{i}$ $\displaystyle=d^{R}_{i}(\tau^{H}_{j})\;\;\forall\
j\;\;\mbox{s.t.}\;\;\tau^{H}_{j}\;\;\mbox{:}\;\;t^{H}_{j}\cap
t^{R}_{i}\neq\emptyset$ $\displaystyle d^{H}_{i}$
$\displaystyle=d^{H}_{i}(\tau^{R}_{j})\;\;\forall\
j\;\;\mbox{s.t.}\;\;\tau^{R}_{j}\;\;\mbox{:}\;\;t^{R}_{j}\cap
t^{H}_{i}\neq\emptyset$

where $\tau^{H}_{j}$ and $\tau^{R}_{j}$ are tasks assigned to humans and
robots. If the coupling between tasks is not modelled, the optimal task (4)
may be unreliable when executed on the real-world system. To overcome this
problem, we extend the formulation given in II-A to include a synergy
coefficient in the task duration.

###  III-A Task Planning Problem with Explicit Task Coupling

For each couple of task indices $(i,j)$, a synergy term is introduced for each
agent and denoted with $s^{R}_{i,j}$ for the robot agent and $s^{H}_{i,j}$ for
the human agent. The synergy term denotes the increment of the duration of
task $\tau_{i}$ when executed by the robot while the human is executing task
$\tau_{j}$. We define this coefficient as:

$s^{R}_{i,j}=\frac{d^{R}_{i,j}}{\hat{d}^{R}_{i}}$ (6)

where $d^{R}_{i,j}$ is the expected duration of task $\tau_{i}$ when the human
executes $\tau_{j}$ and $\hat{d}^{R}_{i}$ is the expected value of the
duration of $\tau_{i}^{R}$ for all concurrent tasks $\tau_{j}^{H}$. Thus, the
duration of a plan $\pi$ becomes:

$\displaystyle d^{R}_{\pi}$
$\displaystyle=\sum_{i}\hat{d}^{R}_{i}a^{R}_{i}\biggl{(}\sum_{j}s^{R}_{i,j}\delta^{R}_{i,j}a^{H}_{j}\biggl{)}$
(7)

where $\delta^{R}_{i,j}$ represents the ratio of the overlapping time between
$\tau^{R}_{i}$ and $\tau^{H}_{j}$ with respect to the duration of the task
$\tau^{R}_{i}$, thus defined as:

$\delta^{R}_{i,j}=\frac{D(t^{H}_{j}\cap t^{R}_{i})}{D(t^{R}_{i})}$ (8)

and $D$ is a function that calculates the duration of a temporal interval
$t=[t^{start},t^{end}]$:

$D(t)=\left\\{\begin{aligned} &t^{end}-t^{start},&&\text{if
$t\neq\emptyset$}\\\ &0,&&\text{if $t=\emptyset$ }\end{aligned}\right.$ (9)

Equations (6), (7), (8) can be rewritten for the human agent as:

$\displaystyle s^{H}_{i,j}$
$\displaystyle=\frac{d^{H}_{i,j}}{\hat{d}^{H}_{i}}$ (10)
$\displaystyle\delta^{H}_{i,j}$ $\displaystyle=\frac{D(t^{R}_{j}\cap
t^{H}_{i})}{D(t^{H}_{i})}$ (11) $\displaystyle d^{H}_{\pi}$
$\displaystyle=\sum_{i}\hat{d}^{H}_{i}a^{H}_{i}\biggl{(}\sum_{j}s^{H}_{i,j}\delta^{H}_{i,j}a^{R}_{j}\biggl{)}$
(12)

Using (7) and (12) in (3), a task planning can minimize the process duration
taking into account the coupling effect between concurrent tasks.

###  III-B Synergy and duration estimation

We estimate both task duration and synergy coefficients from experience. Given
$n$ executions of a task $\tau_{i}$, we approximate its duration with its
expected value:

$\displaystyle\hat{d}^{R}_{i}$
$\displaystyle\approx\mathbb{E}\biggl{[}d^{R}_{i}\rvert_{k}\biggl{]}\quad\forall\
k=\\{1,\dots,n\\}$ (13) $\displaystyle\hat{d}^{H}_{i}$
$\displaystyle\approx\mathbb{E}\biggl{[}d^{H}_{i}\rvert_{k}\biggl{]}\quad\forall\
k=\\{1,\dots,n\\}$

where $d^{R}_{i}\rvert_{k}$ and $d^{H}_{i}\rvert_{k}$ are the task duration
measured in execution $k$.

Furthermore, it is possible to formalize the problem of estimating task
synergy coefficients as a least-square regression problem. For each sample
$k$:

$D(t^{H}_{i})\Big{\rvert}_{k}=\sum_{j=1}^{m}\delta^{H}_{i,j}\Big{\rvert}_{k}s^{R}_{i,j}d^{R}_{i}+T_{\text{idle}}\Big{\rvert}_{k}$
(14)

where $T_{\text{idle}}$ is the time when the robot is not assigned to any task
during $\tau^{H}_{i}$ and $m=|\mathcal{T}^{R}|$ . Equation (14) can be written
in matrix form as:

$\mathbf{D_{i}^{H}}=\mathbf{R^{R}}\mathbf{S^{R}_{i}}+\mathbf{T_{\text{idle}}}$
(15)

where:

  * •

$\mathbf{D_{i}^{H}}\in\mathbb{R}^{n\times 1}$ is a column vector containing
the human execution task duration $D(t^{H}_{i})\Big{\rvert}_{k}$;

  * •

$\mathbf{R^{R}}\in\mathbb{R}^{n\times m}$ is the regression matrix, which
takes the following form:

$\mathbf{R^{R}}=\begin{bmatrix}\delta^{H}_{i,1}\Big{\rvert}_{1}&\cdots&\delta^{H}_{i,m}\Big{\rvert}_{1}\\\
\vdots&\ddots&\vdots\\\
\delta^{H}_{i,1}\Big{\rvert}_{k}&\cdots&\delta^{H}_{i,m}\Big{\rvert}_{k}\\\
\vdots&\ddots&\vdots\\\
\delta^{H}_{i,1}\Big{\rvert}_{n}&\cdots&\delta^{H}_{i,m}\Big{\rvert}_{n}\\\
\end{bmatrix}d^{R}_{i}$ (16)

  * •

$\mathbf{S^{R}_{i}}\in\mathbb{R}^{m\times 1}$ is a column vector containing
the synergy coefficients to be estimated:

$\mathbf{S^{R}_{i}}=\begin{bmatrix}s^{R}_{i,1}&s^{R}_{i,2}&\cdots&s^{R}_{i,m}\end{bmatrix}^{T}$
(17)

  * •

$\mathbf{T_{\text{idle}}}\in\mathbb{R}^{n\times 1}$ is a column vector
containing the robot idle times $T_{\text{idle}}\Big{\rvert}_{k}$.

In conclusion, the solution of the regression problem in (15) is obtained from
the following:

$\mathbf{S^{R}_{i}}=\Bigl{(}\mathbf{R^{R}}^{T}\mathbf{R^{R}}\Bigl{)}^{-1}\mathbf{R^{R}}^{T}\Bigl{[}\mathbf{D_{i}^{H}}-\mathbf{T_{\text{idle}}}\Bigl{]}.$
(18)

and using the same nomenclature for the human agent:

$\mathbf{S^{H}_{i}}=\Bigl{(}\mathbf{R^{H}}^{T}\mathbf{R^{H}}\Bigl{)}^{-1}\mathbf{R^{H}}^{T}\Bigl{[}\mathbf{D_{i}^{R}}-\mathbf{T_{\text{idle}}}\Bigl{]}.$
(19)

##  IV Proposed framework architecture

Figure 1 shows the software architecture developed to measure and estimates
action duration and synergies for real-world deployment of our approach.

Figure 1: Proposed framework architecture. Figure 2: Software Implementation.

The _Task planner_ is at the highest level and has a symbolic knowledge of the
tasks to execute. Based on the temporal constraints between tasks (i.e.,
precedence) and a-priori durations information, it computes the optimized
plan, the scheduling, and the assignment of tasks to agents. This information
is used by the _dispatcher_ , which receives the plan and sends the request at
a lower level; finally, it waits for the result to proceed with the remaining
part of the plan.

The _task planner interface_ module is at the middle layer. The first module
component is the _task manager_ , which receives the task request from the
dispatcher and queries the database to check for properties associated with
the requested symbolic task. In a positive case, the _task manager_ sends the
task execution request to the task executor of the agent specified by the task
planner. When feedback comes from the underlying layer, the _task manager_
module interacts with the _database_ to update the task results. The core of
the _task planner interface_ is a database representing the high-level updated
Knowledge-Base. The database stores information used at different levels: the
definition of high-level task properties, run-time information of duration
execution, statistical information of expected duration, and synergy of
concurrent tasks between agents.

At the lowest level, there is the single-agent _task executor_ , which
converts symbolic information associated with tasks into geometric targets.
First, it interacts with the database to retrieve the action type (e.g., pick,
place, go to) and the symbolic goal associated with the task. This layer has a
geometric knowledge of the system and translates the symbolic goal into a
sequence of robot movements. To do so, it knows the locations of symbolic
goals, integrates a motion planner, and may include an action planner such as
in [13]. Once the task execution is finished, the feedback is sent to the
layer above.

The _task manager_ sends the task request to an _HMI interface_ , which
communicates to the HMI the information of the task to be executed and waits
for the confirmation of execution. Then, it sends the feedback to the _task
manager_.

The advantage of this architecture is that, through the presence of the
database, it is possible to keep the tasks’ representation up-to-date, on
which the high-level planning is based. So, it is possible to update actions’
costs based on the experience of executions. For this purpose, the _task
planner statistics_ module calculates the estimated duration and task synergy
as described in Section III-B.

###  IV-A Software implementation

The functional architecture presented in Section IV was developed as a hybrid
Python/C++ library based on ROS [14].

Software development focused on the task planner interface level. The
crossroads component is the _task service manager_. As shown in Figure 2, this
node exploits the Publisher-Subscriber communication to receive the task
execution request from the dispatcher and then send back the feedback.

The communication with the database occurs through an interface node that acts
as a service server for the task service manager: a service checks for the
existence of task properties associated with the symbolic task, and another
service saves the execution results from the layer below.

Communication with the underlying layer, i.e., _task execution_ nodes, takes
place according to the client-server action model. The task service manager
sends a goal containing a symbolic task name of the task to be executed and
waits for the results, e.g., duration, task type, and agent information.

The database is composed of different collections, each dedicated to storing
different information:

  * •

_task properties_ : for each task, it stores a symbolic identifier, the type
of action associated with it (e.g., pick, place, go to), a textual
description, the agents that can execute it, and the symbolic goals associated
with it.

  * •

_Task results_ : it is updated in real-time with information from the lower
layers about the execution time and the success or failure of the task.

  * •

_Task duration_ : it contains the information of expected durations and
standard deviation associated with each symbolic task.

  * •

_Task synergy_ : for each pair of tasks, it contains the result of the
estimation described in Section III-B.

Notice that both _Task duration_ and _Task synergy_ can be updated
synchronously or asynchronously.

In this framework, MongoDB was used as the database [15]. In addition, PyMongo
library was used for interaction with MongoDB database in the so-called
_Mongo-DB interface node_ shown in Figure 2.

At the lowest level, the _task-execution_ node interacts with the database
interface to retrieve the geometric goal associated with the task. It acts as
a client of the manipulation framework [16], to which it sends requests,
specifying the type of action to be performed. The manipulation framework
defines the sequence of movements required to perform the requested action,
solves the corresponding motion planning problem, and finally executes the
action.

Finally, the _task planner statistics node_ is independent of the framework
information flow and provides ROS-Services to update the _task duration_ and
the _task synergy_ collections and make statistical charts. These services can
be called synchronously or asynchronously concerning the information flow of
the framework.

##  V Experiments

Figure 3: Simulated collaborative workcell

We test our approach in a simulated case study of a collaborative workcell
(Figure 3) composed of a collaborative robot, UR5, mounted on a linear axis
and a human operator. The agents share a portion of the workspace, defining a
collaborative workspace. There are six white cubes on the work table
accessible only by the human operator, six orange cubes on the work table
accessible by the robot, and six blue cubes within the collaborative work area
accessible by both agents. Thus, humans can handle only white and blue
objects, while robots can handle only orange and blue objects.

Each agent has a dedicated object release area, as illustrated in Figure 3.
The process goal is that: (i) the robot picks four orange and two blue boxes
by placing them in its release area; (ii) the human picks four white and two
blue boxes by arranging them in its release area. Therefore, the set of tasks
that the human and the robot can perform is defined by:

$\mathcal{T}^{H}=\\{\text{Pick Orange Box},\text{Place Orange Box},\\\
\text{Pick Blue Box (H)},\text{Place Blue Box (H)}\\}$ (20)
$\mathcal{T}^{R}=\\{\text{Pick White Box},\text{Place White Box},\\\
\text{Pick Blue Box (R)},\text{Place Blue Box (R)}\\}$ (21)

where $\mathcal{T}=\mathcal{T}^{H}\cup\mathcal{T}^{R}$.

The process requires access to the same work area, especially during a
simultaneous execution of the _Pick Blue Box_ and _Place Blue Box_ tasks. The
symbolic goal associated with each task does not refer to a specific pick or
place slot. Thus, the choice of the best slot associated with a symbolic goal
is performed by the action planner of the manipulation framework [16] at
runtime by solving a multi-goal motion planning problem.

The ISO/TS 15066 [17] is applied to ensure safety. Three collaborative areas
are defined. If any human body part enters inside the red area, the robot is
stopped. If the human enters the orange area, the robot moves at 50 % of its
nominal speed. Finally, if the human operator moves into the remaining area,
the robot is free to move at its nominal speed.

###  V-A Results

In order to estimate the synergy coefficients between tasks, 50 plans are
randomly generated respecting the assignment constraints imposed by the
agent’s domains (20) and (21) and precedence constraints between a pick and
place tasks. A simulation run is performed for each recipe, using the
framework proposed in Section IV: tasks are executed, and the results are
saved to the database.

We exploit the _task planner statistics_ node to calculate the expected
duration and standard deviation for each task (grouped by type and agent).
Results are shown in Figure 4.

Figure 4: Task durations grouped by agent, white circles denote the mean.

For each type of task, outliers are removed using an isolation forest
algorithm to reduce noise from the task results before synergy coefficients
estimation [18].

We use the estimated duration of each task and saved task results to apply the
regression described in (15) and obtain the synergy coefficients estimation
between robot and human tasks. We obtain the column vector in (17) by applying
the regression for each task. Repeating this procedure for each task, we
obtain the synergy matrix shown in Figure 5 (a heatmap graphically
highlighting the positive or negative synergy between concurrent tasks). The
heatmap shows the synergy coefficients of the robot agent task versus the
human tasks. The tasks of the robot (21) are placed in the rows of the
heatmap, and on the columns are the human tasks (20). An $s^{R}_{i,j}$ element
of the heatmap reports the duration increment of the robot’s task $\tau_{i}$
when simultaneous with the human task $\tau_{j}$.

If the value $s^{R}_{i,j}$ is greater than 1, it means that task $\tau_{i}$
performed by the robot is penalized when it is in parallel with task
$\tau_{j}$ performed by the human. Vice versa, if $s^{R}_{i,j}$ is smaller
than 1, it means that the coupling between tasks $\tau_{i}$ and $\tau_{j}$ is
advantageous in terms of execution time.

The heatmap analysis shows that all the robot tasks are penalized when the
human simultaneously performs the _Pick Blue Box_ and _Place Blue Box_ tasks.
Indeed, when the human operator executes _Pick Blue Box_ , he accesses a _red
area_ close to the robot, causing the robot to stop. Conversely, the robot
tasks do not slow down or even improve performance when they are executed in
parallel with the human tasks _Pick White Box_ (which does not cause an access
to the collaborative area), and _Place White Box_.

Figure 6 also shows the synergy coefficients and their standard deviation
estimates by the regression to validate the statistical significance of the
estimate.

Figure 5: Robot task Synergy Matrix (each row corresponds to the vector
$S^{R}_{i}$ obtained from (17)) Figure 6: Synergy coefficients and their
standard deviation

##  VI Conclusions and future works

In this paper, a synergy term that takes into account the coupling between
agent tasks is defined. The synergy coefficient is contextualized in a generic
task planning problem by including it in the duration cost function.

We validated the approach in a simple HRC scenario where the positive of
negative couplings were intuitive in order to show that the proposed method
can learn the expected good or bad synergy of pairs of tasks.

Future works will focus on integrating the estimated synergy term in a task
planner to demonstrate the reduction of task plan duration in complex
scenarios.

Alternative strategies will be compared to access the advantages of the
proposed approach.

## References

  * [1] S. Pellegrinelli, A. Orlandini, N. Pedrocchi, A. Umbrico, and T. Tolio, “Motion planning and scheduling for human and industrial-robot collaboration,” _CIRP Annals_ , vol. 66, no. 1, p. 1 : 4, 2017.
  * [2] C. Tonola, M. Faroni, N. Pedrocchi, and M. Beschi, “Anytime informed path re-planning and optimization for human-robot collaboration,” in _Proceedings of the IEEE International Conference on Robot and Human Interactive Communication_ , Vancouver (Canada), 2021.
  * [3] J. Chen, Y. Ding, B. Xin, Q. Yang, and H. Fang, “A unifying framework for human-agent collaborative systems–part i: Element and relation analysis,” _IEEE T. Cybernetics_ , pp. 1–14, 2020.
  * [4] S. Lemaignan, M. Warnier, E. A. Sisbot, A. Clodic, and R. Alami, “Artificial cognition for social human–robot interaction: An implementation,” _Artif. Intell._ , vol. 247, pp. 45–69, 2017.
  * [5] L. De Silva, R. Lallement, and R. Alami, “The hatp hierarchical planner: Formalisation and an initial study of its usability and practicality,” in _IEEE/RSJ Int. Conf. Int. Robots Syst._ , 2015, pp. 6465–6472.
  * [6] G. Milliez, R. Lallement, M. Fiore, and R. Alami, “Using human knowledge awareness to adapt collaborative plan generation, explanation and monitoring,” in _ACM/IEEE Int. Conf. Hum.-Rob. Inter._ , 2016.
  * [7] E. Sebastiani, R. Lallement, R. Alami, and L. Iocchi, “Dealing with on-line human-robot negotiations in hierarchical agent-based task planner,” in _Int. Conf. Automated Plan. Sched._ , 2017.
  * [8] M. Lippi and A. Marino, “A mixed-integer linear programming formulation for human multi-robot task allocation,” in _IEEE International Conference on Robot and Human Interactive Communication_ , 2021, pp. 1017–1023.
  * [9] A. Akbari, Muhayyuddin, and J. Rosell, “Knowledge-oriented task and motion planning for multiple mobile robots,” _J. Exp. Theor. Artif. Intell._ , vol. 31, no. 1, pp. 137–162, 2019.
  * [10] A. Akbari, M. Diab, and J. Rosell, “Contingent task and motion planning under uncertainty for human–robot interactions,” _Applied Sciences_ , vol. 10, no. 5, p. 1665, 2020.
  * [11] M. Cialdea Mayer and A. Orlandini, “An executable semantics of flexible plans in terms of timed game automata,” in _The 22nd International Symposium on Temporal Representation and Reasoning (TIME)_. IEEE, 2015.
  * [12] S. Pellegrinelli, A. Orlandini, N. Pedrocchi, A. Umbrico, and T. Tolio, “Motion planning and scheduling for human and industrial-robot collaboration,” _CIRP Annals_ , vol. 66, pp. 1–4, 2017.
  * [13] M. Faroni, M. Beschi, S. Ghidini, N. Pedrocchi, A. Umbrico, A. Orlandini, and A. Cesta, “A layered control approach to human-aware task and motion planning for human-robot collaboration,” in _Proceedings of the IEEE International Conference on Robot and Human Interactive Communication_ , Naples (Italy), 2020.
  * [14] M. Quigley, K. Conley, B. Gerkey, J. Faust, T. Foote, J. Leibs, R. Wheeler, A. Y. Ng, _et al._ , “Ros: an open-source robot operating system,” in _ICRA workshop on open source software_ , vol. 3, no. 3.2. Kobe, Japan, 2009, p. 5.
  * [15] C. Győrödi, R. Győrödi, G. Pecherle, and A. Olah, “A comparative study: Mongodb vs. mysql,” in _2015 13th International Conference on Engineering of Modern Electric Systems (EMES)_ , 2015, pp. 1–6.
  * [16] E. Villagrossi, N. Pedrocchi, and M. Beschi, “Simplify the robot programming through an action-and-skill manipulation framework,” in _2021 26th IEEE International Conference on Emerging Technologies and Factory Automation (ETFA )_ , 2021, pp. 1–6.
  * [17] “ISO/TS 15066-2016 Robots and robotic devices -: Collaborative robots,” International Organization for Standardization, Geneva, CH, Standard, 2016.
  * [18] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine learning in Python,” _Journal of Machine Learning Research_ , vol. 12, pp. 2825–2830, 2011.


# An Optimal Stochastic Algorithm for Decentralized Nonconvex Finite-sum
Optimization

Luo Luo Haishan Ye

###### Abstract

This paper studies the decentralized nonconvex optimization problem
$\min_{x\in{\mathbb{R}}^{d}}f(x)\triangleq\frac{1}{m}\sum_{i=1}^{m}f_{i}(x)$,
where $f_{i}(x)\triangleq\frac{1}{n}\sum_{j=1}^{n}f_{i,j}(x)$ is the local
function on the $i$-th agent of the network. We propose a novel stochastic
algorithm called DEcentralized probAbilistic Recursive gradiEnt deScenT
(DEAREST), which integrates the techniques of variance reduction, gradient
tracking and multi-consensus. We construct a Lyapunov function that
simultaneously characterizes the function value, the gradient estimation error
and the consensus error for the convergence analysis. Based on this measure,
we provide a concise proof to show DEAREST requires at most
${\mathcal{O}}(mn+\sqrt{mn}L\varepsilon^{-2})$ incremental first-order oracle
(IFO) calls and
${\mathcal{O}}({L\varepsilon^{-2}}/{\sqrt{1-\lambda_{2}(W)}}\,)$ communication
rounds to find an $\varepsilon$-stationary point in expectation, where $L$ is
the smoothness parameter and $\lambda_{2}(W)$ is the second-largest eigenvalue
of the gossip matrix $W$. We can verify both of the IFO complexity and
communication complexity match the lower bounds. To the best of our knowledge,
DEAREST is the first optimal algorithm for decentralized nonconvex finite-sum
optimization.

##  1 Introduction

We study decentralized optimization problem

$\displaystyle\min_{x\in{\mathbb{R}}^{d}}f(x)\triangleq\frac{1}{m}\sum_{i=1}^{m}f_{i}(x),$
(1)

where $m$ agents form a connected network and each $f_{i}(x)$ is the local
function on the $i$-th agent that can be written as the average of $n$
individual functions as follows

$\displaystyle f_{i}(x)\triangleq\frac{1}{n}\sum_{j=1}^{n}f_{i,j}(x).$ (2)

We suppose each individual function $f_{i,j}(x)$ is smooth but possibly
nonconvex. This formulation includes a lot of applications in statistics,
signal processing and machine learning [24, 39, 36, 10, 34, 33, 3, 5]. In
decentralized scenario, all of agents target to collaboratively minimize the
global function $f(x)$ and each of them can only communicate with its
neighbours. This paper focus on designing the efficient algorithm to find an
expected ${\varepsilon}$-stationary point $x_{\rm out}\in{\mathbb{R}}^{d}$
such that ${\mathbb{E}}\left\|\nabla f(x_{\rm out})\right\|\leq{\varepsilon}$.

Incremental first-order oracle (IFO) algorithms [9, 15, 29, 2, 1, 13, 44, 30,
45, 25, 11, 19, 26] are the most popular methods to solve the optimization
problem with finite-sum objective. They evaluate the first-order information
by accessing the individual functions to reduce the computational cost. For
the optimization on single machine, the algorithms with stochastic recursive
gradient estimator [25, 11, 26, 19] achieve the optimal IFO complexity of
${\mathcal{O}}\big{(}N+\sqrt{N}L{\varepsilon}^{-2}\big{)}$, where $N$ is the
total number of individual functions and $L$ is the average-smoothness
parameter. In decentralized setting, the topology of the network leads to the
local machine cannot access the global information at each round and the
communication cost could be expensive. The gradient tracking [23, 27, 28] is a
useful technique to track the average of local gradients and make the first-
order information accurate. Ye et al. [40, 41, 42] showed that combining the
gradient tracking with multi-consensus [17, 4] encourages the analysis of
decentralized algorithm be close to the centralize counterpart and it is
useful to decentralized convex optimization.

Table 1: We compare the order of upper bound complexities for our algorithm with related work. Here we use $\lambda_{2}$ to abbreviate the second-largest eigenvalue of gossip matrix $\lambda_{2}(W)$. The column entitled with “#IFO” presents the IFO complexity on per-agent. Algorithms | #IFO | #Communication
---|---|---
D-GET [32] | $n+\dfrac{\sqrt{n}L{\varepsilon}^{-2}}{(1-\lambda_{2})^{2}}$ | same as #IFO
GT-SARAH [38] | $n+\left(\dfrac{1}{(1-\lambda_{2})^{2}}+\sqrt{\dfrac{n}{m}}+\dfrac{(n/m+1)^{1/3}}{1-\lambda_{2}}\right)L{\varepsilon}^{-2}$ | same as #IFO
DESTRESS [18] | $n+\sqrt{\dfrac{n}{m}}\,L{\varepsilon}^{-2}$ | $\dfrac{\sqrt{mn}+L{\varepsilon}^{-2}}{\sqrt{1-\lambda_{2}}}$
DEAREST (this work) | $n+\sqrt{\dfrac{n}{m}}\,L{\varepsilon}^{-2}$ | $\dfrac{L{\varepsilon}^{-2}}{\sqrt{1-\lambda_{2}}}$

Directly extending (stochastic) gradient descent methods to decentralized
setting [20, 22, 35, 8, 7, 43, 31] does not take the advantage of finite-sum
structure in our local function (2). Sun et al. [32] first applied variance
reduction and gradient tracking to decentralized nonconvex finite-sum
optimization and proposed the algorithm called Decentralized Gradient
Estimation and Tracking (D-GET). Later, Xin et al. [38] proposed GT-SARAH,
which improved the complexity of D-GET in terms of the dependency on $m$ and
$n$. Li et al. [18] further improved the result of GT-SARAH by proposing
DEcentralized STochastic REcurSive gradient methodS (DESTRESS), which requires
${\mathcal{O}}\big{(}n+\sqrt{n/m}L{\varepsilon}^{-2}\big{)}$ per-agent IFO
calls and
${\mathcal{O}}\big{(}\sqrt{mn}+L{\varepsilon}^{-2}/\sqrt{1-\lambda_{2}(W)}\,\big{)}$
rounds of communication to obtain an ${\varepsilon}$-stationary point in
expectation, where $\lambda_{2}(W)$ is the second-largest eigenvalues of the
gossip matrix $W$. Since the formulations (1) and (2) mean there are $mn$
individual functions, the total IFO upper bound of DESTRESS for all agents is
${\mathcal{O}}\big{(}mn+\sqrt{mn}\,L{\varepsilon}^{-2}\big{)}$, matching the
best result for single machine optimization and it is optimal. However, the
lower bound of communication complexity for finding an
${\varepsilon}$-stationary point of smoothed nonconvex function by first-order
algorithm [22, 31] is
$\Omega\big{(}L{\varepsilon}^{-2}/\sqrt{1-\lambda_{2}(W)}\,\big{)}$. The
additional term related to $\sqrt{mn}$ in the communication complexity of
DESTRESS means the algorithm could be communication inefficient when $mn$ is
extremely large.

In this paper, we proposed a novel algorithm named as DEcentralized
probAbilistic Recursive gradiEnt deScenT (DEAREST), which simultaneously
matches the IFO lower bound [11, 19] and communication lower bound [22, 31]
for decentralized nonconvex finite-sum optimization. We compare the
theoretical guarantees of this paper and related work in Table 1. The design
of DEAREST is based on the framework of loopless variance reduced iteration
[14, 12, 16, 19, 41], gradient tracking [23, 27, 28] and multi-consensus [21].
We construct a Lyapunov function to characterize the function value, the
global gradient estimation error, the local gradient estimation error and the
consensus error. By verifying this measure is non-increasing, we provide a
concise analysis to show the tightness for the IFO upper bound and the
communication upper bound of DEAREST. To the best of our knowledge, DEAREST is
the first optimal algorithm for decentralized nonconvex finite-sum
optimization.

##  2 Notations and Assumptions

We use $\|\cdot\|$ to denote the Frobenius norm of a matrix or the
$\ell_{2}$-norm of a vector. We let
${\bf{1}}=[1,\cdots,1]^{\top}\in{\mathbb{R}}^{m}$ and denote $I$ as the
identity matrix. We define aggregate variables for all agents as

$\displaystyle{\bf{x}}=\begin{bmatrix}{\bf{x}}(1)^{\top}\\\ \vdots\\\
{\bf{x}}(m)^{\top}\end{bmatrix}\in{\mathbb{R}}^{m\times d},$

where each ${\bf{x}}(i)\in{\mathbb{R}}^{d}$ are the local variable on the
$i$-th agent. We use the lower case with the bar to represent the mean vector,
such that $\bar{x}=\frac{1}{m}\sum_{i=1}^{m}{\bf{x}}(i)\in{\mathbb{R}}^{d}$.

We also introduce the matrix of aggregate gradients as

$\displaystyle\nabla{\bf{f}}({\bf{x}})=\begin{bmatrix}\nabla
f_{1}({\bf{x}}(1))^{\top}\\\ \vdots\\\ \nabla
f_{m}({\bf{x}}(m))^{\top}\end{bmatrix}\in{\mathbb{R}}^{m\times d}.$

We use the gossip (mixing) matrix $W\in{\mathbb{R}}^{m\times m}$ to
characterize the behavior of agents updating local variables by weighted sum
of information form the neighbors. We suppose the matrix $W$ satisfies the
following assumption.

######  Assumption 1.

We suppose the gossip (mixing) matrix $W=[w_{ij}]\in{\mathbb{R}}^{m\times m}$
is symmetric and satisfies $w_{ij}=0$ if the $i$-th agent and the $j$-th agent
are not connected. Furthermore, we suppose it holds that $W{\bf{1}}={\bf{1}}$.

We assume the objective function in (1) and the local function (2) satisfy the
following assumptions.

######  Assumption 2.

We suppose $f(x)$ is lower bounded. That is, we have

$\displaystyle f^{*}=\inf_{x\in{\mathbb{R}}^{d}}f(x)>-\infty.$

######  Assumption 3.

We suppose the individual functions $\\{f_{i,j}\\}_{j=1}^{n}$ on each agent
are $L$-average smooth for some $L>0$. That is, we have

$\displaystyle\frac{1}{n}\sum_{j=1}^{n}\left\|\nabla f_{i,j}(x)-\nabla
f_{i,j}(x^{\prime})\right\|^{2}\leq L^{2}\left\|x-x^{\prime}\right\|^{2}$

for any $x$ and $x^{\prime}\in{\mathbb{R}}^{d}$.

##  3 The Algorithm and Main Results

Algorithm 1 DEAREST

1: Input: initial parameter ${{\bar{x}}}_{0}\in{\mathbb{R}}^{d}$, stepsize
$\eta>0$, probability $p\in(0,1]$, mini-batch size $b$, numbers of
communication rounds $K_{\rm in}$, $K$ and $\hat{K}$.

2: ${\bf{x}}_{0}={\bf{1}}{{\bar{x}}}_{0}$,
${\bf{g}}_{0}=\nabla{\bf{f}}({\bf{x}}_{0})$

3: ${\bf{s}}_{0}={{\texttt{FastMix}}}({\bf{g}}_{0},K_{\rm in})$

4: for $t=0,\dots,T-1$ do

5: $y_{t}\sim{\rm Bernoulli}(p)$

6: $\displaystyle{K_{t}=\begin{cases}\begin{aligned}
&K,~{}~{}~{}~{}~{}\text{if~{}}y_{t}=1,\\\\[2.84544pt]
&\hat{K},~{}~{}~{}~{}\text{otherwise},\\\ \end{aligned}\end{cases}}$

7: ${\bf{x}}_{t+1}={{\texttt{FastMix}}}({\bf{x}}_{t}-\eta{\bf{s}}_{t},K_{t})$

8: parallel for $i=1,\dots,m$ do

9: $\displaystyle{{\bf{g}}_{t+1}(i)=\begin{cases}\begin{aligned}
&\displaystyle{\nabla
f_{i}({\bf{x}}_{t+1}(i))},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{if~{}}y_{t}=1,\\\\[2.84544pt]
&\displaystyle{{\bf{g}}_{t}(i)+\dfrac{1}{b}\sum_{j=1}^{b}\big{(}\nabla
f_{i,\xi_{j}}({\bf{x}}_{t+1}(i))-\nabla
f_{i,\xi_{j}}({\bf{x}}_{t}(i))\big{)}},~{}~{}~{}~{}\text{otherwise},\\\
\end{aligned}\end{cases}}$ where each $\xi_{i,\xi_{j}}$ is uniformly and
independently sampled from $\\{1,\dots,n\\}$

10: end parallel for

11:
${\bf{s}}_{t+1}={{\texttt{FastMix}}}({\bf{s}}_{t}+{\bf{g}}_{t+1}-{\bf{g}}_{t},K_{t})$

12: end for

13: Output: $x_{\rm out}$ by uniformly sampling from
$\\{{\bf{x}}_{t}(i):t=0,\dots,T-1;~{}i=1,\dots,m\\}$

We propose DEcentralized probAbilistic Recursive gradiEnt deScenT (DEAREST) in
Algorithm 1. It uses the subroutine FastMix (Algorithm 2) and gradient
tracking (line 11) to reduce the consensus error. The update rules of local
variables can be written as

$\displaystyle\begin{cases}{\bf{x}}_{t+1}={{\texttt{FastMix}}}({\bf{x}}_{t}-\eta{\bf{s}}_{t},K_{t}),\\\
{\bf{s}}_{t+1}={{\texttt{FastMix}}}({\bf{s}}_{t}+{\bf{g}}_{t+1}-{\bf{g}}_{t},K_{t}),\end{cases}$
(3)

which has the following properties.

######  Lemma 1 ([40, Lemma 2]).

For Algorithm 1, we have ${{\bar{s}}}_{t}={{\bar{g}}}_{t}$.

######  Lemma 2 ([21, Lemma 2]).

Under Assumption 1, Algorithm 2 holds that

$\displaystyle\frac{1}{m}{\bf{1}}^{\top}{\bf{u}}^{(K)}={{\bar{u}}}^{(0)}\quad\text{and}\quad\big{\|}{\bf{u}}^{(K)}-{\bf{1}}{{\bar{u}}}^{(0)}\big{\|}\leq\left(1-\sqrt{1-\lambda_{2}(W)}\right)^{K}\big{\|}{\bf{u}}^{(0)}-{\bf{1}}{{\bar{u}}}^{(0)}\big{\|},$

where ${{\bar{u}}}^{(0)}=\frac{1}{m}{\bf{1}}^{\top}{\bf{u}}^{(0)}$.

Using Lemma 1, we can analyze the update (3) by focusing on the relationship
of the corresponding mean vectors, that is

$\displaystyle\begin{split}{{\bar{x}}}_{t+1}=&\frac{1}{m}{\bf{1}}^{\top}{{\texttt{FastMix}}}({\bf{x}}_{t}-\eta{\bf{s}}_{t},K_{t})\\\
=&\frac{1}{m}{\bf{1}}^{\top}({\bf{x}}_{t}-\eta{\bf{s}}_{t})\\\
=&{{\bar{x}}}_{t}-\eta{{\bar{s}}}_{t}={{\bar{x}}}_{t}-\eta{{\bar{g}}}_{t}.\end{split}$
(4)

The equation (4) can be viewed as the step of inexact gradient descent on the
mean variable ${{\bar{x}}}_{t}$ with gradient estimator ${{\bar{s}}}_{t}$ and
stepsize $\eta$. Additionally, Lemma 2 leads to each ${\bf{x}}_{t}(i)$ and
${\bf{g}}_{t}(i)$ be close to ${{\bar{x}}}_{t}$ and ${{\bar{g}}}_{t}$
respectively.

The construction of gradient estimator ${\bf{g}}_{t}(i)$ follows the
probabilistic recursive way [19]. Intuitively, it estimates the gradient
recursively by using the mini-batch gradient with high probability $1-p$ and
computing the exact gradient with low probability $p$, where
$p=\Omega(1/\sqrt{mn}\,)$. By taking the mini-batch size
$b={\mathcal{O}}\big{(}\sqrt{n/m}\,\big{)}$, the stepsize $\eta=\Theta(1/L)$
and the iteration number $T={\mathcal{O}}\left(L{\varepsilon}^{-2}\right)$,
the IFO complexity on per-agent is at most
${\mathcal{O}}\big{(}n+\sqrt{n/m}\,L{\varepsilon}^{-2}\big{)}$ in expectation.
Note that each agent uses the same random seed to generate the Bernoulli
distributed variable $y_{t}$, which enforces all agents always share the
identical $y_{t}$.

The communication rounds $K_{t}$ at each iteration is also decided by the
random variable $y_{t}$. Concretely, we let
$K_{t}=K=\Theta\big{(}{\log(n/m)}/{\sqrt{1-\lambda_{2}(W)}}\,\big{)}$ with
probability $p$ and
$K_{t}=\hat{K}={\mathcal{O}}\big{(}1/{\sqrt{1-\lambda_{2}(W)}}\,\big{)}$ with
probability $1-p$. Combining with the initial communication rounds $K_{\rm
in}={\mathcal{O}}\big{(}\log(1/(m{\varepsilon}))/{\sqrt{1-\lambda_{2}(W)}}\,\big{)}$,
the total communication complexity is
${\mathcal{O}}\big{(}{L{\varepsilon}^{-2}}/{\sqrt{1-\lambda_{2}(W)}}\,\big{)}$
in expectation.

Algorithm 2 ${{\texttt{FastMix}}}({\bf{u}}^{(0)},K)$

1: Initialize: ${\bf{u}}^{(-1)}={\bf{u}}^{(0)}$,
$\eta_{u}=\frac{1-\sqrt{1-\lambda_{2}^{2}(W)}}{1+\sqrt{1-\lambda_{2}^{2}(W)}}$.

2: for $k=0,1,\dots,K$ do

3: ${\bf{u}}^{(k+1)}=(1+\eta_{u})W{\bf{u}}^{(k)}-\eta_{u}{\bf{u}}^{(k-1)}$

4: end for

5: Output: ${\bf{u}}^{(K)}$.

Based on above settings, we guarantee the output of DEAREST is an
${\varepsilon}$-stationary point in expectation. We formally describe our main
result in Theorem 1.

######  Theorem 1.

Under Assumption 1, 2 and 3, we run Algorithm 1 with

$\displaystyle\begin{split}&\eta=\frac{1}{2L},\qquad b=\left\lceil
6\sqrt{\frac{n}{m}}\,\right\rceil,\qquad p=\frac{b}{b+n},\qquad
T=\left\lceil\frac{16L(f({{\bar{x}}}_{0})-f^{*})}{{\varepsilon}^{2}}\right\rceil,\\\
&K=\left\lceil\frac{\max\big{\\{}12,\log((\sqrt{mn}+6)/(24m))\big{\\}}}{2\sqrt{1-\lambda_{2}(W)}}\right\rceil,\qquad\hat{K}=\left\lceil\frac{6}{\sqrt{1-\lambda_{2}(W)}}\right\rceil\\\
&\text{and}\quad K_{\rm
in}=\left\lceil\frac{\log\big{(}{\left\|{\bf{g}}_{0}-{\bf{1}}{{\bar{g}}}_{0}\right\|^{2}}/(m{\varepsilon}^{2})\big{)}}{\sqrt{1-\lambda_{2}(W)}}\right\rceil.\end{split}$

Then the output $x_{\rm out}$ satisfies ${\mathbb{E}}\left\|\nabla f(x_{\rm
out})\right\|\leq{\varepsilon}$. Furthermore, the algorithm requires at most
${\mathcal{O}}\big{(}n+\sqrt{n/m}\,L{\varepsilon}^{-2}\big{)}$ IFO calls on
per-agent and
${\mathcal{O}}\big{(}{L{\varepsilon}^{-2}}/{\sqrt{1-\lambda_{2}(W)}}\,\big{)}$
rounds of communication in expectation.

Now we give a brief sketch for the proof of Theorem 1 and the details are
shown in Section 4. We first introduce the following quantities:

  * •

the global gradient estimation error:
$U_{t}=\left\|\frac{1}{m}\sum_{i=1}^{m}\left({\bf{g}}_{t}(i)-\nabla
f_{i}({\bf{x}}_{t}(i))\right)\right\|^{2}$;

  * •

the local gradient estimation error:
$V_{t}=\frac{1}{m}{\mathbb{E}}\left\|{\bf{g}}_{t}-\nabla{\bf{f}}({\bf{x}}_{t})\right\|^{2}$;

  * •

the consensus error:
$C_{t}=\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}+\eta^{2}\left\|{\bf{s}}_{t}-{\bf{1}}{{\bar{s}}}_{t}\right\|^{2}$.

Then we define the Lyapunov function:

$\displaystyle\Phi_{t}\triangleq
f({{\bar{x}}}_{t})+\frac{\eta}{p}U_{t}+\frac{\eta}{p}V_{t}+\frac{1}{m\eta}C_{t}.$

By establishing the recursion for each terms of $\Phi_{t}$, we can achieve

$\displaystyle\frac{\eta}{2}{\mathbb{E}}\left\|\nabla
f({{\bar{x}}}_{t})\right\|^{2}\leq{\mathbb{E}}\left[\Phi_{t}-\Phi_{t+1}-\frac{1}{m\eta}C_{t}\right].$
(5)

Taking the average on (5) over $t=0,\dots,T-1$, we have

$\displaystyle\frac{\eta}{2T}\sum_{t=0}^{T}{\mathbb{E}}\left\|\nabla
f({{\bar{x}}}_{t})\right\|^{2}\leq{\mathbb{E}}\left[\frac{\Phi_{0}-\Phi_{T}}{T}-\frac{1}{m\eta
T}\sum_{t=0}^{T}C_{t}\right].$ (6)

On the other hand, the multi-consensus steps (FastMix) in DEAREST guarantee
the difference between ${\mathbb{E}}\|\nabla f(x_{\rm out})\|^{2}$ and
$\frac{1}{T}\sum_{t=0}^{T-1}{\mathbb{E}}\left\|\nabla
f({{\bar{x}}}_{t})\right\|^{2}$ could be controlled by the term related to
$\sum_{t=0}^{T}C_{t}$ in (6). As a result, we can prove the output $x_{\rm
out}$ is an ${\varepsilon}$-stationary point in expectation by taking
$\eta=\Omega(1/L)$ and $T={\mathcal{O}}(L{\varepsilon}^{-2})$.

##  4 Complexity Analysis

We first provide some useful lemmas.

######  Lemma 3 ([40, Lemma 3]).

For any ${\bf{z}}\in{\mathbb{R}}^{m\times d}$, we have
$\left\|{\bf{z}}-{\bf{1}}{{\bar{z}}}\right\|\leq\left\|{\bf{z}}\right\|$ where
${{\bar{z}}}=\frac{1}{m}{\bf{1}}^{\top}{\bf{z}}$.

######  Lemma 4.

Under Assumption 3, we have
$\left\|\nabla{\bf{f}}({\bf{x}})-\nabla{\bf{f}}({\bf{x}}^{\prime})\right\|\leq
L\left\|{\bf{x}}-{\bf{x}}^{\prime}\right\|$ for any
${\bf{x}},{\bf{x}}^{\prime}\in{\mathbb{R}}^{m\times d}$.

###### Proof.

Assumption 3 means each $f_{i}(\cdot)$ is $L$-Lipschitz continuous. Then we
have

$\displaystyle\left\|\nabla{\bf{f}}({\bf{x}})-\nabla{\bf{f}}({\bf{x}}^{\prime})\right\|$
$\displaystyle=$ $\displaystyle\sum_{i=1}^{m}\left\|\nabla
f_{i}({\bf{x}}(i))-\nabla f_{i}({\bf{x}}^{\prime}(i))\right\|^{2}$
$\displaystyle\leq$ $\displaystyle
L^{2}\sum_{i=1}^{m}\left\|{\bf{x}}(i)-{\bf{x}}^{\prime}(i)\right\|^{2}$
$\displaystyle=$ $\displaystyle
L^{2}\left\|{\bf{x}}-{\bf{x}}^{\prime}\right\|^{2}.$

∎

We describe the decease of function value in following lemma.

######  Lemma 5.

Algorithm 1 holds that

$\displaystyle f({{\bar{x}}}_{t+1})\leq
f({{\bar{x}}}_{t})-\frac{\eta}{2}\left\|\nabla
f({{\bar{x}}}_{t})\right\|^{2}+\eta
U_{t}+\frac{L^{2}\eta}{m}C_{t}-\left(\frac{1}{2\eta}-\frac{L}{2}\right)\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}.$
(7)

###### Proof.

Using Lemma 2 of Li et al. [19] and the relation (4), we have

$\displaystyle\begin{split}f({{\bar{x}}}_{t+1})\leq
f({{\bar{x}}}_{t})-\frac{\eta}{2}\left\|\nabla
f({{\bar{x}}}_{t})\right\|^{2}-\left(\frac{1}{2\eta}-\frac{L}{2}\right)\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}+\frac{\eta}{2}\left\|{{\bar{s}}}_{t}-\nabla
f({{\bar{x}}}_{t})\right\|^{2}.\end{split}$ (8)

We bound the last term of (8) as follows

$\displaystyle\begin{split}&\left\|{{\bar{s}}}_{t}-\nabla
f({{\bar{x}}}_{t})\right\|^{2}\\\
=&\left\|\frac{1}{m}\sum_{i=1}^{m}\big{(}{\bf{g}}_{t}(i)-\nabla
f_{i}({{\bar{x}}}_{t})\big{)}\right\|^{2}\\\
\leq&2\left\|\frac{1}{m}\sum_{i=1}^{m}\big{(}{\bf{g}}_{t}(i)-\nabla
f_{i}({\bf{x}}_{t}(i))\big{)}\right\|^{2}+2\left\|\frac{1}{m}\sum_{i=1}^{m}\big{(}\nabla
f_{i}({\bf{x}}_{t}(i))-\nabla f_{i}({{\bar{x}}}_{t})\big{)}\right\|^{2}\\\
\leq&2\left\|\frac{1}{m}\sum_{i=1}^{m}\big{(}{\bf{g}}_{t}(i)-\nabla
f_{i}({\bf{x}}_{t}(i))\big{)}\right\|^{2}+\frac{2}{m}\sum_{i=1}^{m}\left\|\nabla
f_{i}({\bf{x}}_{t}(i))-\nabla f_{i}({{\bar{x}}}_{t})\right\|^{2}\\\
\leq&2\left\|\frac{1}{m}\sum_{i=1}^{m}\big{(}{\bf{g}}_{t}(i)-\nabla
f_{i}({\bf{x}}_{t}(i))\big{)}\right\|^{2}+\frac{2L^{2}}{m}\sum_{i=1}^{m}\left\|{\bf{x}}_{t}(i)-{{\bar{x}}}_{t}\right\|^{2}\\\
=&2\left\|\frac{1}{m}\sum_{i=1}^{m}\big{(}{\bf{g}}_{t}(i)-\nabla
f_{i}({\bf{x}}_{t}(i))\big{)}\right\|^{2}+\frac{2L^{2}}{m}\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}\\\
\leq&2U_{t}+\frac{2L^{2}}{m}C_{t},\end{split}$ (9)

where the equality uses Lemma 1; the first inequality uses Young’s inequality;
the second inequality uses the fact
$\left\|\frac{1}{m}\sum_{i=1}^{m}a_{i}\right\|^{2}\leq\frac{1}{m}\sum_{i=1}^{m}\left\|a_{i}\right\|^{2}$
for $a_{1},\dots,a_{m}\in{\mathbb{R}}^{d}$; the third inequality is based on
the Assumption 3; the last step is based on the definitions of $U_{t}$ and
$C_{t}$. We finish the proof by combining the results of (8) and (9). ∎

We define the variables
$\rho_{t}=\big{(}1-\sqrt{1-\lambda_{2}(W)}\,\big{)}^{K_{t}}$,
$\rho=\big{(}1-\sqrt{1-\lambda_{2}(W)}\,\big{)}^{K}$ and
$\hat{\rho}=\big{(}1-\sqrt{1-\lambda_{2}(W)}\,\big{)}^{\hat{K}}$ to
characterize the effect of FastMix (Algorithm 2). Note that the setting of
$K_{t}$ in Theorem 1 means

$\displaystyle\rho\leq\rho_{t}\leq\hat{\rho}<\frac{1}{\sqrt{296}}.$ (10)

Then we provide the recursion of consensus error as follows.

######  Lemma 6.

Under the setting of Theorem 1, we have

$\displaystyle{\mathbb{E}}\left[C_{t+1}\right]\leq{\mathbb{E}}\left[\frac{1}{4}C_{t}+4\rho^{2}p\eta^{2}mV_{t}+6\left(\dfrac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\right)m\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}\right].$

###### Proof.

The relation of (4) means

$\displaystyle\left\|{\bf{x}}_{t+1}-{\bf{1}}{{\bar{x}}}_{t+1}\right\|$
$\displaystyle=$
$\displaystyle\left\|{{\texttt{FastMix}}}({\bf{x}}_{t}-\eta{\bf{s}}_{t},K_{t})-\frac{1}{m}{\bf{1}}{\bf{1}}^{\top}{{\texttt{FastMix}}}({\bf{x}}_{t}-\eta{\bf{s}}_{t},K_{t})\right\|$
$\displaystyle\leq$
$\displaystyle\rho_{t}\left\|({\bf{x}}_{t}-\eta{\bf{s}}_{t})-\frac{1}{m}{\bf{1}}{\bf{1}}^{\top}\left({\bf{x}}_{t}-\eta{\bf{s}}_{t}\right)\right\|$
$\displaystyle=$
$\displaystyle\rho_{t}\left\|{\bf{x}}_{t}-\eta{\bf{s}}_{t}-{\bf{1}}({{\bar{x}}}_{t}-\eta{{\bar{s}}}_{t})\right\|$
$\displaystyle\leq$
$\displaystyle\rho_{t}\big{(}\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|+\eta\left\|{\bf{s}}_{t}-{\bf{1}}{{\bar{s}}}_{t}\right\|\big{)},$

where the first inequality is based on Lemma 2 and the last step is due to
triangle inequality.

Consequently, we use Young’s inequality to obtain

$\displaystyle\begin{split}&\left\|{\bf{x}}_{t+1}-{\bf{1}}{{\bar{x}}}_{t+1}\right\|^{2}\leq
2\rho_{t}^{2}\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}+2\rho_{t}^{2}\eta^{2}\left\|{\bf{s}}_{t}-{\bf{1}}{{\bar{s}}}_{t}\right\|^{2}=2\rho_{t}^{2}C_{t}\leq
2\hat{\rho}^{2}C_{t}.\end{split}$ (11)

The update of ${\bf{g}}_{t+1}(i)$ means

$\displaystyle\begin{split}&{\mathbb{E}}\left[\rho_{t}^{2}\left\|{\bf{g}}_{t+1}(i)-{\bf{g}}_{t}(i)\right\|^{2}\right]\\\
\leq&\rho^{2}p{\mathbb{E}}\left\|\nabla
f_{i}({\bf{x}}_{t+1}(i))-{\bf{g}}_{t}(i)\right\|^{2}+\frac{\hat{\rho}^{2}(1-p)}{b}{\mathbb{E}}\left\|\nabla
f_{i,\xi_{j}}({\bf{x}}_{t+1}(i))-\nabla
f_{i,\xi_{j}}({\bf{x}}_{t}(i))\right\|^{2}\\\
\leq&2\rho^{2}p{\mathbb{E}}\left\|\nabla f_{i}({\bf{x}}_{t+1}(i))-\nabla
f_{i}({\bf{x}}_{t}(i))\right\|^{2}+2\rho^{2}p{\mathbb{E}}\left\|\nabla
f_{i}({\bf{x}}_{t}(i))-{\bf{g}}_{t}(i)\right\|^{2}\\\
&+\frac{\hat{\rho}^{2}(1-p)L^{2}}{b}{\mathbb{E}}\left\|{\bf{x}}_{t+1}(i)-{\bf{x}}_{t}(i)\right\|^{2}\\\
\leq&2\rho^{2}pL^{2}{\mathbb{E}}\left\|{\bf{x}}_{t+1}(i)-{\bf{x}}_{t}(i)\right\|^{2}+2\rho^{2}p{\mathbb{E}}\left\|\nabla
f_{i}({\bf{x}}_{t}(i))-{\bf{g}}_{t}(i)\right\|^{2}\\\
&+\frac{\hat{\rho}^{2}(1-p)L^{2}}{b}{\mathbb{E}}\left\|{\bf{x}}_{t+1}(i)-{\bf{x}}_{t}(i)\right\|^{2}\\\
=&2\rho^{2}p{\mathbb{E}}\left\|\nabla
f_{i}({\bf{x}}_{t}(i))-{\bf{g}}_{t}(i)\right\|^{2}+\left(\frac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\right)L^{2}{\mathbb{E}}\left\|{\bf{x}}_{t+1}(i)-{\bf{x}}_{t}(i)\right\|^{2},\end{split}$
(12)

where the first inequality is based on the update rule; the second inequality
is based on Young’s inequality and the last inequality is due to Assumption 3
and (10).

Summing over (12) above result over $i=1,\dots,m$, we obtain

$\displaystyle\begin{split}&{\mathbb{E}}\left[\rho_{t}^{2}\left\|{\bf{g}}_{t+1}-{\bf{g}}_{t}\right\|^{2}\right]\\\
\leq&2\rho^{2}p{\mathbb{E}}\left\|\nabla{\bf{f}}({\bf{x}}_{t})-{\bf{g}}_{t}\right\|^{2}+\left(\frac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\right)L^{2}{\mathbb{E}}\left\|{\bf{x}}_{t+1}-{\bf{x}}_{t}\right\|^{2}\\\
\leq&2\rho^{2}p{\mathbb{E}}\left\|\nabla{\bf{f}}({\bf{x}}_{t})-{\bf{g}}_{t}\right\|^{2}\\\
&+3\left(\frac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\right)L^{2}{\mathbb{E}}\left[\left\|{\bf{x}}_{t+1}-{\bf{1}}{{\bar{x}}}_{t+1}\right\|^{2}+\left\|{\bf{1}}{{\bar{x}}}_{t+1}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}+\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}\right]\\\
\leq&2\rho^{2}p{\mathbb{E}}\left\|\nabla{\bf{f}}({\bf{x}}_{t})-{\bf{g}}_{t}\right\|^{2}\\\
&+3\left(\frac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\right)L^{2}\left(2\rho^{2}\left(\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}+\eta^{2}\left\|{\bf{s}}_{t}-{\bf{1}}{{\bar{s}}}_{t}\right\|^{2}\right)\right)\\\
&+3\left(\frac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\right)\left(m{\mathbb{E}}\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}+{\mathbb{E}}\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}\right)\\\
=&2\rho^{2}p{\mathbb{E}}\left\|\nabla{\bf{f}}({\bf{x}}_{t})-{\bf{g}}_{t}\right\|^{2}+9\left(\frac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\right)L^{2}\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}\\\
&+6\left(\frac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\right)L^{2}\eta^{2}\left\|{\bf{s}}_{t}-{\bf{1}}{{\bar{s}}}_{t}\right\|^{2}\\\
&+3\left(\frac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\right)L^{2}m{\mathbb{E}}\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2},\end{split}$
(13)

where the second inequality is based on Young’s inequality and the third
inequality uses the result of (11). Using Lemma 2, we have

$\displaystyle\begin{split}&\left\|{\bf{s}}_{t+1}-{\bf{1}}{{\bar{s}}}_{t+1}\right\|^{2}\\\
=&\Big{\|}{{\texttt{FastMix}}}({\bf{s}}_{t}+{\bf{g}}_{t+1}-{\bf{g}}_{t},K)-\frac{1}{m}{\bf{1}}{\bf{1}}^{\top}{{\texttt{FastMix}}}({\bf{s}}_{t}+{\bf{g}}_{t+1}-{\bf{g}}_{t},K)\Big{\|}^{2}\\\
\leq&\rho_{t}^{2}\Big{\|}{\bf{s}}_{t}+{\bf{g}}_{t+1}-{\bf{g}}_{t}-\frac{1}{m}{\bf{1}}{\bf{1}}^{\top}({\bf{s}}_{t}+{\bf{g}}_{t+1}-{\bf{g}}_{t})\Big{\|}^{2}\\\
\leq&2\rho_{t}^{2}\left\|{\bf{s}}_{t}-{\bf{1}}{{\bar{s}}}_{t}\right\|^{2}+2\rho_{t}^{2}\Big{\|}{\bf{g}}_{t+1}-{\bf{g}}_{t}-\frac{1}{m}{\bf{1}}{\bf{1}}^{\top}({\bf{g}}_{t+1}-{\bf{g}}_{t})\Big{\|}^{2}\\\
\leq&2\rho_{t}^{2}\left\|{\bf{s}}_{t}-{\bf{1}}{{\bar{s}}}_{t}\right\|^{2}+2\rho_{t}^{2}\left\|{\bf{g}}_{t+1}-{\bf{g}}_{t}\right\|^{2},\end{split}$
(14)

where the second inequality is based on triangle inequality and the last step
uses Lemma 3. Combining the results of (13) and (14), we have

$\displaystyle\begin{split}&{\mathbb{E}}\left\|{\bf{s}}_{t+1}-{\bf{1}}{{\bar{s}}}_{t+1}\right\|^{2}\\\
\leq&{\mathbb{E}}\Big{[}2\rho_{t}^{2}\left\|{\bf{s}}_{t}-{\bf{1}}{{\bar{s}}}_{t}\right\|^{2}+4\rho^{2}p\left\|\nabla{\bf{f}}({\bf{x}}_{t})-{\bf{g}}_{t}\right\|^{2}\\\
&+18\Big{(}\frac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\Big{)}L^{2}\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}+12\Big{(}\frac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\Big{)}L^{2}\eta^{2}\left\|{\bf{s}}_{t}-{\bf{1}}{{\bar{s}}}_{t}\right\|^{2}\\\
&+6\Big{(}\frac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\Big{)}L^{2}m\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}\Big{]}\\\
\leq&{\mathbb{E}}\Big{[}18\Big{(}\frac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\Big{)}L^{2}\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}+2\Big{(}\frac{6\hat{\rho}^{2}(1-p)}{b}+12\rho^{2}p+\rho_{t}^{2}\Big{)}\left\|{\bf{s}}_{t}-{\bf{1}}{{\bar{s}}}_{t}\right\|^{2}\\\
&+4\rho^{2}p\left\|\nabla{\bf{f}}({\bf{x}}_{t})-{\bf{g}}_{t}\right\|^{2}+6\Big{(}\frac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\Big{)}L^{2}m\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}\Big{]}\\\
\leq&{\mathbb{E}}\Big{[}18\Big{(}\frac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\Big{)}L^{2}\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}+38\hat{\rho}^{2}\left\|{\bf{s}}_{t}-{\bf{1}}{{\bar{s}}}_{t}\right\|^{2}\\\
&+4\rho^{2}p{\mathbb{E}}\left\|\nabla{\bf{f}}({\bf{x}}_{t})-{\bf{g}}_{t}\right\|^{2}+6\Big{(}\frac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\Big{)}L^{2}m\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}\Big{]}.\end{split}$
(15)

We introduce the vector
$z_{t}=[\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2};\eta^{2}\left\|{\bf{s}}_{t}-{\bf{1}}{{\bar{s}}}_{t}\right\|^{2}]$.
Connecting inequalities (11), (15) and the definitions of $\rho$ and
$\hat{\rho}$, we achieve

$\displaystyle\begin{split}{\mathbb{E}}[z_{t+1}]\leq&2\hat{\rho}^{2}\begin{bmatrix}1&1\\\\[2.84544pt]
9\left(\dfrac{1-p}{b}+2p\right)&25\end{bmatrix}{\mathbb{E}}[z_{t}]\\\
&+{\mathbb{E}}\begin{bmatrix}0\\\\[2.84544pt]
4\rho^{2}p\eta^{2}{\mathbb{E}}\left\|\nabla{\bf{f}}({\bf{x}}_{t})-{\bf{g}}_{t}\right\|^{2}+6\left(\dfrac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\right)m\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}\end{bmatrix}\\\
\leq&\frac{1}{4}{\mathbb{E}}[z_{t}]+{\mathbb{E}}\begin{bmatrix}0\\\\[2.84544pt]
4\rho^{2}p\eta^{2}{\mathbb{E}}\left\|\nabla{\bf{f}}({\bf{x}}_{t})-{\bf{g}}_{t}\right\|^{2}+6\left(\dfrac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\right)m\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}\end{bmatrix},\end{split}$

where we use $9\left(\dfrac{1-p}{b}+2p\right)\leq 27$,
$\left\|\begin{bmatrix}1&1\\\ 27&25\end{bmatrix}\right\|_{2}\leq 37$ and (10).
Hence, we have

$\displaystyle C_{t+1}\leq$
$\displaystyle{\mathbb{E}}\Big{[}\frac{1}{4}\big{(}\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}+\eta^{2}\left\|{\bf{s}}_{t}-{\bf{1}}{{\bar{s}}}_{t}\right\|^{2}\big{)}$
$\displaystyle+4\rho^{2}p\eta^{2}{\mathbb{E}}\left\|\nabla{\bf{f}}({\bf{x}}_{t})-{\bf{g}}_{t}\right\|^{2}+6\Big{(}\dfrac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\Big{)}m{\mathbb{E}}\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}\Big{]}$
$\displaystyle=$
$\displaystyle{\mathbb{E}}\Big{[}\frac{1}{4}C_{t}+4\rho^{2}p\eta^{2}mV_{t}+6\Big{(}\dfrac{\hat{\rho}^{2}(1-p)}{b}+2\rho^{2}p\Big{)}m\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}\Big{]}.$

∎

Now we provide the recursion for local and global error of gradient
estimation.

######  Lemma 7.

Under the setting of Theorem 1, we have

$\displaystyle
V_{t+1}\leq{\mathbb{E}}\left[(1-p)V_{t}+\dfrac{9(1-p)L^{2}}{mb}C_{t}+\dfrac{3(1-p)L^{2}}{b}\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}\right].$
(16)

###### Proof.

The update of ${\bf{g}}_{t+1}(i)$ means

$\displaystyle\begin{split}&{\mathbb{E}}\left\|{\bf{g}}_{t+1}(i)-\nabla
f_{i}({\bf{x}}_{t+1}(i))\right\|^{2}\\\ =&p{\mathbb{E}}\left\|\nabla
f_{i}({\bf{x}}_{t+1}(i))-\nabla f_{i}({\bf{x}}_{t+1}(i))\right\|^{2}\\\
&+(1-p){\mathbb{E}}\left\|{\bf{g}}_{t}(i)+\dfrac{1}{b}\sum_{j=1}^{b}\big{(}\nabla
f_{i,\xi_{j}}({\bf{x}}_{t+1}(i))-\nabla
f_{i,\xi_{j}}({\bf{x}}_{t}(i))\big{)}-\nabla
f_{i}({\bf{x}}_{t+1}(i))\right\|^{2}\\\
=&(1-p){\mathbb{E}}\left\|{\bf{g}}_{t}(i)-\nabla
f_{i}({\bf{x}}_{t}(i))\right\|^{2}\\\
&+(1-p){\mathbb{E}}\left\|\dfrac{1}{b}\sum_{j=1}^{b}\big{(}\nabla
f_{i,\xi_{j}}({\bf{x}}_{t+1}(i))-\nabla f_{i,\xi_{j}}({\bf{x}}_{t}(i))-\nabla
f_{i}({\bf{x}}_{t+1}(i))+\nabla f_{i}({\bf{x}}_{t}(i))\big{)}\right\|^{2}\\\
\leq&(1-p){\mathbb{E}}\left\|{\bf{g}}_{t}(i)-\nabla
f_{i}({\bf{x}}_{t}(i))\right\|^{2}+\dfrac{1-p}{b}{\mathbb{E}}\left\|\nabla
f_{i,\xi_{j}}({\bf{x}}_{t+1}(i))-\nabla
f_{i,\xi_{j}}({\bf{x}}_{t}(i))\right\|^{2}\\\
\leq&(1-p){\mathbb{E}}\left\|{\bf{g}}_{t}(i)-\nabla
f_{i}({\bf{x}}_{t}(i))\right\|^{2}+\dfrac{(1-p)L^{2}}{b}{\mathbb{E}}\left\|{\bf{x}}_{t+1}(i)-{\bf{x}}_{t}(i)\right\|^{2},\end{split}$

where the first equality is based on the update rule; the second equality uses
the property of Martingale [11, Proposition 1]; the first inequality use the
property of variance and independence of $\xi_{1},\dots,\xi_{b}$; the last
step is based on Assumption 3.

Taking the average over on above result over $i=1,\dots,m$, we obtain

$\displaystyle\begin{split}&{\mathbb{E}}[V_{t+1}]\\\
=&\frac{1}{m}{\mathbb{E}}\left\|{\bf{g}}_{t+1}-\nabla{\bf{f}}({\bf{x}}_{t+1})\right\|^{2}\\\
\leq&\frac{1-p}{m}{\mathbb{E}}\Big{[}\left\|{\bf{g}}_{t}-\nabla{\bf{f}}({\bf{x}}_{t})\right\|^{2}+\dfrac{(1-p)L^{2}}{mb}\left\|{\bf{x}}_{t+1}-{\bf{x}}_{t}\right\|^{2}\Big{]}\\\
\leq&\frac{1-p}{m}{\mathbb{E}}\left\|{\bf{g}}_{t}-\nabla{\bf{f}}({\bf{x}}_{t})\right\|^{2}\\\
&+\dfrac{3(1-p)L^{2}}{mb}{\mathbb{E}}\Big{[}\left\|{\bf{x}}_{t+1}-{\bf{1}}{{\bar{x}}}_{t+1}\right\|^{2}+{\mathbb{E}}\left\|{\bf{1}}{{\bar{x}}}_{t+1}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}+{\mathbb{E}}\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}\Big{]}\\\
\leq&(1-p){\mathbb{E}}\Big{[}V_{t}+\dfrac{9(1-p)L^{2}}{mb}C_{t}+\dfrac{3(1-p)L^{2}}{b}\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}\Big{]},\end{split}$

where the last step use inequality (11). ∎

######  Lemma 8.

Under the setting of Theorem 1, we have

$\displaystyle{\mathbb{E}}[U_{t+1}]\leq{\mathbb{E}}\left[(1-p)U_{t}+\dfrac{9(1-p)L^{2}}{m^{2}b}C_{t}+\dfrac{3(1-p)L^{2}}{mb}\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}\right].$

###### Proof.

The update of ${\bf{g}}_{t+1}(i)$ means

$\displaystyle{\mathbb{E}}[U_{t+1}]$ $\displaystyle=$ $\displaystyle
p{\mathbb{E}}\left\|\frac{1}{m}\sum_{i=1}^{m}\left(\nabla
f_{i}({\bf{x}}_{t+1}(i))-\nabla f_{i}({\bf{x}}_{t+1}(i))\right)\right\|^{2}$
$\displaystyle+(1-p){\mathbb{E}}\left\|\frac{1}{m}\sum_{i=1}^{m}\left({\bf{g}}_{t}(i)+\dfrac{1}{b}\sum_{j=1}^{b}\big{(}\nabla
f_{i,\xi_{j}}({\bf{x}}_{t+1}(i))-\nabla
f_{i,\xi_{j}}({\bf{x}}_{t}(i))\big{)}-\nabla
f_{i}({\bf{x}}_{t+1}(i))\right)\right\|^{2}$ $\displaystyle=$
$\displaystyle(1-p){\mathbb{E}}\left\|\frac{1}{m}\sum_{i=1}^{m}\left({\bf{g}}_{t}(i)-\nabla
f_{i}({\bf{x}}_{t}(i))\right)\right\|^{2}$
$\displaystyle+(1-p){\mathbb{E}}\left\|\dfrac{1}{mb}\sum_{i=1}^{m}\sum_{j=1}^{b}\big{(}\nabla
f_{i,\xi_{j}}({\bf{x}}_{t+1}(i))-\nabla f_{i,\xi_{j}}({\bf{x}}_{t}(i))-\nabla
f_{i}({\bf{x}}_{t+1}(i))+\nabla f_{i}({\bf{x}}_{t}(i))\big{)}\right\|^{2}$
$\displaystyle\leq$
$\displaystyle(1-p){\mathbb{E}}[U_{t}]+\dfrac{1-p}{m^{2}b^{2}}\sum_{i=1}^{m}\sum_{j=1}^{b}{\mathbb{E}}\left\|\nabla
f_{i,\xi_{j}}({\bf{x}}_{t+1}(i))-\nabla
f_{i,\xi_{j}}({\bf{x}}_{t}(i))\right\|^{2}$ $\displaystyle\leq$
$\displaystyle(1-p){\mathbb{E}}[U_{t}]+\dfrac{(1-p)L^{2}}{m^{2}b^{2}}\sum_{i=1}^{m}\sum_{j=1}^{b}{\mathbb{E}}\left\|{\bf{x}}_{t+1}(i)-{\bf{x}}_{t}(i)\right\|^{2}$
$\displaystyle=$
$\displaystyle(1-p){\mathbb{E}}[U_{t}]+\dfrac{(1-p)L^{2}}{m^{2}b}{\mathbb{E}}\left\|{\bf{x}}_{t+1}-{\bf{x}}_{t}\right\|^{2}$
$\displaystyle\leq$
$\displaystyle(1-p){\mathbb{E}}[U_{t}]+\dfrac{3(1-p)L^{2}}{m^{2}b}{\mathbb{E}}\big{[}\left\|{\bf{x}}_{t+1}-{\bf{1}}{{\bar{x}}}_{t+1}\right\|^{2}+\left\|{\bf{1}}{{\bar{x}}}_{t+1}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}+\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}\big{]}$
$\displaystyle\leq$
$\displaystyle(1-p){\mathbb{E}}\Big{[}U_{t}+\dfrac{6\rho^{2}(1-p)L^{2}}{m^{2}b}C_{t}+\dfrac{3(1-p)L^{2}}{mb}\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}+\dfrac{3(1-p)L^{2}}{m^{2}b}\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}\Big{]}$
$\displaystyle\leq$
$\displaystyle(1-p){\mathbb{E}}\Big{[}U_{t}+\dfrac{9(1-p)L^{2}}{m^{2}b}C_{t}+\dfrac{3(1-p)L^{2}}{mb}\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}\Big{]},$

where the second equality use the property of Martingale [11, Proposition 1];
the first inequality is based on the property of variance; the second
inequality is based on Assumption 3; the third inequality use Young’s
inequality; the last two steps use the upper bound (11) and the definition of
$C_{t}$. ∎

Finally, we give the proof of Theorem 1.

###### Proof.

Combining results of Lemma 5, 6, 7 and 8, we have

$\displaystyle\begin{split}&{\mathbb{E}}[\Phi_{t+1}]\\\
=&{\mathbb{E}}\left[f({{\bar{x}}}_{t+1})+\frac{\eta}{p}U_{t+1}+\frac{\eta}{mp}V_{t+1}+\frac{1}{m\eta}C_{t+1}\right]\\\
\leq&{\mathbb{E}}\Big{[}f({{\bar{x}}}_{t})-\frac{\eta}{2}\left\|\nabla
f({{\bar{x}}}_{t})\right\|^{2}+\eta
U_{t}+\frac{L^{2}\eta}{m}C_{t}-\left(\frac{1}{2\eta}-\frac{L}{2}\right)\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}\\\
&\quad+\frac{\eta(1-p)}{p}U_{t}+\dfrac{9(1-p)\eta
L^{2}}{m^{2}bp}C_{t}+\dfrac{3(1-p)L^{2}\eta}{mbp}\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}\\\
&\quad+\frac{(1-p)\eta}{mp}V_{t}+\dfrac{9(1-p)L^{2}\eta}{m^{2}bp}C_{t}+\dfrac{3(1-p)L^{2}\eta}{mbp}\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}\\\
&\quad+\frac{1}{4m\eta}C_{t}+\frac{4\rho^{2}p\eta^{2}}{\eta}V_{t}+\frac{6}{\eta}\left(\dfrac{\hat{\rho}^{2}(1-p)}{b}+3\rho^{2}p\right)\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}\Big{]}\\\
=&{\mathbb{E}}\Bigg{[}f({{\bar{x}}}_{t})-\frac{\eta}{2}\left\|\nabla
f({{\bar{x}}}_{t})\right\|^{2}+\left(\eta+\frac{\eta(1-p)}{p}\right)U_{t}+\left(\frac{(1-p)\eta}{mp}+4\rho^{2}p\eta\right)V_{t}\\\
&\quad+\left(\frac{L^{2}\eta}{m}+\dfrac{9(1-p)\eta
L^{2}}{m^{2}bp}+\dfrac{9(1-p)L^{2}\eta}{m^{2}bp}+\frac{1}{4m\eta}\right)C_{t}\\\
&-\left(\frac{1}{2\eta}-\frac{L}{2}-\dfrac{3(1-p)L^{2}\eta}{mbp}-\dfrac{3(1-p)L^{2}\eta}{mbp}-\frac{6}{\eta}\left(\dfrac{\hat{\rho}^{2}(1-p)}{b}+3\rho^{2}p\right)\right)\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}\Bigg{]}\\\
\leq&{\mathbb{E}}\left[\Phi_{t}-\frac{\eta}{2}\left\|\nabla
f({{\bar{x}}}_{t})\right\|^{2}-\frac{3}{8m\eta}C_{t}\right],\end{split}$

where the last step is based on the settings of $\eta$, $p$, $b$, $K$ and
$\hat{K}$ in the theorem.

Taking the average on above inequality over $t=0,\dots,T-1$, we have

$\displaystyle\begin{split}&{\mathbb{E}}\left[\frac{1}{T}\sum_{t=0}^{T-1}\left\|\nabla
f({{\bar{x}}}_{t})\right\|^{2}\right]\\\ \leq&\frac{2}{\eta
T}\sum_{t=0}^{T-1}{\mathbb{E}}\Bigg{[}\Phi_{t}-\Phi_{t+1}-\frac{3}{8m\eta}C_{t}\Bigg{]}\\\
\leq&\frac{2}{\eta
T}\sum_{t=0}^{T-1}{\mathbb{E}}\Bigg{[}\Phi_{t}-\Phi_{t+1}-\frac{3}{8m\eta}\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}\Bigg{]}\\\
=&\frac{2{\mathbb{E}}[\Phi_{0}-\Phi_{T}]}{\eta
T}-\frac{3}{4m\eta^{2}T}\sum_{t=0}^{T-1}\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}.\end{split}$
(17)

Based on Assumption 2, we have

$\displaystyle\begin{split}&\Phi_{0}-\Phi_{T}\\\
=&f({{\bar{x}}}_{0})+\frac{\eta}{p}U_{0}+\frac{\eta}{p}V_{0}+\frac{1}{m\eta}C_{0}-\left(f({{\bar{x}}}_{T})+\frac{\eta}{p}U_{T}+\frac{\eta}{p}V_{T}+\frac{1}{m\eta}C_{T}\right)\\\
=&f({{\bar{x}}}_{0})-f({{\bar{x}}}_{T})+\frac{\eta}{p}\left\|\frac{1}{m}\sum_{i=1}^{m}\left({\bf{g}}_{0}(i)-\nabla
f_{i}({\bf{x}}_{0}(i))\right)\right\|^{2}+\frac{\eta}{mp}\left\|{\bf{g}}_{0}-\nabla{\bf{f}}({\bf{x}}_{0})\right\|^{2}+\frac{\eta}{m}\left\|{\bf{s}}_{0}-{\bf{1}}{{\bar{s}}}_{0}\right\|^{2}\\\
\leq&f({{\bar{x}}}_{0})-f^{*}+\frac{\eta}{m}\left\|{\bf{s}}_{0}-{\bf{1}}{{\bar{s}}}_{0}\right\|^{2}.\end{split}$
(18)

Connecting the inequalities (17) and (18), we have

$\displaystyle\begin{split}{\mathbb{E}}\left[\frac{1}{T}\sum_{t=0}^{T-1}\left\|\nabla
f({{\bar{x}}}_{t})\right\|^{2}\right]\leq\frac{2(f({{\bar{x}}}_{0})-f^{*})}{\eta
T}+\frac{2}{mT}\left\|{\bf{s}}_{0}-{\bf{1}}{{\bar{s}}}_{0}\right\|^{2}-\frac{2}{m\eta^{2}T}\sum_{t=0}^{T-1}{\mathbb{E}}\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}.\end{split}$
(19)

Hence, the output $x_{\rm out}$ satisfies that

$\displaystyle\begin{split}&{\mathbb{E}}\left\|\nabla f(x_{\rm
out})\right\|^{2}\\\
=&\frac{1}{mT}\sum_{i=1}^{m}\sum_{t=0}^{T-1}{\mathbb{E}}\left\|\nabla
f({\bf{x}}_{t}(i))\right\|^{2}\\\
\leq&\frac{2}{mT}\sum_{i=1}^{m}\sum_{t=0}^{T-1}{\mathbb{E}}\left[\left\|\nabla
f(\bar{x}_{t})\right\|^{2}+\left\|\nabla f({\bf{x}}_{t}(i))-\nabla
f({{\bar{x}}}_{t})\right\|^{2}\right]\\\
\leq&\frac{2}{mT}\sum_{i=1}^{m}\sum_{t=0}^{T-1}{\mathbb{E}}\left[\left\|\nabla
f(\bar{x}_{t})\right\|^{2}+L^{2}\left\|{\bf{x}}_{t}(i)-{{\bar{x}}}_{t}\right\|^{2}\right]\\\
=&\frac{2}{T}\sum_{t=0}^{T-1}{\mathbb{E}}\left\|\nabla
f(\bar{x}_{t})\right\|^{2}+\frac{2L^{2}}{mT}\sum_{t=0}^{T-1}{\mathbb{E}}\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}\\\
\leq&\frac{4(f({{\bar{x}}}_{0})-f^{*})}{\eta
T}+\frac{4}{mT}\left\|{\bf{s}}_{0}-{\bf{1}}{{\bar{s}}}_{0}\right\|^{2}-\left(\frac{6L^{2}}{mT}-\frac{2L^{2}}{mT}\right)\left\|{\bf{x}}_{t}-{\bf{1}}{{\bar{x}}}_{t}\right\|^{2}\\\
\leq&\frac{4(f({{\bar{x}}}_{0})-f^{*})}{\eta
T}+\frac{4}{m}\left\|{\bf{s}}_{0}-{\bf{1}}{{\bar{s}}}_{0}\right\|^{2}\leq{\varepsilon}^{2},\end{split}$
(20)

where the first step uses Young’s inequality; the second inequality is based
on Assumption 3; the third inequality uses inequality (19); the other steps
are based on the parameter settings. Using Jensen’s inequality, we obtain

$\displaystyle{\mathbb{E}}\left\|\nabla f(x_{\rm
out})\right\|\leq\sqrt{{\mathbb{E}}\left\|\nabla f(x_{\rm
out})\right\|^{2}\,}\leq{\varepsilon}.$

Furthermore, the total expected IFO calls on per-agent is

$\displaystyle{\mathcal{O}}\big{(}n+T(pn+(1-p)b)\big{)}={\mathcal{O}}\left(n+\sqrt{\frac{n}{m}}\,L{\varepsilon}^{-2}\right)$

and the total expected number of communication rounds is

$\displaystyle{\mathcal{O}}\big{(}K_{\rm
in}+T(pK+(1-p)\hat{K})\big{)}={\mathcal{O}}\left(\frac{L{\varepsilon}^{-2}}{\sqrt{1-\lambda_{2}(W)}}\right).$

∎

######  Remark 1.

The terms of $\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}$ are
crucial to making the analysis be concise. In the proof of Theorem 1, we take
the advantage of the negative coefficient before
$\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}$ in (7) to cancel the
term $\left\|{{\bar{x}}}_{t+1}-{{\bar{x}}}_{t}\right\|^{2}$ appears in the
recursion of $C_{t}$, $U_{t}$ and $V_{t}$, which directly leads to the
decrease of $\Phi_{t}$.

##  5 The Optimality

We first consider the IFO complexity. Theorem 1 means the total IFO complexity
of DEAREST for all of $m$ agents is
${\mathcal{O}}\left(mn+\sqrt{mn}\,L{\varepsilon}^{-2}\right)$, which matches
the lower bound of finding ${\varepsilon}$-stationary point in nonconvex
finite-sum optimization with $mn$ individual functions by IFO algorithm on
single machine [11, 19]. Since there are $mn$ individual functions in our
settings (1)-(2), the IFO complexity of DEAREST is optimal.

Then we consider the communication complexity. Lu et al. [22] study the lower
bound for solving decentralized nonconvex problem by using unbiased stochastic
first-order oracle (SFO) with bounded variance. Concretely, they considered
the decentralized optimization problem

$\displaystyle\min_{x\in{\mathbb{R}}^{d}}\hat{f}(x)\triangleq\frac{1}{m}\sum_{i=1}^{m}\hat{f}_{i}(x),$

where $\hat{f}_{i}(x)$ is the $L$-smooth local function on the $i$-th agent
and there are $m$ agents form a network with gossip matrix $\hat{W}$. It is
shown that there exist some $\\{\hat{f}_{i}(x)\\}_{i=1}^{m}$ and $\hat{W}$
such that finding an ${\varepsilon}$-stationary point of $\hat{f}(x)$ by SFOs
$\\{\tilde{g}_{i}(x)\\}$ satisfying

$\displaystyle{\mathbb{E}}[\tilde{g}_{i}(x)]=\nabla\hat{f}_{i}(x)\qquad\text{and}\qquad{\mathbb{E}}\big{\|}\nabla\hat{f}_{i}(x)-\tilde{g}_{i}(x)\big{\|}^{2}\leq\sigma^{2}$
(21)

for all $x\in{\mathbb{R}}^{d}$ and gossip protocol associate with $\hat{W}$
requires at least

$\displaystyle\Omega\left(\frac{L\sigma^{2}{\varepsilon}^{-4}}{m}+\frac{L{\varepsilon}^{-2}}{\sqrt{1-\lambda_{2}(\hat{W})}}\right)$
(22)

rounds of communication [22, Corollary 1].

For the setting of decentralized finite-sum optimization in this paper,
consider the formulations (1)-(2) with $f_{ij}(x)\triangleq\hat{f}_{i}(x)$ and
network topology with gossip matrix $\hat{W}$, where $\hat{f}_{i}(x)$ and
$\hat{W}$ correspond to the ones we mentioned in last paragraph. We let
$\sigma=0$ for (21), then the SFOs $\\{\tilde{g}_{i}(x)\\}$ satisfying (21)
are equivalent to our IFOs $\\{\nabla f_{ij}(x)\\}$. Hence, the lower bound
(22) directly means there exists some $\\{\hat{f}_{i}(x)\\}_{i=1}^{m}$ and
$\hat{W}$ such that using IFO algorithm to find an ${\varepsilon}$-stationary
point for finite-sum problem (1) with
$f_{i}(x)\triangleq\frac{1}{n}\sum_{j=1}^{n}f_{ij}(x)$ and
$f_{ij}(x)\triangleq\hat{f}_{i}(x)$ requires at least

$\displaystyle\Omega\left(\frac{L{\varepsilon}^{-2}}{\sqrt{1-\lambda_{2}\big{(}{\hat{W}}\big{)}}}\right)$

rounds of communication, which matches the upper bound of DEAREST shown in
Theorem 1.

##  6 Related Work and Discussion

We summarize the results of our algorithm and related work in Table 1. The
proposed algorithm DEAREST is the only one with both optimal IFO complexity
and optimal communication complexity. D-GET [32] and GT-SARAH [37] are also
based on variance reduction and gradient tracking but they only take single
consensus step when communication takes place. DESTRESS includes the multi-
consensus steps like DEAREST while its procedure has double-loop structure,
leading to the unnecessary complicated implementation and more communication
cost. In contrast, our DEAREST enjoys the single-loop framework, which results
to the easier implementation, the simpler convergence analysis and the optimal
IFO/communication complexity.

##  7 Numerical Experiments

|  |
---|---|---
(a) a9a | (b) w8a | (c) rcv1
Figure 1: The comparison for the number of IFO calls vs. $\left\|\nabla f(x)\right\|$ on random graph. |  |
---|---|---
(a) a9a | (b) w8a | (c) rcv1
Figure 2: The comparison for the number of communication rounds vs.
$\left\|\nabla f(x)\right\|$ on random graph.

In this section, we provide the empirical study on regularized logistic
regression for binary classification [3] The individual function for the model
is defined as

$\displaystyle
f_{ij}(x)\triangleq\log\left(1+\exp(-b_{ij}a_{ij}^{\top}x)\right)+\lambda\sum_{k=1}^{d}\frac{x_{k}^{2}}{1+x_{k}^{2}},$

where $a_{j}\in{\mathbb{R}}^{d}$ and $b_{j}\in\\{-1,1\\}$ are the feature and
label of the $j$-th sample on the $i$-th agent; $\lambda>0$ is the
regularization parameter; $x=[x_{1},\dots,x_{d}]^{\top}\in{\mathbb{R}}^{d}$ is
the vector of the classifier.

We compare the proposed DEAREST with baseline algorithms GT-SARAH [37] and
DESTRESS [18] on three real-world datasets “a9a” ($mn=32,560$, $d=123$), “w8a”
($mn=49,740$, $d=300$) and “rcv1” ($mn=20,420$, $d=47,236$). All of the
datasets can be downloaded from LIBSVM repository [6]. We set the number of
agents $m=20$ and split the dataset uniformly at random to all of the agents.
We conduct the experiments on two communication graphs: a) the random graph
where each pair of agents has a connection with the probability of $0.15$; b)
the circle graph. We set $W=I-{L}/{\lambda_{1}(L)}$, where $L$ is the
Laplacian matrix associated with the adjacency matrix and $\lambda_{1}(L)$ is
the largest eigenvalue of $L$. We observe that $1-\lambda_{2}(W)=0.0482$ for
the random graph and and $1-\lambda_{2}(W)=0.0245$ for the circle graph. We
set $\lambda=10^{-4}$ in experiments and the parameters of all algorithms are
well-tuned.

|  |
---|---|---
(a) a9a | (b) w8a | (c) rcv1
Figure 3: The comparison for the number of IFO calls vs. $\left\|\nabla f(x)\right\|$ on circle graph. |  |
---|---|---
(a) a9a | (b) w8a | (c) rcv1
Figure 4: The comparison for the number of communication rounds vs.
$\left\|\nabla f(x)\right\|$ on circle graph.

We report the results on random graph and circle graph in Figure 1-2 and
Figure 3-4 respectively. We observe that DEAREST outperforms the other two
algorithms is all of the settings. Especially, the communication cost of
DEAREST is apparently less than the ones of baselines. This phenomenon matches
our theoretical result that the upper bound of communication complexity for
DEAREST is tighter than the ones of GT-SARAH and DESTRESS. The results also
show that both DEAREST and DESTRESS achieve better performance than GT-SARAH.
The reason is GT-SARAH lacks the multi-consensus step, which leads to more
computation cost and communication cost when $1-\lambda_{2}(W)$ is close to
zero. Additionally, DEAREST and DESTRESS has the similar computational cost in
most of the cases, which also validates the theoretical results that both of
these algorithms has the optimal IFO complexity.

##  8 Conclusion

In this paper, we have proposed a novel stochastic algorithm DEcentralized
probAbilistic Recursive gradiEnt deScenT (DEAREST) for decentralized nonconvex
finite-sum optimization. The algorithm has single-loop structure and easy to
be implemented. We have proved DEAREST achieves both optimal incremental
first-order oracle (IFO) complexity and optimal communication complexity. Our
convergence analysis depends on constructing the Lyapunov function consists of
function value, global/local gradient estimation error and consensus error,
which is intuitive and easy to follow. The results of numerical experiments
also support our theoretical analysis. We believe the idea of DEAREST is not
limited to the IFO algorithm for nonconvex optimization and it also can be
extended to solve other types of decentralized optimization problems.

## References

  * Allen-Zhu [2017] Zeyuan Allen-Zhu.  Katyusha: The first direct acceleration of stochastic gradient methods.  _Journal of Machine Learning Research_ , 18(1):8194–8244, 2017.
  * Allen-Zhu and Hazan [2016] Zeyuan Allen-Zhu and Elad Hazan.  Variance reduction for faster non-convex optimization.  In _ICML_ , 2016.
  * Antoniadis et al. [2011] Anestis Antoniadis, Irène Gijbels, and Mila Nikolova.  Penalized likelihood regression for generalized linear models with non-quadratic penalties.  _Annals of the Institute of Statistical Mathematics_ , 63(3):585–615, 2011.
  * Arioli and Scott [2014] Mario Arioli and Jennifer Scott.  Chebyshev acceleration of iterative refinement.  _Numerical Algorithms_ , 66(3):591–608, 2014\.
  * Bottou et al. [2018] Léon Bottou, Frank E. Curtis, and Jorge Nocedal.  Optimization methods for large-scale machine learning.  _Siam Review_ , 60(2):223–311, 2018.
  * Chang and Lin [2011] Chih-Chung Chang and Chih-Jen Lin.  LIBSVM: a library for support vector machines.  _ACM transactions on intelligent systems and technology_ , 2(3):1–27, 2011.  URL https://www.csie.ntu.edu.tw/~cjlin/libsvm/.
  * Chen and Sayed [2012] Jianshu Chen and Ali H. Sayed.  Diffusion adaptation strategies for distributed optimization and learning over networks.  _IEEE Transactions on Signal Processing_ , 60(8):4289–4305, 2012.
  * Chen and Sayed [2015] Jianshu Chen and Ali H. Sayed.  On the learning behavior of adaptive networks—part I: Transient analysis.  _IEEE Transactions on Information Theory_ , 61(6):3487–3517, 2015.
  * Defazio et al. [2014] Aaron Defazio, Francis Bach, and Simon Lacoste-Julien.  SAGA: A fast incremental gradient method with support for non-strongly convex composite objectives.  In _NIPS_ , 2014.
  * Facchinei et al. [2015] Francisco Facchinei, Gesualdo Scutari, and Simone Sagratella.  Parallel selective algorithms for nonconvex big data optimization.  _IEEE Transactions on Signal Processing_ , 63(7):1874–1889, 2015.
  * Fang et al. [2018] Cong Fang, Chris Junchi Li, Zhouchen Lin, and Tong Zhang.  SPIDER: Near-optimal non-convex optimization via stochastic path-integrated differential estimator.  In _NeurIPS_ , 2018.
  * Hofmann et al. [2015] Thomas Hofmann, Aurelien Lucchi, Simon Lacoste-Julien, and Brian McWilliams.  Variance reduced stochastic gradient descent with neighbors.  In _NIPS_ , 2015.
  * Johnson and Zhang [2013] Rie Johnson and Tong Zhang.  Accelerating stochastic gradient descent using predictive variance reduction.  In _NIPS_ , 2013.
  * Kovalev et al. [2020] Dmitry Kovalev, Samuel Horváth, and Peter Richtárik.  Don’t jump through hoops and remove those loops: SVRG and Katyusha are better without the outer loop.  In _ALT_ , 2020.
  * Lan and Zhou [2018] Guanghui Lan and Yi Zhou.  An optimal randomized incremental gradient method.  _Mathematical programming_ , 171(1):167–215, 2018\.
  * Lei et al. [2017] Lihua Lei, Cheng Ju, Jianbo Chen, and Michael I. Jordan.  Non-convex finite-sum optimization via SCSG methods.  In _NIPS_ , 2017.
  * Li et al. [2020] Boyue Li, Shicong Cen, Yuxin Chen, and Yuejie Chi.  Communication-efficient distributed optimization in networks with gradient tracking and variance reduction.  _Journal of Machine Learning Research_ , 21:1–51, 2020\.
  * Li et al. [2022] Boyue Li, Zhize Li, and Yuejie Chi.  DESTRESS: Computation-optimal and communication-efficient decentralized nonconvex finite-sum optimization.  _SIAM Journal on Mathematics of Data Science_ , 4(3):1031–1051, 2022.
  * Li et al. [2021] Zhize Li, Hongyan Bao, Xiangliang Zhang, and Peter Richtárik.  PAGE: A simple and optimal probabilistic gradient estimator for nonconvex optimization.  In _ICML_ , 2021.
  * Lian et al. [2017] Xiangru Lian, Ce Zhang, Huan Zhang, Cho-Jui Hsieh, Wei Zhang, and Ji Liu.  Can decentralized algorithms outperform centralized algorithms? a case study for decentralized parallel stochastic gradient descent.  _NIPS_ , 2017.
  * Liu and Morse [2011] Ji Liu and A. Stephen Morse.  Accelerated linear iterations for distributed averaging.  _Annual Reviews in Control_ , 35(2):160–165, 2011\.
  * Lu and De Sa [2021] Yucheng Lu and Christopher De Sa.  Optimal complexity in decentralized training.  In _ICML_ , 2021.
  * Nedic et al. [2017] Angelia Nedic, Alex Olshevsky, and Wei Shi.  Achieving geometric convergence for distributed optimization over time-varying graphs.  _SIAM Journal on Optimization_ , 27(4):2597–2633, 2017.
  * Nedić et al. [2018] Angelia Nedić, Alex Olshevsky, and Michael G. Rabbat.  Network topology and communication-computation tradeoffs in decentralized optimization.  _Proceedings of the IEEE_ , 106(5):953–976, 2018\.
  * Nguyen et al. [2017] Lam M. Nguyen, Jie Liu, Katya Scheinberg, and Martin Takáč.  SARAH: A novel method for machine learning problems using stochastic recursive gradient.  In _ICML_ , 2017.
  * Pham et al. [2020] Nhan H Pham, Lam M Nguyen, Dzung T Phan, and Quoc Tran-Dinh.  ProxSARAH: An efficient algorithmic framework for stochastic composite nonconvex optimization.  _Journal of Machine Learning Research_ , 21(110):1–48, 2020.
  * Pu and Nedić [2021] Shi Pu and Angelia Nedić.  Distributed stochastic gradient tracking methods.  _Mathematical Programming_ , 187(1):409–457, 2021\.
  * Qu and Li [2019] Guannan Qu and Na Li.  Accelerated distributed nesterov gradient descent.  _IEEE Transactions on Automatic Control_ , 65(6):2566–2581, 2019.
  * Reddi et al. [2016] Sashank J. Reddi, Ahmed Hefny, Suvrit Sra, Barnabas Poczos, and Alex Smola.  Stochastic variance reduction for nonconvex optimization.  In _ICML_ , 2016.
  * Schmidt et al. [2017] Mark Schmidt, Nicolas Le Roux, and Francis Bach.  Minimizing finite sums with the stochastic average gradient.  _Mathematical Programming_ , 162(1-2):83–112, 2017.
  * Sun and Hong [2019] Haoran Sun and Mingyi Hong.  Distributed non-convex first-order optimization and information processing: Lower complexity bounds and rate optimal algorithms.  _IEEE Transactions on Signal processing_ , 67(22):5912–5928, 2019.
  * Sun et al. [2020] Haoran Sun, Songtao Lu, and Mingyi Hong.  Improving the sample and communication complexity for decentralized non-convex optimization: Joint gradient estimation and tracking.  In _ICML_ , 2020.
  * Sun et al. [2016] Ju Sun, Qing Qu, and John Wright.  Complete dictionary recovery over the sphere I: Overview and the geometric picture.  _IEEE Transactions on Information Theory_ , 63(2):853–884, 2016.
  * Sun et al. [2018] Ju Sun, Qing Qu, and John Wright.  A geometric analysis of phase retrieval.  _Foundations of Computational Mathematics_ , 18(5):1131–1198, 2018.
  * Sundhar Ram et al. [2010] S. Sundhar Ram, Angelia Nedić, and Venugopal V. Veeravalli.  Distributed stochastic subgradient projection algorithms for convex optimization.  _Journal of optimization theory and applications_ , 147(3):516–545, 2010.
  * Wai et al. [2017] Hoi-To Wai, Jean Lafond, Anna Scaglione, and Eric Moulines.  Decentralized Frank–Wolfe algorithm for convex and nonconvex problems.  _IEEE Transactions on Automatic Control_ , 62(11):5522–5537, 2017.
  * Xian et al. [2021] Wenhan Xian, Feihu Huang, Yanfu Zhang, and Heng Huang.  A faster decentralized algorithm for nonconvex minimax problems.  _NeurIPS_ , 2021.
  * Xin et al. [2022] Ran Xin, Usman A. Khan, and Soummya Kar.  Fast decentralized nonconvex finite-sum optimization with recursive variance reduction.  _SIAM Journal on Optimization_ , 32(1):1–28, 2022\.
  * Yang et al. [2019] Tao Yang, Xinlei Yi, Junfeng Wu, Ye Yuan, Di Wu, Ziyang Meng, Yiguang Hong, Hong Wang, Zongli Lin, and Karl H Johansson.  A survey of distributed optimization.  _Annual Reviews in Control_ , 47:278–305, 2019.
  * Ye et al. [2020a] Haishan Ye, Luo Luo, Ziang Zhou, and Tong Zhang.  Multi-consensus decentralized accelerated gradient descent.  _arXiv preprint arXiv:2005.00797_ , 2020a.
  * Ye et al. [2020b] Haishan Ye, Wei Xiong, and Tong Zhang.  PMGT-VR: A decentralized proximal-gradient algorithmic framework with variance reduction.  _arXiv preprint arXiv:2012.15010_ , 2020b.
  * Ye et al. [2020c] Haishan Ye, Ziang Zhou, Luo Luo, and Tong Zhang.  Decentralized accelerated proximal gradient descent.  In _NeurIPS_ , 2020c.
  * Yuan et al. [2016] Kun Yuan, Qing Ling, and Wotao Yin.  On the convergence of decentralized gradient descent.  _SIAM Journal on Optimization_ , 26(3):1835–1854, 2016.
  * Zhang et al. [2013] Lijun Zhang, Mehrdad Mahdavi, and Rong Jin.  Linear convergence with condition number independent access of full gradients.  In _NIPS_ , 2013.
  * Zhou et al. [2018] Dongruo Zhou, Pan Xu, and Quanquan Gu.  Stochastic nested variance reduction for nonconvex optimization.  In _NeurIPS_ , 2018.