# Communication-Efficient Modeling with Penalized Quantile Regression for Distributed Data

Authors: Aijun Hu; Chujin Li; Jing Wu
Journal: Complexity
 (2021)

---

## Abstract

In order to deal with high-dimensional distributed data, this article develops a novel and communication-efficient approach for sparse and high-dimensional data with the penalized quantile regression. In each round, the proposed method only requires the master machine to deal with a sparse penalized quantile regression which could be realized fastly by proximal alternating direction method of multipliers (ADMM) algorithm and the other worker machines to compute the subgradient on local data. The advantage of the proximal ADMM algorithm is that it could make every parameter of iteration to have closed formula even in high-dimensional case, which greatly improves the speed of calculation. As for the communication efficiency, the proposed method does not sacrifice any statistical accuracy and provably improves the estimation error obtained by centralized method, provided the penalty levels are chosen properly. Moreover, the asymptotic properties of the proposed estimation and the convergence of the algorithm are convincible. Especially, it presents extensive experiments on both the numerical simulations and the HIV drug resistance data analysis, which all confirm the significant efficiency of our proposed method in quantile regression for distributed data by comparative and empirical analysis.


---

## Body

## 1. Introduction

The quick developments of modern science and technology bring us to the distributed data which are characterized as not only the multiple scales and high dimension, but also the abundant diversity and variability. Distributed data promise and often deliver big dividends and also make diverse challenges to the statistic analysis. One special challenge is that the storage and analysis of such data cannot be performed on a single machine.Since the data are stored in multiple machines, it is natural to consider methods that split the dataset across multiple machines and conduct statistical inference in a distributed manner [1–4]. One popular consideration is the divide-and-conquer (DC) method [5–8], noting the main idea of DC is to take the simple average and only use one round of communication between different data segmentations.However, these averaging-based one-shot communication methods tend to suffer from several drawbacks. Firstly, when the number of machines for data segmentations is large, the local sample size in each machine will be reduced, which could result in the unsatisfactory estimation of each local part, and the final estimation after average would not be accurate. Secondly, when the underlying data model is nonlinear, empirical studies show that the average estimator can only be improved slightly on accuracy as for the local estimators. In order to effectively overcome the above deficiencies, Wang et al. [9] and Jordan et al. [10] proposed the communication-efficient surrogate likelihood (CSL) procedure to solve distributed statistical learning problems. Especially, Jordan et al. [10] pointed out that CSL can replace the total likelihood function to do some statistical inferences. Therefore, it only needs to exchange a number of local data gradients which could effectively reduce the transmission cost of the distributed data. Moreover, the estimation of the method can reach the same convergence rate as the global likelihood-based estimation.The aforementioned studies of DC are mainly focused on linear regression, which only consider the central trend of the conditional distribution of the response variable. As an alternative model, quantile regression proposed by Koenker and Bassett [11] is to analyze the impact of regressors on the conditional distribution of response. Buchinsky [12] captured the heterogeneous impact of regressors on different parts of the distribution, and Koenker [13] exhibited great robustness to outliers by measuring the effect of the regression variable in the upper or lower tails of the distribution. Especially when the error does not follow normal distribution, the quantile regression estimator is more effective and interesting than the least square estimator. Also, quantile regression for high-dimensional data has recently been systematically studied under the sparsity assumption, i.e., only a small number of predictors correlating with the response [14–20]. The quantile regression model selection procedure has been proved to have oracle property under some appropriate penalties [15, 16, 21]. Yu and Lin [22] and Yu et al. [23] used the popular ADMM algorithm to solve the calculation problem of large-scale penalized quantile regression. Recently Chen et al. [24] studied certain high-dimensional distributed quantile regression. The idea is to transform the quantile estimate into a least squares estimate between the transformed response variable and the covariate. However, the calculation needs to estimate the density function of the error and the covariance matrix of the covariate. The main difference from our method is that we approximate the total loss function of the sample with a surrogate loss function, which greatly reduces the amount of calculation.Here, we propose a communication-efficient method with penalized quantile regression on distributed data in this article. The observation data are supposed to be stored randomly in multiple machines. Inspired by Jordan et al. [10] at each round, we only solve the problem on the master machine and calculate the subgradient on the other local machines. After one round of communication, the estimators are passed to each local machine to update the subgradient, and then they are transferred to the master machine for the next round of solution.Importantly, noting that the loss function of quantile regression does not have strong convexity and the subderivative does not satisfy Lipschitz-continuity, the theoretical results by Jordan et al. [10] and Wang et al. [9] cannot be directly applied. Here, under more general conditions, we prove the estimates obtained by our proposed method have oracle properties as presented by Fan and Li [25] and Zou [26]. In terms of algorithms, because the loss function of quantile regression is not smooth, the usual calculation method may not be very accurate. In addition, the traditional ADMM algorithm on iteration of some parameters in high-dimensional quantile regression has no closed-form solution, which could affect the calculation speed and efficiency. Gu et al. [27] proposed the proximal ADMM algorithm to effectively resolve the calculation problem of sparse penalized quantile regression.Hence, we combine it with CSL and propose the proximal ADMM algorithm based on distributed computing framework, which could effectively solve the computational problems of quantile regression on distributed data. Since our method has a display solution for each parameter in each iteration and could be performed point by point, it is very suitable for parallel computing. In each round of communication, only the subgradient calculated by the local data of each machine needs to be transmitted to the host, which greatly saves the cost of data transmission and improves the calculation efficiency.Moreover, the convergence property of the proposed algorithm is still convincible. Considering the extensive numerical simulations and the analysis on HIV Drug Resistance Database, with our high-dimensional penalized quantile regression model, the estimation errors and variable selection results obtained by our proposed approach are compared with those obtained by the centralized method which put all the data on one supercomputer. The proposed method can not only effectively solve the problem by dealing with the heterogeneity of the distributed data but also reduce the cost of storage and transmission. Most importantly, the efficiency of the communication-efficient method with penalized quantile regression on distributed data could match or transcend the centralized method.The rest of this article is organized as follows. We propose the communication-efficient method with penalized quantile regression in the next section and present the proximal ADMM algorithm for solving the communication-efficient penalized sparse quantile regression in Section3. In Section 4, we prove the oracle properties of the proposed estimator. We conduct simulation studies to evaluate the finite-sample performance of our proposed method in Section 5. In Section 6, we demonstrate our method with application to one real data example. Some remarks are concluded in Section 7. The proofs are outlined in Appendix.


## 2. Penalized Distributed Quantile Regression

Letyi,xii=1N be N sample observations, where y is the response variable and x=x1,…,xpT is the p-dimensional covariate; here, p and N are both very large. Assume the data are stored across k machines, and let yji,xji:i=1,…,n denote the subsample stored in the j-th machine ℳj for j=1,…,k, where N=nk denoted the global number of sample observations.Given the covariatex=x1,…,xpT, considering the linear model Y=xTβ+ε, the τth 0<τ<1 quantile of response variable Y is a linear function of the covariate x=x1,…,xpT, that is,(1)QτY|x=xTβ0τ,with Pε≤0|x=τ. Here, the conditional distributed function of Y is written as PY≤y|xi=FYy|xi=Fiy, and so, QτY|xi=Fi−1τ|xi≡ξiτ.Hence, the quantile regression model is to resolve the following optimization problem:(2)β^=argminβ∑i=1Nρτyi−xiTβ,where ρτt=tτ−1t≤0 is the asymmetric absolute deviation function. For ease of the exposition, we simplify the argument τ in β0τ and ξiτ and denote them by β0, ξi, respectively.Define the local and global quantile regression loss function as(3)ℒjβ=1n∑i=1nρτyji−xjiTβ,j=1,…,k,ℒNβ=1k∑j=1kℒjβ=1N∑j=1k∑i=1nρτyji−xjiTβ,where ℒjβ is calculated at β, using the local data stored in j-th machine ℳj.We now apply communication-efficient distributed approach to quantile regression, by proposing the surrogate loss function as(4)ℒ˜β≔ℒ1β−β,∇ℒ1β0−∇ℒNβ0,where β0 is any initial estimator of β, ⋅,⋅ denotes the inner product, and ∇ℒjβ denotes the subgradient of ℒjβ with respect to β.For high-dimensional quantile regression, we consider the following weightedL1-penalized surrogate estimate of quantile regression:(5)minβℒ˜β+∑j=1ppλβj,as a surrogate estimate for β^τ in (2), where λ>0 is the regularization parameter.In the following, we use a convex ALasso (adaptive lasso) penalty and nonconvex SCAD (smoothly clipped absolute deviation) penalty to show the significance of our model. For ALasso penalty,pλβj=wjβj, where w=w1,…,wpT is the vector of nonnegative weight, wj≥0,j=1,…,p. The typical choice of wj is wj=β^jlasso+1/n−ν,j=1,…,p for some appropriately chosen ν>0, where β^lasso=β^jlasso,j=1,…,pT denotes the quantile lasso estimator.Note that (5) is a convex optimization problem, and the calculation is relatively easy. However, if the penalty is nonconvex, such as SCAD [28] case, we have(6)pλβ=λβI0≤β<λ+aλβ−β2+λ2/2a−1⋅Iλ≤β<aλ+a+1λ22Iβ>aλ,for some a>2. A typical choice is a=3.7 as suggested by Fan and Li [25]. Although the nonconvex penalty has a promising theoretical property, the singularity and nonconvexity of the penalty function have led us to encounter several challenges in computing. Zou and Li [29] proposed to replace the nonconvex penalized function by the local linear approximation (LLA).Then, the penalized procedure with both convex penalty and nonconvex penalty can be transformed as the following weighted optimization problem:(7)minβℒ˜β+λw∘β1,where w∘β1=∑j=1pwjβj=∑j=1pwjβj. In the LLA case, we usually take wj=λ−1pλ′β^js−1,j=1,…,p, where β^s−1=β^js−1,j=1,…,pT is the estimation of the s−1th step iteration.


## 3. Algorithm Analysis

Now, we set up to develop a proximal communication-efficient surrogate likelihood ADMM (PCA) algorithm for solving (7). We firstly solve the quantile on the master machine and calculate the subgradient of the local data on other machines. Secondly, we transfer these values back to the master machine and combine them by formula (7), then use Algorithm 1 to calculate the result, which constitutes one communication. Using Algorithm 2 and after only a few rounds of communications T, we can obtain excellent results comparable to the centralized method. The advantage of our algorithm is that each step of Algorithms 1 and 2 has a simple closed-form updata formula, which makes the calculation speed being very fast, especially for the high-dimensional case. Furthermore, Algorithm 2 combined with the CSL idea can achieve very accurate results.Definez=y−Xβ and ℚτz=1/n∑i=1nρτzi, where X=x11,…,x1nT, y=y11,…,y1nT, and zii=1,…,n is the i-th element of z. Hence, we consider the following optimization problem:(8)minβ,zℚτz+gTβ+λw∘β1s.t.Xβ+z=y,where g is a given subgradient vector with g=∇ℒNβ0−∇ℒ1β0, and β0 is any initial estimate of β, and the subgradient of loss function at zero can take any value between two slopes.By direct calculations, we have(9)∇ℒ1β=−1n∑i=1nx1iψτy1i−xi1Tβ,∇ℒNβ=−1N∑j=1k∑i=1nxjiψτyji−xjiTβ,where ψτu=∇ρτu=τIu>0+τ−1Iu<0+ξIu=0,ξ∈τ−1,τ.According to the idea of [10], we calculate the subgradient vectors on each local machine by (9), then calculate (7) for the master machine. Since ℒ˜β here is a nonsmooth function and the direct calculation efficiency is very low, we use ADMM [30] to approximate ℒ˜β as follows. For fixed σ>0, the augmented Lagrangian function of (8) is defined as(10)Lσβ,z,θ=ℚτz+gTβ+λw∘β1−θ,Xβ+z−y+σ2Xβ+z−y22,where θ∈ℝn is the Lagrangian multiplier and ⋅,⋅ and ⋅2 denote the inner product and L2-norm in the Euclidean space, respectively. Denoting βm,zm,θm as the mth iteration of the algorithm for m≥0, the standard ADMM algorithm is given by(11)βstep:βm+1=argminβLσβ,zm,θm=argminβgTβ+λw∘β1−θm,Xβ+σ2Xβ+zm−y22,zstep:zm+1=argminzLσβm+1,z,θm=argminzℚτz−θm,z+σ2Xβm+1+z−y22,θstep:θm+1=θm−γσXβm+1+zm+1−y,where γ is a constant controlling the step length for the θ step.Note thatβ step does not have closed-form formula with a general design matrix X, usually just a numerical solution. Gu et al. [27] proposed a proximal ADMM algorithm for computing the sparse penalized quantile regression. We use this idea to improve our algorithm so that the iteration of each parameter has a closed-form expression, which speeds up the calculation obviously, and the details are as follows. We consider to adding a proximal term to the objective function in the β step and modify the β step in (11) with the following augmented β step.Augmentedβ step:(12)βm+1≔argminβgTβ+λw∘β1−θm,Xβ+σ2Xβ+zm−y22+12β−βmS2,where S is a positive semidefinite matrix.LetS=σηIp−XTX with η≥ΛmaxXTX, where Λmax⋅ denotes the largest eigenvalue of a real symmetric matrix. Here, vS2≔v,Sv is the seminorm induced by the semiinner product via S. In the augmented β step, the update of β can also be carried out component-wisely,(13)βm+1≔argminβgTβ+λw∘β1−θm,Xβ+ση2β−σηβm−g+XTθm+σy−σXβm−σzmση22=Shrinkβjm−gjση+1σηΔj,λwjση1≤j≤p,where Δj=XjTθm+σy−σXβm−σzm,Shinku,α=sgnumaxu−α,0 denotes the soft shrinkage operator with sgn⋅ being the sign function and Xj denotes the jth column of X,j=1,…,p, βjm,andgjj=1,…,p are the jth components of βm and g, respectively.The update ofzk has a closed-form solution and can be implemented component-wisely too. Namely, for i=1,…,n, we have(14)zim+1≔argminzi1nρτzi−θimzi+σ2x1iTβm+1+zi−y1i2=argminziρτzi+nσ2zi−y1i−x1iTβm+1+1σθim2=Proxρτy1i−x1iTβm+1+1σθim,nσ,where the proximal mapping operator Proxρτ⋅,⋅ is given by Gu et al. [27] as follows:(15)Proxρτξ,α≔argminu∈ℝρτu+α2u−ξ2=maxξ−τa,0−max−ξ−1−τa,0.In the following, we present the algorithms. Algorithm1 is to solve the problem on the master machine and the final communication-efficient distributed estimation can be obtained by Algorithm 2.Algorithm 1: pADMM-proximal ADMM algorithm for solving the weighted L1-penalized quantile regression.
Input: data yi,xii=1N. Constants σ>0 and τ>0. Set the maximum number of iterations M.(1)
Form=0,1,2,…,M−1do(2)
Updateβm+1 via (13).(3)
Updatezm+1 via (14).(4)
Updateθm+1 via (11).(5)
End forOutputβM,zM,θM.Algorithm 2: (PCA)-proximal communication-efficient surrogate likelihood ADMM algorithm for solving the weighted L1-penalized quantile regression.
Input: data yi,xii=1N on machine ℳj,j=1,2,…,k. Constants σ>0 and τ>0. Initialize the algorithm with β0,z0,θ0=0,0,0 and g=0.(1)
Setℳ1 as the master machine and obtain the initial value β0,z0,θ0=βM,zM,θM on ℳ1 using Algorithm 1.(2)
Fort=0,1,2,…,T−1do(3)
Master machineℳ1 broadcasts βt to each worker (i.e., ℳ2,…,ℳk).(4)
Each machineℳm calculates subgradient ∇ℒmβt,m=1,…,k and reduces them to master machine ℳ1 to form gt=∇ℒNβt−∇ℒ1βt. Update (8) by replacing g with gt.(5)
Master machineℳ1 updates βt,zt,θt via Algorithm 1.(6)
End forOutput:βT.


## 4. Asymptotic Properties

In this section, for establishing the theoretical properties of the penalized communication-efficient quantile regression, we assume the dataxi,yi,i=1,…,N consist of N observations with the quantile regression model(16)yi=xiTβ0+εi=xi1Tβ10+xi2Tβ20+εi,i=1,…,N,with Pεi≤0|x=τ. Here, xi=xi1T,xi2TT,β0=β10T,β20TT,xi1∈ℝs,xi2∈ℝp−s. The model is sparse, and the true regression coefficient is β10 with each component being nonzero, and β20=0 (as a result β0=β10T,0TT). This means the first s regressors are important while the remaining p−s regressors are noise variables.In the following, we impose the conditions as follows:C1. The distributed functions Fi are absolutely continuous, with continuous densities fiξi uniformly bounded away from 0 and ∞ at the points ξi,i=1,2,….C2. For j=1,…,k, there exists positive definite matrices Σ0 and Σ1τ such that(i)
limn⟶∞n−1∑i=1nxjixjiT=Σ0,(ii)
n−1∑i=1nfiξixjixjiT=Σ1τ,(iii)
maxi=1,…,nxji2/n⟶0.C3.λn⟶0 and nλn⟶∞ as n⟶∞ (SCAD penalty condition).C4.nλn⟶0 and nν+1/2λn⟶∞ for some appropriately chosen ν>0 (ALasso penalty condition).For ease of the exposition, we simplify the argumentτ in Σ1τ and denote it by Σ1.Remark 1.
The condition C1 here is a common condition in general quantile regression. The condition C2 is the technical condition being used in [31]. The conditions C3 and C4 here are the classic conditions for tuning parameters in the SCAD and ALasso penalties, respectively.
Now, we present the asymptotic properties of the communication-efficient distributed quantile regression for the convex penalty (ALasso) and the nonconvex penalty (SCAD).Theorem 1 (consistency).
Consider a samplexi,yi,i=1,…,N from model (16). For the SCAD penalty with conditions C1, C2, and C3, or for the ALasso penalty with conditions C1, C2, and C4, there exists a local minimizer β˜ such that β˜−β0=Opn−1/2.The following theorem gives the asymptotic oracle properties of the estimate.Theorem 2 (oracle).
Consider a samplexi,yi,i=1,…,N from model (16). For the SCAD penalty if conditions C1, C2, and C3 are satisfied, or for the ALasso penalty if conditions C1, C2 and C4 are satisfied, then we have(a)
Sparsity:β˜2=0;(b)
Asymptotic normality:(17)Nβ˜1−β10⟶dN0,DΣ1−1Σ01Σ1−1,where ⟶d means convergence in terms of distribution and Σ01 is defined as the top-left s-by-s submatrix of Σ0, D=kτ1−τ+k−1Varψτε¯−2Covψτε¯,ψτε, ε¯=y−xTβ0, ε=y−xTβ0.Note the initial valueβ0 satisfies β0−β02=Opn−1/2 and the regression error εi is independently and identically distributed, with the τth quantile zero and a continuous density f⋅ in a neighborhood of zero.We could obtain by Theorem2 that D=τ1−τ, and Nβ˜1−β10⟶dN0,τ1−τΣ01−1/f20. In this way, the covariance of the proposed method is exactly same as that of the centralized method [16]. Commonly, since the initial value β0 can be chosen from the first machine satisfying β0−β02=Opn−1/2, this means that our proposed estimate is unbiased and applicable.In the following, we establish the convergence of Algorithm1 and Algorithm 2. Note that g is a constant; by Theorem 1 of Gu et al. [27], the convergence of Algorithm 1 is direct. Also, by choosing the appropriate step length γ, PCA Algorithm 2 can own a sequence βt,zt,t=1,2,… that converges to the global minimizer of problem (7).Theorem 3.
Givenλ>0,σ>0,0<τ<1,0<γ<5+1/2 and a nonnegative weight vector w, let βt,zt,θt be generated by PCA Algorithm 2. Then, the sequence βt,zt,t=0,1,… converges to an optimal solution β∗,z∗ to (8), and θt,t=0,1,… converges to an optimal solution θ∗ to the dual problem of (8). Equivalently, βt,t=0,1,… converges to a global minimizer of problem (7). Moreover, when γ=1, the sequence of norms βt−β∗S2+σzt−z∗22+σ−1θt−θ∗22,t≥0 is nonincreasing and satisfies βt−β∗S2+σzt−z∗22+σ−1θt−θ∗22=Op1/t as t⟶∞.


## 5. Numerical Analysis

In this section, we demonstrate simulation studies to evaluate the finite-sample performance of the communication-efficient distributed quantile regression approach. Since least squares regression has no advantage in dealing with heterogeneous data compared with quantile regression, we only focus on distributed quantile regression model in simulation comparison.Now, we consider the heteroscedastic location-scale model as in Wang et al. [18] for comparing the proposed method with the centralized method which is supposed to store all the data together in one supercomputer. Denote the proposed method and the centralized method with ALasso and SCAD penalties by E-ALasso, E-SCAD, C-ALasso, and C-SCAD, respectively.Let the random vectorsZ1,Z2,…,ZpT generate by the multivariate normal distribution Np0,Σ, with Σ=σijp×p and σij=0.5i−j, and X1=ΦZ1 and Xj=Zj for j=2,3,…,p. The scalar response is generated according to the following heteroscedastic model:(18)Y=X6+X12+X15+X20+0.7X1ε,where ε is distributed as the standard normal distribution N0,1, the student’s distribution t2, and the Laplace distribution L0,1 with the density being 1/2e−x, respectively. In that, these three distributions represent the light or heavy tail characteristic.Here, we setn=100,p=300, the machine number k=10 or 50, and the total sample size N=nk. We randomly divide all the data S into k subsets and denote the data stored in the j-th machine as Sj,j=1,…,k. We use the data Sj in one machine as the test set and the data S−Sj of the remaining k−1 machines as the training set to estimate the parameter β˜jλ. By selecting the tuning parameter λ, the k-fold cross-validation is used to minimize the estimate prediction error in our communication-efficient distributed model. As for each λ, the k-fold cross-validation standard is(19)CVλ=∑j=1k∑xj,yj∈Sjρτyj−xjTβ˜jλ,then it is to resolve λ^=argminλCVλ. In the centralized case, we select the regularization parameters according to the usual 5-fold cross-validation.In each case, we present 100 repeated simulations and compare the results obtained by the two methods mentioned above in accordance with the following four criteria:Size: the average number containing nonzero regression coefficients:β^j≠0,j=1,2,…,p.Pa: the proportion including all true important regression variables, namely, β^j≠0, for any j≥1, satisfying βj≠0. This means the percentage of time including X6,X12,X15,X20, and X1 at τ=0.3 and τ=0.7, while X1 does not have to be included at τ=0.5, if the error distribution is symmetric about the zero.P1: the proportion of simulation running if X1 is selected.AE: the average of absolute estimation error defined by∑j=1pβ^j−βj.Time: the average running time of the CPU in 100 runs.For the two cases ofk=10 and k=50, Tables 1 and 2 , respectively, depict the simulation results of the abovementioned five criteria. From Tables 1 and 2, we can see that for three different types of error distributions and two different numbers of storage machines, our proposed method has a very good performance, and the results could match or transcend the centralized situation. Specifically, the results of the two methods are basically the same on the first three criteria. As for the five criteria, the proposed method is even smaller in the majority of the average absolute error than the centralized method, but the standard deviation does not change much. This can be explained as our method is to reuse the data information of each machine in communication, and similar results were also obtained by Wang et al. [9]. For the error distribution of three different tails and the two different machine numbers, the calculation results of the simulated AE are also very close, which explain the advantages of our method in dealing with variable selection and estimation of distributed data.Table 1
Simulation results fork=10,n=100,p=300,andN=nk.
ετMethodSizePa%P1%AETimeN0,10.3E-ALasso5 (0)1001000.38 (0.03)0.18C-ALasso5 (0)1001000.40 (0.03)0.21E-SCAD5 (0)1001000.22 (0.05)0.18C-SCAD5 (0)1001000.33 (0.03)0.210.5E-ALasso4 (0)10000.03 (0.01)0.18C-ALasso4 (0)10000.04 (0.02)0.21E-SCAD4 (0)10000.03 (0.01)0.18C-SCAD4 (0)10000.07 (0.02)0.210.7E-ALasso5 (0)1001000.35 (0.03)0.18C-ALasso5 (0)1001000.40 (0.04)0.21E-SCAD5.01 (0.10)1001000.22 (0.04)0.18C-SCAD5 (0)1001000.31 (0.03)0.21t20.3E-ALasso5 (0)1001000.45 (0.04)0.19C-ALasso5 (0)1001000.47 (0.04)0.21E-SCAD5.02 (0.14)1001000.31 (0.06)0.19C-SCAD5 (0)1001000.34 (0.04)0.210.5E-ALasso4 (0)10000.03 (0.01)0.18C-ALasso4 (0)10000.04 (0.01)0.21E-SCAD4 (0)10000.05 (0.02)0.19C-SCAD4 (0)10000.04 (0.01)0.210.7E-ALasso5 (0)1001000.44 (0.04)0.19C-ALasso5 (0)1001000.46 (0.04)0.23E-SCAD5 (0)1001000.31 (0.06)0.19C-SCAD5 (0)1001000.35 (0.04)0.21L0,10.3E-ALasso5 (0)1001000.37 (0.04)0.19C-ALasso5 (0)1001000.39 (0.04)0.21E-SCAD5.07 (0.26)1001000.27 (0.06)0.18C-SCAD5 (0)1001000.32 (0.04)0.210.5E-ALasso4 (0)10000.03 (0.01)0.17C-ALasso4 (0)10000.04 (0.01)0.19E-SCAD4 (0)10000.04 (0.02)0.18C-SCAD4 (0)10000.05 (0.02)0.210.7E-ALasso5 (0)1001000.37 (0.04)0.18C-ALasso5 (0)1001000.40 (0.04)0.21E-SCAD5 (0)1001000.19 (0.04)0.18C-SCAD5 (0)1001000.30 (0.04)0.21Size, the average number containing nonzero regression coefficients;Pa, the proportion including all true important regression variables; P1, the proportion of simulation run X1 is selected; AE, the average of absolute estimation error, the numbers in parentheses are the corresponding standard deviations; Time, CPU average running time. L0,1, Laplace0,1 distribution.Table 2
Simulation results fork=50,n=100,p=300,andN=nk.
ετMethodSizePa%P1%AETimeN0,10.3E-ALasso5 (0)1001000.36 (0.02)0.37C-ALasso5 (0)1001000.37 (0.02)1.24E-SCAD5.05 (0.26)1001000.27 (0.05)0.37C-SCAD5 (0)1001000.32 (0.02)1.240.5E-ALasso4 (0)10000.03 (0.01)0.36C-ALasso4 (0)10000.03 (0.01)1.24E-SCAD4 (0)10000.03 (0.01)0.37C-SCAD4 (0)10000.01 (0.00)1.240.7E-ALasso5 (0)1001000.35 (0.02)0.36C-ALasso5 (0)1001000.38 (0.02)1.24E-SCAD5 (0)1001000.27 (0.04)0.36C-SCAD5 (0)1001000.32 (0.02)1.23t20.3E-ALasso5 (0)1001000.43 (0.02)0.36C-ALasso5 (0)1001000.46 (0.02)1.21E-SCAD5.03 (0.17)1001000.39 (0.05)0.36C-SCAD5 (0)1001000.39 (0.02)1.210.5E-ALasso4 (0)10000.02 (0.01)0.36C-ALasso4 (0)10000.04 (0.01)1.21E-SCAD4 (0)10000.03 (0.01)0.36C-SCAD4 (0)10000.02 (0.01)1.210.7E-ALasso5 (0)1001000.43 (0.02)0.36C-ALasso5 (0)1001000.45 (0.03)1.19E-SCAD5 (0)1001000.29 (0.04)0.35C-SCAD5 (0)1001000.38 (0.03)1.19L0,10.3E-ALasso5 (0)1001000.34 (0.02)0.34C-ALasso5 (0)1001000.37 (0.02)1.13E-SCAD5.02 (0.14)1001000.25 (0.05)0.36C-SCAD5 (0)1001000.30 (0.02)1.250.5E-ALasso4 (0)10000.03 (0.02)0.33C-ALasso4 (0)10000.03 (0.01)1.10E-SCAD4 (0)10000.03 (0.02)0.33C-SCAD4 (0)10000.01 (0.00)1.110.7E-ALasso5 (0)1001000.36 (0.02)0.34C-ALasso5 (0)1001000.36 (0.02)1.11E-SCAD5.01 (0.10)1001000.25 (0.04)0.37C-SCAD5 (0)1001000.29 (0.02)1.25Size, the average number containing nonzero regression coefficients;Pa, the proportion including all true important regression variables; P1, the proportion of simulation run X1 is selected; AE, the average of absolute estimation error, the numbers in parentheses are the corresponding standard deviations; Time, CPU average running time. L0,1, Laplace0,1 distribution.Tables1 and 2 show the comparison of CPU running time for 10 rounds of communication. It is clear that when the number of machines is 10, the running times of the two methods are similar. However, when the number of machines is increased to 50, the running time of our method is significantly shorter than that of the oracle method, which fully demonstrates the effective communication superiority of our method for dealing with big data.Figure1 plots how the estimation error βt−β∗2 varies for the proposed method with the ALasso penalty when k=50, but the estimation error of the centralized method looks as a horizontal line. In addition, the results of k=10 and SCAD penalty are similar to Figure 1, which we will omit here. Moreover, by Figure 1, the estimation error obtained by our method decreases to be truly competitive with the centralized method within very few rounds of communications. Usually, it just requires less than 5 communications and could achieve smaller prediction error than the centralized method. This phenomenon also appeared in [9], which can be explained as our method is to use data information repeatedly.Figure 1
Comparison of two methods for the estimation error based onk=50; the sparse ALasso penalized quantile regression. L0,1 denote the Laplace0,1 distribution. (a) N0,1,τ=0.3. (b) N0,1,τ=0.5. (c) N0,1,τ=0.7. (d) t2,τ=0.3. (e) t2,τ=0.5. (f) t2,τ=0.7. (g) L0,1,τ=0.3. (h) L0,1,τ=0.5. (i) L0,1,τ=0.7.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 6. Study on HIV Drug Resistance

As an illustration, we apply the communication-efficient penalized quantile regression model on a human immunodeficiency virus (HIV) drug susceptibility dataset from the HIV Drug Resistance Database (http://hivdb.stanford.edu), a public resource for the study of sequence variation and mutations in the molecular targets of HIV drugs [32]. The HIV data are also available from Wang et al. [33]. When a patient begins antiretroviral therapy, infected HIV may form new mutations. Some mutations may not respond to existing drugs, a feature known as drug resistance or reduced drug sensitivity which means that the drug’s role in preventing viral reproduction is diminished, and researchers have estimated that in an untreated HIV-infected person, every possible single-point mutation occurs 104 to 105 times a day [34]. Drug resistance has become a major obstacle to the treatment of HIV. Therefore, understanding the impact of mutations on drug resistance is an important research topic.We analyze the susceptibility data for the drug efavirenz (EFV). After excluding some rare mutations, the dataset includes 1472 HIV isolates and 197 locations of mutations. Although the sample size is not too big, the data background traits conform to the distributed effect extremely. The susceptibility of an HIV sample is defined as the fold decrease in susceptibility of a single virus isolate compared with the susceptibility of a wild type control isolate, that is, the virus that has never been challenged by drugs. We focus on predicting thelog10-susceptibility, denoted by Y, related to the EFV based on Xk,k=1,…,p=197, where Xk indicates the presence of a mutation of interest in the k-th viral sequence position. It is noted that the susceptibility data are often stored in different locations and highly nonnormal, even after logarithmic transformation. Therefore, communication-efficient quantile regression could provide a valuable method for analyzing these data. In practice, the location of virus mutation is usually very scarce, so we would get a sparse solution by adding penalties. Note the analysis of the upper quantile of sensitivity is particularly important, which is related to stronger drug resistance. Therefore, we consider two quantile levels: τ=0.5 and 0.75.In our analysis, we use the ALasso and SCAD penalized quantile regression model with the proposed approach and the centralized approach to select those mutations that are interested. We conduct 20 random partitions. For each partition, we randomly select 1200 HIV isolates as the training data and the other 272 as the testing data. We set 1200 training data on 12 machines and store 100 data on each machine (n=100,k=12). A twelve-fold cross-validation is applied to the training data to select the tuning parameters. But, we use the usual 5-fold cross-validation to select it, while using the centralized method. Table 3 records the average number of nonzero regression coefficients (ave # nonzero) and evaluates the average prediction error using the quantile loss function with τ=0.5 and τ=0.75, where the prediction error is defined as 1/272∑i=1272ρτyi−y^i and the numbers in the parentheses are the corresponding standard errors across 20 random partitions. Table 4 lists the frequency of the important variables selected by the various methods for 20 random partitions. Figure 2 plots how the prediction error varies for the proposed method with the rounds of communication, but a horizontal line for centralized method. Moreover, from Figure 2, we can see that as long as a certain number of communications, the prediction error of our method can be very close to even smaller results than the centralized method, which shows our method can match the performance of the centralized method. This phenomenon is consistent with the results of simulated data.Table 3
Analysis of HIV drug susceptibility dataset. The numbers in parentheses are the corresponding standard deviations.
Methodτ=0.5τ=0.75AnonPEAnonPEE-ALasso24.60 (2.19)0.16 (0.01)28.65 (2.43)0.15 (0.01)C-ALasso24.25 (1.83)0.20 (0.01)27.25 (2.47)0.17 (0.01)E-SCAD27.05 (2.54)0.17 (0.01)30.25 (2.27)0.17 (0.01)C-SCAD25.15 (3.20)0.20 (0.01)29.65 (3.23)0.17 (0.01)Anon and PE represent the number of average nonzero variables and prediction error, respectively.Table 4
Frequency table for the HIV data.
Variableτ=0.5Variableτ=0.75CALEALCSESCALEALESCSX.230L20202020X.230L20202020X.227L20202018X.227L20202020X.225H20202020X.225H20202020X.190A20202020X.190A20202020X.190S20202020X.190S20202020X.188L20202020X.188L20202020X.179D20202020X.179D20202020X.108I20202020X.108I20202020X.103N20202020X.103N20202020X.101H20202019X.102Q20202020X.101Q20202020X.101Q20202020X.101P20202020X.101P20202020X.101E20202020X.101E20202020X.100I20202020X.100I20202020X.90I20202020X.90I20192020X.219N19171817X.181C19191820X.98G.119201920X.221Y18181418X.98G19201920X.101H18201918X.221Y17191720X.219N16131817X.215Y16151412X.219E16171512X.138A13141311X.103R15201817X.181C11201320X.98G.115161718X.219E1051014X.98G15161718X.135L104139X.74V1451719X.172K7111012X.65R1401210X.172K112136X.67G1021112X.189I9131617X.102R8171212X.106I7692Here, CAL, EAL, CS, and ES representC-ALasso, E-ALasso, C-SCAD, and E-SCAD, respectively.Figure 2
Comparison of two methods for the prediction error with the HIV data. (a)ALasso,τ=0.5. (b) ALasso,τ=0.75. (c) Scad,τ=0.5. (d) Scad,τ=0.75.
(a)(b)(c)(d)From Table3, we can see the average number of variables selected is very close to 25 when τ=0.5; however, when τ=0.75, it is close to 30. This shows that the dataset has a certain heteroscedasticity. From the perspective of the prediction error, the proposed method yields a smaller prediction error, which further demonstrates the good performance of our method in the HIV data.As we can see from Table4, the variables selected using our method are very consistent with those selected using the centralized method. Similar phenomena also appear for the case τ=0.75 in Table 4. But, the ALasso penalized communication-efficient distributed method performs poor in variables X.172K, X.65R, and X.67G; for other variables, the abovementioned two methods are very perfect. However, the variables selected by the SCAD penalized communication-efficient distributed method are highly consistent with the variables selected by the centralized method. These results show that our method performs well in variable selection in the HIV drug susceptibility dataset.


## 7. Discussion

We propose a new communication-efficient distributed method to solve sparse penalized quantile regression on massive data, with the parameter estimation obtained by the proposed method having oracle properties. In terms of computation, with a proximal ADMM algorithm under the distributed framework, which makes every parameter of iteration have closed formulas, we also establish the convergence of the algorithm which owns nice efficiency and accuracy. The computational efficiency and accuracy of the proposed method are verified by extensive numerical simulation. The simulation results show that our method usually takes only a few communications to achieve priority performance compared with the centralized approach even having the heteroscedasticity, and the real data example demonstrates that our proposed method has fine performance in parameter estimation and variable selection in high-dimensional quantile regression, especially for multiple and distributed data with heterogeneity and outliers.It is noted that our method and the PCA algorithm can be modified to the quantile regression with the elastic net penalty, and the actual implementation of our proposed PCA algorithm can be further improved to make it more efficient in communication. For example, in each iteration of (13), (14), and (11), we have to update all the data on each machine to complete the iteration before it is done. In addition, due to the synchronization in (13), (14), and (11), the total computation speed must be limited by the slowest computing machine. Zhang and Kwok [35] proposed an asynchronous ADMM algorithm to resolve the problem. Moreover, if the distributed data are ultrahigh-dimensional, one necessary choice is to divide the data along the p direction, which may be of further research interest.



---

*Source: 6341707-2021-01-16.xml*