Least Squares Approximation for a Distributed System

TL;DR

This paper introduces a distributed least squares approximation method that efficiently solves various regression problems on large datasets with minimal communication, ensuring statistical efficiency and computational speed.

Contribution

It proposes a novel distributed regression approach using local quadratic approximations, achieving statistical efficiency with only one communication round and integrating adaptive Lasso for variable selection.

Findings

Achieves statistical efficiency comparable to global estimators.

Requires only one round of communication in distributed settings.

Demonstrates computational efficiency on large-scale datasets.

Abstract

In this work, we develop a distributed least squares approximation (DLSA) method that is able to solve a large family of regression problems (e.g., linear regression, logistic regression, and Cox's model) on a distributed system. By approximating the local objective function using a local quadratic form, we are able to obtain a combined estimator by taking a weighted average of local estimators. The resulting estimator is proved to be statistically as efficient as the global estimator. Moreover, it requires only one round of communication. We further conduct a shrinkage estimation based on the DLSA estimation using an adaptive Lasso approach. The solution can be easily obtained by using the LARS algorithm on the master node. It is theoretically shown that the resulting estimator possesses the oracle property and is selection consistent by using a newly designed distributed Bayesian information criterion (DBIC). The finite sample performance and computational efficiency are further illustrated by an extensive numerical study and an airline dataset. The airline dataset is 52 GB in size. The entire methodology has been implemented in Python for a {\it de-facto} standard Spark system. The proposed DLSA algorithm on the Spark system takes 26 minutes to obtain a logistic regression estimator, which is more efficient and memory friendly than conventional methods.

Tables

Table 1. Table 1: Simulation results for Example 1 (Setting 2) with 500 replications. The numerical performances are evaluated for different sample sizes N 𝑁 N ( × 10 3 absent superscript 10 3 \times 10^{3} ) and numbers of workers K 𝐾 K . The REE is reported for all estimators. For the CSL, DLSA, TDLSA method, the CP is further reported in parentheses. Finally, the MS and CM are reported for the SDLSA method to evaluate the variable selection accuracy.

Est.	${\theta }_1$	${\theta }_2$	${\theta }_3$	${\theta }_4$	${\theta }_5$	${\theta }_6$	${\theta }_7$	${\theta }_8$	MS	CM
$N=20$, $K=10$
OS	0.98	1.00	0.99	1.00	1.00	0.99	0.99	0.98
CSL	0.94	0.93	0.94	0.94	0.95	0.96	0.93	0.96
	(94.2)	(93.6)	(95.4)	(92.2)	(92.8)	(92.6)	(91.0)	(91.2)
DLSA	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	(96.0)	(94.8)	(96.8)	(94.0)	(94.6)	(94.0)	(94.0)	(93.0)
TDLSA	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	(96.0)	(94.8)	(96.8)	(94.0)	(94.6)	(94.0)	(94.0)	(93.0)
SDLSA	1.08	-	-	1.13	-	-	1.19	-	3.01	1.00
$N=100$, $K=20$
OS	1.00	0.99	1.00	1.00	1.00	0.99	1.00	0.99
CSL	0.93	0.91	0.92	0.91	0.93	0.92	0.92	0.92
	(93.2)	(91.0)	(92.0)	(92.6)	(92.0)	(91.6)	(92.4)	(90.6)
DLSA	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	(95.4)	(95.4)	(94.6)	(94.6)	(95.6)	(93.8)	(94.4)	(94.6)
TDLSA	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	(95.4)	(95.4)	(94.6)	(94.6)	(95.6)	(93.8)	(94.4)	(94.6)
SDLSA	1.10	-	-	1.13	-	-	1.16	-	3.01	1.00

Table 2. Table 2: Simulation results for Example 2 (Setting 2) with 500 replications. The numerical performances are evaluated for different sample sizes N 𝑁 N ( × 10 3 absent superscript 10 3 \times 10^{3} ) and numbers of workers K 𝐾 K . The REE is reported for all estimators. For the CSL, DLSA, TDLSA method, the CP is further reported in parentheses. Finally, the MS and CM are reported for the SDLSA method to evaluate the variable selection accuracy.

Est.	${\theta }_1$	${\theta }_2$	${\theta }_3$	${\theta }_4$	${\theta }_5$	${\theta }_6$	${\theta }_7$	${\theta }_8$	MS	CM
$N=20$, $K=10$
OS	0.76	0.87	0.87	0.79	0.89	0.91	0.78	0.88
CSL	0.15	0.17	0.18	0.16	0.19	0.17	0.16	0.17
	(40.8)	(43.0)	(42.6)	(42.0)	(44.8)	(42.0)	(43.6)	(40.0)
DLSA	0.88	1.01	1.01	0.94	1.01	1.01	0.90	1.01
	(90.0)	(95.8)	(95.8)	(92.6)	(96.0)	(96.0)	(92.4)	(95.0)
TDLSA	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	(95.6)	(96.0)	(95.6)	(95.0)	(95.6)	(95.4)	(95.0)	(95.0)
SDLSA	1.02	-	-	1.09	-	-	1.05	-	3.00	1.00
$N=100$, $K=20$
OS	0.77	0.95	0.95	0.84	0.95	0.93	0.82	0.93
CSL	0.12	0.13	0.13	0.12	0.14	0.13	0.12	0.12
	(31.8)	(29.0)	(29.2)	(29.8)	(30.0)	(26.2)	(29.4)	(29.2)
DLSA	0.92	1.00	1.01	0.94	1.00	1.01	0.92	1.01
	(92.6)	(94.6)	(94.2)	(93.2)	(93.2)	(96.0)	(92.8)	(95.4)
TDLSA	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	(95.6)	(94.8)	(94.2)	(94.0)	(93.4)	(95.6)	(94.6)	(95.4)
SDLSA	1.02	-	-	1.11	-	-	1.07	-	3.00	1.00

Table 3. Table 3: Simulation results for Example 3 (Setting 2) with 500 replications. The numerical performances are evaluated for different sample sizes N 𝑁 N ( × 10 3 absent superscript 10 3 \times 10^{3} ) and numbers of workers K 𝐾 K . The REE is reported for all estimators. For the CSL, DLSA, TDLSA method, the CP is further reported in parentheses. Finally, the MS and CM are reported for the SDLSA method to evaluate the variable selection accuracy.

Est.	${\theta }_1$	${\theta }_2$	${\theta }_3$	${\theta }_4$	${\theta }_5$	${\theta }_6$	${\theta }_7$	${\theta }_8$	MS	CM
$N=20$, $K=10$
OS	0.76	0.87	0.87	0.79	0.89	0.91	0.78	0.88
CSL	0.15	0.17	0.18	0.16	0.19	0.17	0.16	0.17
	(40.8)	(43.0)	(42.6)	(42.0)	(44.8)	(42.0)	(43.6)	(40.0)
DLSA	0.88	1.01	1.01	0.94	1.01	1.01	0.90	1.01
	(90.0)	(95.8)	(95.8)	(92.6)	(96.0)	(96.0)	(92.4)	(95.0)
TDLSA	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	(95.6)	(96.0)	(95.6)	(95.0)	(95.6)	(95.4)	(95.0)	(95.0)
SDLSA	1.02	-	-	1.09	-	-	1.05	-	3.00	1.00
$N=100$, $K=20$
OS	0.77	0.95	0.95	0.84	0.95	0.93	0.82	0.93
CSL	0.12	0.13	0.13	0.12	0.14	0.13	0.12	0.12
	(31.8)	(29.0)	(29.2)	(29.8)	(30.0)	(26.2)	(29.4)	(29.2)
DLSA	0.92	1.00	1.01	0.94	1.00	1.01	0.92	1.01
	(92.6)	(94.6)	(94.2)	(93.2)	(93.2)	(96.0)	(92.8)	(95.4)
TDLSA	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	(95.6)	(94.8)	(94.2)	(94.0)	(93.4)	(95.6)	(94.6)	(95.4)
SDLSA	1.02	-	-	1.11	-	-	1.07	-	3.00	1.00

Table 4. Table 4: Variable description for the U.S. airline data. Non-categorical numerical variables are standardized to have mean zero and variance one.

Variable	Description	Variable used in the model
Delayed	Whether the flight is delayed, 1 for Yes; 0 for No.	Used as the response variable
Year	Year between 1987 and 2008	Used as numerical variable
Month	Which month of the year	Converted to $11$ dummies
DayofMonth	Which day of the month	Used as numerical variable
DayofWeek	Which day of the week	Converted to $6$ dummies
DepTime	Actual departure time	Used as numerical variable
CRSDepTime	Scheduled departure time	Used as numerical variable
CRSArrTime	Scheduled arrival time	Used as numerical variable
ElapsedTime	Actual elapsed time	Used as numerical variable
Distance	Distance between the origin and destination in miles	Used as numerical variable
Carrier	Flight carrier code for $29$ carriers	Top $7$ carries converted to $7$ dummies
Destination	Destination of the flight (total $348$ categories)	Top $75$ destination cities converted to $75$ dummies
Origin	Departing origin (total $343$ categories)	Top $75$ origin cities converted to $75$ dummies

Table 5. Table 5: Log likelihood, computational time and communication time comparison with the airline data. Note that the inevitable overhead for initializing the spark task is also included in the computational time. It is difficult to monitor the exact communication time per code interaction in Spark due to its lazy evaluation technique used. The overall communication time is estimated from the computational time by subtracting the computational time used on the master machine and the average computational time used on the workers.

Algorithm	Log likelihood	Computational time (min)	Communication time (min)	Worker’s RAM	No. of workers
DLSA	$-1.62\times {10}^8$	26.2	3.8	64GB	2
TDLSA	$-1.62\times {10}^8$	28.7	4.4	64GB	2
OS	$-1.65\times {10}^8$	25.3	2.3	64GB	2
Spark ML	$-1.63\times {10}^8$	54.7	18.6	128GB	2
Serialized SGD	$-1.71\times {10}^8$	913.2	—	64GB	1

Table 6. Table 6: Coefficients for the logistic model estimated with SDLSA using the BIC. The carrier and airport abbreviations are assigned by the International Air Transport Association. We denote “Airport”, “Origin” and “Destination” by “A”, “O” and “D”, respectively. The variances for all coefficients are all smaller than 0.001 0.001 0.001 . The notation *** indicates 0.001 level of significance.

	Intercept	Year	CRSArrTime	Distance	CRSDepTime	DayofMonth	ElapsedTime	DepTime
	$-0.30$	${6.12}^{\ast \ast \ast }$	$-0.13$	$-{5.68}^{\ast \ast \ast }$	$-0.12$	$-0.06$	$0.08$	$0.15$
Month	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
	$-0.13$	$-0.01$	$0.03$	${6.03}^{\ast \ast \ast }$	$-0.28$	$-0.09$	$-0.02$	$-0.07$	${0.8}^{\ast \ast \ast }$	$0.4$	$0.19$
Day of Week	Tue	Wed	Thu	Fri	Sat	Sun
	${0.6}^{\ast \ast \ast }$	${0.85}^{\ast \ast \ast }$	$-0.57$	$0.39$	$0.22$	$0.26$
Carrier	AA	CO	DL	NW	UA	US	WN
	${0.49}^{\ast \ast \ast }$	$-0.16$	$0.18$	$0.39$	${0.7}^{\ast \ast \ast }$	$-0.43$	$-0.6$
A	ABQ	ANC	ATL	AUS	BDL	BHM	BNA	BOS	BUF	BUR	BWI
O	$-0.46$	$-0.14$	$-0.02$	$0.39$	$0.38$	$0.13$	$-0.17$	$0.13$	$-0.55$	$-0.01$	$0.54$
D	$-0.61$	$0.78$	$-0.84$	$-0.73$	$-0.74$	$-0.46$	$-0.28$	$0.51$	$-0.50$	$-0.90$	$0.69$
A	CLE	CLT	CMH	CVG	DAL	DAY	DCA	DEN	DFW	DTW	ELP
O	$0.49$	$-0.87$	$-0.80$	$0.29$	$-0.08$	$0.40$	$-0.07$	$-0.22$	$0.53$	$0.32$	$-0.56$
D	$0.14$	$-0.46$	$0.64$	$-1.60$	$-0.72$	$-0.79$	${1.07}^{\ast \ast \ast }$	$0.37$	$-1.04$	$-0.19$	$-0.54$
A	EWR	FLL	GSO	HNL	HOU	IAD	IAH	IND	JAX	JFK	LAS
O	$0.21$	$-0.99$	$0.11$	$0.19$	$-0.51$	$-0.29$	${0.82}^{\ast \ast \ast }$	$0.46$	$-0.71$	$0.07$	$-0.94$
D	$-0.44$	$-1.70$	$0.77$	$0.79$	$-1.62$	$0.00$	$-0.54$	$-0.27$	$-0.12$	$-0.70$	$-0.29$
A	LAX	LGA	MCI	MCO	MDW	MEM	MIA	MKE	MSP	MSY	OAK
O	$0.38$	${0.87}^{\ast \ast \ast }$	$0.61$	$0.52$	$-0.48$	$0.71$	$0.39$	$-0.64$	$-0.10$	$0.21$	$-1.04$
D	$-0.20$	$0.13$	${1.00}^{\ast \ast \ast }$	$0.30$	$0.10$	${1.16}^{\ast \ast \ast }$	$-0.11$	$-0.83$	$-0.61$	$0.97$	$-1.47$
A	OKC	OMA	ONT	ORD	ORF	PBI	PDX	PHL	PHX	PIT	PVD
O	$-0.71$	$-0.85$	${1.58}^{\ast \ast \ast }$	$-0.41$	$-0.67$	$-1.12$	$-0.13$	${0.94}^{\ast \ast \ast }$	$0.37$	$0.18$	$-0.42$
D	$0.69$	$-0.74$	$-1.49$	$-0.72$	$-0.79$	${0.99}^{\ast \ast \ast }$	$-0.64$	$-0.54$	$0.83$	$0.33$	$-0.40$
A	RDU	RIC	RNO	ROC	RSW	SAN	SAT	SDF	SEA	SFO	SJC
O	$-5.80$	$0.02$	$-0.06$	$0.20$	$0.05$	$0.00$	$-0.04$	$0.00$	$0.24$	$0.24$	$0.14$
D	$0.97$	$0.83$	$0.71$	$-1.07$	$0.37$	${1.15}^{\ast \ast \ast }$	$0.78$	$-0.35$	$-0.79$	$-0.51$	$0.21$
A	SJU	SLC	SMF	SNA	STL	SYR	TPA	TUL	TUS
O	$-0.14$	$0.59$	${1.59}^{\ast \ast \ast }$	$-0.15$	$0.14$	$-1.21$	$-0.28$	$-0.36$	$-1.51$
D	$0.84$	$0.09$	$-0.03$	$0.02$	$0.38$	$-0.06$	$0.31$	$0.43$	$0.03$

Table 7. Table 7: Simulation results for Example 1 (Setting 1) with 500 replications. The numerical performances are evaluated for different sample sizes N 𝑁 N ( × 10 3 absent superscript 10 3 \times 10^{3} ) and numbers of workers K 𝐾 K . The REE is reported for all estimators. For the CSL, DLSA, TDLSA method, the CP is further reported in parentheses. Finally, the MS and CM are reported for the SDLSA method to evaluate the variable selection accuracy.

Est.	${\theta }_1$	${\theta }_2$	${\theta }_3$	${\theta }_4$	${\theta }_5$	${\theta }_6$	${\theta }_7$	${\theta }_8$	MS	CM
$N=20$, $K=10$
OS	1.00	1.00	1.00	1.00	1.00	1.00	1.00	0.99
CSL	0.96	0.93	0.95	0.93	0.95	0.96	0.93	0.96
	(94.0)	(94.0)	(94.8)	(92.8)	(92.8)	(93.2)	(90.8)	(90.4)
DLSA	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	(96.4)	(96.0)	(96.4)	(94.6)	(94.6)	(94.4)	(94.0)	(93.8)
TDLSA	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	(96.4)	(96.0)	(96.4)	(94.6)	(94.6)	(94.4)	(94.0)	(93.8)
SDLSA	1.03	-	-	1.07	-	-	1.08	-	3.01	1.00
$N=100$, $K=20$
OS	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
CSL	0.95	0.92	0.93	0.94	0.95	0.93	0.92	0.96
	(93.8)	(92.2)	(91.4)	(92.6)	(92.2)	(92.2)	(92.0)	(94.6)
DLSA	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	(95.4)	(94.8)	(94.8)	(95.6)	(95.2)	(94.0)	(94.6)	(95.2)
TDLSA	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	(95.4)	(94.8)	(94.8)	(95.6)	(95.2)	(94.0)	(94.6)	(95.2)
SDLSA	1.03	-	-	1.05	-	-	1.08	-	3.01	1.00

Table 8. Table 8: Simulation results for Example 2 (Setting 1) with 500 replications. The numerical performances are evaluated for different sample sizes N 𝑁 N ( × 10 3 absent superscript 10 3 \times 10^{3} ) and numbers of workers K 𝐾 K . The REE is reported for all estimators. For the CSL, DLSA, TDLSA method, the CP is further reported in parentheses. Finally, the MS and CM are reported for the SDLSA method to evaluate the variable selection accuracy.

Est.	${\theta }_1$	${\theta }_2$	${\theta }_3$	${\theta }_4$	${\theta }_5$	${\theta }_6$	${\theta }_7$	${\theta }_8$	MS	CM
$N=20$, $K=10$
OS	0.84	0.98	0.99	0.90	0.98	0.99	0.88	0.99
CSL	0.88	0.98	0.94	0.91	0.95	0.94	0.86	0.96
	(91.0)	(93.6)	(93.0)	(93.8)	(92.0)	(92.4)	(90.2)	(92.0)
DLSA	0.93	1.01	1.01	0.96	1.01	1.01	0.94	1.01
	(92.2)	(95.6)	(94.8)	(94.6)	(94.2)	(95.0)	(92.8)	(93.2)
TDLSA	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	(94.8)	(95.2)	(95.0)	(95.4)	(93.8)	(95.0)	(95.2)	(93.4)
SDLSA	1.00	-	-	1.01	-	-	1.03	-	3.01	1.00
$N=100$, $K=20$
OS	0.87	1.00	1.00	0.91	1.00	1.00	0.92	0.99
CSL	0.88	0.98	0.91	0.93	0.94	0.93	0.90	0.96
	(90.6)	(94.0)	(91.8)	(92.4)	(94.8)	(93.2)	(91.8)	(94.8)
DLSA	0.92	1.00	1.00	0.98	1.00	1.00	0.92	1.00
	(92.0)	(95.4)	(96.8)	(93.8)	(97.2)	(95.8)	(92.2)	(95.6)
TDLSA	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
	(95.0)	(94.8)	(96.8)	(94.4)	(97.0)	(95.6)	(95.4)	(95.6)
SDLSA	1.01	-	-	1.04	-	-	1.03	-	3.00	1.00

Table 9. Table 9: Simulation results (with 100 replications) for logistic regression model by interleaving the bootstrap methods. The REE of one-shot (OS), DLSA, and bootstrap estimators are reported.

$N$	$K$	Est.	${\theta }_1$	${\theta }_2$	${\theta }_3$	${\theta }_4$	${\theta }_5$	${\theta }_6$	${\theta }_7$	${\theta }_8$
			Setting 1: i.i.d Covariates
5	8	OS	0.70	0.97	0.97	0.78	0.94	0.97	0.73	0.96
		DLSA	0.90	1.03	1.02	0.90	1.03	1.02	0.95	1.02
		Bootstrap	0.68	1.02	1.04	0.71	1.01	1.04	0.76	1.01
20	10	OS	0.90	0.99	0.99	0.91	0.98	1.00	0.92	0.99
		DLSA	0.89	1.00	1.01	0.97	1.00	1.00	0.92	1.01
		Bootstrap	0.86	0.97	1.00	0.95	1.00	0.99	0.90	0.99
			Setting 2: Heterogenous Covariates
5	8	OS	0.54	0.85	0.82	0.63	0.78	0.87	0.62	0.77
		DLSA	0.77	1.02	1.02	0.86	1.04	1.04	0.84	1.04
		Bootstrap	0.52	0.97	1.00	0.65	0.99	1.01	0.61	1.05
20	10	OS	0.68	0.88	0.99	0.79	0.91	0.91	0.75	0.89
		DLSA	0.90	1.02	1.01	1.01	1.00	1.01	0.92	1.02
		Bootstrap	0.83	1.04	0.99	0.98	0.98	0.99	0.86	1.00

Table 10. Table 10: Simulation results for screening result with 100 replications. The CRs of 4 screening methods are reported.

$N$	$K$	SIS	SIRS	KC	DC		SIS	SIRS	KC	DC
		Linear Regression					Logistic Regression
		Setting 1: i.i.d Covariates
2000	5	1.00	0.71	0.87	0.91		0.99	0.48	0.73	0.89
10000	8	1.00	0.92	1.00	0.99		1.00	0.81	0.73	0.99
		Setting 2: Heterogenous Covariates
2000	5	0.98	0.73	0.79	0.84		0.95	0.55	0.73	0.85
10000	8	1.00	0.86	0.91	0.95		1.00	0.77	0.73	0.94

Table 11. Table 11: Simulation results for linear and logistic regression models after screening procedure (using SIS) with 100 replications. The numerical performances are evaluated for different sample sizes N 𝑁 N ( × 10 3 absent superscript 10 3 \times 10^{3} ) and numbers of workers K 𝐾 K . The CR screen , REE (of five modelling methods for all covariates), CR select , CM and are reported.

$N$	$K$	CRscreen	OS	CSL	DLSA	TDLSA	SDLSA	CRselect	CM
		Linear Regression
		Setting 1: i.i.d Covariates
2	5	1.00	0.97	0.92	1.00	1.00	1.02	0.97	0.97
10	8	1.00	0.98	0.92	1.00	1.00	1.02	1.00	1.00
		Setting 2: Heterogenous Covariates
2	5	0.98	1.00	1.00	1.00	1.00	1.00	0.74	0.74
10	8	1.00	0.98	0.87	1.00	1.00	1.01	1.00	1.00
		Logistic Regression
		Setting 1: i.i.d Covariates
2	5	0.99	0.18	0.03	0.45	1.09	0.96	0.87	0.82
10	8	1.00	0.37	0.30	0.58	1.07	1.08	1.00	0.98
		Setting 2: Heterogenous Covariates
2	5	0.95	0.55	0.05	0.71	0.95	0.91	0.42	0.39
10	8	1.00	0.44	0.32	0.58	1.02	0.98	0.99	0.98

Table 12. Table 12: Simulation results for comparison of the shrinkage estimation with the hypothesis testing approach with 500 replications. The coverage rate (CR) and percentage of identifying the true model (CM) are presented.

		Setting 1: i.i.d Case						Setting 2: Heterogenous Case
		Testing			SDLSA			Testing			SDLSA
$N$	$K$	CR	CM		CR	CM		CR	CM		CR	CM
20	5	1.00	0.47		1.00	0.94		1.00	0.47		1.00	0.90
100	10	1.00	0.49		1.00	0.99		1.00	0.46		1.00	0.99

Equations

$$\sqrt{N}(\hat{\theta }-{\theta }_0){\to }_{d}N(0,\Sigma )$$

$$\sqrt{n_k}({\hat{\theta }}_k-{\theta }_0){\to }_{d}N(0,{\Sigma }_k)$$

$$\mathcal{L}(\theta )$$

$$\approx N^{-1}\sum _{k=1}^K\sum _{i\in {\mathcal{S}}_k}(\theta -{\hat{\theta }}_k{)}^{\top }\ddot{\mathcal{L}}({\hat{\theta }}_k;Z_i)(\theta -{\hat{\theta }}_k)+C_2,$$

$$\tilde{\mathcal{L}}(\theta )$$

$$\stackrel{\mathrm{def}}{=}\sum _k(\theta -{\hat{\theta }}_k{)}^{\top }{\alpha }_k{\hat{\Sigma }}_k^{-1}(\theta -{\hat{\theta }}_k),$$

$$\widetilde{\theta}=\arg\min_{\theta}\widetilde{\mathcal{L}}(\theta)=\Big{(}\sum_{k}\alpha_{k}\widehat{\Sigma}_{k}^{-1}\Big{)}^{-1}\Big{(}\sum_{k}\alpha_{k}\widehat{\Sigma}_{k}^{-1}\widehat{\theta}_{k}\Big{)}.$$

$$\displaystyle\lim_{n\rightarrow\infty}\inf P\Big{\{}\inf_{\|\theta^{*}-\theta_{0}\|\geq\delta,1\leq k\leq K}\big{(}\mathcal{L}_{k}(\theta^{*})-\mathcal{L}_{k}(\theta_{0})\big{)}\geq\epsilon\Big{\}}=1,$$

$$\displaystyle~{}~{}\mbox{and}~{}~{}E\Big{\{}\frac{\partial\mathcal{L}_{k}(\theta)}{\partial\theta}\Big{|}_{\theta=\theta_{0}}\Big{\}}=\mathbf{0}.$$

$$\Omega_{k}(\theta)=E\Big{\{}\frac{\partial\mathcal{L}(\theta;Z_{i})}{\partial\theta}\frac{\partial\mathcal{L}(\theta;Z_{i})}{\partial\theta^{\top}}\Big{|}i\in\mathcal{S}_{k}\Big{\}}=E\Big{\{}\frac{\partial^{2}\mathcal{L}_{k}(\theta;Z_{i})}{\partial\theta\partial\theta^{\top}}\Big{|}i\in\mathcal{S}_{k}\Big{\}}$$

$$\Big{|}\frac{\partial^{3}}{\partial\theta_{i}\partial\theta_{j}\partial\theta_{l}}\mathcal{L}(\theta^{*};Z)\Big{|}\leq M_{ijl}(Z),~{}~{}~{}\mbox{for all}~{}~{}\theta^{*}\in B(\delta),$$

$$\sqrt{N}(\tilde{\theta }-{\theta }_0)=V({\theta }_0)+B({\theta }_0)$$

$$\widehat{\theta}_{k}^{(2)}=\widetilde{\theta}-\Big{(}\frac{\partial^{2}\mathcal{L}_{k}(\theta)}{\partial\theta\partial\theta^{\top}}\Big{|}_{\theta=\widetilde{\theta}}\Big{)}^{-1}\frac{\partial\mathcal{L}_{k}(\theta)}{\partial\theta}\Big{|}_{\theta=\widetilde{\theta}}.$$

$$\widetilde{\theta}^{(2)}=\Big{(}\sum_{k}\alpha_{k}\widehat{\Sigma}_{k}^{(2)-1}\Big{)}^{-1}\Big{(}\sum_{k}\alpha_{k}\widehat{\Sigma}_{k}^{(2)-1}\widehat{\theta}_{k}^{(2)}\Big{)}.$$

$$Q_{\lambda }(\theta )=\tilde{\mathcal{L}}(\theta )+\sum _j{\lambda }_j|{\theta }_j|.$$

$$P({\tilde{\theta }}_{\lambda }^{(-{\mathcal{M}}_T)}=0)\to 1$$

$$P(Y_i=1|X_i)\stackrel{\mathrm{def}}{=}p(X_i^{({\mathcal{M}}_T)})=\frac{\exp ({\theta }_0^{({\mathcal{M}}_T)\top }{X}_i^{({\mathcal{M}}_T)})}{1+\exp ({\theta }_0^{({\mathcal{M}}_T)\top }{X}_i^{({\mathcal{M}}_T)})}.$$

$$\sqrt{N}\Big{(}\widetilde{\theta}_{\lambda}^{(\mathcal{M}_{T})}-\theta^{(\mathcal{M}_{T})}\Big{)}\rightarrow_{d}N\Big{(}0,\Sigma_{\mathcal{M}_{T}}\Big{)}.$$

$$\text{mbox}{DBIC}_{\lambda }=({\tilde{\theta }}_{\lambda }-\tilde{\theta }{)}^{\top }{\hat{\Sigma }}^{-1}({\tilde{\theta }}_{\lambda }-\tilde{\theta })+\log N\times d{f}_{\lambda },$$

$${\mathbb{R}}_{-}=\{$$

$${\mathbb{R}}_{+}=\{\lambda \in {\mathbb{R}}^p:{\mathcal{M}}_{\lambda }\supset {\mathcal{M}}_T,{\mathcal{M}}_{\lambda }\ne {\mathcal{M}}_T\},$$

$$P\big{(}\inf_{\lambda\in\mathbb{R}_{-}\cup\mathbb{R}_{+}}\mbox{DBIC}_{\lambda}>\mbox{DBIC}_{\lambda_{N}}\big{)}\rightarrow 1.$$

$$Y_i=X_i^{\top }{\theta }_0+{ϵ}_i,$$

$$P(Y_i=1|X_i)=\frac{\exp (X_i^{\top }{\theta }_0)}{1+\exp (X_i^{\top }{\theta }_0)}.$$

$$P(Y|X_i,{\theta }_0)=\frac{{\lambda }^{Y_i}}{Y_i!}\exp (-\lambda ),$$

$${\omega }_j=g({\nu }_1,\cdots ,{\nu }_{j5})=\left|\frac{E(X_{ij}{Y}_i)-E(X_{ij})E(Y_i)}{\sqrt{(E{X}_{ij}^2-E^2(X_{ij}))(E{Y}_i^2-E^2(Y_i))}}\right|,$$

$${\omega }_j=g({\nu }_1,\cdots ,{\nu }_8)=\frac{{\nu }_{j1}+{\nu }_{j2}{\nu }_{j3}-2{\nu }_{j4}}{\sqrt{({\nu }_{j5}+{\nu }_{j2}^2-2{\nu }_{j6})({\nu }_{j7}+{\nu }_{j3}^2-2{\nu }_{j8})}},$$

$${\hat{\nu }}_{j1}^{(k)}=\frac{1}{n_k(n_k-1)}\sum _{i_1\ne i_2,i_1,i_2\in {\mathcal{S}}_k}|Y_{i1}-Y_{i2}|\cdot |X_{i_{1}j}-X_{i_{2}j}|.$$

$$\widetilde{\theta}-\theta_{0}=\big{\{}\sum_{k}\alpha_{k}\widehat{\Sigma}_{k}^{-1}\big{\}}^{-1}\big{\{}\sum_{k}\alpha_{k}\widehat{\Sigma}_{k}^{-1}(\widehat{\theta}_{k}-\theta_{0})\big{\}}.$$

$$\sum _k{\alpha }_k\{{\hat{\Sigma }}_k({\hat{\theta }}_k){\}}^{-1}{\to }_p{\Sigma }^{-1},$$

$$\displaystyle\sqrt{N}\Big{[}\sum_{k}\alpha_{k}\{\widehat{\Sigma}_{k}(\widehat{\theta}_{k})\}^{-1}(\widehat{\theta}_{k}-\theta_{0})\Big{]}=V^{*}(\theta_{0})+B^{*}(\theta_{0}),$$

$${\hat{\Sigma }}_k^{-1}({\hat{\theta }}_k)-{\Sigma }_k^{-1}={\hat{\Sigma }}_k^{-1}({\hat{\theta }}_k)-{\hat{\Sigma }}_k^{-1}({\theta }_0)+{\hat{\Sigma }}_k^{-1}({\theta }_0)-{\Sigma }_k^{-1}$$

$$\displaystyle=\sum_{j}\frac{\partial^{3}\mathcal{L}_{k}(\theta)}{\partial\theta_{j}\partial\theta\partial\theta^{\top}}\Big{|}_{\theta=\theta^{*}}(\theta_{j}^{*}-\theta_{j})+\widehat{\Sigma}_{k}^{-1}(\theta_{0})-\Sigma_{k}^{-1}$$

Full text

**KEY WORDS: ** Distributed System; Least Squares Approximation; Shrinkage Estimation; Distributed BIC.

11footnotetext: Feng Li is the corresponding author. Xuening Zhu is supported by the National Natural Science Foundation of China (nos. 11901105, 71991472, U1811461), the Shanghai Sailing Program for Youth Science and Technology Excellence (19YF1402700), and funds from Fudan University. Feng Li’s research is supported by the National Natural Science Foundation of China (no. 11501587). Hansheng Wang’s research is partially supported by National Natural Science Foundation of China (No. 11831008).

1 Introduction

Modern data analysis often needs to address huge datasets. In many cases, the size of the dataset could be too large to be conveniently handled by a single computer. Consequently, the dataset must be stored and processed on many connected computer nodes, which thereafter are referred to as a distributed system. More precisely, a distributed system refers to a large cluster of computers, which are typically connected with each other via wire protocols such as RPC and HTTP (Zaharia et al., 2012). Consequently, they are able to communicate with each other and accomplish the intended data analysis tasks at huge scales in a collective manner.

By using a distributed system, we are able to break a large-scale computation problem into many small pieces and then solve them in a distributed manner. A key challenge faced by statistical computation on a distributed system is the communication cost. The communication cost refers to the wall-clock time cost needed for data communication between different computer nodes, which could be expensive in distributed systems (Zhang et al., 2013; Shamir et al., 2014; Jordan et al., 2019). In this work, we consider a “master-and-worker”-type distributed system with strong workers. We assume that the workers are strong in the sense they are modern computers with reasonable storage and computing capacities. For example, a worker with 32 CPU cores, 128 GB RAM, and 512 GB SSD hard disk could be a very strong worker. As one will see later, the most widely used distributed environment, Hadoop (Apache Software Foundation, 2019a, version 2.7.2) and Spark (Apache Software Foundation, 2019b, version 2.3.1), belong to this category. Typically, the workers do not communicate with each other directly. However, they should be connected to a common master node, which is another computer with outstanding capacities. Consequently, most data should be distributed to workers, and most computations should be conducted by the workers. This enables us to solve a large-scale computation problem in a distributed manner. In contrast, the master should take the responsibility to coordinate with different workers.

For this “master-and-worker”-type distributed system, the communication cost is mostly between the master and workers. One can easily verify that good algorithms for some simple moment estimates (e.g., sample mean) can be easily developed using this type of distributed system. For example, to compute the sample mean on a distributed system, one can first compute the sample mean on each worker, which is known as a map process. Then, each worker reports to the master the resulting sample mean and the associated sample size. Thereafter, the master can compute the overall sample mean by a weighted average of the sample means from each worker, which is known as a reduce process. Such a “MapReduce” algorithm requires only one “master-and-worker” communication for each worker. It requires no direct communication between workers. Because most computations are accomplished by the workers, it also makes good use of the strong worker capacities. As a result, the algorithm can be considered effective. Unfortunately, cases such as the sample mean are rather rare in statistical analysis. Most statistical algorithms do not have an analytical solution (e.g., the maximum likelihood estimation of a logistic regression model) and thus require multiple iterations (e.g., Newton-Raphson iteration or stochastic-gradient-descent-type algorithms). These iterations unfortunately lead to substantial “master-and-worker” communication, which is communicationally expensive. Therefore, developing algorithms that are highly efficient computationally, communicationally and statistically for distributed systems has become a problem of great interest.

In the literature, the common wisdom for addressing a distributed statistical problem can be classified into two categories. The first category is the “one-shot” (OS) or “embarrassingly parallel” approach, which requires only one round of communication. Specifically, the local worker computes the estimators in parallel and then communicate to the master to obtain an average global estimator (Zhang et al., 2013; Liu and Ihler, 2014; Lee et al., 2015; Battey et al., 2015; Fan et al., 2017; Chang et al., 2017b, a). Although this approach is highly efficient in terms of communication, it might not achieve the best efficiency in statistical estimation in most occasions (Shamir et al., 2014; Jordan et al., 2019). The second approach includes iterative algorithms, which require multiple rounds of communication between the master and the workers. This approach typically requires multiple iterations to be taken so that the estimation efficiency can be refined to match the global (or centralized) estimator (Shamir et al., 2014; Wang et al., 2017a, b; Jordan et al., 2019). In addition, see Yang et al. (2016); Heinze et al. (2016); Smith et al. (2018); Li et al. (2019) for distributed statistical modelling methods when the data are distributed according to features rather than samples.

The aforementioned two approaches are also studied for the sparse learning problem using $\ell_{1}$ shrinkage estimation. For the first approach, Lee et al. (2015) investigated the distributed high-dimensional sparse regression using the OS approach by combining local debiased $\ell_{1}$ estimates. Battey et al. (2015) revisited the same problem but further considered distributed testing and estimation methods in a unified likelihood framework, in which a refitted estimation is used to obtain an oracle convergence rate. For the second approach, both Wang et al. (2017a) and Jordan et al. (2019) have developed iterative algorithms to solve the sparse estimation problem, and they theoretically proved that the error bounds match the centralized estimator. Beyond the $\ell_{1}$ shrinkage estimation, Chen and Xie (2014) studied a penalized likelihood estimator with more general penalty function forms in a high-dimensional setting. However, to the best of our knowledge, there are no guarantees that simultaneously ensure the model selection consistency (Fan and Li, 2001) and establish a criterion for consistent tuning parameter selection (Wang et al., 2007). In addition, all of the above methods assume independent and identical samples stored by each worker, which is questionable in practice because the distributed dataset might experience great heterogeneity from worker to worker. We would like to remark that the heterogeneity cannot be avoided because it is mainly due to the practical need to record data across time or space (for example). Ignoring the heterogeneity will lead to suboptimal or even biased estimation result.

In this work, we aim to develop a novel methodology to address a sparse estimation problem with low dimensions ($p<n$, where $n$ is the local sample size). Under the high dimensional setting, we refer to the recent work of Li et al. (2020) for a distributed pre-feature screening procedure, which consumes a communication cost of $O(p)$. Hence, the feature dimension can be greatly reduced and then we could launch our estimation procedure. Specifically, we assume the data possessed by different workers are allowed to be heterogeneous but share the same regression relationship. The proposed method borrows the idea of the least squares approximation (LSA, Wang and Leng, 2007) and can be used to handle a large class of parametric regression models on a distributed system. Specifically, let $Y\in\mathbb{R}$ be the response of interest, let $X$ be the associated predictor, and let $\theta\in\mathbb{R}^{p}$ be the corresponding regression coefficient. The objective is to estimate the regression parameter $\theta$ and conduct a variable selection on a distributed system that has one master and many strong workers. Assume data, denoted by $(Y_{i},X_{i})$, with $1\leq i\leq N$, that are distributed across different workers. Further, assume that the sample size on each worker is sufficiently large and of the same order. Under this setting, we propose a distributed LSA (DLSA) method. The key idea is as follows:

• (1)

First, we estimate the parameter $\theta$ on each worker separately by using local data on distributed workers. This can be done efficiently by using standard statistical estimation methods (e.g., maximum likelihood estimation). By assuming that the sample size of each worker is sufficiently large, the resulting estimator and its asymptotic covariance estimate should be consistent but not statistically efficient, as compared with the global estimates.

• (2)

Next, each worker passes the local estimator of $\theta$ and its asymptotic covariance estimate to the master. Because we do not consider a high-dimensional model setting, the communication cost in this regard should be negligible.

• (3)

Once the master receives all the local estimators from the workers, a weighted least squares-type objective function can be constructed. This can be viewed as a local quadratic approximation of the global log-likelihood functions. As one can expect, the resulting estimator shares the same asymptotic covariance with the full-size MLE method (i.e., the global estimator) under appropriate regularity conditions.

The major steps of the DLSA method are further illustrated in Figure 1. As one can see, the DLSA method reduces communication costs mainly by using only one round of communication and avoids further iterative steps. Given the DLSA objective function on the master node, we can further conduct shrinkage estimation on the master. This is done by formulating an adaptive Lasso-type (Zou, 2006; Zhang and Lu, 2007) objective function. The objective functions can be easily solved by the LARS algorithm (Efron et al., 2004) with minimal computation cost on the master. Thus, no communication is required. Accordingly, a solution path can be obtained on the master node. Thereafter, the best estimator can be selected from the solution path in conjunction with the proposed distributed Bayesian information criterion (DBIC). We theoretically show that the resulting estimation is selection consistent and as efficient as the oracle estimator, which is the global estimator obtained under the true model.

To summarize, we aim to make the following important contributions to the existing literature. First, we propose a master with strong workers (MSW) distributed system framework, which solves a large-scale computation problem in a communication-efficient way. Second, given this MSW system, we propose a novel DLSA method, which easily handles a large class of classical regression models such as linear regression, generalized linear regression, and Cox’s model. Third, due to the simple quadratic form of the objective function, the analytical solution path can be readily obtained using the LARS algorithm on the master. Then, the best model can be easily selected by the DBIC criterion. Finally, but also most importantly, the proposed DLSA method fully takes advantage of the specialty of the MSW system, which pushes the intensive computation to the workers and therefore is as computationally, communicationally and statistically efficient as possible. Furthermore, we would like to make a remark here that although the proposed DLSA is designed for a distributed system, it can also be applied to a single computer when there are memory constraints (Chen et al., 2018; Wang et al., 2018).

The remainder of this article is organized as follows. Section 2 introduces the model setting and the least squares approximation method. Section 3 presents a communication-efficient shrinkage estimation and a distributed BIC criterion. Numerical studies are given in Section 4. An application to U.S. Airline data with datasets greater than 52 GB is illustrated using the DLSA method on the Spark system in Section 5. The article concludes with a brief discussion in Section 6. All technical details are delegated to the Appendix.

2 Statistical modelling on distributed systems

2.1 Model and Notations

Suppose in the distributed system that there are in total $N$ observations, which are indexed as $i=1,\cdots,N$. The $i$th observation is denoted as $Z_{i}=(X_{i}^{\top},Y_{i})^{\top}\in\mathbb{R}^{p+1}$, where $Y_{i}\in\mathbb{R}$ is the response of interest and $X_{i}\in\mathbb{R}^{p}$ is the corresponding covariate vector. Specifically, the observations are distributed across $K$ local workers. Define $\mathcal{S}=\{1,\cdots,N\}$ to be all sample observations. Decompose $\mathcal{S}=\cup_{k=1}^{K}\mathcal{S}_{k}$, where $\mathcal{S}_{k}$ collects the observations distributed to the $k$th worker. Obviously, we should have $\mathcal{S}_{k_{1}}\cap\mathcal{S}_{k_{2}}=\emptyset$ for any $k_{1}\neq k_{2}$. Define $n=N/K$ as the average sample size for each worker. Then, we assume $|\mathcal{S}_{k}|=n_{k}$ and that all $n_{k}$ diverge in the same order $O(n)$. Specifically, $c_{1}\leq\min_{k}n_{k}/n\leq\max_{k}n_{k}/n\leq c_{2}$ for some positive constants $c_{1}$ and $c_{2}$. We know immediately that $N=\sum n_{k}$. In practice, due to the data storing strategy, the data in different workers could be quite heterogeneous, e.g., they might be collected according to spatial regions. Despite the heterogeneity here, we assume they share the same regression relationship, and the parameter of interest is given by $\theta_{0}\in\mathbb{R}^{p}$. We focus on the case in which $p$ is fixed.

Let $\mathcal{L}(\theta;Z)$ be a plausible twice-differentiable loss function. Define the global loss function as $\mathcal{L}(\theta)=N^{-1}\sum_{i=1}^{N}\mathcal{L}(\theta;Z_{i})$, whose global minimizer is $\widehat{\theta}=\arg\min\mathcal{L}(\theta)$ and the true value is $\theta_{0}$. It is assumed that $\widehat{\theta}$ admits the following asymptotic rule

[TABLE]

for some positive definite matrix $\Sigma\in\mathbb{R}^{p\times p}$ as $N\rightarrow\infty$. If $\mathcal{L}(\theta;Z)$ is the negative log-likelihood function, then $\widehat{\theta}$ is the global MLE estimator. Correspondingly, define the local loss function in the $k$th worker as $\mathcal{L}_{k}(\theta)=n_{k}^{-1}\sum_{i\in\mathcal{S}_{k}}\mathcal{L}(\theta;Z_{i})$, whose minimizer is $\widehat{\theta}_{k}=\arg\min_{\theta}\mathcal{L}_{k}(\theta)$. We assume that

[TABLE]

as $n_{k}\rightarrow\infty$ for a positive definite matrix $\Sigma_{k}$. The goal is to conduct statistical analysis based on the data on the local worker and minimize the communication cost as much as possible.

2.2 Least Squares Approximation and Variance Optimality

In this section, we motivate our approach through least squares approximation to the global loss function, which takes a local quadratic form. To motivate this idea, we begin by decomposing and approximating the global loss function using Taylor’s expansion techniques as follows:

[TABLE]

where the last equation uses the fact that $\dot{\mathcal{L}}_{k}(\widehat{\theta}_{k})=0$, and $C_{1}$ and $C_{2}$ are some constants. Typically, the minimizer $\widehat{\theta}_{k}$ will achieve the convergence rate $\sqrt{n_{k}}$. Intuitively, the quadratic form in (2.1) should be a good local approximation of the global loss function (Wang and Leng, 2007). This inspires us to consider the following weighted least squares objective function:

[TABLE]

where $\alpha_{k}=n_{k}/N$. This leads to a weighted least squares estimator (WLSE), which takes an analytical form as follows:

[TABLE]

It is remarkable that the estimator $\tilde{\theta}$ in (2.2) can be easily computed on a distributed system. Specifically, the local worker sends $\widehat{\theta}_{k}$ and $\widehat{\Sigma}_{k}$ to the master node, and then, the master node produces the WLSE by (2.2). As a result, the above WLSE requires only one round of communication. Hence, it is highly efficient in terms of communication.

Note that instead of taking a simple average of local estimators $\widehat{\theta}_{k}$ in the literature, the analytical solution in (2.2) takes a weighted average of $\widehat{\theta}_{k}$ using weights $\widehat{\Sigma}_{k}^{-1}$. This will result in a higher statistical efficiency if the data are stored heterogeneously. To investigate the asymptotic properties of the WLSE, we assume the following conditions.

• (C1)

(Parameter Space) The parameter space $\Theta$ is a compact and convex subset of $\mathbb{R}^{p}$. In addition, the true value $\theta_{0}$ lies in the interior of $\Theta$.

• (C2)

(Covariates Distribution) Assume the covariates $X_{i}$ ($i\in\mathcal{S}_{k}$) from the $k$th worker are independently and identically distributed from the distribution $F_{k}(x)$.

• (C3)

(Identifiability) For any $\delta>0$, there exists $\epsilon>0$ such that

[TABLE]

• (C4)

(Local Convexity) Define

[TABLE]

Assume that $\Omega_{k}(\theta)$ is nonsingular at the true value $\theta_{0}$. In addition, let $\Sigma_{k}=\{\Omega_{k}(\theta_{0})\}^{-1}$, and $\Sigma=\{\sum_{k}\alpha_{k}\Omega(\theta_{0})\}^{-1}$.

• (C5)

(Smoothness) Define $B(\delta)=\{\theta^{*}\in\Theta|\|\theta^{*}-\theta_{0}\|\leq\delta\}$ as a ball around the true value $\theta_{0}$ with radius $\delta>0$. Assume for almost all $Z\in\mathbb{R}^{p}$ that the loss function $\mathcal{L}(\theta;Z)$ admits all third derivatives $\partial^{3}\mathcal{L}(\theta;Z)/(\partial\theta_{i}\partial\theta_{j}\partial\theta_{l})$ for all $\theta\in B(\delta)$. In addition, assume that there exist functions $M_{ijl}(Z)$ and $\delta>0$ such that

[TABLE]

where $E\{M_{ijl}(Z_{m})|m\in\mathcal{S}_{k}\}<\infty$ for all $1\leq i,j,l\leq p$ and $1\leq k\leq K$.

The above conditions are standard conditions to establish the asymptotic properties for $M$-estimators. First, Condition (C1) assumes the parameter space to be convex (Jordan et al., 2019). Next, Condition (C2) concerns the distribution of the covariates $\{X_{i}:i\in\mathcal{S}_{k}\}$. Specifically, there are different $F_{k}(x)$ for different $1\leq k\leq K$. In particular, it allows for the heterogeneous distribution of covariates across workers. We would like to remark that heterogeneity is a common phenomenon in distributed systems, and it has been ignored in much of the literature. Condition (C3) assures the identifiability of the local loss functions across all workers. Finally, Conditions (C4) and (C5) are standard regularity conditions of the loss functions, which require certain degrees of local convexity and smoothness of the loss functions. These conditions are widely assumed in the literature to guarantee asymptotic convergence of the estimators (Fan and Li, 2001; Lehmann and Casella, 2006; Jordan et al., 2019).

Given the conditions, we can establish the asymptotic properties of WLSE in the following Proposition 1 and Theorem 1.

Proposition 1**.**

Assume Conditions (C1)–(C5). Then, we have

[TABLE]

with $\mbox{cov}\{V(\theta_{0})\}=\Sigma$ and $B(\theta_{0})=O_{p}(K/\sqrt{N})$, where $\Sigma=(\sum_{k=1}^{K}\alpha_{k}\Sigma_{k}^{-1})^{-1}.$

The proof of Proposition 1 is given in Appendix A.1. Proposition 1 separates $\sqrt{N}(\widetilde{\theta}-\theta_{0})$ into two parts, namely, the variance part and the bias part. Particularly, one should note that the variance order is the same as the global estimator $\widehat{\theta}$, which is $O(N^{-1})$, while the bias order is related to the number of local workers $K$. In addition, the covariance of $V(\theta_{0})$ is exactly the same with the global estimator (i.e., $\Sigma$). Consequently, if the local sample size is sufficiently large, the bias should be sufficiently small, and thus, the estimation efficiency will be the same as the global estimator. We state this result in the following theorem.

Theorem 1**.**

(Global Asymptotic Normality)* Assume conditions (C1)–(C5) and further assume $n/N^{1/2}\rightarrow\infty$. Then, we have $\sqrt{N}(\widetilde{\theta}-\theta_{0})\rightarrow_{d}N(0,\Sigma)$, which achieves the same asymptotic normality as the global estimator $\widehat{\theta}$.*

The proof of Theorem 1 is given in Appendix A.2. It can be concluded that we should require the local sample size to be of order larger than $\sqrt{N}$, which is easy to satisfy in practice. Otherwise, we should have $N/n^{2}=K/n\rightarrow\infty$. This implies that the number of workers is even larger than the average local sample size $n$. This is obviously not the case in practice if we have limited computational resources.

2.3 Two-Step Estimation

In the previous analysis, to guarantee a $\sqrt{N}$-consistency rate of the WLSE, we require that $n_{k}/\sqrt{N}\rightarrow\infty$. The assumption might not be suitable if we have sufficient computational resources (e.g., a large number of workers). As a result, it could violate the condition that $K/n\rightarrow 0$ if $K$ is much larger than $n$. To allow the case with large $K$ and small $n$, we further propose a two-step estimator to refine $\widetilde{\theta}$. The details are given as follows.

In the first step, we obtain WLSE $\widetilde{\theta}$ using one round of communication. Next, the WLSE is broadcasted from the master node to the workers. In the second step, on the $k$th work, we perform a one-step iteration using $\widetilde{\theta}$ as the initial value to obtain

[TABLE]

Subsequently, the local estimator $\widehat{\theta}_{k}^{(2)}$ as well as $\widehat{\Sigma}_{k}^{(2)}=\{\partial^{2}\mathcal{L}_{k}(\widetilde{\theta})/\partial\theta\partial\theta^{\top}\}^{-1}$ are then transmitted from the worker to the master. Similar as the design of WLSE, we obtain a two-step WLSE (TWLSE) on the master as

[TABLE]

As we will show in the theoretical analysis, TWLSE enjoys smaller bias rate than WLSE. It is mainly because that it takes a global estimator as an initial value with a further one-step iteration. Hence it is capable to borrow the power of the global estimator and further reduce the bias from $O(n^{-1})$ to $O(n^{-2})$. We state the result rigorously in the following Theorem.

Theorem 2**.**

Assume Condition (C1)–(C5). Then we have $\sqrt{N}(\widetilde{\theta}^{(2)}-\theta_{0})=V_{2}(\theta_{0})+B_{2}(\theta_{0})$ with $\mbox{cov}\{V_{2}(\theta_{0})\}=\Sigma$ and $B_{2}(\theta_{0})=O_{p}(n^{-2}N^{1/2})$. In addition, we have $\sqrt{N}(\widetilde{\theta}^{(2)}-\theta_{0})\rightarrow_{d}N(0,\Sigma)$ if we assume $n/N^{1/4}\rightarrow\infty$.

The proof of Theorem 2 is given in Appendix A.3. The assumption $n/N^{1/4}\rightarrow\infty$ is much milder than WLSE and in the meanwhile it allows $K=O(N^{3/4})$. We would like to remark that the same theoretical result can be achieved if we use one-shot estimator in our first step to reduce the communication cost.

One could see that our TWLSE can be easily extended to $m$-WLSE with $m$ steps of iterations. Specifically, we could repeat step (2.5) and (2.6) for $m$ rounds to obtain $\widehat{\theta}_{k}^{(m)}$ and $\widetilde{\theta}^{(m)}$ respectively. Particularly, if we let $m/\log N\rightarrow\infty$, we could allow for an extreme case that $n=O(1)$. This is suitable for the case when we have memory constraint of local workers. The theoretical result is given in the following Proposition.

Proposition 2**.**

Assume Condition (C1)–(C5) and $m/\log N\rightarrow\infty$. Then we have $\sqrt{N}(\widetilde{\theta}^{(m)}-\theta_{0})\rightarrow_{d}N(0,\Sigma)$.

The proof of Proposition 2 is given in Appendix A.4. To achieve the desired property, we trade off the communication cost to reduce local computational burden with a large number of workers.

3 Communication-efficient shrinkage estimation

3.1 Distributed Adaptive Lasso Estimation and Oracle Property

Variable selection is a classical but critically important problem. That is because in practice, the number of available covariates is typically large, but only a small number of covariates are related to the response. Given an appropriate variable selection technique, one can discover the important variables with high probability. In recent decades, various variable selection techniques have been well studied (Tibshirani, 1996; Fan and Li, 2001; Zou, 2006; Wang and Leng, 2007; Zhang, 2010). However, how to conduct variable selection on a distributed system has not been sufficiently investigated. Existing approaches mostly focus on the $\ell_{1}$ shrinkage estimation and develop corresponding algorithms (Lee et al., 2015; Battey et al., 2015; Wang et al., 2017a; Jordan et al., 2019). However, to the best of our knowledge, there are three problems that remain unsolved on a distributed system: (a) most works do not establish the oracle properties of the shrinkage estimators, (b) no consistent tuning parameter selection criterion is given or investigated, and (c) the computation will be heavy if one needs to conduct estimation and select the tuning parameters simultaneously.

To solve the above problems, we first define some notations. Without loss of generality, we assume the first $d_{0}$ ($0<d_{0}<p$) to be nonzero, i.e., $\theta_{j}\neq 0$ for $1\leq j\leq d_{0}$ and $\theta_{j}=0$ for $j>d_{0}$. Correspondingly, we denote $\mathcal{M}_{T}=\{1,\cdots,d_{0}\}$ to be true model. In addition, let $\mathcal{M}=\{i_{1},\cdots,i_{d}\}$ be an arbitrary candidate model with size $|\mathcal{M}|=d$. In addition, for an arbitrary vector $v=(v_{j}:1\leq j\leq p)^{\top}$, define $v^{(\mathcal{M})}=(v_{i}:i\in\mathcal{M})^{\top}\in\mathbb{R}^{|\mathcal{M}|}$ and $v^{(-\mathcal{M})}=(v_{i}:i\not\in\mathcal{M})^{\top}\in\mathbb{R}^{p-|\mathcal{M}|}$. For an arbitrary matrix $M=(m_{ij})$, define $M^{(\mathcal{M})}=(m_{j_{1}j_{2}}:j_{1},j_{2}\in\mathcal{M})\in\mathbb{R}^{|\mathcal{M}|\times|\mathcal{M}|}$.

For simultaneous variable selection and parameter estimation, we follow the idea of Wang and Leng (2007) and consider the adaptive Lasso objective function on the master (Zou, 2006; Zhang and Lu, 2007),

[TABLE]

By the adaptive Lasso method, different amounts of shrinkage $\lambda_{j}$ are imposed on each estimator to improve the estimation efficiency (Zou, 2006; Zou and Li, 2008). Compared to the LSA approach of Wang and Leng (2007), we have the following key differences. First, $\widetilde{\theta}$ is the combined WLSE from local workers. Second, $\widehat{\Sigma}$ is constructed by the local asymptotic covariance estimators $\widehat{\Sigma}_{k}$. Consequently, to achieve a global convergence rate, one needs to carefully balance the local convergence rate of $\widehat{\theta}_{k}$ and $\widehat{\Sigma}_{k}$ with that of the global ones.

Remark 1**.**

Under the high dimensional setting with large $p$, the communication of $\widehat{\Sigma}_{k}$ is costly and inefficient. Under such a case, we recommend to conduct a pre-feature screening procedure (Li et al., 2020) to screening important features by using one round of communication. That consumes a communication cost of $O(p)$. This enables us to reduce the feature dimension $p$ greatly to a feasible size. Next, we could communicate the covariance estimator of important features to the master to complete shrinkage estimation using (3.1), which is communication friendly under the high dimensional setting. We give the details of the pre-feature screening procedure in Section 4.2.

Define $\widetilde{\theta}_{\lambda}=\arg\min_{\theta}Q_{\lambda}(\theta)$. Then, we can establish the $\sqrt{N}$-consistency as well as the selection consistency result of $\widetilde{\theta}_{\lambda}$ under certain conditions of $\lambda=(\lambda_{1},\cdots,\lambda_{p})^{\top}$.

Theorem 3**.**

*Assume the conditions (C1)–(C5) and $n/N^{1/2}\rightarrow\infty$. Let $a_{\lambda}=\max\{\lambda_{j},j\leq d_{0}\}$ and $b_{\lambda}=\min\{\lambda_{j},j>d_{0}\}$. Then, the following result holds.

a. ($\sqrt{N}$-consistency). If $\sqrt{N}a_{\lambda}\rightarrow_{p}0$, then $\widetilde{\theta}_{\lambda}-\theta=O_{p}(N^{-1/2})$.

b. (Selection Consistency). If $\sqrt{N}a_{\lambda}\rightarrow_{p}0$ and $\sqrt{N}b_{\lambda}\rightarrow_{p}\infty$, then*

[TABLE]

The proof of Theorem 3 is given in Appendix B.1. Note here that $a_{\lambda}$ controls the largest amount of penalty on the true nonzero parameters. Consequently, this amount cannot be too large; otherwise, it will result in a highly biased estimator. In contrast, $b_{\lambda}$ is responsible for producing sparse solutions of irrelevant covariates. Therefore, $b_{\lambda}$ should be sufficiently large to produce an effective amount of shrinkage.

Remark 2**.**

Note that the objective function of (3.1) constitutes of two parts. The first part is the weighted least squares type objective function, and the second part is a penalty term. Hence the Theorem 3 and the following theoretical properties are established based on Theorem 1. If we could allow for more steps of iterations (with TWLSE or $m$-WLSE), we could further relax the condition about $n$ here as in Theorem 2 and Proposition 2.

By Theorem 3, we know that with probability tending to one, we have $\widetilde{\theta}_{\lambda}^{(-\mathcal{M}_{T})}=\mathbf{0}$. Meanwhile, $\widetilde{\theta}_{\lambda}^{(\mathcal{M}_{T})}-\theta_{0}^{(\mathcal{M}_{T})}=O_{p}(N^{-1/2})$. It is then natural to ask whether the statistical efficiency of $\widetilde{\theta}_{\lambda}^{(\mathcal{M}_{T})}$ can be as good as the oracle estimator, which is the global estimator obtained under the true model, i.e., $\widehat{\theta}_{oracle}=\arg\min_{\theta\in\mathbb{R}^{p},\theta_{j}=0,\forall j\not\in\mathcal{M}_{T}}\mathcal{L}(\theta)$. To this end, we require the following technical condition.

• (C6)

(Covariance Assumption) Define the global unpenalized estimator as $\widehat{\theta}_{\mathcal{M}}=\arg\min_{\{\theta\in\mathbb{R}^{p}:\theta_{j}=0,\forall j\not\in\mathcal{M}\}}\mathcal{L}(\theta).$ Assume for the global estimator $\widehat{\theta}_{\mathcal{M}}$ with $\mathcal{M}\supset\mathcal{M}_{T}$ that $\sqrt{N}(\widehat{\theta}_{\mathcal{M}}^{(\mathcal{M})}-\theta_{0}^{(\mathcal{M})})\rightarrow_{d}N(0,\Sigma_{\mathcal{M}})=N(0,\Omega_{\mathcal{M}}^{-1})$, where $\Sigma_{\mathcal{M}}\in\mathbb{R}^{|\mathcal{M}|\times|\mathcal{M}|}$ is a positive-definite matrix, and $\widehat{\theta}_{\mathcal{M}}^{(\mathcal{M})}=(\widehat{\theta}_{\mathcal{M},j}:j\in\mathcal{M})^{\top}\in\mathbb{R}^{|\mathcal{M}|}$. Further assume for any $\mathcal{M}\supset\mathcal{M}_{T}$ that $\Omega_{\mathcal{M}}=\Omega_{\mathcal{M}_{F}}^{(\mathcal{M})}$ holds, where $\mathcal{M}_{F}=\{1,\cdots,p\}$ denotes the whole set.

Condition (C6) does not seem very intuitive. Nevertheless, it is a condition that is well satisfied by most maximum likelihood estimators. For instance, consider a logistic regression model with response $Y_{i}\in\{0,1\}$, i.e.,

[TABLE]

For an overfitted model, the inverse asymptotic covariance of $\sqrt{N}\widehat{\theta}_{\mathcal{M}}^{(\mathcal{M})}$ is $\Omega_{\mathcal{M}}=\sum_{k=1}^{K}\alpha_{k}E\{$ $p(X_{i}^{(\mathcal{M}_{T})})(1-p(X_{i}^{(\mathcal{M}_{T})}))X_{i}^{(\mathcal{M})}X_{i}^{(\mathcal{M})\top}\}$, which is a submatrix of $\Omega_{\mathcal{M}_{F}}=\sum_{k=1}^{K}\alpha_{k}E\{p(X_{i}^{(\mathcal{M}_{T})})$ $(1-p(X_{i}^{(\mathcal{M}_{T})}))X_{i}^{(\mathcal{M}_{F})}X_{i}^{(\mathcal{M}_{F})\top}\}$. Hence Condition (C6) holds. By similar argument, one can show that for any regression model with likelihood function of the form $\prod_{i}f(Y_{i},X_{i}^{\top}\theta)$, the Condition (C6) is satisfied. A more detailed discussion has been provided by Wang and Leng (2007). We then have the oracle property in the following theorem.

Theorem 4**.**

(Oracle Property)* Assume Conditions (C1)–(C6) and $n/N^{1/2}\rightarrow\infty$. Let $\sqrt{N}a_{\lambda}\rightarrow_{p}0$ and $\sqrt{N}b_{\lambda}\rightarrow_{p}\infty$; then, it holds that*

[TABLE]

By Theorem 3 and 4, we know that as long as the tuning parameters are approximately selected, the resulting estimator is selection consistent and as efficient as the oracle estimator. It is remarkable that tuning a total of $p$ parameters simultaneously is not feasible in practice. To fix this problem, we follow the tradition of Zou (2006) and Wang et al. (2007) to specify $\lambda_{j}=\lambda_{0}|\widetilde{\theta}_{j}|^{-1}$. Since $\widetilde{\theta}_{j}$ is $\sqrt{N}$-consistent, then as long as as $\lambda_{0}$ satisfies the condition $\lambda_{0}\sqrt{N}\rightarrow 0$ and $\lambda_{0}N\rightarrow\infty$, then the conditions $\sqrt{N}a_{\lambda}\rightarrow_{p}0$ and $\sqrt{N}b_{\lambda}\rightarrow_{p}\infty$ are satisfied. Thereafter, the original problem of tuning parameter selection for $\lambda$ can be replaced by selection for $\lambda_{0}$.

3.2 The Distributed Bayes Information Criterion

Although it has been shown that asymptotically, the oracle property can be guaranteed as long as the tuning parameters are approximately selected, it is still unclear how to conduct variable selection in practice. That motivates us to design a BIC-type criterion that can select the true model consistently in a completely data-driven manner (Zhang and Lu, 2007; Chen and Chen, 2008; Zou and Zhang, 2009; Wang et al., 2013). Specifically, to consistently recover the sparsity pattern, we consider a distributed Bayesian information criterion (DBIC)-based criterion as follows:

[TABLE]

where $\widehat{\Sigma}=(\sum_{k=1}^{K}\alpha_{k}\widehat{\Sigma}_{k}^{-1})^{-1}$ and $df_{\lambda}$ is the number of nonzero elements in $\widetilde{\theta}_{\lambda}$.

The design of the DBIC criterion is in the spirit of the BIC criterion used in Wang and Leng (2007). The difference is that the DBIC uses the WLSE estimator $\widetilde{\theta}$ and the average of distributed covariance estimators $\widehat{\Sigma}$ to construct the least squares objective function. Intuitively, if $\widetilde{\theta}$ and $\widehat{\Sigma}$ approximate the global estimator $\widehat{\theta}=\arg\max\mathcal{L}(\theta)$ and asymptotic covariance very well, then the DBIC criterion should be able to facilitate consistent tuning parameter selection. Specifically, the resulting model should be selection consistent (Shao, 1997).

To formally investigate the theoretical performance of DBIC, we first define some notations. First, we define the set of nonzero elements of $\widehat{\theta}_{\lambda}$ by $\mathcal{M}_{\lambda}$. Given a tuning parameter $\lambda$, $\mathcal{M}_{\lambda}$ could be underfitted, correctly fitted or overfitted. We could then have the following partition:

[TABLE]

where $\mathbb{R}_{-}$ denotes the underfitted model, and $\mathbb{R}_{+}$ denotes an overfitted model. We show in the following Theorem that the DBIC can consistently identify the true model.

Theorem 5**.**

Assume Conditions (C1)–(C6) and $n/N^{1/2}\rightarrow\infty$. Define a reference tuning parameter sequence $\{\lambda_{N}\in\mathbb{R}^{p}\}$, where the first $d_{0}$ elements of $\lambda_{N}$ are $1/N$ and the remaining elements are $\log N/\sqrt{N}$. Then, we have

[TABLE]

By Theorem 3 and 4, we know that with probability tending to one, we should have $\mathcal{M}_{\lambda_{N}}=\mathcal{M}_{T}$. Consequently, the sequence $\lambda_{N}$ here plays a role as a reference sequence that leads to the true model. Accordingly, Theorem 5 implies that the optimal $\lambda$ selected by the DBIC will consistently identify the true model. This is because any $\lambda$ leading to an inconsistent model selection result should perform worse than $\lambda_{N}$ in terms of DBIC values. Namely, denote the estimated set selected by DBIC is $\widehat{\mathcal{M}}_{DBIC}$, then Theorem 5 implies that $\widehat{\mathcal{M}}_{DBIC}=\mathcal{M}_{T}$ with probability tending to 1.

4 Numerical studies

4.1 Simulation Models and Settings

To demonstrate the finite sample performance of the DLSA method, we conduct a number of simulation studies in this section. Five classical regression models are presented, and the corresponding DLSA algorithms are implemented. For each model, we consider two typical settings to verify the numerical performance of the proposed method. They represent two different data storing strategies together with competing methods. The first strategy is to distribute data in a complete random manner. Thus, the covariates on different workers are independent and identically distributed (i.i.d). In contrast, the second strategy allows for covariate distribution on different workers to be heterogeneous. The estimation efficiency as well as the as the variable selection accuracy are evaluated. Examples are given as follows.

Example 1. (Linear Regression). We first consider one of the most popular regression analysis tools, i.e., linear regression. In particular, we generate the continuous response $Y_{i}$ by a linear relationship with the covariates $X_{i}$ as follows:

[TABLE]

where the noise term $\varepsilon_{i}$ is independently generated using a standard normal distribution $N(0,1)$. Following Fan and Li (2001), the true parameter is set as $\theta_{0}=(3,1.5,0,0,2,0,0,0)^{\top}$.

Example 2. (Logistic Regression). The logistic regression is a classical model that addresses binary responses (Hosmer Jr et al., 2013). In this example, we generate the response $Y_{i}$ independently by the Bernoulli distribution given the covariate $X_{i}$ as

[TABLE]

We follow Wang and Leng (2007) to set the true parameter $\theta_{0}=(3,0,0,1.5,0,0,2,0)^{\top}$.

Example 3. (Poisson Regression). In this example, we consider the Poisson regression, which is used to model counted responses (Cameron and Trivedi, 2013). The responses are generated according to the Poisson distribution as

[TABLE]

where $\lambda=\exp(X_{i}^{\top}\theta_{0})$. The true parameter $\theta_{0}$ is set to $(0.8,0,0,1,0,0,-0.4,0,0)^{\top}$.

For each example, two different data storage strategies are considered. They lead to different covariate distributions $F_{x}(x)$. Specifically, the following two settings are investigated.

• •

Setting 1 (i.i.d Covariates). We first consider the setting in which the data are distributed independently and identically across the workers. Specifically, the covariates $X_{ij}$ ($1\leq i\leq N,1\leq j\leq p$) are sampled from the standard normal distribution $N(0,1)$.

• •

Setting 2 (Heterogeneous Covariates). Next, we look at the case whereby the covariates distributed across each worker are heterogeneous. This is a common case in practice. Specifically, on the $k$th worker, the covariates are sampled from the multivariate normal distribution $N(\mu_{k},\Sigma_{k})$, where $\mu_{k}=k/K$, and $\Sigma_{k}=(\sigma_{k,ij})=(\rho_{k}^{|j_{1}-j_{2}|})$ with $\rho_{k}$ sampled from $U[0.2,0.3]$.

The i.i.d setting is consistent with the existing methods such as the one-shot estimation method. The second heterogenous data setting is not widely considered in literature but more realistic in practice. We evaluate the numerical performances of the proposed method and existing methods under both settings.

4.2 Pre-Feature Screening for High Dimensional Data

For high dimensional data, transmitting covariance estimators $\widehat{\Sigma}_{k}\in\mathbb{R}^{p\times p}$ could be costly. To reduce the feature dimension, we use a pre-feature screening procedure to screening important features before we conduct model estimation. As a consequence, the communication cost can be well controlled. In a recent work of Li et al. (2020), they design a distributed feature screening method using componentwise debiasing approach. The method first expresses specific correlation measures using a non-linear function of several component parameters, and then conduct distributed unbiased estimation of each component parameter. The method is shown to have better performance than simple average of correlation measures among workers.

Specifically, we express the correlation measure between the response $Y$ and the $j$th feature as $\omega_{j}=g(\nu_{j1},\cdots,\nu_{js})$, where $g$ is pre-specified function and $\nu_{j1},\cdots,\nu_{js}$ are $s$ components. Suppose we obtain the estimation of the $\nu_{jm}$ on the $j$th worker as $\widehat{\nu}_{jm}^{(k)}$. Then by distributed estimation we can obtain $\overline{\nu}_{jm}=K^{-1}\sum_{k=1}^{K}\widehat{\nu}_{jm}^{(k)}$ and this enables us to resemble $\widehat{\omega}_{j}=g(\overline{\nu}_{j1},\cdots,\overline{\nu}_{js})$ on the master. As a result, the pre-feature screening procedure will use one round of communication with communication cost of $O(p)$. After the screening procedure, we broadcast the indexes of screened important features to workers for further estimation. For illustration purpose, we give the following two examples to explain that how the correlation measure is calculated. We also include the Kendall $\tau$ rank correlation (Li et al., 2012a) and SIRS correlation (Zhu et al., 2011) for comparison. See Li et al. (2020) for more details.

1. Pearson Correlation.

The Pearson correlation is used by Fan and Lv (2008) in the sure independence screening (SIS) procedure. In this case we have

[TABLE]

where $\nu_{j1}=E(X_{ij}Y_{i})$, $\nu_{j2}=E(X_{ij})$, $\nu_{j3}=E(Y_{i})$, $\nu_{j4}=E(X_{ij}^{2})$, and $\nu_{j5}=E(Y_{i}^{2})$. Using simple moments estimation we could estimate the above five components. For instance we have $\widehat{\nu}_{j1}^{(k)}=n_{k}^{-1}\sum_{i\in\mathcal{S}_{k}}X_{ij}Y_{i}$. As a result, the estimated $\widehat{\omega}_{j}$ is exactly the same with using the whole dataset in this case.

2. Distance Correlation.

The distance correlation (DC) is used by Li et al. (2012b) as a model-free screening index. In this case $\omega_{j}$ can be expressed as

[TABLE]

where $\nu_{j1}=E\{|Y_{i}-Y_{i}^{\prime}|\cdot|X_{ij}-X_{ij}^{\prime}|\}$, $\nu_{j2}=E\{|Y_{i}-Y_{i}^{\prime}|\}$, $\nu_{j3}=E\{|X_{ij}-X_{ij}^{\prime}|\}$, $\nu_{j4}=E\{E(|Y_{i}-Y_{i}^{\prime}|~{}|Y_{i})E(|X_{ij}-X_{ij}^{\prime}|~{}|X_{ij})\}$, $\nu_{j5}=E\{|Y_{i}-Y_{i}^{\prime}|^{2}\}$, $\nu_{j6}=E\{E^{2}(|Y_{i}-Y_{i}^{\prime}|^{2}|Y_{i})\}$, $\nu_{j7}=E\{|X_{ij}-X_{ij}^{\prime}|^{2}\}$, $\nu_{j8}=E\{E^{2}(|X_{ij}-X_{ij}^{\prime}|~{}|X_{ij})\}$, and $(Y_{i}^{\prime},X_{ij}^{\prime})$ is an independent copy of $(Y_{i},X_{ij})$. On the $k$th worker, we could estimate $\nu_{j1}$ to $\nu_{j8}$ using local dataset. For instance, we can estimate $\nu_{j1}$ as follows,

[TABLE]

We refer to Li et al. (2012b) for detailed estimations for $\nu_{j2}$ to $\nu_{j8}$.

Under the high dimensional setting we set $(N,p,K)=(2000,10^{3},5)$ and $(10^{4},10^{4},8)$ respectively. Hence the second case is more challenging with $N=p$. Correspondingly, the non-zero coefficients are set as $\theta_{0j}=U_{1}\cdot\mbox{sign}(U_{2})$ for $1\leq j\leq 15$, where $U_{1}\sim U[0.5,2]$ and $U_{2}\sim U[-0.2,0.8]$. The numerical performances are evaluated for linear regression (Example 1) and logistic regression (Example 2) models respectively. After the pre-screening procedure, we select top $p_{0}=$40 features with highest screening measure values ($\widehat{\omega}_{j}$) as our candidate important feature set $\widehat{\mathcal{M}}_{screen}$. Based on $\widehat{\mathcal{M}}_{screen}$, we further conduct post-estimation with our DLSA method. We would like to further remark that if the post-inference is needed, the pre-feature screening procedure may introduce post-selective bias (Berk et al., 2013; Lee and Taylor, 2014; Taylor and Tibshirani, 2015; Lee et al., 2016). Such bias may not be a major issue for estimation, but is critical for statistical inference. Hence, selection-adjusted statistical inference should be investigated and devised. We leave this problem as a future research topic.

4.3 Performance Measurements

In this section, we give performance measurements and summarize simulation results with respect to the estimation efficiency as well as the variable selection accuracy. The sample sizes are set as $N=(20,100)\times 10^{3}$. Correspondingly, the number of workers is set to $K=(10,20)$.

For a reliable evaluation, we repeat the experiment $R=500$ times. Take the WLSE for example. For the $r$th replication, denote $\widehat{\theta}^{(r)}$ and $\widetilde{\theta}^{(r)}$ as the global estimator and WLSE, respectively. To measure the estimation efficiency, we calculate the root mean square error (RMSE) for the $j$th estimator as $\mbox{RMSE}_{\widetilde{\theta},j}=\{R^{-1}\sum_{r}\|\widetilde{\theta}_{j}^{(r)}-\theta_{0j}\|^{2}\}^{1/2}$. The RMSE for the global estimator $\widehat{\theta}^{(r)}$ can be defined similarly. Then, the relative estimation efficiency (REE) with respect to the global estimator is given by $\mbox{REE}_{j}=\mbox{RMSE}_{\widehat{\theta},j}/\mbox{RMSE}_{\widetilde{\theta},j}$ for $j=1,\cdots,p$. Next, for each $1\leq j\leq p$, we could construct a 95% confidence interval for $\theta_{0j}$ as CI${}_{j}^{(r)}=(\widetilde{\theta}_{j}^{(r)}-z_{0.975}\widetilde{\mbox{SE}}_{j}^{(r)},\widetilde{\theta}_{j}^{(r)}+z_{0.975}\widetilde{\mbox{SE}}_{j}^{(r)})$, where $\widetilde{\mbox{SE}}_{j}^{(r)}$ is the $j$th diagonal element of $\widehat{\Sigma}=(\sum_{k}\alpha_{k}\widehat{\Sigma}_{k}^{-1})^{-1}$, and $z_{\alpha}$ is the $\alpha$th quantile of a standard normal distribution. Then the coverage probability (CP) is computed as CP${}_{j}=R^{-1}\sum_{r=1}^{R}I(\theta_{j}\in\mbox{CI}_{j}^{(r)})$.

Next, based on the TWLSE, we further conduct shrinkage estimation on the master node. Let $\widehat{\mathcal{M}}^{(r)}$ be the set of selected variables in the $r$th replication using the DBIC. Correspondingly, $\widetilde{\theta}_{\lambda}^{(r)}$ is the shrinkage estimator. To measure the sparse discovery accuracy, we calculate the average model size (MS) as MS = $R^{-1}\sum_{r}|\widehat{\mathcal{M}}^{(r)}|$. Next, the percentage of the true model being correctly identified is given by CM$=R^{-1}\sum_{r}I(\widehat{\mathcal{M}}^{(r)}=\mathcal{M}_{T})$. In addition, to further investigate the estimation accuracy, we calculate the REE of the shrinkage estimation with respect to the global estimator as REE${}_{j}^{s}=\mbox{RMSE}_{\widehat{\theta},j}/\mbox{RMSE}_{\widetilde{\theta}_{\lambda},j}$ for $j\in\mathcal{M}_{T}$.

Lastly, for the high dimensional setting, we further evaluate the performances of four screening measures and post-estimation result. Specifically, we select top 40 features with highest screening measure values ($\widehat{\omega}_{j}$), which is denoted as $\widehat{\mathcal{M}}_{screen}^{(r)}$ in the $r$th replication. The coverage rate is defined as CR${}_{screen}=R^{-1}\sum_{r=1}^{R}I(\mathcal{M}_{T}\subset\widehat{\mathcal{M}}_{screen}^{(r)})$. Similarly, one could define CRselect for the post shrinkage estimation result. Due to the high dimensionality, we define an overall REE$=(\sum_{j}\mbox{RMSE}_{\widehat{\theta},j})/(\sum_{j}\mbox{RMSE}_{\widetilde{\theta},j})$ for the estimator $\widetilde{\theta}$.

4.4 Simulation Results

We compare the proposed DLSA method with (a) the OS estimator (Zhang et al., 2013), and (b) the CSL estimator (Jordan et al., 2019). Specifically, we denote the DLSA method with two-step estimation (TWLSE) as TDLSA method. The simulation results for the i.i.d setting is presented in Table 7–8 (in the Appendix). In addition, the results for the heterogenous setting are summarized in Table 1–3. First, in the i.i.d case, one can observe that all three methods are as efficient as the global estimator when $N$ is increased. For example, for the linear regression (i.e., Table 7), all the methods could achieve $\mbox{REE}\approx 1$ when $N=100,000$ and $K=20$. However, in the heterogeneous setting (i.e., Setting 2), the finite sample performances of the competing methods are quite different. The proposed DLSA and TDLSA method achieve higher efficiency than the other methods, which is also asymptotically efficient as the global estimator. Furthermore, the TDLSA method is more efficient than the DLSA method. For instance, the REEs of the TDLSA estimation for the logistic regression (i.e., Table 2) is near 1 in the second setting with $N=100,000$ and $K=20$, while the REE for $\theta_{1}$ of the OS, CSL, DLSA methods are approximately 0.77, 0.12, and 0.92, respectively. Although the OS estimator is less efficient, it is still consistent as $N$ increases. The CSL method performs worst under this situation. That is because it only uses the local Hessian matrix; this could result in a highly biased estimator.

Lastly, the CPs for both DLSA and TDLSA methods are all around 95%, while the CP for the CSL method is largely underestimated especially under the heterogenous setting. With respect to the shrinkage estimation, one can observe that the adaptive Lasso estimator is able to achieve higher estimation efficiency in the second setting than the global estimator. For example, the REE for the shrinkage DLSA (SDLSA) method could be even higher than 1 for Poisson model in Setting 2 (i.e., Table 3). In addition, we observe that the proposed DBIC does a great job at identifying the nonzero variables with high accuracy.

Next, we summarize the simulation results for the high dimensional setting in Table 10–11 for linear and logistic regression models. Specifically, the screening accuracies of four screening methods are presented in Table 10 and the post-estimation result is given in Table 11. Among all screening methods, the SIS and DC are able to achieve higher screening coverage rate. For the post-estimation result, both TDLSA and SDLSA methods have better estimation accuracy than other methods.

5 Application to airline data

For illustration purposes, we study a large real-world dataset. Specifically, the dataset considered here is the U.S. Airline Dataset. The dataset is available at http://stat-computing.org/dataexpo/2009. It contains detailed flight information about U.S. airlines from 1987 to 2008. The task is to predict the delayed status of a flight given all other flight information with a logistic regression model. Each sample in the data corresponds to one flight record, which consists of a binary response variable for delayed status (Delayed), and departure time, arrival time, distance of the flight, flight date, delay status at departures, carrier information, origin and destination as regressors. The complete variable details are described in Table 4. The data contain six continuous variables and five categorical variables. The categorical variables are converted to dummies with appropriate dimensions. We treat the Year and DayofMonth variables as numerical to capture the time effects. To capture possible seasonal patterns, we also convert the time variables Month and DayofWeek to dummies. Ultimately, a total of 181 variables are used in the model. The total sample size is 113.9 million observations. This leads to the raw dataset being 12 GB on a hard drive. After the dummy transformation described in Table 4, the overall in-memory size is over 52 GB, even if all the dummies are stored in a sparse matrix format. Thus, this dataset can hardly be handled by a single computer. All the numerical variables are standardized to have a mean of zero and a variance of one.

5.1 The Spark System and MLE

To demonstrate our method, we set up a Spark-on-YARN cluster on the Alibaba cloud server (https://www.alibabacloud.com/products/emapreduce). This is a standard industrial-level architecture setup for a distributed system.

The system consists of one master node and two worker nodes. Each node contains 32 virtual cores, 128 GB of RAM and two 80 GB SSD local hard drives. The dataset is stored on the Hadoop data file system (HDFS). We use such hardware specification to prevent the failure of Spark’s default algorithm due to the well-known out-of-memory issue for comparison purposes. It is worth mentioning that our DLSA algorithm works well on workers with less than 64GB RAM.

Because the system’s RAM is larger than the raw data size, one may wonder whether the logistic regression task can be run on a single node. Unfortunately, this is infeasible in practice. This is because much more memory (typically $>128$ GB with a double-precision floating-point format) is needed for operating on matrices of such a huge size. Even for the Spark system, a task of this magnitude cannot be directly performed using an existing algorithm library (e.g., Spark ML) unless each worker node is equipped with at least 128G RAM. This is because Spark is a very memory-intensive system. For example, to compute a single distributed matrix with a size of approximately 1 GB in memory, one might need each worker to have 2-5 GB of memory in practice. This overhead memory consumption grows significantly as the size of the data matrix increases. For discussions, see the Spark documentation at https://spark.apache.org/docs/latest/tuning.html#memory-tuning and Chen and Guestrin (2016).

We compared our approach with Spark’s builtin distributed algorithm (implemented in the spark.ml.classification.LogisticRegression module). If one insists on computing the traditional MLE based on the entire dataset with a single node, we employ a stochastic gradient descent (SGD) algorithm (Zhang, 2004) in the Python scikit-learn module (Pedregosa et al., 2011) to allow for a memory-constraint situation.

Fortunately, both the proposed DLSA and OS methods allow us to develop a user-friendly Spark algorithm with minimal computer resources. Our DLSA works well for a Spark system with only 64G RAM. As the algorithm is designed in a batch manner, it is highly efficient under memory constraints. The algorithm is developed with the Spark Python API (PySpark) and run on a Spark system (version $>$ 2.3) (see Algorithm 1 for details). It can be freely downloaded from https://github.com/feng-li/dlsa. We also provide the implementation of the aforementioned algorithms in Table 5 in the repository. We then use our algorithm to fit a logistic regression model for modelling the delayed status of a flight.

To this end, the entire dataset is randomly partitioned into 1139 subgroups. The sample size for each subgroup is approximately $100,000$. Next, for each subgroup of data, we create a virtual worker (i.e., an executor in the Spark system) so that the computation for each worker can be conducted in a parallel manner. By doing so, the computation power of the entire Spark system can be maximized for both DLSA and OS methods. The corresponding log-likelihood values, computational time as well as the estimated communication time are reported in Table 5, respectively. We remark that the computing time for the MLE includes the data shuffling time with our DLSA and OS algorithms. This serves as an important benchmark to gauge the performance of the other competing methods (e.g., DLSA and OS methods). Comparing these results with that of the traditional MLE, we find that the traditional MLE is extremely difficult to compute without a distributed system.

Table 5 also depicts that the log-likelihood value of the DLSA is the best with an affordable RAM consumption. The Spark’s ML module takes longer time than DLSA because DLSA only requires one round communication and the parameter optimization is done within each executor. A serialized SGD takes more that $15$ hours and obtains an inferior result (i.e., smaller log-likelihood value).

5.2 Variable Selection Results with BIC

We next apply the proposed shrinkage DLSA method (referred to as SDLSA) with the BIC criterion to conduct variable selection. It is remarkable that this can be fully conducted on the master, and no further communication is needed. It takes only $0.2$ seconds to accomplish the task. After the shrinkage estimation, we are able to reduce the 181 variables to 157 variables.

The detailed results are summarized in Table 6. First, with respect to time effects, both yearly and seasonal trends are found. The coefficient for the Year is 6.12, which implies that as the year proceeds, the airline delays become more severe. Next, the passengers are likely to encounter delays in May, October, November, and December (with coefficients of 6.03, 0.8, 0.4, and 0.19, respectively). In terms of days of the week, more delays are expected for certain working days, i.e., Tuesday, Wednesday and Friday (with coefficients of 0.6, 0.85 and 0.39, respectively) compared to weekends. Finally, within a given day, we find the coefficients for both the scheduled departure time (CRSDepTime, $-0.12$) and the scheduled arrival time (CRSArrTime, $-0.13$) are negative, indicating that late departing and arriving flights suffer less so from delays.

Next, with respect to the airline carriers, AA (American Airlines), DL (Delta Air Lines), NW (Northwest Airlines), and UA (United Airlines) have estimated coefficients of 0.49, 0.18, 0.39, and 0.70, which indicates how more likely they are to be delayed compared with other airlines. In addition, one may be interested in with which airports are more likely to have delayed flights. According to our estimation results, the top five origin airports that cause delays are IAH (George Bush Intercontinental Airport), LGA (LaGuardia Airport), PHL (Philadelphia International Airport), RDU (Raleigh-Durham International Airport), ONT (Ontario International Airport), and SMF (Sacramento International Airport), with coefficients of $0.82$, $0.87$, $0.94$, $1.58$, and $1.59$, respectively. The top five destination airports that cause delays are PBI (Palm Beach International Airport), MCI (Kansas City International Airport), DCA (Ronald Reagan Washington National Airport), SAN (San Diego International Airport), and MEM (Memphis International Airport), with coefficients of $0.99$, $1.00$, $1.07$, $1.15$, and $1.16$, respectively.

[FIGURE:]

6 Concluding remarks

In this article, we develop a novel DLSA algorithm that is able to perform large-scale statistical estimation and inference on a distributed system. The DLSA method can be applied to a large family of regression models (e.g., logistic regression, Poisson regression, and Cox’s model). First, it is shown that the DLSA estimator is as statistically optimal as the global estimator. Moreover, it is computationally efficient and only requires one round of communication.

Furthermore, we develop the corresponding shrinkage estimation by using an adaptive Lasso approach. The oracle property is theoretically proven. A new DBIC measure for distributed variable selection, which only needs to be performed on the master and requires no further communication, is designed. We prove the DBIC measure to be selection consistent. Finally, numerical studies are conducted with five classical regression examples. In addition, a Spark toolbox is developed, which is shown to be computationally efficient both through simulation and in airline data analysis.

To facilitate future research, we now discuss several interesting topics. First, the DLSA method requires the objective function to have continuous second-order derivatives. This assumption might be restrictive and cannot be satisfied for certain regression models, e.g., the quantile regression. Consequently, the relaxation of this assumption can be investigated, and corresponding distributed algorithms should be designed for such regression models. Second, the dimension considered in our framework is finite. As a natural extension, one could study the shrinkage estimation properties in high-dimensional settings. Furthermore, as we commented on our pre-feature screening procedure, further post-selection statistical inference should be investigated (Berk et al., 2013; Taylor and Tibshirani, 2015; Lee et al., 2016). Third, the algorithm is designed for independent data. In practice, dependent data (e.g., time series data and network data) are frequently encountered. It is thus interesting to develop corresponding algorithms by considering the dependency structure. Lastly, note that the penalty term given in (3.1) can be added to other surrogate loss functions as in Jordan et al. (2019). Hence the theoretical analysis in Section 3 as a great potential to extend to other methodologies.

7 Acknowledgements

The authors thank the editor, associate editor, and two referees for their insightful comments that have led to significant improvement of this article. We also thank Dr. Zhihui Jin for his constructive suggestions.

Appendix A

Appendix A.1: Proof of Proposition 1

Note that $\widetilde{\theta}-\theta_{0}$ takes the form

[TABLE]

Define $\widehat{\Sigma}_{k}(\theta)=\{\partial^{2}\mathcal{L}_{k}(\theta)/\partial\theta\partial\theta^{\top}\}^{-1}$. In the following section, we denote $\widehat{\Sigma}_{k}$ by $\widehat{\Sigma}_{k}(\widehat{\theta}_{k})$ to make it clearer. By Slutsky’s Theorem, to prove (2.4), it suffices to verify that

[TABLE]

where $\mbox{cov}\{V^{*}(\theta_{0})\}=\Sigma^{-1}$ and $B^{*}(\theta_{0})=O_{p}(K/\sqrt{N})$. We prove them in the following.

1. Proof of (A.1). Recall that $\widehat{\theta}_{k}$ is a $\sqrt{n}$-consistent estimator of $\theta_{0}$. This enables us to conduct a Taylor’s expansion of $\widehat{\Sigma}_{k}^{-1}(\widehat{\theta}_{k})$ at $\theta_{0}$, which yields

[TABLE]

where $\theta^{*}$ lies on the line joining $\theta_{0}$ and $\widehat{\theta}_{k}$. By Condition (C5), we have \frac{\partial^{3}\mathcal{L}_{k}(\theta)}{\partial\theta_{j}\partial\theta\partial\theta^{\top}}\Big{|}_{\theta=\theta^{*}}=O_{p}(1). Therefore, the order of the first term is $O_{p}(1/\sqrt{n}_{k})$. In addition, we have $\widehat{\Sigma}_{k}^{-1}(\theta_{0})-\Sigma_{k}^{-1}=\widehat{\Sigma}_{k}^{-1}(\theta_{0})-E\{\widehat{\Sigma}_{k}^{-1}(\theta_{0})\}=O_{p}(n_{k}^{-1/2})$. Consequently, it can be derived that $\widehat{\Sigma}_{k}^{-1}(\widehat{\theta}_{k})-\Sigma_{k}^{-1}=O_{p}(n_{k}^{-1/2})$. Further note that $\alpha_{k}=n_{k}/N$ and $\sum_{k}\alpha_{k}=1$. Then, we have

[TABLE]

Hence (A.1) is proven.

2. Proof of (A.2). Recall that $\widehat{\theta}_{k}$ is the local minimizer of $\mathcal{L}_{k}(\theta)$. Therefore, it holds that

[TABLE]

where $\theta^{*}$ lies between $\theta_{0}$ and $\widehat{\theta}_{k}$. By standard arguments,

[TABLE]

Further note that $\widehat{\theta}_{k}-\theta_{0}=O_{p}(n_{k}^{-1/2})$. Then, we have

[TABLE]

where $B_{k}(\theta_{0})=O(1)$ is the bias term. Then, it holds that

[TABLE]

where the second equation is implied by $\sqrt{N}\sum_{k}\alpha_{k}\{\widehat{\Sigma}_{k}^{-1}(\widehat{\theta}_{k})-\Sigma_{k}^{-1}\}(\widehat{\theta}_{k}-\theta_{0})=O_{p}(K/\sqrt{N})$. By condition (C4), it can be concluded that $\mbox{cov}\{\frac{1}{\sqrt{N}}\sum_{k}n_{k}\frac{\partial\mathcal{L}_{k}(\theta)}{\partial\theta}|_{\theta=\theta_{0}}\}=\Sigma^{-1}$. Consequently, (A.2) can be proven.

Appendix A.2: Proof of Theorem 1

By Slutsky’s Theorem, the asymptotic normality is directly implied by (A.1) and (A.2) with $V^{*}(\theta_{0})\rightarrow_{d}N(0,\Sigma^{-1})$ and $B^{*}(\theta_{0})=o_{p}(1)$. First, by Condition (C6) and the Lyapunov central limit theorem, we have $V^{*}(\theta_{0})\rightarrow_{d}N(0,\Sigma^{-1})$. Next, by the condition that $n\gg\sqrt{N}$, we have $K\gg\sqrt{N}$, and thus, $B^{*}(\theta_{0})=o_{p}(1)$.

Appendix A.3: Proof of Theorem 2

Note that we have

[TABLE]

It suffices to show that,

[TABLE]

where $\mbox{cov}\{V_{2}^{*}(\theta_{0})\}=\Sigma^{-1}$ and $B_{2}^{*}(\theta_{0})=O_{p}(1/(n^{2}\sqrt{N}))$. We prove them in the following three steps.

Step 1. $\widehat{\theta}_{k}^{(2)}-\theta_{0}=O_{p}(n^{-1/2})$.

Note that we have

[TABLE]

By performing a Taylor’s expansion of $\partial\mathcal{L}_{k}(\widetilde{\theta})/\partial\theta$ at $\theta_{0}$, we could obtain,

[TABLE]

where $\overline{\theta}$ lies on the line joining $\widetilde{\theta}$ and $\theta_{0}$. As a result, we have

[TABLE]

Using Proposition 1 we have $\Delta_{k1}=O_{p}(n^{-1/2})$. Then it suffices to derive the rate of $\Delta_{k2}$.

By using Taylor’s expansion again we have

[TABLE]

where $\overline{\theta}^{(2)}$ lies on the line joining $\overline{\theta}$ and and $\widetilde{\theta}$. By Proposition 1 we can obtain ${\partial^{3}\mathcal{L}_{k}(\overline{\theta}^{(2)})}/{\partial\theta\partial\theta^{\top}\partial\theta_{j}}\rightarrow_{p}E\{{\partial^{3}\mathcal{L}_{k}(\theta_{0})}/{\partial\theta\partial\theta^{\top}\partial\theta_{j}}\}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathcal{L}_{kj}^{(3)}$ and $\partial^{2}\mathcal{L}_{k}(\widetilde{\theta})/{\partial\theta\partial\theta^{\top}}\rightarrow_{p}\Sigma_{k}^{-1}$. By Proposition 1, we have

[TABLE]

where $V^{*(j)}(\theta_{0})=\{V(\theta_{0})B_{j}(\theta_{0})+B(\theta_{0})V_{j}(\theta_{0})\}/(n\sqrt{N})=O_{p}(1/(n\sqrt{N}))$, and $B^{*(j)}(\theta_{0})=B(\theta_{0})B_{j}(\theta_{0})=O_{p}(n^{-2})$.

Step 2. Proof of (A.5).

Following the same procedure of proving (A.1) and using $\widehat{\theta}_{k}^{(2)}-\theta_{0}=O_{p}(n^{-1/2})$ proved in the first step we could obtain the result.

Step 3. Proof of (A.6).

By (A.8) and (A.9) we have

[TABLE]

Define

[TABLE]

Then we have $\mbox{cov}\{V_{2}^{*}(\theta_{0})\}=\Sigma^{-1}$, $V_{22}(\theta_{0})=o_{p}(1)$, and $B_{2}^{*}(\theta_{0})=O_{p}(n^{-2}N^{1/2})$ by Step 2.

Step 4. Proof of (A.7).

by Condition (C6) and the Lyapunov central limit theorem, we have $V_{2}^{*}(\theta_{0})\rightarrow_{d}N(0,\Sigma^{-1})$. Under the condition that $n^{2}/N^{1/2}\rightarrow\infty$, we conclude that $B_{2}^{*}(\theta_{0})=o_{p}(1)$. Using the Slutsky’s Theorem we could obtain the result.

Appendix A.4: Proof of Proposition 2

The result can be proved by induction method. By Theorem 2, the bias term of $\widetilde{\theta}^{(2)}$ reduces from $O_{p}(n^{-1})$ to $O_{p}(n^{-2})$ when we take one more round of iteration. Applying the proof of Theorem 2 sequentially for $m$ times, we can show that the bias term of $\widetilde{\theta}^{(m)}$ is in the order of $O_{p}(n^{-m})$. If $n^{-m}\ll N^{-1/2}$, i.e., $m\gg\log N/\log n$, then the bias term can be ignored. As a result, following the same proof procedure of Theorem 2, we can obtain the result.

Appendix B

Appendix B.1: Proof of Theorem 3

1. Proof of $\sqrt{N}$-consistency.

Note that the objective function $Q_{\lambda}(\theta)$ in (3.1) is a strictly convex function. Then, the local minimizer is also a global minimizer. To establish $\sqrt{N}$-consistency results, it suffices to verify the following result (Fan and Li, 2001), i.e., for an arbitrarily small $\epsilon>0$, there exists a sufficiently large constant $C$ such that

[TABLE]

Let $u=(u_{1},\cdots,u_{p})^{\top}$ and \widehat{\Delta}_{N}=\sum_{k}\alpha_{k}\widehat{\Sigma}_{k}^{-1}(\widehat{\theta}_{k})\big{\{}\sqrt{N}(\theta_{0}-\widehat{\theta}_{k})\big{\}}, Then, we have

[TABLE]

where the second equality holds because we assume $\theta_{0j}=0$ for $j>d_{0}$. Further note that we assume that $\sqrt{N}a_{N}\rightarrow_{p}0.$ Consequently, the last term (B.2) is $o_{p}(1)$. Next, note that the first term of (B.2) is lower bounded by $\lambda_{\max}^{-1}(\widehat{\Sigma})C^{2}$ because $\|u\|=C$. By (A.3), we have $\lambda_{\max}(\widehat{\Sigma})\rightarrow_{p}\lambda_{\max}(\Sigma)$. Consequently, with probability tending to 1, we have the first term uniformly larger than $0.5\lambda_{\max}^{-1}(\Sigma)C^{2}$, which is positive due to Condition (C4). In addition, by $K/\sqrt{N}\rightarrow 0$, we have $\widehat{\Delta}_{N}=O_{p}(1)$. Consequently, as long as $C$ is sufficiently large, the first term will dominate the last two terms. Then, the result of (B.1) is proven.

2. Proof of Selection Consistency.

It suffices to verify that $P(\widetilde{\theta}_{\lambda,j}=0)\rightarrow 1$ for any $d_{0}<j\leq p$. Note that $Q_{\lambda}(\theta)$ can be rewritten as

[TABLE]

where $C$ is a constant. Define $\widehat{\Omega}=\widehat{\Sigma}^{-1}$, and $\widehat{\Omega}^{(j)}$ denotes the $j$th row of the matrix $\widehat{\Omega}$. If $\widetilde{\theta}_{\lambda,j}\neq 0$ for some $j>d_{0}$, then the partial derivative can be calculated as

[TABLE]

Note that $\widehat{\Omega}\rightarrow_{p}\Sigma^{-1}$ and $\sqrt{N}(\widetilde{\theta}_{\lambda}-\widetilde{\theta})=\sqrt{N}(\widetilde{\theta}_{\lambda}-\theta_{0})-\sqrt{N}(\widetilde{\theta}-\theta_{0})=O_{p}(1)$, by (A.1), Theorem 1, and Theorem 3 (a). Consequently, the first term (B.3) is $O_{p}(1)$. Next, by this condition, we know that $\sqrt{N}\lambda_{j}\geq\sqrt{N}b_{\lambda}\rightarrow\infty$ for $j>d_{0}$. Because $\widetilde{\theta}_{\lambda,j}\neq 0$, we have $\mbox{sign}(\widetilde{\theta}_{\lambda,j})=1$ or -1; thus, the second term (B.3) goes to infinity. Obviously, the equation will not be equal to zero. This implies $P(\widetilde{\theta}_{\lambda,j}=0)\rightarrow 1$ as a result.

Appendix B.2: Proof of Theorem 4

We first rewrite the asymptotic covariance $\Sigma$ into the following block matrix form:

[TABLE]

where $\Sigma_{11}\in\mathbb{R}^{d_{0}\times d_{0}}$. Similarly, we partition its inverse matrix $\Omega$ into 4 corresponding parts, $\Omega=(\Omega_{11},\Omega_{12};\Omega_{21},\Omega_{22})$. By Theorem 3, with probability tending to 1, we have $\widetilde{\theta}_{\lambda}^{(-\mathcal{M}_{T})}=0$. Therefore, $\widetilde{\theta}_{\lambda}^{(\mathcal{M}_{T})}$ should be the global minimizer of the objective function,

[TABLE]

By Theorem 3, it can be concluded that with probability tending to 1, $\widetilde{\theta}_{\lambda}^{(\mathcal{M}_{T})}$ should be nonzero (otherwise, the $\sqrt{N}$-consistency result in Theorem 3 will not hold). As a result, The partial derivative $\partial Q_{\lambda}(\theta)/\partial\theta_{j}$ should exist for $1\leq j\leq d_{0}$, which yields

[TABLE]

where $D(\widetilde{\theta}_{\lambda}^{(\mathcal{M}_{T})})$ is a $d_{0}$-dimensional vector, with its $j$th component given by $0.5\lambda_{j}\mbox{sign}(\widetilde{\theta}_{\lambda,j})$. By (B.4), it can be derived that

[TABLE]

where the second equality follows because $\sqrt{N}\widetilde{\theta}^{(-\mathcal{M}_{T})}=O_{p}(1)$ by Theorem 1, $\widehat{\Omega}_{11}\rightarrow_{p}\Omega_{11}$ and $\widehat{\Omega}_{12}\rightarrow_{p}\Omega_{12}$ by (A.1), and $\sqrt{N}\lambda_{j}=o_{p}(1)$ ($1\leq j\leq d_{0}$) by Theorem 3. Furthermore, by the matrix inverse formula, it holds that $\Omega_{11}^{-1}\Omega_{12}=-\Sigma_{12}\Sigma_{22}^{-1}$. Consequently, we have

[TABLE]

By Theorem 1, we have that the above is asymptotically normal with a mean of 0, and the inverse asymptotic covariance matrix is given by $(\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21})^{-1}=\Omega_{11}$. By condition (C6), we have $\Omega_{11}=\Omega_{\mathcal{M}_{T}}$. Consequently, the estimator $\widetilde{\theta}_{\lambda}^{(\mathcal{M})}$ shares the same asymptotic distribution with the oracle estimator $\widehat{\theta}_{\mathcal{M}_{T}}^{(\mathcal{M}_{T})}$.

Appendix B.3: Proof of Theorem 5

To establish the selection consistency property of the DBIC, we consider the following two cases for any $\mathcal{M}_{\lambda}\neq\mathcal{M}_{T}$. The first case is the underfitted case, and the second case is the overfitted case.

1. Underfitted Model. Note that $\lambda_{N}$ satisfies the condition in Theorem 3. Consequently, we have that $\widetilde{\theta}_{\lambda_{N}}$ is $\sqrt{N}$-consistent. We thus also have $\mbox{DBIC}_{\lambda_{N}}=o_{p}(1)$. For $\mathcal{M}\not\supset\mathcal{M}_{T}$, it can be derived that

[TABLE]

Define $\widetilde{\theta}_{\mathcal{M}}=\arg\min_{\theta\in\mathbb{R}^{p}:\theta_{j}=0,\forall j\not\in\mathcal{M}}(\theta-\widetilde{\theta})^{\top}\widehat{\Sigma}^{-1}(\theta-\widetilde{\theta})$ as the unpenalized estimator given the model identified by $\mathcal{M}$. Consequently, by definition, we should have

[TABLE]

where the last convergence is due to (A.1). Because $\Sigma$ is positive definite by Condition (C4), we have $\min_{\mathcal{M}\not\supset\mathcal{M}_{T}}(\widetilde{\theta}_{\mathcal{M}}-\widetilde{\theta})^{\top}\Sigma^{-1}(\widetilde{\theta}_{\mathcal{M}}-\widetilde{\theta})>0$ with probability tending to 1. One could then conclude immediately that $P(\inf_{\lambda\in\mathbb{R}_{-}}\mbox{DBIC}_{\lambda}>\mbox{DBIC}_{\lambda_{N}})\rightarrow 1$.

2. Overfitted Model. We next consider the overfitted model. In contrast, let $\lambda$ be an arbitrary tuning parameter that over selects the parameters. We then have $df_{\lambda}-d_{0}\geq 1$. It can then be concluded that $N(\mbox{DBIC}_{\lambda}-\mbox{DBIC}_{\lambda_{N}})=$

[TABLE]

First note that for $\mathcal{M}\supset\mathcal{M}_{T}$, we have that $\widetilde{\theta}_{\mathcal{M}}$ is $\sqrt{N}$-consistent. As a result, the first term of (B.6) is $O_{p}(1)$. Similarly, by Theorem 3, $\widetilde{\theta}_{\lambda_{N}}$ is $\sqrt{N}$-consistent, Thus, the second term (B.6) is also $O_{p}(1)$. As a result, (B.6) diverges to infinity as $N\rightarrow\infty$. This implies that $P(\inf_{\lambda\in\mathbb{R}_{+}}\mbox{DBIC}_{\lambda}>\mbox{DBIC}_{\lambda_{N}})\rightarrow 1$. This completes the proof.

Appendix B.4: Comparison of the Shrinkage Estimations

One could consider an alternative approach to conduct variable selection. Based on the asymptotic theory given in Theorem 1, we could perform simultaneous hypothesis testing for the variable selections. Specifically, for each $1\leq j\leq p$, we could construct a 95% confidence interval for $\theta_{0j}$ as CI${}_{j}=(\widetilde{\theta}_{j}-z_{0.975}\widetilde{\mbox{SE}}_{j},\widetilde{\theta}_{j}+z_{0.975}\widetilde{\mbox{SE}}_{j})$, where $\widetilde{\mbox{SE}}_{j}$ is the $j$th diagonal element of $\widehat{\Sigma}=(\sum_{k}\alpha_{k}\widehat{\Sigma}_{k}^{-1})^{-1}$, and $z_{\alpha}$ is the $\alpha$th quantile of a standard normal distribution. As a result, if $0\not\in\mbox{CI}_{j}$, we could identify the $j$th variable as an important variable.

We compare the finite sample performance of the above variable selection method with the proposed shrinkage estimation method for the logistic regression model. We set $p=30$, and for $1\leq j\leq 15$, we specify $\theta_{0j}=U_{1}\mbox{sign}(U_{2})$, where $U_{1}$ and $U_{2}$ are generated from uniform distribution $U[0.5,2]$ and $U[-0.2,0.8]$ respectively. In addition, for $15<j\leq p$, we set $\theta_{0j}=0$. We also denote the percentage for covering true model and identifying the true model as CR and CM respectively. As shown in Table 12, both methods are able to cover the true model with high accuracy. However, the hypothesis testing method tends to over select the important features while the SDLSA method could still guarantee the variable selection consistency property.

Appendix B.5: Interleave with Bootstrap Methods

An alternative way to deal with statistical inference for massive dataset is to use bootstrap methods. Recently, a surge of researches have investigated related methodologies. For instance, Kleiner et al. (2014) and Sengupta et al. (2016) propose to use subsample and resample methods to conduct statistical inference to reduce computational cost of traditional bootstrap methods. Volgushev et al. (2019) and Yu et al. (2020) further discuss how to incorporate the bootstrap methods into distributed learning algorithms. Particularly, we note that if we have a large number of workers, we could adopt the idea of Volgushev et al. (2019) to transmit estimators $\widehat{\theta}_{k}$ to the master and estimate the asymptotic covariance $\widehat{\Sigma}$ as $\widehat{\Sigma}=K^{-1}\sum_{k=1}^{K}(\widehat{\theta}_{k}-\overline{\theta})(\widehat{\theta}_{k}-\overline{\theta})^{\top}$, where $\overline{\theta}=K^{-1}\sum_{k=1}^{K}\widehat{\theta}_{k}$. The advantage here is to avoid the communication of the second order derivative of $\mathcal{L}_{k}(\theta)$. However, the estimation requires a large number of workers to be valid and in addition, the estimated asymptotic covariance is for the one-shot estimator. Hence, the methodology cannot be directly applied especially for our heterogenous data setting.

As an alternative, we can interleave the bootstrap methods on workers for the estimation of the covariance $\widehat{\Sigma}_{k}$. The method can be useful when it is hard to obtain an analytical form of the asymptotic covariance $\widehat{\Sigma}_{k}$. To show that it is a possible direction for future study, we conduct a simulation study for the logistic regression model. Specifically, we use the subsampled double bootstrap (SDB) estimation method on workers to obtain the covariance estimator $\widehat{\Sigma}_{k}^{boot}$. Next, we substitute $\widehat{\Sigma}_{k}^{boot}$ to $\widehat{\Sigma}_{k}$ in (2.2) to obtain the bootstraped-WLSE (BWLSE) estimator. We set the bootstrap sample size $b=n^{\gamma}$ with $\gamma=0.9$ and the resampling times $R=n_{k}$. The estimation REE with respect to global estimator is presented in Table 9. Under the i.i.d data setting, all the three methods are comparable when $N=20,000$ and $K=10$. However, under the heterogenous data setting, one could observe that the bootstrap method is better than the OS method and is comparable to the DLSA method.

Appendix B.6: Supplementary results for simulation studies