# Fast networked data selection via distributed smoothed quantile estimation

Xu Zhang and Marcos M. Vasconcelos X. Zhang was supported by the Postdoctoral
Fellowship Program of CPSF under Grant Number GZC20232038. M. M. Vasconcelos
was partially supported by the Commonwealth Cyber Initiative.X. Zhang is with
the School of Artificial Intelligence, Xidian University, Xi’an, China
(e-mail: zhang.xu@xidian.edu.cn).M. M. Vasconcelos is with the Department of
Electrical Engineering at the FAMU-FSU College of Engineering, Florida State
University, USA (e-mail: mm22eo@fsu.edu).

###### Abstract

Collecting the most informative data from a large dataset distributed over a
network is a fundamental problem in many fields, including control, signal
processing and machine learning. In this paper, we establish a connection
between selecting the most informative data and finding the top-$k$ elements
of a multiset. The top-$k$ selection in a network can be formulated as a
distributed nonsmooth convex optimization problem known as quantile
estimation. Unfortunately, the lack of smoothness in the local objective
functions leads to extremely slow convergence and poor scalability with
respect to the network size. To overcome the deficiency, we propose an
accelerated method that employs smoothing techniques. Leveraging the piecewise
linearity of the local objective functions in quantile estimation, we
characterize the iteration complexity required to achieve top-$k$ selection, a
challenging task due to the lack of strong convexity. Several numerical
results are provided to validate the effectiveness of the algorithm and the
correctness of the theory.

##  I Introduction

Multi-agent networks have been widely used to model many applications such as
robotic, sensor and social networks, as well as client-server architectures
for distributed machine learning. With inexpensive sensing devices and storage
now readily available, there has been an exponential increase in data
generation, leading to the production of vast amounts of data. However, data
processing and wireless communication consume much more power than sensing.
Therefore, selecting and transmitting only the most valuable information from
a potentially very large collection of random data becomes a fundamental
problem in many control, signal processing and machine learning applications
[1]. In this paper, we relate this problem to choosing the largest entries in
a multiset. While such a selection problem can be easily solved through
sorting in centralized systems, a significant challenge arises when dealing
with decentralized systems where agents are locally connected over a network
through peer-to-peer communication links (see Fig. 1).

Figure 1: System architecture for the top-$k$ distributed sensor selection
problem, where a multi-robot network employs sensors to gather observations,
but only the most informative top-$k$ data can be relayed to the remote
station via wireless links.

We address the following problem setting: a dataset is distributed among a
potentially large network of $n$ agents interconnected by a local and
incomplete communication network, where agents can only communicate with
neighbors through peer-to-peer (P2P) links. Each agent computes an
informativeness score for each data point in its local dataset. For example,
in sensor selection based on linear and quadratic estimation models [2, 3],
informativeness scores are functions of the observation matrices and vectors
(c.f. Section II); in the context of polling clients in Federated Learning,
informativeness scores can be a criterion for selection, such as choosing
clients with the top-$k$ scores based on the required data upload time to the
server [4], the gradient norms of the local loss functions[5], or the current
value of the local loss functions [6]. The objective of the agents is to
collaborate locally to identify which nodes in the network possess one or more
of the top-$k$ most informativeness scores. Subsequently, the selected agents
transmit the relevant data points over long-range wireless links to a remote
central station for further processing and decision-making.

###  I-A Motivation

We would like to solve the top-$k$ selection problem as fast and reliably as
possible, and develop a scalable algorithm that accommodates a growing number
of agents in the network. Despite its apparent simplicity, identifying the
top-$k$ largest numbers across a network of interconnected agents poses
significant challenges. Numerous studies explore distributed settings in which
agents iteratively communicate with a central server to collectively determine
the top-$k$ highest informativeness scores [7, 8, 9, 10, 11]. However, relying
on a server compromises the robustness and scalability of distributed
networks. Moreover, these works entail actual data transmission and require
significant storage and communication expenses. To address these limitations,
distributed sub-gradient methods [12, 13, 14] have been proposed to achieve
top-$k$ selection by formulating the problem as a quantile estimation problem.
These methods are fully distributed and require smaller storage and
communication costs. Nevertheless, they encounter challenges such as slow
convergence and lack of analysis regarding iteration complexity. These issues
stem from the non-smoothness and non-strong convexity of the objective
function. The former renders the application of momentum methods for
acceleration, while the latter complicates the analysis of iteration
complexity for the iterative variable.

###  I-B Contributions

To address the above challenges, this paper proposes an efficient and scalable
algorithm for distributed top-$k$ selection based on smoothed distributed
quantile estimation. The main contributions are listed below.

  1. 1.

By applying smoothing techniques to the local objective functions [15, 16], we
propose an efficient and scalable distributed top-$k$ selection method. This
approach is complemented by integrating it with one of the advanced
distributed smooth optimization methods, namely the EXTRA algorithm [17, 18].
Each agent is required to store only two units of memory and transmit just one
unit to each agent at every iteration. Moreover, the storage and communication
costs remain independent of the number of agents $n$, and the number $k$ for
each agent at each iteration, thereby presenting a highly desirable feature in
practical applications.

  2. 2.

We characterize the iteration complexity of the distributed top-$k$ selection
problem for our method, whose expression captures the existing trade-offs
between the smoothing parameter, the number $k$, the number of agents, the
graph connectivity and the resolution. To give the iteration complexity, the
connection between the sequence of objective values and the sequence of
iterative variables is established by making full use of the piecewise linear
property of the loss function and the uniqueness of the solution.

  3. 3.

Extensive numerical results substantiate the effectiveness of the proposed
method compared to traditional sub-gradient and vanilla message passing
methods. Before our current work, distributed sub-gradient algorithms for
top-$k$ selection were confined to a small number of nodes. In contrast, our
algorithm can scale to very large-scale scenarios and significantly reduce the
number of iterations.

###  I-C Related Work

The problem of top-$k$ selection has been an active research topic since the
pioneering work of Blum et al. in 1973 [19]. However, it remains a relevant
topic, in which researchers continuously develop algorithms and
implementations that harness GPUs to efficiently perform top-$k$ selection on
immense datasets [20, 21, 22].

We highlight that the terminology top-$k$ may correspond to multiple classes
of problems in the literature. One such class is called top-$k$ selective
gossip [23, 24, 25], where each agent has a vector and engages to reach
consensus on the $k$ largest entries of the average of all initial vectors.
Despite the apparent similarity in name, top-$k$ selective gossip cannot
perform top-$k$ selection, whose task is to identify within the network the
$k$ largest numbers in a multiset. The distributed implementations of top-$k$
selection have been the focus of many fundamental works, including top-$k$
queries [26, 9, 10, 27, 28], top-$k$ monitoring [29, 30, 11], gossip-based
algorithms [31, 32, 33], and other message passing algorithms[34, 35]. The
work reported herein contributes to the state-of-the-art in the top-$k$
selection in networked datasets.

The overwhelming majority of the works in top-$k$ selection, considers the
case of aggregating data over a spanning tree for the underlying network.
Using message passing algorithms, the top-$k$ data can be consolidated at the
root and collected or transmitted to a server. This approach has three
drawbacks: (1) It requires the construction or availability of a spanning
tree, which may not be available, or the generation may be difficult in very
large-scale ad-hoc networks; (2) It requires that the actual data are
communicated between nodes in the network, and finally aggregated at the root
of the spanning tree, which violates privacy and poses security threats in
sensitive applications; (3) It requires communication between any two agents
in the network to be noiseless, which is an unrealistic assumption as most
modern networks are wireless, and therefore are subject to several forms of
communication imperfections such as packet drops and additive noise [36, 37,
38, 39, 40, 41].

Our work uses a different approach to solve the top-$k$ selection problem
based on distributed convex optimization. A classic result from statistics is
that selecting the top-$k$ numbers in a dataset is equivalent to estimating
the quantile by minimizing the corresponding pinball loss [42], which has
recently been used to solve the distributed top-$k$ problem [12]. In [13], the
authors noticed that quantile estimation could be solved using distributed
optimization. The combination of these results was used in the context of
distributed estimation in [14], where the connection between the top-$k$ and
quantile estimation was formalized, and solved by using a standard distributed
subgradient method [43]. One of the big advantages of distributed optimization
is that it works under minimal assumptions on the communication network and
does not require the construction of a spanning tree. Additionally, there are
simple ways of modifying distributed convex optimization algorithms by
introducing an extra time scale to handle several types of communication
imperfections [44, 45, 46]. Finally, since our goal is to compute the $k$-th
quantile, the actual data are never communicated to any other node. The data
remain local, and hence privacy is preserved. One of the drawbacks of current
optimization-based algorithms is that they are slow due to the non-smooth
nature of the objective function, where diminishing step size is required to
guarantee convergence to an exact quantile.

###  I-D Paper Organization

The rest of this paper is organized as follows. Section II introduces two
typical applications in sensor subset selection and formulates them as top-$k$
selection problems. Subsequently, Section III provides the problem setup by
exploring the fundamentals of quantile estimation and modeling top-$k$
selection problem into a distributed quantile estimation problem. Section IV
relates the function error of the smoothed objective function and the variable
error to the solution in the optimal solution interval. Section V presents two
typical smoothing techniques. Section VI proposes an accelerated distributed
quantile algorithm via EXTRA, and provides the iteration complexity.
Simulations in Section VII demonstrate the effectiveness of the proposed
method. Conclusion and future work are provided in Section VIII.

##  II Prelude – Sensor Subset Selection

Sensor subset selection is an NP-hard problem traditionally formulated as an
integer programming problem [47, 2, 48]. Recent work in this area casts the
sensor subset selection as a submodular optimization problem, and by using
greedy algorithms, it is possible to obtain good solutions within reasonable
time [2, 3]. In what follows, we illustrate that the problem of subset
selection in linear models and quadratic models is equivalent to top-$k$
selection under the so-called T-optimality criterion.

###  II-A Subset Selection for Linear Models

Consider a linear model $y_{i}=\bm{a}_{i}^{T}\bm{x}+v_{i}$, where
$\\{y_{i}\\}_{i=1}^{n}\in\mathbb{R}$ are the observations,
$\\{\bm{a_{i}}\\}_{i=1}^{n}\in\mathbb{R}^{m}$ are the observation vectors,
$\bm{x}\in\mathbb{R}^{m}$ is the unknown vector and $\\{v_{i}\\}$ are i.i.d.
Gaussian noise variables satisfying $v_{i}\sim\mathcal{N}(0,\sigma_{i}^{2})$
and $\sigma_{i}$ denotes the standard deviation of $v_{i}$. Consider a
scenario where communication limitations dictate that at most $k$ out of $n$
pairs $\\{(\bm{a}_{i},y_{i})\\}_{i\in\mathcal{S}}$ can be used to estimate the
desired vector $\bm{x}$, where $\mathcal{S}\subset\\{1,\ldots,n\\}$ with
cardinality $|\mathcal{S}|\leq k$. Assume that $\bm{x}$ has a Gaussian prior
distribution, i.e., $\boldsymbol{x}\sim\mathcal{N}(\bm{0},\mathbf{P})$, where
$\mathbf{P}$ denotes its covariance matrix. Following [3], the task of
selecting the top-$k$ most informative data points such as to maximize the
T-optimality criterion is:

$\displaystyle\max_{\mathcal{S}}~{}~{}\text{Tr}(\mathbf{M}_{\mathcal{S}}^{-1})-\text{Tr}(\mathbf{P}^{-1})~{}~{}\text{
s.t. }\ \ |\mathcal{S}|\leq k,$ (1)

where $\mathbf{M}_{\mathcal{S}}$ denotes the error covariance matrix [49]

$\mathbf{M}_{\mathcal{S}}\operatorname{\overset{def}{=}}\left(\mathbf{P}^{-1}+\sum_{i\in\mathcal{S}}\frac{1}{\sigma_{i}^{2}}\bm{a}_{i}\bm{a}_{i}^{T}\right)^{-1}.$
(2)

Incorporating Eq. 2 into Eq. 1, the problem is reformulated as

$\displaystyle\max_{\mathcal{S}}~{}~{}\sum_{i\in\mathcal{S}}\frac{\left\lVert\bm{a}_{i}\right\rVert_{2}^{2}}{\sigma_{i}^{2}}~{}~{}\text{
s.t. }\ \ |\mathcal{S}|\leq k,$ (3)

which is equivalent to finding the top-$k$ scores in $\\{s_{i}\\}_{i=1}^{n}$,
where

$s_{i}\operatorname{\overset{def}{=}}\frac{\left\lVert\bm{a}_{i}\right\rVert_{2}^{2}}{\sigma_{i}^{2}}.$
(4)

###  II-B Subset Selection for Quadratic Models

Consider a quadratic measurement model

$y_{i}=\frac{1}{2}\bm{x}^{T}\mathbf{A}_{i}\bm{x}+\bm{b}_{i}^{T}\bm{x}+v_{i},$
(5)

where $\\{y_{i}\\}_{i=1}^{n}\in\mathbb{R}$ are the observations,
$\\{\mathbf{A}_{i}\\}_{i=1}^{n}\in\mathbb{R}^{m\times m}$ and
$\\{\bm{b}_{i}\\}_{i=1}^{n}\in\mathbb{R}^{m}$ are known observation matrices
and vectors, respectively, $\bm{x}\in\mathbb{R}^{m}$ is the unknown vector to
be estimated and $\\{v_{i}\\}$ are i.i.d. Gaussian noise variables satisfying
$v_{i}\sim\mathcal{N}(0,\sigma_{i}^{2})$. We would like to select at most $k$
out of $n$ tuples $\\{(\mathbf{A}_{i},\bm{b}_{i},y_{i})\\}_{i\in\mathcal{S}}$
to estimate the unknown vector $\bm{x}$. Suppose that
$\bm{x}\sim\mathcal{N}(\bm{0},\mathbf{P})$. Following [3], one way to
formulate this problem is to solve:

$\displaystyle\max_{\mathcal{S}}~{}~{}\text{Tr}(\mathbf{B}_{\mathcal{S}}^{-1})-\text{Tr}(\bm{\Lambda})~{}~{}\text{
s.t. }\ \ |\mathcal{S}|\leq k,$ (6)

where $\mathbf{B}_{\mathcal{S}}$ denotes Bayesian Cramér-Rao lower bound
according to Theorem 2 in [3]

$\mathbf{B}_{\mathcal{S}}=\left(\bm{\Lambda}+\sum_{i\in\mathcal{S}}\frac{1}{\sigma_{i}^{2}}(\mathbf{A}_{i}\mathbf{P}\mathbf{A}_{i}^{T}+\bm{b_{i}}\bm{b_{i}}^{T})\right)^{-1},$
(7)

and $\bm{\Lambda}$ denotes the Fisher information matrix. The optimization
problem is reformulated as

$\displaystyle\max_{\mathcal{S}}~{}~{}\sum_{i\in\mathcal{S}}\frac{1}{\sigma_{i}^{2}}\left(\text{Tr}(\mathbf{A}_{i}\mathbf{P}\mathbf{A}_{i}^{T})+\left\lVert\bm{b_{i}}\right\rVert_{2}^{2}\right)~{}~{}\text{
s.t. }\ \ |\mathcal{S}|\leq k,$ (8)

which is also equivalent to the top-$k$ selection problem with

$s_{i}\operatorname{\overset{def}{=}}\text{Tr}(\mathbf{A}_{i}\mathbf{P}\mathbf{A}_{i}^{T})+\left\lVert\bm{b_{i}}\right\rVert_{2}^{2}.$
(9)

Sensor subset selection is one of the possible applications of the results
developed herein. The underlying assumption is that the data is scored using a
certain rule. After the scores are obtained, we invoke our algorithms to find
in a distributed manner the top-$k$ scores.

##  III Problem Setup

Consider a distributed system with $n$ agents, which interact locally via a
connected undirected graph $\mathcal{G}=([n],\mathcal{E})$. Consider a dataset
with $N\geq n$ points $\mathcal{D}=\\{(x_{i},y_{i})\\}_{i=1}^{N}$ distributed
over the agents in the network. Without loss of generality, we assume that
$N=n$, and that each agent has a single data point. The analysis put forth
here can easily be extended to a local dataset with more than one point at
each node111The extension to $N\geq n$ is discussed in Appendix B.. The $i$-th
agent ascribes to its data-point $(x_{i},y_{i})$ an informativeness score
$s_{i}$, according to an application appropriate metric (c.f. Section II). We
highlight that the focus of this work is not on scoring the data, but rather
on ranking the scores in a distributed manner.

######  Definition 1 (The $k$-th largest score)

Let $\\{s_{i}\\}_{i=1}^{n}$ denote the collection of all the scores. We
arrange the $n$ scores $\\{s_{i}\\}_{i=1}^{n}$ in a new sequence
$\\{\theta_{i}\\}_{i=1}^{n}$ in descending order such that
$\theta_{1}\geq\theta_{2}\geq\cdots\geq\theta_{n}$. Define $\theta_{k}$ as the
$k$-th largest score. We allow for the possibility of repeated scores by
denoting $m$ as the number of data points with scores equal to $\theta_{k}$,
and $\overline{m}$ as the number of scores equal to $\theta_{k}$ whose index
in $\\{\theta_{i}\\}_{i=1}^{n}$ is less than or equal to $k$, and
$\underline{m}$ as the number of scores equal to $\theta_{k}$ whose index in
$\\{\theta_{i}\\}_{i=1}^{n}$ is strictly larger than $k$.

###### Example 1

Consider the multiset $\\{2,2,5,1,2\\}$, we obtain the ordered list
$\\{5,2,2,2,1\\}$, the 2nd, 3rd and 4th largest scores are $2$. Furthermore,
for $k=3$, we have $\theta_{3}=2$, $m=3$, $\overline{m}=2$ and
$\underline{m}=1$.

Figure 2: Pinball loss function used in quantile estimation. Notice it is
neither smooth nor strongly convex.

Our goal is to determine the top-$k$ largest scores via local communication.
We achieve this by obtaining a decentralized algorithm to compute a threshold
$T$ between top-$k$ largest score $\theta_{k}$ and the largest score smaller
than $\theta_{k}$. Once the $i$-th agent has the threshold $T$, it compares
$T$ with its score $s_{i}$. Then the agent knows whether it is holding one of
the top-$k$ data points or not. The data point is then transmitted to a remote
access point (c.f. Fig. 1).

###  III-A Background on Quantile Estimation

We proceed by relating the computation of the $k$-th largest score
$\theta_{k}$ with sample quantile estimation [42, 12, 14]. Let $F(x)$ be the
empirical cumulative distribution function (CDF) of the scores
$\\{s_{i}\\}_{i=1}^{n}$, i.e,

$F(x)=\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}(s_{i}\leq x),$ (10)

and $\omega_{p}$ be the $p$-th sample quantile of $\\{s_{i}\\}_{i=1}^{n}$,
i.e.,

$\omega_{p}\operatorname{\overset{def}{=}}\inf\Big{\\{}x\in\mathbb{R}\
\Big{|}\ F(x)\geq p\Big{\\}},$ (11)

where $\mathbf{1}(\cdot)$ denotes the indicator function. The following lemma
establishes the equivalence between the $k$-largest score $\theta_{k}$ with
the $p$-th sample quantile $\omega_{p}$ when
$p\in(\frac{n-k}{n},\frac{n-k+1}{n})$, which extends the result in Proposition
of [14] from continuous random variables to arbitrary random variables.

###### Lemma 1

Let $\\{s_{i}\\}_{i=1}^{n}$ be a sequence of scores, then if
$p\in(\frac{n-k}{n},\frac{n-k+1}{n})$, we have

$\theta_{k}=\omega_{p}.$ (12)

###### Proof:

The proof can be found in Appendix A-A. ∎

From to [42, Section 1.3], if $np$ is not an integer, then the $p$-th sample
quantile $\omega_{p}$ is the unique solution of the following quantile
estimation problem

$\omega_{p}=\arg\min_{x\in\mathbb{R}}\
\sum_{i=1}^{n}\rho_{p}\big{(}s_{i}-x\big{)},$ (13)

where

$\rho_{p}(x)\operatorname{\overset{def}{=}}\begin{cases}p\cdot x,&\text{if}\
x\geq 0.\\\ -(1-p)\cdot x,&\text{otherwise},\end{cases}$ (14)

is the so-called pinball loss function (c.f. Fig. 2).

Lemma 1 implies that if $p\in(\frac{n-k}{n},\frac{n-k+1}{n})$, the $k$-th
largest score can be computed as the solution of the quantile estimation
problem:

$\theta_{k}=\arg\min_{x\in\mathbb{R}}\
\sum_{i=1}^{n}\rho_{p}\big{(}s_{i}-x\big{)}\operatorname{\overset{def}{=}}\sum_{i=1}^{n}f_{i}(x)\operatorname{\overset{def}{=}}f(x),$
(15)

where $f_{i}(x)=\rho_{p}\big{(}s_{i}-x\big{)}$.

Notice that $f(x)$ is a piecewise linear convex function with a unique
minimizer. As shown in Fig. 3, we present the empirical CDF $F(x)$, and its
corresponding aggregate pinball loss function $f(x)$ with $n=15$ and $k=5$.
When $p=\frac{2n-2k+1}{2n}=0.7$, the $p$-th sample quantile corresponds to the
minimizer of $f(x)$.

Figure 3: Empirical CDF $F(x)$, and its corresponding aggregate pinball loss
function $f(x)$ with $n=15$ and $k=5$. The horizontal dotted line denotes the
choice of quantile $p$.

###  III-B Quantile Estimation via Distributed Optimization

The above problem can be solved using distributed algorithms for non-smooth
convex optimization such as the distributed subgradient method [43],
distributed dual averaging [50] and decentralized alternating direction method
of multipliers [51], to name a few. However, the aforementioned algorithms
tend to be slow due to the non-smoothness of the objective function, which
requires either an increased number of communication rounds to achieve
convergence or a higher computational cost to address a subproblem in each
iteration. Moreover, these algorithms are highly sensitive to noise in the
communication links between the agents [52, 46].

We are interested in characterizing the iteration complexity of top-$k$
selection for a list of $n$ real numbers. Intuitively, the number of
iterations required to identify the $k$ largest elements of a list should
depend on how close or far apart the elements of the list are, i.e., the
search is easier to perform if the elements are distinctively separated.
Therefore, the iteration complexity should increase if the minimum gap between
the $k$-th largest score and other scores in the list decrease. To precisely
characterize this dependence, we introduce the gap parameter
$\Delta(\theta_{k})$ to capture the tradeoff between the iteration complexity
of our algorithm and the estimation precision, which is defined as follows.

######  Definition 2 (Minimum gap from $\theta_{k}$)

The minimum gap from $\theta_{k}$, denoted by $\Delta(\theta_{k})$, is defined
as the minimum absolute difference between $k$-th largest score and other
scores in the list:

$\Delta(\theta_{k})=\min\Big{\\{}|s_{i}-\theta_{k}|\ \Big{|}\
s_{i}\neq\theta_{k},\ i=1,\ldots,n\Big{\\}}.$ (16)

Using Definition 2, we proceed by defining the optimal solution interval and
optimal threshold interval as follows.

######  Definition 3 (Optimal solution interval)

Let the optimal solution interval, $\mathcal{I}(\theta_{k})\subset\mathbb{R}$,
be defined as the open interval centered at $\theta_{k}$ with radius
$\Delta(\theta_{k})/2$, i.e.,

$\mathcal{I}(\theta_{k})\operatorname{\overset{def}{=}}\Big{(}\theta_{k}-\frac{\Delta(\theta_{k})}{2},\theta_{k}+\frac{\Delta(\theta_{k})}{2}\Big{)}.$
(17)

######  Definition 4 (Optimal threshold interval)

Let the optimal threshold interval,
$\mathcal{T}(\theta_{k})\subset\mathbb{R}$, be defined as the following open
interval:

$\mathcal{T}(\theta_{k})\operatorname{\overset{def}{=}}\big{(}\theta_{k}-\Delta(\theta_{k}),\theta_{k}\big{)}.$
(18)

Notice that for any $s\in\mathcal{I}(\theta_{k})$, we can compute a threshold
$T$ such that:

$T\operatorname{\overset{def}{=}}s-\frac{\Delta(\theta_{k})}{2}\in\mathcal{T}(\theta_{k}).$
(19)

Any threshold $T$ within the interval $\mathcal{T}(\theta_{k})$ is equally
suitable for selecting the top-$k$ elements. This is because all scores below
$\theta_{k}$ are strictly less than $T$, whereas those equal to or exceeding
$\theta_{k}$ are strictly larger than $T$. Therefore, if an iterative
algorithm produces a sequence of points $\\{w^{t}\\}_{t=0}^{\infty}$ that
converges to a value $s$ in the optimal solution interval
$\mathcal{I}(\theta_{k})$, then we can obtain a desired threshold $T$ using
Eq. 19 that guarantees a correct top-$k$ selection.

One key aspect of our analysis is that the definition of a minimum gap
$\Delta(\theta_{k})$ allows us to converge to any number within the optimal
solution interval, $\mathcal{I}(\theta_{k})$ instead of the exact optimal
solution $\theta_{k}$, which is a less stringent convergence condition. In
this paper, we design an accelerated algorithm to achieve top-$k$ selection by
exploiting the piecewise linearity of the objective functions and the optimal
solution interval. We accomplish this task by solving a smoothed distributed
quantile estimation problem. Specifically, we will initially establish a
correspondence between the (smoothed) optimal value error and the optimization
variable error within the optimal solution interval. Subsequently, we will
devise appropriate smoothed functions and employ accelerated smooth
distributed algorithms to attain our objective. From hereon, we will suppress
the dependence of $\theta_{k}$ from $\Delta(\theta_{k})$, to simplify the
notation.

Figure 4: Piecewise linear function $f(x)$, its corresponding linear functions
$f(\theta_{k})+g_{l}(x-\theta_{k})$ and $f(\theta_{k})+g_{r}(x-\theta_{k})$
around the optimal solution $\theta_{k}$. The shaded blue area denotes the
optimal solution interval $\mathcal{I}(\theta_{k})$ and the shaded purple area
denotes the optimal threshold interval $\mathcal{T}(\theta_{k})$.

##  IV Analysis

In this section, we initially establish a connection between the objective
function and variable errors for Problem (15) . Subsequently, we employ a
smooth approximation technique to refine the non-smooth objective function.
Next, we elucidate the relationship between the smoothed function and variable
errors. Based on the result, we identify the requirements for the smooth
approximation and the distributed algorithm, which will be designed in
subsequent sections.

###  IV-A Properties of the Original Objective Function

The first step in our analysis of the distributed quantile estimation problem
in Eq. 15 is to establish a correspondence between the function and the
variable errors. Let $\theta_{k}$ be the $k$-th largest score among all the
nodes in the network, which is found by solving the optimization problem in
Eq. 15. Let the function error be defined as $f(x)-f(\theta_{k})$ and variable
error be defined as $|x-\theta_{k}|$. In general, for non-smooth functions, we
can only obtain the iteration complexity to achieve a predetermined function
error, and cannot guarantee the iteration complexity to achieve a
predetermined variable error [53]. Fortunately, as shown in Fig. 4, using the
piecewise linearity of the objective function and uniqueness of the optimal
solution, we can establish a connection between the function error and
variable error in the optimal solution interval as follows.

###### Lemma 2

Let $g_{r}$ and $g_{l}$ are defined as the right-hand derivative and left-hand
derivative of $f$ at $\theta_{k}$, respectively, i.e.,

$\displaystyle g_{r}=$ $\displaystyle\lim_{\delta\to
0^{+}}\frac{f(\theta_{k}+\delta)-f(\theta_{k})}{\delta},$ (20) $\displaystyle
g_{l}=$ $\displaystyle\lim_{\delta\to
0^{-}}\frac{f(\theta_{k}+\delta)-f(\theta_{k})}{\delta}.$ (21)

Suppose $n\cdot p$ is not an integer. If there exists $x$ such that

$f(x)-f(\theta_{k})\leq\min\\{g_{r},-g_{l}\\}\cdot\frac{\Delta}{2},$ (22)

then

$|x-\theta_{k}|\leq\frac{\Delta}{2}.$ (23)

###### Proof:

The proof can be found in Appendix A-B. ∎

###### Remark 1

This lemma also applies to other piecewise linear convex local objective
functions and is not exclusive to the top-$k$ problem.

Let $p=\frac{n-k}{n}+\frac{1}{2n}$, then the exact expression for $g_{r}$ and
$g_{l}$ are given by the following lemma.

###### Lemma 3

Choosing $p=\frac{n-k}{n}+\frac{1}{2n}$, we have

$\displaystyle
g_{r}=\overline{m}-\frac{1}{2}\text{~{}~{}and~{}~{}}g_{l}=-\underline{m}-\frac{1}{2}.$
(24)

###### Proof:

The proof can be found in Appendix A-C. ∎

Combining Lemmas 2 and 3, we have the following corollary.

###### Corollary 1

Let $p=\frac{n-k}{n}+\frac{1}{2n}$. If there exists an $x\in\mathbb{R}$ such
that

$f(x)-f(\theta_{k})\leq\frac{g_{m}\Delta}{2},$ (25)

then

$|x-\theta_{k}|\leq\frac{\Delta}{2},$ (26)

where

$g_{m}\operatorname{\overset{def}{=}}\min\big{\\{}\overline{m}-\frac{1}{2},\underline{m}+\frac{1}{2}\big{\\}}.$
(27)

###  IV-B Properties of the Smoothed Objective Function

Due to the lack of smoothness of the objective function $f(x)$, existing sub-
gradient algorithms for quantile estimation suffer from very slow convergence.
To improve convergence speed, one effective strategy is to use smoothing,
which requires approximating the non-smooth function with a smooth version and
subsequently optimizing the resulting function.

Let function $\tilde{f}_{i}^{h}$ be a convex smooth approximation of $f_{i}$
indexed by a smoothing parameter $h$ and let

$\tilde{f}_{h}(x)\operatorname{\overset{def}{=}}\sum_{i=1}^{n}\tilde{f}_{i}^{h}(x),\
\ x\in\mathbb{R}.$ (28)

We say that $\tilde{f}_{h}$ is a convex smooth approximation of $f$. Let
$\theta_{k}^{h}$ be the minimizer of $\tilde{f}_{h}(x)$, i.e.,

$\theta_{k}^{h}\operatorname{\overset{def}{=}}\arg\min_{x\in\mathbb{R}}\tilde{f}_{h}(x).$
(29)

In this section, we will establish the relationship between the function error
of the smooth approximation function error, denoted as
$\tilde{f}_{h}(x)-\tilde{f}_{h}(\theta_{k}^{h})$ and variable error
$|x-\theta_{k}|$. Notice that we are interested in minimizing the smooth
approximation, but obtaining a guarantee on the variable error of the original
non-smooth function.

Before stating our result, we must introduce the following Assumption 1 on the
smooth approximation. The two typical approximations obtained using Nesterov’s
and convolution smoothing techniques satisfy Assumption 1 and will be
presented in Section V-A and V-B, respectively.

###### Assumption 1

Let the function $\tilde{f}_{h}$ denote a convex smooth approximation of $f$
with smoothing parameter $h$. The function $\tilde{f}_{h}$ is
$L_{h}$-Lipschitz continuous, i.e.,

$|\tilde{f}_{h}(x)-\tilde{f}_{h}(y)|\leq L_{h}|x-y|,\ \ x,y\in\mathbb{R},$
(30)

$M_{h}$-smooth, i.e.,

$\displaystyle|\nabla\tilde{f}_{h}(x)-\nabla\tilde{f}_{h}(y)|\leq M_{h}|x-y|,\
\ x,y\in\mathbb{R}$ (31)

and the approximation error is uniformly bounded, i.e.,

$\displaystyle|f(x)-\tilde{f}_{h}(x)|\leq U_{h},\ \ x\in\mathbb{R}.$ (32)

Furthermore, the function $\tilde{f}_{i}^{h}$ is $\frac{M_{h}}{n}$-smooth,
i.e.,

$|\nabla\tilde{f}_{i}^{h}(x)-\nabla\tilde{f}_{i}^{h}(y)|\leq\frac{M_{h}}{n}|x-y|,\
\ x,y\in\mathbb{R}.$ (33)

Theorem 1 provides the theoretical connection between
$\tilde{f}_{h}(x)-\tilde{f}_{h}(\theta_{k}^{h})$ and $|x-\theta_{k}|$ within
the optimal solution interval, which means that if we can find $x$ satisfying
$\tilde{f}_{h}(x)-\tilde{f}_{h}(\theta_{k}^{h})<g_{m}\Delta/4$, then we can
obtain the optimal solution interval and the $k$-th largest score.

######  Theorem 1 (Connection between smoothed function and variable errors)

Suppose Assumption 1 holds, and let $p=\frac{n-k}{n}+\frac{1}{2n}$. If there
exists $x$ such that the smooth approximation $\tilde{f}_{h}$ satisfies

$\tilde{f}_{h}(x)-\tilde{f}_{h}(\theta_{k}^{h})<\frac{g_{m}\Delta}{4},$ (34)

and a smoothing parameter $h$ such that

$U_{h}\leq\frac{g_{m}\Delta}{8},$ (35)

then we have

$|x-\theta_{k}|\leq\frac{\Delta}{2}.$ (36)

###### Proof:

Using the triangle inequality, we have

$|f(x)-f(\theta_{k})|\leq|f(x)-\tilde{f}_{h}(x)|+|\tilde{f}_{h}(x)-\tilde{f}_{h}(\theta_{k}^{h})|\\\
+|f(\theta_{k})-\tilde{f}_{h}(\theta_{k}^{h})|.$ (37)

From Assumption 1, the first term satisfies

$|f(x)-\tilde{f}_{h}(x)|\leq U_{h}.$ (38)

Combining Assumption 1 and the optimality of $\theta_{k}$ and
$\theta_{k}^{h}$, we have

$\displaystyle f(\theta_{k})\leq
f(\theta_{k}^{h})\leq\tilde{f}_{h}(\theta_{k}^{h})+U_{h},$ (39)

and

$\displaystyle
f(\theta_{k})\geq\tilde{f}_{h}(\theta_{k})-U_{h}\geq\tilde{f}_{h}(\theta_{k}^{h})-U_{h}.$
(40)

Therefore, the third term satisfies

$\displaystyle|f(\theta_{k})-\tilde{f}_{h}(\theta_{k}^{h})|\leq U_{h}.$ (41)

Since there exists $h$ such that $U_{h}\leq\frac{g_{m}\Delta}{8}$ and

$\tilde{f}_{h}(x)-\tilde{f}_{h}(\theta_{k}^{h})<\frac{g_{m}\Delta}{4},$ (42)

we get

$|f(x)-f(\theta_{k})|\leq\frac{g_{m}\Delta}{8}+\frac{g_{m}\Delta}{4}+\frac{g_{m}\Delta}{8}=\frac{g_{m}\Delta}{2}.$
(43)

Invoking Corollary 1, we obtain then

$|x-\theta_{k}|\leq\frac{\Delta}{2}.$ (44)

∎

Based on Theorem 1, we can design an accelerated distributed top-$k$ algorithm
by accomplishing two tasks:

  * (1)

Find a smooth approximation $\tilde{f}_{h}$ that satisfies Assumption 1 with

$U_{h}\leq\frac{g_{m}\Delta}{8};$ (45)

  * (2)

Find a distributed algorithm such that

$\tilde{f}_{h}(w_{i}^{t})-\tilde{f}_{h}(\theta_{k}^{h})\leq\frac{g_{m}\Delta}{4},\
\ i=1,\ldots,n$ (46)

with $t$ as small as possible, where $w_{i}^{t}$ is the local estimate of
$\theta_{k}^{h}$ computed by agent $i$ at the $t$-th iteration.

In the subsequent sections, we will design smoothing approximations in Section
V and provide fast distributed algorithms with corresponding iteration
complexity guarantees in Section VI.

##  V Smooth approximation

In this section, we consider two smoothing techniques to approximate $f$: The
Nesterov’s and convolution smoothing approaches [15, 53, 16, 54]. The
rationale here is to introduce the approximations for the pinball loss
function in terms of a certain smoothing parameter, and obtain sufficient
conditions that guarantee the validity of 1.

###  V-A Nesterov’s smoothing

Nesterov’s smoothing [15, 53] is a systematic way to approximate non-smooth
objective functions $f$ by a function with Lipschitz-continuous gradient. Let
$h$ be a positive smoothness parameter. The following function is Nesterov’s
smooth approximation of $f$:

$\tilde{f}_{h}^{\mathrm{nest}}(x)\operatorname{\overset{def}{=}}\sum_{i=1}^{n}\rho_{p}^{h}\big{(}s_{i}-x\big{)},$
(47)

where $\rho_{p}^{h}(x)$ is the smooth approximation of $\rho_{p}(x)$

$\rho_{p}^{h}(x)\operatorname{\overset{def}{=}}\max_{z\in\mathbb{R}}\Big{\\{}zx-\phi(z)-\frac{h}{2}z^{2}\Big{\\}},$
(48)

and $\phi(z)$ is the conjugate function of $\rho_{p}(x)$

$\phi(z)\operatorname{\overset{def}{=}}\max_{x\in\mathbb{R}}\\{zx-\rho_{p}(x)\\}.$
(49)

Nesterov’s smooth approximation of the score function, $\rho_{p}^{h}(x)$,
admits a closed-form expression, which is stated as the following result.

###### Lemma 4

The smooth approximation of $\rho_{p}(x)$ under Nesterov’s smoothing, denoted
as $\rho_{p}^{h}(x)$, is given by

$\rho_{p}^{h}(x)\operatorname{\overset{def}{=}}\begin{cases}px-\frac{h}{2}p^{2}&\text{if}\
\ x>hp\\\ \frac{x^{2}}{2h}&\text{if}\ \ h(p-1)<x\leq hp\\\
(p-1)x-\frac{h}{2}(p-1)^{2}&\text{if}\ \ x\leq h(p-1).\\\ \end{cases}$ (50)

###### Proof:

The proof can be found in Appendix A-D. ∎

Using Nesterov’s smoothing, the optimization problem in Eq. 15 can be
approximated by

$\theta^{\mathrm{nest}}_{h}\operatorname{\overset{def}{=}}\arg\min_{x\in\mathbb{R}}\tilde{f}_{h}^{\mathrm{nest}}(x)=\sum_{i=1}^{n}\rho_{p}^{h}\big{(}s_{i}-x\big{)}.$
(51)

The next Lemma establishes a sufficient condition on the smoothness parameter
$h$ such that the Nesterov’s approximation satisfies all the conditions in 1.

###### Lemma 5

If the smoothing parameter $h$ satisfies

$h\leq\frac{g_{m}\Delta}{4n\max\\{p^{2},(1-p)^{2}\\}},$ (52)

then the Nesterov’s smooth approximation
$\tilde{f}_{h}^{\mathrm{nest}}(\cdot)$ satisfies Assumption 1 with

$\displaystyle L_{h}$ $\displaystyle=$ $\displaystyle n\max\\{p,1-p\\},$ (53)
$\displaystyle M_{h}$ $\displaystyle=$ $\displaystyle\frac{n}{h},$ (54)
$\displaystyle U_{h}$ $\displaystyle=$
$\displaystyle\frac{nh}{2}\max\\{p^{2},(1-p)^{2}\\}\leq\frac{g_{m}\Delta}{8}.$
(55)

###### Proof:

The proof can be found in Appendix A-E. ∎

Figure 5: Examples of smooth approximation for different smoothing parameter
$h$ with $n=10$, $k=4$, and $p=0.65$: (a) Nesterov’s smoothing, (b)
Convolution smoothing. Here, $f(x)$ denotes the original piecewise linear
function, $f_{h}^{\mathrm{nest}}(x)$ denotes the Nesterov’s smoothed function,
$f_{h}^{\mathrm{conv}}(x)$ denotes the convolution smoothed function, and the
marker $\times$ denotes the minimizer of a function.

###  V-B Convolution smoothing

Another way to obtain an approximation with a tunable smoothness parameter use
convolutions. By using the so-called convolution smoothing approach [16], the
optimization problem in Eq. 15 becomes

$\theta_{h}^{\mathrm{conv}}\operatorname{\overset{def}{=}}\arg\min_{x\in\mathbb{R}}\tilde{f}_{h}^{\mathrm{conv}}(x)=\sum_{i=1}^{n}l_{p}^{h}\big{(}s_{i}-x\big{)},$
(56)

where

$l_{p}^{h}(x)\operatorname{\overset{def}{=}}(\rho_{p}*K_{h})(x)=\int_{-\infty}^{\infty}\rho_{p}(y)K_{h}(y-x){\rm{d}}y,$
(57)

where $*$ denotes the convolution operation. The function $K_{h}(x)$ is a
scaling of a convolution kernel, $K(x)$, [54]. That is

$K_{h}(x)=\frac{1}{h}K\left(\frac{x}{h}\right).$ (58)

######  Definition 5 (Convolution Kernel)

A convolution kernel function $K(x)$ is a symmetric, nonnegative and bounded
function that integrates to one, i.e.,

$\displaystyle K(x)=K(-x),\ x\in\mathbb{R},$ $\displaystyle K(x)\geq 0,\
x\in\mathbb{R},$ (59)
$\displaystyle\overline{K}=\sup_{x\in\mathbb{R}}K(x)<\infty,$
$\displaystyle\int_{-\infty}^{\infty}K(x){\rm d}x=1.$ (60)

###### Example 2

When $K(x)$ is the uniform kernel with $K(x)=\frac{1}{2}\mathbf{1}(|x|\leq
1)$, the smooth approximation of the score function is:

$l_{p}^{h}(x)=\begin{cases}px&\text{if}\ \ x>h\\\
\frac{(x-h)^{2}}{4h}+px&\text{if}\ \ -h<x\leq h\\\ (p-1)x&\text{if}\ \
x\leq-h.\\\ \end{cases}$ (61)

The following lemma provides as sufficient condition on the smoothing
parameter $h$ such that $\tilde{f}_{h}^{\mathrm{conv}}(\cdot)$ satisfies
Assumption 1 with $U_{h}\leq\frac{g_{m}\Delta}{8}$.

###### Lemma 6

For any convolution kernel defined according to Definition 5, if the smoothing
parameter $h$ satisfies

$h\leq\frac{g_{m}\Delta}{8n\max\\{p,1-p\\}\int_{-\infty}^{\infty}|t|K(t){\rm{d}}t},$
(62)

then the smooth approximation $\tilde{f}_{h}^{\mathrm{conv}}(\cdot)$ satisfies
Assumption 1 with

$\displaystyle L_{h}$ $\displaystyle=$ $\displaystyle n\max\\{p,1-p\\},$ (63)
$\displaystyle M_{h}$ $\displaystyle=$ $\displaystyle\frac{n\overline{K}}{h},$
(64) $\displaystyle U_{h}$ $\displaystyle=$ $\displaystyle
nh\max\\{p,1-p\\}\int_{-\infty}^{\infty}|t|K(t){\rm{d}}t\leq\frac{g_{m}\Delta}{8}.$
(65)

###### Proof:

The proof can be found in Appendix A-F. ∎

Figure 5 shows the curves for $f(x)$ along with its Nesterov’s smooth
approximation $\tilde{f}_{h}^{\mathrm{nest}}(x)$ and its convolution smooth
approximation $\tilde{f}_{h}^{\mathrm{conv}}(x)$ for randomly generated data
points for different parameters $h$ with $n=10,k=4$ and $p=0.65$. The
convolution kernel function is chosen as $K(x)=\frac{1}{2}\mathbf{1}(|x|\leq
1)$. First, notice that $\tilde{f}_{h}^{\mathrm{nest}}(x)$ approximates $f(x)$
from below and $\tilde{f}_{h}^{\mathrm{conv}}(x)$ approximates $f(x)$ from
above. Both approximations converge to the original function as $h$ converges
to zero. The main distinction between the two smoothing techniques is the
location of the minimizer with respect to the changes in $h$. The minimizers
of Nestrov’s smooth approximation are close to the solution of quantile
estimation, however the distance between the minimizers and the solution does
not necessarily decrease as $h$ decreases; On the other hand, the distance
between the minimizers of the convolution smooth approximation and the
solution of the original function is monotone decreasing as $h$ decreases.

Algorithm 1 Fast distributed quantile estimation via EXTRA

Input: $\tilde{f}_{1:n}^{h},T,w_{1:n}^{0},v_{1:n}^{0}$

for $t=0,1,\cdots,T$ do

for $i=1,2,\ldots,n$ do

$w_{i}^{t+1}=w_{i}^{t}-\alpha\Big{(}\nabla\tilde{f}_{i}^{h}(w_{i}^{t})+v_{i}^{t}+\frac{\beta}{2}\Big{(}w_{i}^{t}-\sum_{j\in\mathcal{N}_{i}}W_{i,j}w_{j}^{t}\Big{)}\Big{)}$

$v_{i}^{t+1}=v_{i}^{t}+\frac{\beta}{2}\Big{(}w_{i}^{t+1}-\sum_{j\in\mathcal{N}_{i}}W_{i,j}w_{j}^{t+1}\Big{)}$

end for

end for

Output: $\bar{w}_{1:n}^{T+1}$

##  VI Complexity of distributed top-$k$ selection

Combining the state-of-the-art algorithm EXTRA [17, 18] and the auxiliary
results on the smoothed optimization problem obtained in the previous
sections, we are equipped to obtain the iteration complexity of distributed
top-$k$ selection. EXTRA is a sophisticated decentralized optimization
algorithm for distributed smooth convex problems, which uses the differences
of gradients and achieves convergence with a constant step-size.

Let $w_{i}^{t}$ denote the local estimate of $\theta_{k}^{h}$ computed by node
$i$ at the $t$-th iteration and $\bar{w}_{i}^{T}$ be a value generated from
$\\{w_{i}^{t}\\}_{t=1}^{T}$. Our goal is to obtain an iteration complexity
result, i.e., the minimum value of $T$ such that

$\tilde{f}_{h}(\bar{w}_{i}^{T})-\tilde{f}_{h}(\theta_{k}^{h})<\frac{g_{m}\Delta}{4},\
\ i=1,\ldots,n.$ (66)

We adopt the EXTRA algorithm and the improved analysis from [18], where the
convergence is given via the running local average

$\bar{w}_{i}^{T}\operatorname{\overset{def}{=}}\frac{1}{T}\sum_{t=1}^{T}w_{i}^{t}.$
(67)

Let
$\bar{\bm{w}}\operatorname{\overset{def}{=}}[\bar{w}_{1},\ldots,\bar{w}_{n}]^{T}$.
The algorithm is presented in Algorithm 1, which is written in a primal-dual
framework. The primal variable is updated by

$w_{i}^{t+1}=w_{i}^{t}-\alpha\Big{(}\nabla\tilde{f}_{i}^{h}(w_{i}^{t})+v_{i}^{t}+\frac{\beta}{2}\Big{(}w_{i}^{t}-\sum_{j\in\mathcal{N}_{i}}W_{i,j}w_{j}^{t}\Big{)}\Big{)},$
(68)

and the dual variable is updated by

$v_{i}^{t+1}=v_{i}^{t}+\frac{\beta}{2}\Big{(}w_{i}^{t+1}-\sum_{j\in\mathcal{N}_{i}}W_{i,j}w_{j}^{t+1}\Big{)},$
(69)

where $\alpha,\beta>0$ are constants, $\mathcal{N}_{i}$ denotes the set of
neighbors that can communicate locally with node $i$ and $W_{i,j}$ is the
$(i,j)$-th entry of the weight matrix $\bm{W}$. The readers are referred to
[55] for properties and design of weight matrices.

Let $\gamma(\bm{x})\operatorname{\overset{def}{=}}(1/n)\sum_{i=1}^{n}{x}_{i}$.
Following [18, Lemma 3], using the fact that each $\tilde{f}_{i}^{h}(x)$ is
$M_{h}/n$-smooth, we may rewrite the convergence of function error
$\tilde{f}_{h}\left(\gamma\left(\bar{\bm{w}}^{T}\right)\right)-\tilde{f}_{h}\left(\theta_{k}^{h}\right)$
and variable error
$(1/n)\sum_{i=1}^{n}\left|\bar{{w}}_{i}^{T}-\gamma\left(\bar{\bm{w}}^{T}\right)\right|^{2}$.

###### Lemma 7

Suppose that Assumption 1 holds and each $\tilde{f}_{i}^{h}(x)$ is
$M_{h}/n$-smooth. Let
$\alpha\operatorname{\overset{def}{=}}\frac{1}{2(M_{h}/n+\beta)}$,
$\beta\operatorname{\overset{def}{=}}\frac{M_{h}}{n\sqrt{1-\sigma_{2}(\bm{W})}}$,
and
$\bar{\bm{w}}^{T}\operatorname{\overset{def}{=}}\frac{1}{T}\sum_{t=1}^{T}\bm{w}^{t}$.
If $T\geq\frac{1}{\sqrt{1-\sigma_{2}(\bm{W})}}$, then the following
inequalities hold:

$\tilde{f}_{h}\left(\gamma\left(\bar{\bm{w}}^{T}\right)\right)-\tilde{f}_{h}\left(\theta_{k}^{h}\right)\leq\frac{34}{T\sqrt{1-\sigma_{2}(\bm{W})}}\left(R_{1}M_{h}+\frac{R_{2}n^{2}}{M_{h}}\right),$
(70)

and

$\frac{1}{n}\sum_{i=1}^{n}\left|\bar{{w}}_{i}^{T}-\gamma\left(\bar{\bm{w}}^{T}\right)\right|^{2}\leq\frac{16}{T^{2}\left(1-\sigma_{2}(\bm{W})\right)}\left(R_{1}+\frac{R_{2}n^{2}}{M_{h}^{2}}\right),$
(71)

where $\sigma_{2}(\bm{W})$ is the second largest singular value of $\bm{W}$,
$R_{1}$ is a constant that satisfies $|w_{i}^{0}-\theta_{k}^{h}|^{2}\leq
R_{1}$ and $|\theta_{k}^{h}|^{2}\leq R_{1}$ and
$R_{2}=\max\\{p^{2},(1-p)^{2}\\}$.

###### Proof:

Eq. (70) is derived by assigning $L=\frac{M_{h}}{n}$ and then multiplying both
sides of the first inequality by $n$ in Lemma 3 from [18]. Similarly, Eq. (71)
is obtained by setting $L=\frac{M_{h}}{n}$ in the second inequality of Lemma 3
from [18].

∎

We provide the exact number of iterations required to enter the optimal
solution interval.

######  Theorem 2 (Exact Complexity of Distributed Top-$k$ Selection)

Suppose that Assumption 1 holds, and $h$ is chosen such that
$U_{h}\leq\frac{g_{m}\Delta}{8}$. Then, we need

$T\geq\frac{1}{g_{m}\Delta\sqrt{1-\sigma_{2}(\bm{W})}}\max\Bigg{\\{}272\left(R_{1}M_{h}+\frac{R_{2}n^{2}}{M_{h}}\right),\\\
32L_{h}\sqrt{n}\left(\sqrt{R_{1}}+\sqrt{R_{2}}\frac{n}{M_{h}}\right)\Bigg{\\}}$
(72)

to reach the optimal solution interval using Algorithm 1.

###### Proof:

From Lemma 7 and Assumption 1, we obtain the following inequalities hold:

$\displaystyle|$ $\displaystyle\tilde{f}_{h}$
$\displaystyle(\bar{w}_{i}^{T})-\tilde{f}_{h}(\theta_{k}^{h})|$ (73)
$\displaystyle\leq$
$\displaystyle|\tilde{f}_{h}(\bar{w}_{i}^{T})-\tilde{f}_{h}(\gamma(\bar{\bm{w}}^{T}))|+|\tilde{f}_{h}(\gamma\left(\bar{\bm{w}}^{T}\right))-\tilde{f}_{h}(\theta_{k}^{h})|$
$\displaystyle\leq$ $\displaystyle
L_{h}|\bar{w}_{i}^{T}-\gamma(\bar{\bm{w}}^{T})|+\frac{34}{T\sqrt{1-\sigma_{2}(\bm{W})}}\left(R_{1}M_{h}+\frac{R_{2}n^{2}}{M_{h}}\right)$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}$
$\displaystyle\frac{4L_{h}\sqrt{n}}{T\sqrt{1-\sigma_{2}(\bm{W})}}\left(\sqrt{R_{1}}+\sqrt{R_{2}}\frac{n}{M_{h}}\right)$
$\displaystyle\qquad\qquad+\frac{34}{T\sqrt{1-\sigma_{2}(\bm{W})}}\left(R_{1}M_{h}+\frac{R_{2}n^{2}}{M_{h}}\right),$

where inequality $(a)$ follows from Eq. 71 and $x^{2}+y^{2}\leq(x+y)^{2}$,
when $x,y\geq 0$.

Thus, if

$T\geq\frac{1}{g_{m}\Delta\sqrt{1-\sigma_{2}(\bm{W})}}\max\Bigg{\\{}272\left(R_{1}M_{h}+\frac{R_{2}n^{2}}{M_{h}}\right),\\\
32L_{h}\sqrt{n}\left(\sqrt{R_{1}}+\sqrt{R_{2}}\frac{n}{M_{h}}\right)\Bigg{\\}},$
(74)

we have

$\displaystyle|\tilde{f}_{h}(\bar{w}_{i}^{T})-\tilde{f}_{h}(\theta_{k}^{h})|\leq\frac{g_{m}\Delta}{4}.$
(75)

Therefore, invoking Theorem 1, we have

$|\bar{w}_{i}^{T}-\theta_{k}|\leq\frac{\Delta}{2},~{}~{}i=1,\ldots,n.$ (76)

∎

###  VI-A Order of Complexity of Distributed Top-k Selection

Incorporating the definition of $L_{h}$, $M_{h}$ and the largest $h$
satisfying $U_{h}\leq\frac{g_{m}\Delta}{8}$ in Lemmas 5 and 6, Nesterov’s
smoothing approach leads to the following iteration complexity to reach the
optimal solution interval:

$\mathcal{O}\bigg{(}\frac{1}{\sqrt{1-\sigma_{2}(\bm{W})}}\max\bigg{\\{}\frac{n^{2}}{g_{m}^{2}\Delta^{2}},\frac{n^{\frac{3}{2}}}{g_{m}\Delta},\sqrt{n}\bigg{\\}}\bigg{)},$
(77)

Similarly, the convolution smoothing approach leads to the following iteration
complexity to reach the optimal solution interval:

$\mathcal{O}\bigg{(}\frac{1}{\sqrt{1-\sigma_{2}(\bm{W})}}\max\bigg{\\{}\frac{n^{2}\overline{K}\int_{-\infty}^{\infty}|t|K(t){\rm{d}}t}{g_{m}^{2}\Delta^{2}},\\\
\frac{n^{\frac{3}{2}}}{g_{m}\Delta},\frac{\sqrt{n}}{\overline{K}\int_{-\infty}^{\infty}|t|K(t){\rm{d}}t}\bigg{\\}}\bigg{)}.$
(78)

The expressions above indicate that with the increase of resolution, $\Delta$,
and the multiplicity parameter, $g_{m}$ (c.f. Eq. 27), the iteration
complexity initially decreases and then stabilizes when $g_{m}\Delta$ reaches
a sufficiently large value. However, increasing the resolution $\Delta$
diminishes the estimation precision of the $k$-th largest value. Moreover, as
the network connectivity increases, the parameter $\sigma_{2}(\bm{W})$
diminishes, resulting in a reduction of the iteration complexity.

Increasing the number of agents contributes to an increase in iteration
complexity, primarily because the communication network and its associated
diameter expand with the growing number of agents.

Finally, we discuss the choice of convolution Kernel in Eq. 78. If the
resolution parameter $\Delta$ is small, the first term in the max function,
that is,

$\frac{n^{2}\overline{K}\int_{-\infty}^{\infty}|t|K(t){\rm{d}}t}{g_{m}^{2}\Delta^{2}},$
(79)

will be the dominant term. Therefore, we should choose a convolution Kernel
such that $\int_{-\infty}^{\infty}|t|K(t){\rm{d}}t$ is minimized. According to
Definition 5, we know $K(t)\leq\overline{K}$ for all $t\in\mathbb{R}$.
Therefore, the best choice is to place most mass around $t=0$ with magnitude
$\overline{K}$, that is,

$K^{\star}(x)=\overline{K}\,\mathbf{1}\left(|x|\leq\frac{1}{2\overline{K}}\right).$
(80)

Figure 6: Variable error as a function of the number of iterations with
$k=10^{3}$ and $\Delta=0.1$, (a) $N=n=10^{4}$; (b) $N=10^{5}$, $n=10^{3}$. The
black horizontal line denotes the function
$||\bm{w}^{t}-\theta_{k}\mathbf{1}||_{\infty}=\frac{\Delta}{2}$.

The implication is that for convolution smoothing the kernel function that
minimizes the iteration complexity is the uniform kernel.

###  VI-B Fixing the resolution via score quantization

One of the key parameters in the iteration complexity analysis is the
resolution parameter, $\Delta$. We have previously mentioned that this
parameter depends on the desired quantile, $\theta_{k}$, which is unknown.
Therefore, the precise iteration complexity varies with each dataset. However,
we can adopt a scoring mechanism that guarantees that any two nonidentical
consecutive scores are at least $\Delta$ apart. This is accomplished using the
following data preconditioning operation.

We process the list of original scores $\\{z_{i}\\}_{i=1}^{n}$ by rounding
each element using $\Delta>0$ as follows

$s_{i}=\texttt{round}\left(\frac{z_{i}}{\Delta}\right)\times\Delta,~{}~{}i=1,\ldots,n,$
(81)

where $\texttt{round}(\cdot)$ means rounding to the nearest integer. It is
important to notice that by quantizing the scores, information about the true
quantile with respect to the original score list is lost. Consequently, it is
important to select an appropriate resolution $\Delta$ that strikes a good
balance between iteration complexity and quantile estimation accuracy.

##  VII Numerical Results

In this section, we provide numerical simulations comparing our top-$k$
selection scheme via smoothed quantile estimation via EXTRA with traditional
nonsmooth quantile estimation via distributed subgradient descent (DGD) [43].
We randomly generate scores using the quantization scheme to guarantee a
minimum gap $\Delta$ for any $\theta_{k}$. Here, we exclusively used
convolution smooth approximations. Although the numerical results obtained via
Nesterov’s smooth approximations are highly comparable, we have omitted them
for brevity. Nevertheless, these results can be readily accessed by utilizing
the provided code 222The code is available at Github
https://github.com/connorzhangxu/DistributedFastTopKSelection..

The connected graph is an Erdős-Rényi graph, which is generated randomly with
$|\mathcal{E}|$ edges. The weight matrix is chosen as
$\bm{W}=\bm{I}-\varrho\bm{L}$ according to [55], where
$\varrho=\frac{2}{\lambda_{1}(\bm{L})+\lambda_{n-1}(\bm{L})}$, $\bm{L}$ is
Laplacian matrix of the graph and $\lambda_{i}(\bm{L})$ denotes the $i$-th
largest eigenvalue of a $\bm{L}$. The convolution kernel function is chosen as
$K(x)=\frac{1}{2}\mathbf{1}(|x|\leq 1)$. The scores are randomly generated
Gaussian variables with variance $\sigma^{2}$, and then rounded using the
scheme in Section VI-B. The constraint on the smoothing parameter in Lemma 6
directly affects the convergence rate. To accelerate the convergence of the
algorithm, we manually adjust the parameter $h$. Following the setting in
[18], we set $\alpha=h$ and $\beta=\frac{1}{h}$. Additionally, we opt for a
constant step size manually to accelerate the DGD algorithm towards the
optimal solution interval. To guarantee all nodes reach the optimal solution
interval, we compare the maximum variable error
$||\bm{w}^{t}-\theta_{k}\mathbf{1}||_{\infty}$ with $\Delta/2$ to find the
required iterations satisfying
$||\bm{w}^{t}-\theta_{k}\mathbf{1}||_{\infty}<\Delta/2$.

###  VII-A Convergence rate

Figure 6 shows the convergence of distributed top-$k$ selection by plotting
the variable error as a function of the number of iterations of Algorithm 1.
In Fig. 6 (a), the agents are collectively estimating the quantile
corresponding to the top-$10^{3}$ scores. The network is sampled from a random
graph ensemble with $n=10^{4}$ nodes and $|\mathcal{E}|=5n$ edges. In this
simulation, each agent has a single score, i.e., $N=n=10^{4}$, which has been
randomly drawn from a Gaussian distribution with variance $\sigma^{2}=10$ and
then quantized such that the minimum gap from $\theta_{k}$ (resolution) is
constant, $\Delta=0.1$. Fig. 6 (a) demonstrates that the number of iterations
required to identify the top-$k$ for our algorithm is $230$ while for
distributed gradient descent (DGD) algorithm is much larger at $12693$. In
this case, our algorithm performs approximately $55$ times faster than DGD.
This result indicates that we can run the algorithm efficiently in large-scale
settings.

In Fig. 6 (b), we consider a slightly modified setting. The network is sampled
from a random graph ensemble with $n=10^{3}$ nodes and $|\mathcal{E}|=3n$
edges. Each agent has a random local dataset with at least one score and the
total number of scores across all agents is $N=10^{5}$. Figure 6 (b) shows
that the number of iterations required to distributedly identify the
top-$1000$ scores using our algorithm is $211$, which is approximately $17$
times faster than DGD. Our implementation for top-$k$ selection in this case
is very efficient with respect to communication requirements because each
agent only communicates a single variable with their neighbors.

Figure 7: Required iterations to reach the optimal solution interval as a
function of $k$ with $n=1000$ and $\Delta=0.1$. The blue line with a square
marker is the average of 100 trials, and the shaded area denotes the standard
error.

###  VII-B Iteration complexity as a function of desired quantile

We show the required number of iterations to reach the optimal solution
interval as a function of $k$ for $N=n=10^{3}$ connected by a random graph
with $|\mathcal{E}|=3n$ in Fig. 7. Here we generated a random dataset from a
Gaussian distribution with variance $\sigma^{2}=10$ and quantized them to
obtain scores with a resolution of $\Delta=0.1$. We conducted 100 Monte Carlo
trials and calculated the average number of required iterations to enter the
optimal solution interval. It can be seen that the required iterations are
small when $k$ is approximately $n/2$, while they become large when $k$ is
close to $1$ and $n$. Therefore, finding the minimum and the maximum elements
over the network using Algorithm 1 is harder than finding the median, for
example.

Figure 8: Required iterations to reach the optimal solution interval as a
function of $\Delta$ with $n=1000$ and $k=100$.  Figure 9: Required iterations
to reach the optimal solution interval as a function of $n$ with ${k}/{n}=0.1$
and $\Delta=0.1$.

###  VII-C Iteration complexity as a function of the quantization gap

Figure 8 shows the dependence of the total number of iterations required to
reach the optimal solution interval and the number of declared top-$k$ scores
as a function of $\Delta$ for $n=10^{3}$ and $k=10^{2}$. We performed 100
Monte Carlo simulations to obtain the results. For each Monte Carlo
simulation, we generated random scores $\\{s_{i}\\}_{i=1}^{n}$ and quantized
them using different values of $\Delta$. The results show that by increasing
$\Delta$, the iterations to reach the optimal solution interval decrease,
which coincides with our theoretical results for iteration complexity. In
addition, the number of declared top-$k$ scores increases as $\Delta$
increases. This is due to the fact that the number of scores equal to the
$k$-th largest value increases as the quantization gap $\Delta$ increases.

###  VII-D Iteration complexity as a function of the number of agents

We also compared the total number of iterations required to reach the optimal
solution interval as a function of $n$ in Fig. 9, where the bold blue line
represents the average of 100 trials, and the shaded area denotes the standard
deviation. Each trial corresponds to a different set of scores generated
randomly sampled from a Gaussian distribution with variance $\sigma^{2}=100$.
To isolate the possible influence of the connectivity graph, we adopted a ring
graph with $n$ nodes in this simulation. It is evident that the required
iterations increase as the number of nodes increases. However, the average
iteration complexity for these randomly generated score lists seems to scale
linearly with $n$, which is a desirable feature from the practical
perspective.

Figure 10: Communication cost to achieve top-$k$ selection as a function of
$k$ with $N=1000$, $n=100$ and $\Delta=0.1$.

###  VII-E Communication cost

Finally, we compare our distributed top-$k$ selection method against a
straightforward and intuitive message-passing algorithm outlined in [46]. The
simple strategy operates as follows: Each node maintains a list of at most $k$
scores with their corresponding indices in its memory. At iteration $t$, each
node transmits this list to its neighboring nodes. At iteration $t+1$, each
node updates its list by selecting the top-$k$ scores and discarding the rest.
Subsequently, each node sorts its list and repeats this process iteratively.
To operate more efficiently, when $k\leq n/2$, each node selects the top-$k$
scores; when $k>n/2$, each node chooses the bottom-$k$ scores. We refer to
this scheme as Simple Top-$k$ (STop-$k$) method. Figure 10 shows the
communication cost is simulated for a fixed graph with $|\mathcal{E}|=3n$
edges using the STop-$k$ and our method using EXTRA. Here, we have used
$N=10^{4}$, $n=10^{2}$ and $\Delta=0.1$. The communication cost denotes the
number of transmitted scalars, without taking into account the cost for
transmitting the indices in the STop-$k$ method. The results are averaged over
100 trials. From this experiment, it is evident that EXTRA outperforms
STop-$k$ for $k$ between 100 and 900, specially when $k$ is approximately
$500$. This difference arises due to the fact that the communication cost of
STop-$k$ scales with $k$, whereas our method using EXTRA relatively
insensitive to $k$ except when used to compute extrema (e.g. minima and
maxima). Therefore, using our algorithm is advantageous for a wide range of
applications, in which we are interested in determining the top-$k$ scores
when $k$ is not small compared to the total number of scores, $N$.

##  VIII Conclusion and Future Work

Top-$k$ selection algorithms have been studied in many forms for decades and
remain relevant in technology, specially in distributed systems. Our
distributed top-$k$ algorithm is based on distributed optimization and offers
a versatile solution to determine the top-$k$ largest entries in a networked
dataset. This framework is applicable in various fields such as wireless
sensor networks, signal processing, and machine learning. Unlike existing
methods relying on spanning trees, our approach employs distributed
optimization and is adaptable to handle noise, packet drops, and other
communication imperfections. Leveraging the properties of our local objective
functions, we introduced an accelerated algorithm based on smoothing
techniques, and expressions for its iteration complexity. As a side product,
our method promotes privacy by avoiding data transmission, as nodes estimate a
single common threshold to determine whether they are holding or not a
top-$k$data-point. Simulations demonstrated its superiority over Distributed
Gradient Descent and a simple aggregation-based method called STop-$k$. Future
work involves addressing communication imperfections and enabling asynchronous
implementation using Gossip mechanisms.

## References

  * [1] M. M. Vasconcelos and U. Mitra, “Extremum information transfer over networks for remote estimation and distributed learning,” _Frontiers in Complex Systems_ , vol. 2, 2024. [Online]. Available: https://www.frontiersin.org/articles/10.3389/fcpxs.2024.1322785
  * [2] A. Krause and D. Golovin, “Submodular function maximization.” _Tractability_ , vol. 3, no. 71-104, p. 3, 2014.
  * [3] A. Hashemi, M. Ghasemi, H. Vikalo, and U. Topcu, “Submodular observation selection and information gathering for quadratic models,” in _International Conference on Machine Learning_. PMLR, 2019, pp. 2653–2662.
  * [4] T. Nishio and R. Yonetani, “Client selection for federated learning with heterogeneous resources in mobile edge,” in _IEEE international conference on communications (ICC)_ , 2019, pp. 1–7.
  * [5] Y. J. Cho, J. Wang, and G. Joshi, “Client selection in federated learning: Convergence analysis and power-of-choice selection strategies,” _arXiv preprint arXiv:2010.01243_ , 2020.
  * [6] O. Marnissi, H. E. Hammouti, and E. H. Bergou, “Client selection in federated learning based on gradients importance,” _arXiv preprint arXiv:2111.11204_ , 2021.
  * [7] L. Hübschle-Schneider and P. Sanders, “Communication efficient algorithms for top-$k$ selection problems,” in _IEEE International Parallel and Distributed Processing Symposium (IPDPS)_ , 2016, pp. 659–668.
  * [8] D. Durfee and R. M. Rogers, “Practical differentially private top-$k$ selection with pay-what-you-get composition,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * [9] J. Shiraishi, H. Yomo, K. Huang, Č. Stefanović, and P. Popovski, “Content-based wake-up for top-$k$ query in wireless sensor networks,” _IEEE Transactions on Green Communications and Networking_ , vol. 5, no. 1, pp. 362–377, 2020.
  * [10] B. Malhotra, M. A. Nascimento, and I. Nikolaidis, “Exact top-$k$ queries in wireless sensor networks,” _IEEE Transactions on Knowledge and Data Engineering_ , vol. 23, no. 10, pp. 1513–1525, 2010.
  * [11] Z. Liu, A. Manousis, G. Vorsanger, V. Sekar, and V. Braverman, “One sketch to rule them all: Rethinking network flow monitoring with univmon,” in _Proceedings of the 2016 ACM SIGCOMM Conference_ , 2016, pp. 101–114.
  * [12] J. Lee, C. Tepedelenlioglu, A. Spanias, and G. Muniraju, “Distributed quantiles estimation of sensor network measurements,” _International Journal of Smart Security Technologies (IJSST)_ , vol. 7, no. 2, pp. 38–61, 2020\.
  * [13] H. Wang and C. Li, “Distributed quantile regression over sensor networks,” _IEEE Transactions on Signal and Information Processing over Networks_ , vol. 4, no. 2, pp. 338–348, 2018.
  * [14] X. Zhang, M. M. Vasconcelos, W. Cui, and U. Mitra, “Distributed remote estimation over the collision channel with and without local communication,” _IEEE Transactions on Control of Network Systems_ , vol. 9, no. 1, pp. 282–294, 2021.
  * [15] Y. Nesterov, “Smooth minimization of non-smooth functions,” _Mathematical programming_ , vol. 103, no. 1, pp. 127–152, 2005.
  * [16] M. Fernandes, E. Guerre, and E. Horta, “Smoothing quantile regressions,” _Journal of Business & Economic Statistics_, vol. 39, no. 1, pp. 338–357, 2021.
  * [17] W. Shi, Q. Ling, G. Wu, and W. Yin, “EXTRA: An exact first-order algorithm for decentralized consensus optimization,” _SIAM Journal on Optimization_ , vol. 25, no. 2, 2014.
  * [18] H. Li and Z. Lin, “Revisiting EXTRA for smooth distributed optimization,” _SIAM Journal on Optimization_ , vol. 30, no. 3, pp. 1795–1821, 2020.
  * [19] M. Blum, R. W. Floyd, V. R. Pratt, R. L. Rivest, R. E. Tarjan _et al._ , “Time bounds for selection,” _Journal of Computer and System Sciences_ , vol. 7, no. 4, pp. 448–461, 1973.
  * [20] V. Zois, V. J. Tsotras, and W. A. Najjar, “GPU accelerated top-$k$ selection with efficient early stopping,” in _Workshop on Accelerating Analytics and Data Management Systems Using Modern Processor and Storage Architectures (ADMS)_ , 2019.
  * [21] J. Zhang, A. Naruse, X. Li, and Y. Wang, “Parallel top-$k$ algorithms on GPU: A comprehensive study and new methods,” in _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_. New York, NY, USA: Association for Computing Machinery, 2023.
  * [22] Y. Li, B. Zhou, J. Zhang, X. Wei, Y. Li, and Y. Chen, “RadiK: Scalable radix top-$k$ selection on GPUs,” in _Proceedings of the 29th ACM SIGPLAN Annual Symposium on Principles and Practice of Parallel Programming_. New York, NY, USA: Association for Computing Machinery, 2024, pp. 472–474.
  * [23] D. Üstebay, R. Castro, and M. Rabbat, “Efficient decentralized approximation via selective gossip,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 5, no. 4, pp. 805–816, 2011.
  * [24] D. Üstebay and M. Rabbat, “Efficiently reaching consensus on the largest entries of a vector,” in _51st IEEE Conference on Decision and Control (CDC)_ , 2012, pp. 56–61.
  * [25] ——, “Top-$k$ selective gossip,” in _IEEE 13th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)_ , 2012, pp. 505–509.
  * [26] I. F. Ilyas, G. Beskales, and M. A. Soliman, “A survey of top-$k$ query processing techniques in relational database systems,” _ACM Computing Surveys (CSUR)_ , vol. 40, no. 4, pp. 1–58, 2008.
  * [27] K. V. Jonsson, K. Palmskog, and Y. Vigfusson, “Secure distributed top-$k$ aggregation,” in _IEEE International Conference on Communications (ICC)_ , 2012, pp. 804–809.
  * [28] J. Liang, C. Jiang, X. Ma, G. Wang, and X. Kui, “Secure data aggregation for top-$k$ queries in tiered wireless sensor networks.” _Adhoc & Sensor Wireless Networks_, vol. 32, 2016.
  * [29] M. Wu, J. Xu, X. Tang, and W.-C. Lee, “Top-$k$ monitoring in wireless sensor networks,” _IEEE Transactions on Knowledge and Data Engineering_ , vol. 19, no. 7, pp. 962–976, 2007.
  * [30] A. Abbasi, A. Khonsari, and N. Farri, “Mote: efficient monitoring of top-$k$ set in sensor networks,” in _IEEE Symposium on Computers and Communications_ , 2008, pp. 957–962.
  * [31] D. Kempe, A. Dobra, and J. Gehrke, “Gossip-based computation of aggregate information,” in _44th Annual IEEE Symposium on Foundations of Computer Science, 2003. Proceedings._ IEEE, 2003, pp. 482–491.
  * [32] S. Kashyap, S. Deb, K. Naidu, R. Rastogi, and A. Srinivasan, “Efficient gossip-based aggregate computation,” in _Proceedings of the twenty-fifth ACM SIGMOD-SIGACT-SIGART symposium on Principles of database systems_ , 2006, pp. 308–317.
  * [33] B. Haeupler, J. Mohapatra, and H.-H. Su, “Optimal gossip algorithms for exact and approximate quantile computations,” in _Proceedings of the 2018 ACM Symposium on Principles of Distributed Computing_ , 2018, pp. 179–188.
  * [34] R. Rajagopal, M. J. Wainwright, and P. Varaiya, “Universal quantile estimation with feedback in the communication-constrained setting,” in _2006 IEEE International Symposium on Information Theory_. IEEE, 2006, pp. 836–840.
  * [35] F. Kuhn, T. Locher, and R. Wattenhofer, “Tight bounds for distributed selection,” in _Proceedings of the nineteenth annual ACM symposium on Parallel algorithms and architectures_ , 2007, pp. 145–153.
  * [36] S. Liu, L. Xie, and H. Zhang, “Distributed consensus for multi-agent systems with delays and noises in transmission channels,” _Automatica_ , vol. 47, no. 5, pp. 920–934, 2011.
  * [37] G. Chen, C. Chen, G. Yin _et al._ , “Critical connectivity and fastest convergence rates of distributed consensus with switching topologies and additive noises,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 12, pp. 6152–6167, 2017.
  * [38] J. Zhou, G. Gu, and X. Chen, “Distributed kalman filtering over wireless sensor networks in the presence of data packet drops,” _IEEE Transactions on Automatic Control_ , vol. 64, no. 4, pp. 1603–1610, 2018.
  * [39] Q. Yang, G. Chen, and T. Wang, “Admm-based distributed algorithm for economic dispatch in power systems with both packet drops and communication delays,” _IEEE/CAA Journal of Automatica Sinica_ , vol. 7, no. 3, pp. 842–852, 2020\.
  * [40] Q. Li, B. Kailkhura, R. Goldhahn, P. Ray, and P. K. Varshney, “Robust decentralized learning using ADMM with unreliable agents,” _IEEE Transactions on Signal Processing_ , vol. 70, pp. 2743–2757, 2022.
  * [41] V. Khatana and M. V. Salapaka, “Noise resilient distributed average consensus over directed graphs,” _IEEE Transactions on Signal and Information Processing over Networks_ , 2023.
  * [42] R. Koenker, _Quantile Regression_ , ser. Econometric Society Monographs. Cambridge University Press, 2005.
  * [43] A. Nedic and A. Ozdaglar, “Distributed subgradient methods for multi-agent optimization,” _IEEE Transactions on Automatic Control_ , vol. 54, no. 1, pp. 48–61, 2009.
  * [44] T. T. Doan, S. T. Maguluri, and J. Romberg, “Convergence rates of distributed gradient methods under random quantization: A stochastic approximation approach,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 10, pp. 4469–4484, 2020.
  * [45] M. M. Vasconcelos, T. T. Doan, and U. Mitra, “Improved convergence rate for a distributed two-time-scale gradient method under random quantization,” in _2021 60th IEEE Conference on Decision and Control (CDC)_. IEEE, 2021, pp. 3117–3122.
  * [46] X. Zhang and M. M. Vasconcelos, “Top-$k$ data selection via distributed sample quantile inference,” in _Learning for Dynamics and Control Conference_. PMLR, 2023, pp. 813–824.
  * [47] D. P. Williamson and D. B. Shmoys, _The design of approximation algorithms_. Cambridge university press, 2011.
  * [48] J. Moon and T. Başar, “Static optimal sensor selection via linear integer programming: The orthogonal case,” _IEEE Signal Processing Letters_ , vol. 24, no. 7, pp. 953–957, 2017.
  * [49] S. M. Kay, _Fundamentals of statistical signal processing: Practical algorithm development_. Pearson Education, 2013, vol. 3.
  * [50] J. C. Duchi, A. Agarwal, and M. J. Wainwright, “Dual averaging for distributed optimization: Convergence analysis and network scaling,” _IEEE Transactions on Automatic control_ , vol. 57, no. 3, pp. 592–606, 2011.
  * [51] A. Makhdoumi and A. Ozdaglar, “Convergence rate of distributed admm over networks,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 10, pp. 5082–5095, 2017.
  * [52] G. Qu and N. Li, “Harnessing smoothness to accelerate distributed optimization,” _IEEE Transactions on Control of Network Systems_ , vol. 5, no. 3, pp. 1245–1260, 2017.
  * [53] A. Beck and M. Teboulle, “Smoothing and first order methods: A unified framework,” _SIAM Journal on Optimization_ , vol. 22, no. 2, pp. 557–580, 2012.
  * [54] X. He, X. Pan, K. M. Tan, and W.-X. Zhou, “Smoothed quantile regression with large-scale inference,” _Journal of Econometrics_ , 2021.
  * [55] L. Xiao and S. Boyd, “Fast linear iterations for distributed averaging,” _Systems & Control Letters_, vol. 53, no. 1, pp. 65–78, 2004.

##  Appendix A Proof of Lemmas

###  A-A Proof of Lemma 1

According to Definition 1, the number of scores equal to $\theta_{k}$ is $m$,
the number of scores larger than $\theta_{k}$ is $k-\overline{m}$ and the
number of scores smaller than $\theta_{k}$ is $n-k-\underline{m}$, where
$\overline{m}+\underline{m}=m$, $\overline{m}\geq 1$ and $\underline{m}\geq
0$. So we have

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}(s_{i}<\theta_{k})=\frac{n-k-\underline{m}}{n}<p$
(82)

and

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}(s_{i}\leq\theta_{k})=\frac{n-k+\overline{m}}{n}>p.$
(83)

Together with the definition of $p$-th sample quantile in (11), we obtain
$\theta_{k}=\omega_{p}$.

###  A-B Proof of Lemma 2

According to Section 1.3 in [42], when $np$ is not an integer, the solution of
Problem (15) is unique. Together with the fact that $f(x)$ is a piecewise
linear function, we can show that if

$\displaystyle|f(x)-f(\theta_{k})|\leq\min\\{g_{r},-g_{l}\\}\cdot\frac{\Delta}{2},$
(84)

then we can guarantee $|x-\theta_{k}|\leq\frac{\Delta}{2}$.

###  A-C Proof of Lemma 3

By direct superposition the left-hand and right-hand derivatives of
$\\{\rho_{p}\big{(}s_{i}-x\big{)}\\}_{i=1}^{n}$ at $\theta_{k}$, we obtain

$\displaystyle g_{r}=$
$\displaystyle\sum_{i=1}^{n}(-p\mathbf{1}(s_{i}>\theta_{k})+(1-p)\mathbf{1}(s_{i}\leq\theta_{k}))$
(85) $\displaystyle=$
$\displaystyle-(k-\overline{m})p+(n-k+\overline{m})(1-p)$ (86)
$\displaystyle=$ $\displaystyle n-k+\overline{m}-np$ (87)

and

$\displaystyle g_{l}=$
$\displaystyle\sum_{i=1}^{n}(-p\mathbf{1}(s_{i}\geq\theta_{k})+(1-p)\mathbf{1}(s_{i}<\theta_{k}))$
(88) $\displaystyle=$
$\displaystyle[-(k+\underline{m})p+(n-k-\underline{m})(1-p)]$ (89)
$\displaystyle=$ $\displaystyle n-k-\underline{m}-np.$ (90)

Incorporating $p=\frac{n-k}{n}+\frac{1}{2n}$ obtains

$g_{r}=\overline{m}-\frac{1}{2},~{}~{}g_{l}=-\underline{m}-\frac{1}{2}.$ (91)

###  A-D Proof of Lemma 4

The conjugate function $\phi$ of $\rho_{p}$ is equivalent to

$\phi(z)=\max_{x\in\mathbb{R}}\\{zx-x(p-\mathbf{1}(x<0))\\}.$ (92)

By simple calculation, we obtain

$\phi(z)=0,\,p-1\leq z\leq p.$ (93)

Notice that

$\rho_{p}^{h}(x)=\max_{z\in\mathbb{R}}\Big{\\{}zx-\phi(z)-\frac{h}{2}z^{2}\Big{\\}}$
(94)

is the conjugate of $\phi(z)+\frac{h}{2}z^{2}$. Since
$\phi(z)+\frac{h}{2}z^{2}$ is $h$-strongly convex, we have $\rho_{p}^{h}(x)$
is $\frac{1}{h}$-smooth. Using (93) gets

$\displaystyle\rho_{p}^{h}(x)=\max_{p-1\leq z\leq
p}\Big{\\{}zx-\frac{h}{2}z^{2}\Big{\\}}$ (95) $\displaystyle=$
$\displaystyle\begin{cases}px-\frac{h}{2}p^{2}&\text{if}\ \ x>hp\\\
\frac{x^{2}}{2h}&\text{if}\ \ h(p-1)\leq x\leq hp\\\
(p-1)x-\frac{h}{2}(p-1)^{2}&\text{if}\ \ x<h(p-1).\\\ \end{cases}$ (96)

###  A-E Proof of Lemma 5

By taking the first-order and second-order derivative of $\rho_{p}^{h}(x)$, we
get

$\displaystyle\nabla\rho_{p}^{h}(x)=\begin{cases}p&\text{if}\ \ x>hp\\\
\frac{x}{h}&\text{if}\ \ h(p-1)\leq x\leq hp\\\ (p-1)&\text{if}\ \ x<h(p-1)\\\
\end{cases}$ (97)

and

$\displaystyle\nabla^{2}\rho_{p}^{h}(x)=\begin{cases}\frac{1}{h}&\text{if}\ \
h(p-1)\leq x\leq hp\\\ 0&\text{otherwise.}\end{cases}$ (98)

For all $x\in\mathbb{R}$, we obtain

$\displaystyle|\nabla\rho_{p}^{h}(x)|\leq\max\\{p,1-p\\}\text{~{}~{}and~{}~{}}0\leq\nabla^{2}\rho_{p}^{h}(x)\leq\frac{1}{h}.$
(99)

Using the fact that
$\tilde{f}_{h}^{\mathrm{nest}}(x)=\sum_{i=1}^{n}\rho_{p}^{h}\big{(}s_{i}-x\big{)}$,
we have

$\displaystyle|\nabla\tilde{f}_{h}^{\mathrm{nest}}(x))|\leq
n\max\\{p,1-p\\}\text{~{}and~{}}0\leq\nabla^{2}\tilde{f}_{h}^{\mathrm{nest}}(x)\leq\frac{n}{h}.$
(100)

So $\tilde{f}_{h}^{\mathrm{nest}}(\cdot)$ is a convex,
$L_{h}=n\max\\{p,1-p\\}$-Lipschitz continuous, $M_{h}=\frac{n}{h}$-smooth
function. By basic calculations, we have

$\displaystyle\rho_{p}^{h}(x)-\rho_{p}(x)=\begin{cases}\frac{h}{2}p^{2}&\text{if}\
\ x>hp\\\ px-\frac{x^{2}}{2h}&\text{if}\ \ 0<x\leq hp\\\
(p-1)x-\frac{x^{2}}{2h}&\text{if}\ \ h(p-1)\leq x\leq 0\\\
\frac{h}{2}(p-1)^{2}&\text{if}\ \ x<h(p-1),\\\ \end{cases}$ (101)

and

$0\leq\rho_{p}^{h}(x)-\rho_{p}(x)\leq\frac{h}{2}\max\\{p^{2},(1-p)^{2}\\}.$
(102)

Using the fact that
$\tilde{f}_{h}^{\mathrm{nest}}(x)=\sum_{i=1}^{n}\rho_{p}^{h}\big{(}s_{i}-x\big{)}$
again yields

$\tilde{f}_{h}^{\mathrm{nest}}(x)\leq
f(x)\leq\tilde{f}_{h}^{\mathrm{nest}}(x)+\frac{nh}{2}\max\\{p^{2},(1-p)^{2}\\},$
(103)

Therefore, $|f(x)-\tilde{f}_{h}^{\mathrm{nest}}(x)|\leq
U_{h}\operatorname{\overset{def}{=}}\frac{nh}{2}\max\\{p^{2},(1-p)^{2}\\}$.

###  A-F Proof of Lemma 6

The first-order and second-order derivative of
$\tilde{f}_{h}^{\mathrm{conv}}(x)$ are

$\displaystyle\nabla\tilde{f}_{h}^{\mathrm{conv}}(x)$
$\displaystyle=\sum_{i=1}^{n}\left[\mathcal{K}_{h}(x-s_{i})-p\right],$ (104)
$\displaystyle\nabla^{2}\tilde{f}_{h}^{\mathrm{conv}}(x)$
$\displaystyle=\sum_{i=1}^{n}K_{h}(x-s_{i}),$ (105)

where

$\displaystyle\mathcal{K}(x)$
$\displaystyle=\int_{-\infty}^{x}K(y){\rm{d}}y\text{~{}~{}and~{}~{}}\mathcal{K}_{h}(x)=\mathcal{K}\left(\frac{x}{h}\right).$
(106)

According to Definition 5, we have

$\displaystyle 0\leq$ $\displaystyle\mathcal{K}_{h}(x)\leq 1,~{}\forall
x\in\mathbb{R},$ (107) $\displaystyle 0\leq$ $\displaystyle
K_{h}(x)\leq\frac{\bar{K}}{h},~{}\forall x\in\mathbb{R}$ (108)

Together with Eqs. (104) and (105), we get

$\displaystyle|\nabla\tilde{f}_{h}^{\mathrm{conv}}(x)|$ $\displaystyle\leq
n\max\\{p,(1-p)\\}\operatorname{\overset{def}{=}}L_{h},$ (109) $\displaystyle
0\leq\nabla^{2}\tilde{f}_{h}^{\mathrm{conv}}(x)$
$\displaystyle\leq\frac{n\overline{K}}{h}\operatorname{\overset{def}{=}}M_{h}.$
(110)

So $\tilde{f}_{h}^{\mathrm{conv}}(\cdot)$ is convex, $L_{h}$-Lipschitz
continuous, $M_{h}$-smooth. Next, let’s give the bound of
$|f(x)-\tilde{f}_{h}^{\mathrm{conv}}(x)|$. First of all, using variable
substitution and Jenson’s inequality yield

$\displaystyle l_{p}^{h}(x)$
$\displaystyle=\int_{-\infty}^{\infty}\rho_{p}(y)K_{h}(y-x){\rm{d}}y$ (111)
$\displaystyle=\int_{-\infty}^{\infty}\rho_{p}(x+t)K_{h}(t){\rm{d}}t$ (112)
$\displaystyle\geq\rho_{p}\left(x+\int_{-\infty}^{\infty}tK_{h}(t){\rm{d}}t\right)$
(113) $\displaystyle=\rho_{p}(x),\forall x\in\mathbb{R},$ (114)

where the last equality uses the even function property of $K_{h}(t)$. So we
have

$\displaystyle\tilde{f}_{h}^{\mathrm{conv}}(x)-f(x)=\sum_{i=1}^{n}\left[l_{p}^{h}\big{(}s_{i}-x\big{)}-\rho_{p}\big{(}s_{i}-x\big{)}\right]\geq
0.$ (115)

Using $\int_{-\infty}^{\infty}K_{h}(t){\rm{d}}t=1$ yields

$\displaystyle l_{p}^{h}(x)-\rho_{p}(x)=$
$\displaystyle\int_{-\infty}^{\infty}\big{(}\rho_{p}(x+t)-\rho_{p}(x)\big{)}K_{h}(t){\rm{d}}t,$
(116) $\displaystyle\leq$
$\displaystyle\int_{-\infty}^{\infty}\big{|}\rho_{p}(x+t)-\rho_{p}(x)\big{|}K_{h}(t){\rm{d}}t$
(117) $\displaystyle\leq$
$\displaystyle\max\\{p,1-p\\}\int_{-\infty}^{\infty}|t|K_{h}(t){\rm{d}}t$
(118) $\displaystyle=$ $\displaystyle
h\max\\{p,1-p\\}\int_{-\infty}^{\infty}|t|K(t){\rm{d}}t,$ (119)

where the first inequality applies that $K_{h}(t)\geq 0$ and
$\rho_{p}(x+t)-\rho_{p}(x)\leq\big{|}\rho_{p}(x+t)-\rho_{p}(x)\big{|}$ for all
$t$, the second inequality uses the property of piecewise linear function
$\big{|}\rho_{p}(x+t)-\rho_{p}(x)\big{|}\leq\max\\{p,1-p\\}|t|$, and the final
equality uses $K_{h}(x)=\frac{1}{h}K\left(\frac{x}{h}\right)$ and variable
substitution. So we have

$\tilde{f}_{h}^{\mathrm{conv}}(x)-f(x)=\sum_{i=1}^{n}\left[l_{p}^{h}\big{(}s_{i}-x\big{)}-\rho_{p}\big{(}s_{i}-x\big{)}\right]\\\
\leq nh\max\\{p,1-p\\}\int_{-\infty}^{\infty}|t|K(t){\rm{d}}t,\forall
x\in\mathbb{R}.$ (120)

Combining Eqs. (115) and (120), we get
$U_{h}=nh\max\\{p,1-p\\}\int_{-\infty}^{\infty}|t|K(t){\rm{d}}t$.

##  Appendix B Extension to multiple data points within each node

Considering the case where the $i$-th node has $n_{i}\geq 1$ scores
$\\{s_{i,j}\\}_{j=1}^{n_{i}}$, our goal is to select the top-$k$ data points
from $N=\sum_{i=1}^{n}n_{i}$, $i=1,\ldots,n$. The objective function of this
problem becomes

$\theta_{k}=\arg\min_{x\in\mathbb{R}}f(x)=\sum_{i=1}^{n}f_{i}(x)=\sum_{i=1}^{n}\sum_{j=1}^{n_{i}}\rho_{p}\big{(}s_{i,j}-x\big{)},$
(121)

where $f_{i}(x)=\sum_{j=1}^{n_{i}}\rho_{p}\big{(}s_{i,j}-x\big{)}$ and
$p\in(\frac{N-k}{N},\frac{N-k+1}{N})$. Actually, this problem with $n_{i}\geq
1$ is equivalent to the problem with $n_{i}=1$, which only rearranges the data
points and increases the number of data points from $n$ to $d$. Therefore, the
analysis and algorithm are also similar with $n_{i}=1$, $i=1,\ldots,n$. Using
the same smoothing technique as described in the manuscript, we have

$\theta_{k}=\arg\min_{x\in\mathbb{R}}\tilde{f}_{h}(x)=\sum_{i=1}^{n}\tilde{f}_{i}^{h}(x)=\sum_{i=1}^{n}\sum_{j=1}^{n_{i}}\rho_{p}^{h}\big{(}s_{i,j}-x\big{)},$
(122)

where
$\tilde{f}_{i}^{h}(x)=\sum_{j=1}^{n_{i}}\rho_{p}^{h}\big{(}s_{i,j}-x\big{)}$.
The Algorithm 1 also works for this case. Besides, we note that for fixed $n$,
the number of transmissions in each round doesn’t change since
$\nabla\tilde{f}_{i}^{h}(x)$ is computed locally. The analysis is without loss
of generality, but the performance of our algorithm will improve, as the
agents engage in communication by sharing only quantile estimates rather than
raw data.

|  Xu Zhang is an assistant professor at the School of Artificial
Intelligence, Xidian University. Previously, he was a Postdoctoral Researcher
at the Academy of Mathematics and Systems Science, Chinese Academy of
Sciences, from 2021 to 2023. He received his B.S. degree and Ph.D. degree in
Electronics Engineering from the School of Information and Electronics,
Beijing Institute of Technology, Beijing, China, in 2015 and 2021. He was a
visiting student in the Ming Hsieh Department of Electrical Engineering at the
University of Southern California from 2018 to 2019. His research interests
include federated learning, distributed optimization, and inverse problems.
---|---
|  Marcos M. Vasconcelos is an Assistant Professor with the Department of
Electrical Engineering at the FAMU-FSU College of Engineering, Florida State
University. He received his Ph.D. from the University of Maryland, College
Park, in 2016. He was a Research Assistant Professor at the Commonwealth Cyber
Initiative and the Bradley Department of Electrical and Computer Engineering
at Virginia Tech from 2021 to 2022. From 2016 to 2020, he was a Postdoctoral
Research Associate in the Ming Hsieh Department of Electrical Engineering at
the University of Southern California. His research interests include
networked control and estimation, multi-scale robotic networks, game theory,
distributed optimization, distributed machine learning, and systems biology.
---|---