In this section, we present a brief introduction to quantile regression. Let S be a scalar random variable, B a L -dimensional random vector, and F S ( s ∣ b ) = P ( S ≤ s ∣ B = b ) represent the conditional cumulative distribution function. The conditional quantile τ is defined as follows: q τ S b = inf s : F S s ∣ b ≥ τ for τ ∈ ( 0,1 ) . A linear model is given by s = b ⊤ w + ϵ , where b is the L × 1 input data vector, w is the L × 1 deterministic unknown parameter vector of interest, and ϵ is the observed noise following a certain distribution. Unlike standard regression methods, which focus on estimating the mean of s , quantile regression provides a more comprehensive analysis by modeling different points (quantiles) in the distribution of s .

In the quantile regression model, q τ S ( b ) is assumed to be linearly related to b as follows: q τ S ( b ) = b ⊤ w + q τ ϵ , where q τ ϵ ∈ R represents the τ -th quantile of the noise. The τ -th quantile of the noise, q τ ϵ , is not necessarily 0 (e.g., for asymmetric noise distributions). Explicitly including q τ ϵ avoids the restrictive assumption that q τ ϵ = 0 , thereby allowing the model to flexibly adapt to the true noise distribution. Omitting q τ ϵ would force the conditional quantile q τ S ( b ) to pass through the origin, which is often inappropriate in real-world applications. Since both q τ ϵ and w are unknown parameters that must be estimated, the optimization problem for estimating w and q τ ϵ can be expressed as follows: w , q τ ϵ = arg min w , q τ ϵ E ρ τ s − b ⊤ w − q τ ϵ . (1) Here, ρ τ is the non-differentiable quantile loss function, also known as the check function, defined as follows: ρ τ u = τ u , if  u ≥ 0 , τ − 1 u , if  u < 0 . This function adjusts the loss asymmetrically depending on whether the residual u is positive or negative, allowing the model to estimate conditional quantiles for different τ values.

Consider a network consisting of K nodes distributed over a certain geographic region. Assume that the network is strongly connected, that is, there is no isolated node in the network. Every node k ∈ { 1,2 , … , K } has access to the realization of zero-mean random data { s k , i , b k , i } at every time instant i and is allowed to communicate only with its neighbors N k , where s k , i is a scalar measurement, b k , i is an L × 1 measurement vector, and N k denotes a set of nodes in the neighborhood of node k including itself. Moreover, { s k , i , b k , i , k = 1 , … , K , i = 1 , … , M } satisfy a standard lineal regression model s k , i = b k , i ⊤ w 0 + ϵ k , i , (2) where ϵ k , i is the measurement noise and w 0 is a deterministic sparse vector of dimension L .

The aim of this study is to develop a distributed quantile regression algorithm to estimate w 0 using the dataset { s k , i , b k , i , k = 1 , … , K , i = 1 , … , M } . Recalling Equation 1, the sparsity-penalized quantile regression estimate of w 0 is obtained by minimizing the global cost function for quantile regression across the network, formulated as follows: min ϖ J glob ϖ ≜ min ϖ 1 K ∑ k = 1 K J k loc ϖ , (3) where J loc ϖ ≜ 1 M ∑ i = 1 M ∑ l ∈ N k c l , k ρ τ s i , l − b i , l ⊤ w − q τ ϵ + λ ‖ w ‖ 1 = 1 M ∑ i = 1 M ∑ l ∈ N k c l , k ρ τ s i , l − a i , l ⊤ ϖ + λ ‖ w ‖ 1 . (4) Here, J glob ( ϖ ) represents the global cost function over the network, and the local cost functions J k loc ( ϖ ) , which reflect the costs at each node k , are aggregated to form the global objective. The first term in Equation 4 captures the quantile regression residuals across all nodes k and data points i , while the second term imposes ℓ 1 -norm regularization on w , encouraging sparsity. In this formulation, a k , i = [ b k , i ⊤ , 1 ] ⊤ is the augmented input data vector, ϖ = [ w ⊤ , q τ ϵ ] ⊤ is the augmented parameter vector to be estimated, and C is a K × K weighting matrix with individual entries { c l , k } . The coefficients c l , k ( l , k = 1,2 , … , K ) are non-negative weighting factors, satisfying c l , k = 0 if l ∉ N k and ∑ l = 1 K c l , k = 1 . (5) This condition (Equation 5) is explicitly satisfied by widely used weight rules in the distributed optimization literature (Cattivelli and Sayed, 2010; Tu and Sayed, 2011). In addition, λ is a regularization parameter controlling the sparsity of the solution, and ‖ w ‖ 1 denotes the ℓ 1 -norm which encourages sparse solutions.

We assume that the underlying network operates under ideal and stable communication conditions. Transient link or node failures—well-studied in the distributed network literature (Gao et al., 2022; Swain et al., 2018)—are effectively handled using established engineering solutions, such as fault-tolerant protocols, redundancy, and consensus mechanisms, thereby ensuring system reliability.

Since the main task of quantile regression is to estimate the parameter vector w 0 , we introduce a diffusion-based distributed estimation framework. In this framework, each node solves its local optimization problem using the subsequently proposed algorithm. The nodes then share their intermediate results with neighboring nodes to collaboratively solve the global quantile regression problem.

By definition in Equation 4, the local cost function of node k can be further expressed as a combination of the local cost functions of its neighboring nodes, i.e., J k loc ( ϖ ) = ∑ l ∈ N k c l , k J ̃ l loc ( ϖ ) , with J ̃ l loc ( ϖ ) = ∑ i = 1 M ρ τ ( s i , l − a i , l ⊤ ϖ ) + λ M w 1 . Therefore, each node obtains its estimate ϖ k by combining its newly generated estimate φ k with the estimates received from its neighboring nodes: ϖ k = ∑ l ∈ N k c l k x l ∈ R L + 1 with x l = arg min x J ̃ l loc x . (6)

Based on this idea, we will use a distributed strategy in the following subsections to solve the problem (Equation 3).

This study focuses on a diffusion-based framework for decentralized quantile regression. The integration of consensus-based strategies with primal–dual hybrid gradient methods for distributed quantile regression, which presents significant algorithmic challenges, is left for future investigation.

By definition, we split J ̃ l loc ( x ) into J ̃ l loc ( x ) = g l ( x ) + λ M f ( x ) with g l ( x ) = ∑ i = 1 M ρ τ ( s i , l − a i , l ⊤ x ) and f ( x ) = ∑ i = 1 L | x i | . For the local optimization problem min x J ̃ l loc ( x ) , the check function ρ τ ( v ) is not differentiable at the origin. To address this, we adopt its conjugate function ρ τ * ( v ) , which allows us to express the problem equivalently. The conjugate function is defined as ρ τ * v = sup u u v − ρ τ u = 0 if  τ − 1 ≤ v < τ , ∞ otherwise . Using the conjugate function ρ τ * ( v ) , we can express ρ τ ( v ) as ρ τ v = sup u u v − ρ τ * u = sup τ − 1 ≤ u < τ u v . This suggests that g ( x ) can be expressed as g l x = max y l ⟨ y l , s l − A l x ⟩ − I τ y l , (7) where A l = a 1 , l , … , a M , l ⊤ , s l = [ s 1 , l , … , s M , l ] ⊤ , y l = [ y 1 , l , … , y M , l ] ⊤ , and I τ ( y ) is an indicator function given by I τ y = 0 , τ − 1 ≤ y i < τ , ∀ i ∞ , else . We can obtain a dual problem associated with the original function minimization problem min x J ̃ l loc ( x ) , which can be expressed as follows: min x k max y k ⟨ y k , s k − A k x k ⟩ − I τ y k + M λ f x k . (8)

This indicates that the global optimization problem (Equation 3) can be formulated as min x k max y k ∑ k = 1 K ∑ l ∈ N k c l , k ⟨ y l , s l − A l x l ⟩ − I τ y l + K M λ f x l , (9) which has a standard saddle-point optimization problem expression with the primal variable x k and the dual variable y k , where k = 1 , … , K . Note that we are primarily concerned with the first L elements of the vector x , which can be interpreted as the estimate of w 0 . Moreover, f ( x ) is defined as f ( x ) = ∑ i = 1 L | x i | .

This section presents the algorithmic principles and derivations of the proposed dPDHG to solve Equation 9. Its basic framework includes the following steps (a)–(d), which involve iterative updates of the primal and dual variables { x k ( n ) , y k ( n ) , k = 1 , … , K , n = 1 , …   } :

(a) Update the primal variable x k ( n ) , k = 0,1 , … , K :

The algorithm continues iterating until the stopping criterion is met, typically based on the difference between successive iterates or the duality gap. Finally, we summarize the proposed dPDHG algorithm in Table 1.

In the above steps, { μ , η } are the primal–dual step sizes chosen to ensure the convergence of the dPDHG algorithm. For the standard form of a convex–concave saddle-point optimization problem min X max Y ⟨ Y , A X ⟩ − g ( Y ) + f ( X ) , the step sizes { μ , η } chosen to ensure convergence of the PDHG algorithm are required μ η Ω 2 < 1 (Esser et al., 2010), where Ω ≔ max Z ∈ R L \ { 0 } ‖ A Z ‖ 2 ‖ Z ‖ 2 = ρ ( A ⊤ A ) . For the problem (Equation 9) mentioned in this paper, A , X , Y can be treated as A = [ s l , − A l ] , X = [ 1 , x l ⊤ ] ⊤ , and Y = y l , respectively.

Note that, for each n , choosing a small μ and large η results in small dual residual d k ( n ) but large primal residual p k ( n ) , and vice versa, where p k n = A k ⊤ △ y k n − 1 μ △ x k n d k n = A k △ x k n − 1 η △ y k n , (10) where △ x k ( n ) = x k ( n ) − x k ( n − 1 ) and △ y k ( n ) = y k ( n ) − y k ( n − 1 ) . Therefore, adaptive strategies—such as those proposed in Goldstein et al. (2015); Chambolle et al. (2024)—can also be employed to balance the progress between the primal and dual updates, thereby enhancing the convergence. If the primal residual p k ( n ) is sufficiently large compared to the dual residual d k ( n ) , for example ‖ p k ( n ) ‖ 2 ≥ 2 ‖ d k ( n ) ‖ 2 , we increase the primal step size μ by a factor of ( 1 − α ) − 1 and decrease the dual stepsize η by a factor of 1 − α . If the primal residual is somewhat smaller than the dual residual, we do the opposite. If both residuals are comparable in size, then let the step sizes remain the same on the next iteration. Moreover, when we modify the step size as the iteration continues, we also shrink the adaptivity level to α ← ζ α , for ζ ∈ ( 0,1 ) .