# Two-Stage Robust and Sparse Distributed Statistical Inference for Large-
Scale Data

Emadaldin Mozafari-Majd, Visa Koivunen The authors are with the Department of
Signal Processing and Acoustics, Aalto University, FI-00076 Aalto, Finland
(email<EMAIL_ADDRESS>visa.koivunen@aalto.fi).

###### Abstract

In this paper, we address the problem of conducting statistical inference in
settings involving large-scale data that may be high-dimensional and
contaminated by outliers. The high volume and dimensionality of the data
require distributed processing and storage solutions. We propose a two-stage
distributed and robust statistical inference procedures coping with high-
dimensional models by promoting sparsity. In the first stage, known as model
selection, relevant predictors are locally selected by applying robust Lasso
estimators to the distinct subsets of data. The variable selections from each
computation node are then fused by a voting scheme to find the sparse basis
for the complete data set. It identifies the relevant variables in a robust
manner. In the second stage, the developed statistically robust and
computationally efficient bootstrap methods are employed. The actual inference
constructs confidence intervals, finds parameter estimates and quantifies
standard deviation. Similar to stage 1, the results of local inference are
communicated to the fusion center and combined there. By using analytical
methods, we establish the favorable statistical properties of the robust and
computationally efficient bootstrap methods including consistency for a fixed
number of predictors, and robustness. The proposed two-stage robust and
distributed inference procedures demonstrate reliable performance and
robustness in variable selection, finding confidence intervals and bootstrap
approximations of standard deviations even when data is high-dimensional and
contaminated by outliers.

###### Index Terms:

High-dimensional, large-scale data, big data, sparsity, robust estimator,
Lasso, distributed computation and storage, fixed-point equations, information
fusion.

##  I Introduction

Massive quantities of ubiquitous and heterogeneous data are generated by
social media, smart phones, IoT, environmental monitoring, astronomical
imaging devices and financial markets. Harnessing information from such large-
scale data provides enterprises with meaningful insights into their
performance and offers tremendous business opportunities. However, these
benefits come with formidable challenges in handling storage, processing,
acquisition and privacy concerns of high-speed and high volume data [1]. In
order to remedy the storage and processing issues, distributed computation and
storage solutions are preferred. In a variety of statistical inference
applications, one is dealing with high-dimensional data where the number of
explaining variables $p$ may be comparable or much larger than the number of
observations $n$. Often, high-dimensional problems exhibit a lower dimensional
structure such as sparsity or low-rank. Regularization may be necessary to
address the ill-posed high-dimensional problems and capture the parsimonious
representation. In the context of linear regression, regularization by
$\ell_{1}$-norm is known as the celebrated Lasso and performs simultaneous
parameter estimation and model selection [2]. It improves the prediction
accuracy by introducing some bias while reducing the variance. There are
certain scenarios where regularization by $\ell_{1}$-norm may cause
significant bias in estimated coefficients. Moreover, the exact
characterization of the limiting distribution for the Lasso estimator is a
challenging task. In recent years, several valid statistical inference
procedures have been introduced to characterize the construction of confidence
intervals and hypothesis testing for high-dimensional problems [3], [4], [5].
In general, the state-of-art procedures to address the uncertainty associated
with parameter estimates belong to three main categories. First category
concerns inference procedures based on bootstrapping. However, the
conventional bootstrap methods fail to provide a reliable approximation to the
distribution of Lasso estimator. In order to address this issue, Chatterjee
and Lahiri [3, 6] proposed two alternative solutions, modified residual
bootstrap Lasso and residual bootstrap adaptive Lasso. These solutions
consistently estimate their limiting distributions and provide valid
approximation of confidence intervals. The second category includes post-Lasso
inference methods where the first stage involves model selection using Lasso
and the actual inference is made in the second stage using the selected
variables. This category includes sample splitting [7], bootstrap Lasso-OLS
[5], bootstrap Lasso-partial Ridge [8] and post-selection inference [9, 10]
methods. Herein, we extend the concept of Post-Lasso estimator to conduct
statistical inference for large-scale data sets using distributed computation
and storage. The third category is based on the debiased-Lasso method [4],
[11] where the key idea is to remove the bias introduced by regularization
from Lasso solution. These methods offer a concrete and general framework to
quantify the uncertainty associated with parameters in high-dimensional
settings, e.g. hypothesis testing and construction of confidence intervals.
Other alternatives are covariance test [12], knockoff filter [13], the ridge
projection and bias correction [14].

Another challenge in dealing with large-scale data is that the probability of
observing outliers may increase as the dimensionality ($p$) and sample size
($n$) grows larger. Outliers may become masked and hence difficult to detect.
Outliers may severely deteriorate the performance of ordinary least square and
regularized least square estimators. Robust multivariate statistical
procedures are required to cope with outlying observations and ensure the
veracity of estimation, classification and decision making. In the context of
linear regression, robust regularized estimators are employed to ensure
robustness and find sparse solutions simultaneously [15, 16, 17, 18, 19, 11].

Modern statistical inference procedures need to accommodate distributed
storage and parallel computations to deal with high volume and high-
dimensional data. Bootstrap is a powerful tool for quantifying the uncertainty
of estimates, i.e., confidence intervals and hypothesis tests. However, the
applicability of conventional bootstrap techniques in large-scale settings may
not be feasible because of computational constraints. The bootstrap resamples
may have the same size as the original large-scale data, and repeating
estimation for each bootstrap replicate becomes prohibitively expensive. The
Bag of Little Bootstraps (BLB) [20] offers a scalable and computationally
efficient inference procedure for quantification of the uncertainty associated
with estimates that accommodates distributed and parallel computing
architectures. It subdivides the complete large-scale data into smaller
distinct subsets. It can be considered as sampling without replacement from
the complete data set. Then, bootstrap is applied to each subset and combines
the inference results from each subset. However, the BLB procedure is highly
sensitive to outliers. In order to overcome this issue, a statistically robust
BLFRB approach that extends the idea of fast fixed-point computations of FRB
method in [21] to MM-estimators was introduced in [22]. The modified bootstrap
replicates are calculated by applying a linear correction factor to the one-
step approximation of bootstrap replicates. Despite its robustness and
computational efficiency, BLFRB performs unreliably in high-dimensional
settings.

This paper focuses on statistical inference in large-scale settings with
distributed architecture where data may be high-dimensional and contaminated
by outliers. This problem has not been thoroughly studied in open literature.
In particular, we propose two-stage statistical inference procedures that are
robust to outlying observations and allow for using distributed storage and
processing architectures for scalability. In the first stage, the relevant
predictors are selected in two steps. First robust Lasso estimator is applied
to distinct subsets of data in order to perform local variable selection. The
variables for the whole data set are then found by applying a fusion rule to
the selections from individual nodes at the fusion center or cloud. In the
second stage, we conduct inference by constructing confidence intervals,
finding point estimates of the selected parameters and their standard
deviations based on the large-scale data. The developed distributed and low-
complexity inference procedures that use linearly corrected one-step robust
estimators are employed. In special cases with very high dimensionality
($p/n\approx 1$ or $p\gg b$, $p$ number of predictors, $n$ sample size and $b$
subsample size), one may accelerate the model selection by a preprocessing
stage excluding the majority of irrelevant variables via a robust variable
screening procedure on the distinct subsets of data stored at each node. We
address this issue in Supplemental Materials. The methods proposed in this
paper extend our previous work on two-stage robust and distributed inference
procedures [23], [24] by deriving computationally more efficient estimation
methods and establishing statistical properties of the inference methods using
analytical tools. We emphasize that our asymptotic analysis is restricted to
the classical fixed $p$ setting. In [5], the asymptotic properties were
established for a special case of post-Lasso estimator Lasso+mLS and the valid
bootstrap approximation while allowing $p$ to grow at an exponential rate in
$n$. However, their analysis does not cover the distributed and robust
settings.

The main new contributions of the paper are summarized as follows:

  * •

Two-Stage Robust and Distributed inference employing the class of
$\tau$-estimators called TSRD-$\tau$ is introduced. A robust $\tau$-Lasso
sparsity promoting estimator is employed in the first stage.

  * •

The proposed TSRD-$\tau$ employs a novel Robust and Scalable linearly
corrected One-step Bootstrap procedure using $\boldsymbol{\tau}$-estimator
(RSOB-$\tau$) for performing the actual inference. Its computational
complexity is reduced by efficient estimation of bootstrap replicates.

  * •

Two-Stage Robust and Distributed inference employing the class of MM-
estimators called TSRD-MM is developed. The sparsity promoting estimator in
the first stage is robust MM-Lasso. It extends the BLFRB procedure [22] by a
robust variable selection stage promoting sparsity.

  * •

Analytical results proving the robustness and consistency of the TSRD-$\tau$
method are derived for fixed $p$ and diverging $n$. In order to formally show
the quantitative robustness of $\tau$-Lasso, its finite-sample breakdown is
characterized.

  * •

Extensive simulations are conducted to assess the performance of robust and
distributed two-stage procedures in variable selection, bootstrap estimation
of standard deviation, robustness of confidence intervals and computational
complexity of the RSOB-$\tau$ procedure.

  * •

Analytical results on the consistency and asymptotic normality of the
RSOB-$\tau$ procedure are verified through computer simulations. Furthermore,
the favorable theoretical findings on robustness of bootstrap replications are
confirmed by extensive simulation studies and comparing them to [22].

The two-stage inference solutions presented in this work achieve scalability
by performing inference over smaller distinct subsets of data in parallel,
using multinomial weighting and discarding irrelevant predictors. Considerable
computational gains are achieved by using the low-complexity procedures to
calculate bootstrap replications while retaining the consistency. The model
selection stage facilitates preventing the impact of undesirable bias
introduced by regularization and allows for the inference free of
regularization parameter tuning.

This paper is organized as follows: In section 2, we describe the proposed
two-stage inference methods and the employed data model. In section 3, the
robust model selection methods and the distributed inference procedures are
explained in more detail. In section 4, theoretical characterization of
finite-sample breakdown point of robust $\tau$-Lasso is provided. Moreover,
the details of RSOB-$\tau$ are explained and its consistency and robustness
properties are established. Section 5 studies the performance of the two-stage
robust inference procedures in simulations. Their consistency, robustness of
confidence intervals, computational complexity, variable selection and
standard deviation are considered. Section 6 concludes the paper. Detailed
derivations and explanations of the theorems and their proofs can be found in
the Appendix and Supplemental Material. Moreover, we present a robust variable
screening procedure [25] prior to model selection, aiming to reduce the
computational complexity in very high-dimensional data sets. We refer the
interested readers to Supplemental Materials.

##  II Overview

In this section, we will define the employed data model and briefly describe
the proposed two-stage inference procedures.

###  II-A Data Model

Consider a large-scale data with $n$ independent and identically distributed
(i.i.d.) observations following a linear regression model

$\mathbf{y}=\mathbf{X}\boldsymbol{\beta}_{0}+\mathbf{v}$ (1)

where
$\mathbf{X}=(\mathbf{x}_{[1]},\mathbf{x}_{[2]},\ldots,\mathbf{x}_{[n]})^{T}\in\mathbb{R}^{n\times
p}$ denotes a regression matrix,
$\mathbf{y}=(y_{1},y_{2},\ldots,y_{n})^{T}\in\mathbb{R}^{n}$ is a response
vector, $\mathbf{v}=(v_{1},v_{2},\ldots,v_{n})^{T}\in\mathbb{R}^{n}$ is a
measurement noise and the errors $v_{l}$ are assumed to be independent of
predictors $\mathbf{x}_{[l]}$. $\boldsymbol{\beta}_{0}\in\mathbb{R}^{p}$
denotes a sparse parameter vector with $k_{s}=|\mathcal{S}|$ non-zero entries
and $\mathcal{S}=\\{j:\mathbf{1}([\beta_{0}]_{(j)}\neq 0)\\}$. In order to
deal with sheer volume of data,
$\mathbf{Y}=\big{(}\mathbf{y},\mathbf{X}\big{)}$ is split into $s$ smaller
distinct subsets of data
$\check{\mathbf{Y}}^{(i)}=\big{(}\check{\mathbf{y}}^{(i)},\check{\mathbf{X}}^{(i)}\big{)}\in\mathbb{R}^{b\times(p+1)},i=1,\cdots,s$
that can be stored and processed separately. The subsets may be formed by
resampling without replacement from rows of the complete data set where
$b=\\{\lfloor n^{\gamma}\rfloor|\gamma\in[0.6,1)\\}$. The same situation would
occur if subsets of data are stored on $s$ storage and computing nodes and
each node contains $b$ observations.

###  II-B Assumptions

Suppose the measurement noise or errors $v_{l}$ follow some distribution
$F_{0}$ and the distribution of the observed predictors $\mathbf{x}_{[l]}$ is
$G_{0}$. Then, the joint distribution $H_{0}$ of $(y_{l},\mathbf{x}_{[l]})$ is

$H_{0}(\mathbf{x}_{[l]},y_{l})=G_{0}(\mathbf{x}_{[l]})F_{0}(y_{l}-\mathbf{x}_{[l]}^{T}\boldsymbol{\beta}_{0}).$
(2)

We make the following assumptions on the distribution of errors and
predictors.

  * •

The probability density $f_{0}(u)$ associated with probability distribution
$F_{0}$ of the errors $v_{l}$ has the following properties: even, monotone
decreasing in $|v|$ and strictly decreasing in a neighborhood of $0$.

  * •

$\mathbb{P}(\mathbf{x}_{[l]}^{T}\boldsymbol{\beta}=\mathbf{0})<1-\delta$ for
all non-zero $\boldsymbol{\beta}$ and $\delta$ as defined by equation (5).

  * •

$G_{0}$ has a finite second moment and
$\mathbb{E}_{G_{0}}[\mathbf{x}_{[l]}\mathbf{x}_{[l]}^{T}]$ is non-singular.

The condition 1 generalizes the result established in this work to extremely
heavy-tailed errors by imposing no moment conditions on the residual
distribution $F_{0}$. The condition 2 guarantees the probability that observed
values of explanatory variable are concentrated on a hyperplane does not get
too large. The condition 3 concerns the second moment of explanatory variables
and very common in the asymptotic analysis of regression estimators.

###  II-C Proposed Two-Stage Inference Methods

In order to perform inference on potentially high-dimensional models in the
presence of sparsity and outlying observations, we propose two-stage robust
and distributed statistical inference methods where robust variable selection
is performed in the first stage and then selected variables from each distinct
subset of data are combined by using a fusion rule in the fusion center or
cloud. The actual inference is done in the second stage by using robust and
low-complexity bootstrapping procedures on the variables selected in the first
stage.

The robust variable selection incorporates features of split-and-conquer
approach [26], Bolasso [27] and robust lasso estimators. The Bolasso discusses
that Lasso fails to produce consistent variable selection results under
certain decays of regularization parameter. More specifically, the probability
of selecting irrelevant variables is strictly positive. In order to overcome
this issue, Bolasso recommends to generate sufficient number of bootstrap
samples of the original data set and intersect the support of the parameter
vector estimated based on each bootstrap. Thus, the irrelevant variables would
randomly be selected by Lasso and could be eliminated from the support during
intersection. The majority voting scheme is often preferred over intersection
because if few of the relevant variables are erroneously not selected by one
bootstrap replication, the AND-rule excludes those variables from the model.
In the split-and-conquer approach, the data set is split into $s$ subsets and
distinct subsets are processed separately. In large-scale settings with
distributed storage and processing architecture, multiple distinct subsets of
i.i.d. data are stored on nodes. One can estimate the support of the parameter
vector based on distinct subsets of data stored at each node and combine them
in the fusion center by using a voting scheme. Due to page limitations, we
chose to discuss only TSRD-$\tau$ and its related equations. The TSRD-MM
method basically replaces the class of $\tau$-estimators in TSRD-$\tau$ with
MM-estimators.

###  II-D TSRD-$\tau$

We describe the core components of two-stage robust and distributed inference
procedure employing class of $\tau$-estimators. In the first stage, the robust
variable selection exploits the $\tau$-Lasso method [28], [15] to select the
sparse basis for each of the $s$ distinct data sets of $b$ observations. The
selection results from each node are fused in a cloud or fusion center by
using a percentage-based voting rule to choose the variables for the complete
large data set. The chosen basis is communicated to each distributed computing
and storage node and used in the second stage of the inference. A new robust
and scalable technique using one-step bootstrap approximation of robust
$\tau$-estimators is proposed to find parameter estimates for the selected
basis and their confidence intervals in the presence of outliers. Bootstrap
replicates are computed by using a robust low-dimensional $\tau$-estimator of
regression [29]. The bootstrap percentile method is used to estimate
confidence intervals associated with selected variables. The estimated
confidence intervals from each node are communicated to the cloud or fusion
center for the inference on complete large scale data. The confidence
intervals estimated for each subset of data are combined at the fusion center
by coordinate-wise averaging over the lower and upper bound of confidence
intervals.

In addition, we show that the distribution of the proposed RSOB-$\tau$
asymptotically converges to the same limiting distribution as the sampling
distribution of $\tau$-estimator for distinct subsets of data. Combining this
with the assumption of $\mathbb{P}(\hat{\mathcal{S}}=\mathcal{S})=1$ as
$n\rightarrow\infty$ implies that the asymptotic distribution of bootstrap
$\tau$-estimates for non-zero variables would converge to the same limiting
normal distribution of $\tau$-estimator as if the true non-zero variables were
known a priori. The robustness properties of $\tau$-Lasso and RSOB-$\tau$ are
established by using analytical methods and computer simulations. A simple
schemtaic of the proposed two-stage inference method is presented in Fig. 1.
The details of the proposed TSDR-$\tau$ method are presented in the following
sections.

Figure 1: A simple schematic of the proposed two-stage robust and distributed
inference.

##  III The Proposed Robust Inference Procedure

###  III-A Stage I: Model Selection

The proposed distributed and robust variable selection algorithm employed in
the first stage of the inference procedure is described in detail. The main
steps of the algorithm are summarized in Fig. 2.

Figure 2: The model selection is carried out in two steps: performing model
selection at each node storing a distinct subset of data and combining the
selected models from each node in the fusion center through a voting scheme.

The large-scale data is first split into $s$ distinct subsets via resampling
without replacement. Then, the robust model selection is carried out in two
steps in a distributed manner as follows:

####  III-A1 Variable Selection Using $\tau$-Lasso

At each node, the relevant variables are selected by applying robust
$\tau$-Lasso estimator [28], [15] to each distinct subset of observations
$\check{\mathbf{Y}}^{(i)}=(\check{\mathbf{y}}^{(i)},\check{\mathbf{X}}^{(i)}),i=1,...,s$.
Hence, the robust variable selection is performed by solving a set of $s$
estimation sub-problems defined as follows:

$\hat{\boldsymbol{\beta}}^{(i)}=\operatorname*{argmin}_{\boldsymbol{\beta}}\mathcal{L}(\boldsymbol{\beta})=\operatorname*{argmin}_{\boldsymbol{\beta}}\Big{(}[\hat{\sigma}_{\tau}^{(i)}]^{2}+\lambda\|\boldsymbol{\beta}\|_{\ell_{1}}\Big{)}$
(3)

where $\lambda$ controls the level of sparsity imposed by the
${\ell_{1}}$-norm penalty term. $\hat{\sigma}_{\tau}^{(i)}$ is a shorthand for
$\hat{\sigma}_{\tau}(\check{\mathbf{r}}^{(i)}(\boldsymbol{\beta}))$ an
efficient estimate of scale defined as follows:

$[\hat{\sigma}_{\tau}^{(i)}]^{2}=\frac{[\hat{\sigma}_{b}^{(i)}]^{2}}{b}\sum_{l=1}^{b}\rho_{1}\Big{(}\frac{\check{r}_{l}^{(i)}(\boldsymbol{\beta})}{\hat{\sigma}_{b}^{(i)}}\Big{)}$
(4)

where
$\check{\mathbf{r}}^{(i)}(\boldsymbol{\beta})=\check{\mathbf{y}}^{(i)}-\check{\mathbf{X}}^{(i)}\boldsymbol{\beta}$
and $\hat{\sigma}_{b}^{(i)}$ denotes a shorthand for
$\hat{\sigma}_{\text{M}}(\check{\mathbf{r}}^{(i)}(\boldsymbol{\beta}))$ an
M-scale estimate of residuals $\check{\mathbf{r}}^{(i)}(\boldsymbol{\beta})$
defined as the solution to

$\frac{1}{b}\sum_{l=1}^{b}\rho_{0}\Big{(}\frac{\check{r}_{l}^{(i)}(\boldsymbol{\beta})}{\hat{\sigma}_{b}^{(i)}}\Big{)}=\delta_{1}$
(5)

where $\delta_{1}$ is a tuning constant controlling the asymptotic breakdown
point of the estimator. $\rho_{0}(\cdot)$ and $\rho_{1}(\cdot)$ are even and
bounded functions satisfying the properties of $\rho$-function defined by
Maronna et al. [30]. Tukey’s bisquare $\rho$-function is considered as a
popular choice in robust regression and defined as
$\rho_{i}(t)=1-\big{(}1-(t/c_{i})^{2}\big{)}^{3}\mathbf{1}(|t|\leq
c_{i}),i=0,1$ where $c_{0}$ and $c_{1}$ are chosen so that the desired normal
efficiency $\zeta^{*}$ and breakdown point $\delta^{*}$ are attained for
$\lambda=0$, respectively. This can be achieved by finding $c_{0}$ and $c_{1}$
that satisfy $\mathbb{E}[\rho_{0}(t)]=\delta^{*}$ and
$\big{(}\mathbb{E}[\psi^{{}^{\prime}}(t)]\big{)}^{2}/\mathbb{E}[\psi^{2}(t)]=\zeta^{*}$
under the normality assumption of errors $t\sim\mathcal{N}(0,1)$ when
$\lambda=0$, simultaneously. $\psi(t)$, $\psi_{0}(t)$, $\psi_{1}(t)$ and $W$
are defined as follows:

$\begin{split}\psi(t)=W\psi_{0}(t)+\psi_{1}(t),\\\
\psi_{0}(t)=\partial\rho_{0}(t)/\partial
t,\quad\psi_{1}(t)=\partial\rho_{1}(t)/\partial t,\\\
W=\big{(}2\mathbb{E}[\rho_{1}(t)]-\mathbb{E}[\psi_{1}(t)t]\big{)}/\mathbb{E}[\psi_{0}(t)t].\end{split}$
(6)

Computation: In order to solve the optimization problem given in equation
(38), we employ the generalized gradient to minimize the composite objective
function consisting of a non-convex term and a non-smooth $\ell_{1}$-norm
penalty term. The generalized gradient of the objective function is defined as
$\partial_{\boldsymbol{\beta}}([\hat{\sigma}_{\tau}^{(i)}]^{2}+\lambda\|\boldsymbol{\beta}\|_{\ell_{1}})$
where $\partial_{\boldsymbol{\beta}}[\hat{\sigma}_{\tau}^{(i)}]^{2}$
associated with the smooth, non-convex, continuously differentiable term is
identical to its gradient
$\nabla_{\boldsymbol{\beta}}[\hat{\sigma}_{\tau}^{(i)}]^{2}$ and
$\partial_{\boldsymbol{\beta}}(\lambda\|\boldsymbol{\beta}\|_{\ell_{1}})$
associated with non-smooth, convex term coincides with its subdifferential
[17], [31]. It follows from the local lipschitzity of the composite objective
function, any point $\bar{\boldsymbol{\beta}}\in\mathbb{R}^{p}$ at which
$\mathbf{0}\in\partial_{\boldsymbol{\beta}}([\hat{\sigma}_{\tau}^{(i)}]^{2}+\lambda\|\boldsymbol{\beta}\|_{\ell_{1}})$
is a local minimum of the $\tau$-Lasso estimation problem. Therefore, the
generalized gradient of the objective function wrt $\boldsymbol{\beta}$ may be
leveraged to find the local minima of the given estimation problem. It can be
shown the generalized gradient of the objective function is equivalent to the
sub-gradient of the weighted least square penalized by $\ell_{1}$-norm except
that the weights $w_{l}^{(i)}(\boldsymbol{\beta})$ here depend on the unknown
$\boldsymbol{\beta}$. Hence, the original optimization problem may be
reformulated as follows:

$\begin{split}\hat{\boldsymbol{\beta}}^{(i)}&=\operatorname*{argmin}_{\boldsymbol{\beta}}\|\boldsymbol{\Omega}^{(i)}(\check{\mathbf{y}}^{(i)}-\check{\mathbf{X}}^{(i)}\boldsymbol{\beta})\|_{\ell_{2}}^{2}+\lambda^{{}^{\prime}}\|\boldsymbol{\beta}\|_{\ell_{1}},\end{split}$
(7)

where $\lambda^{{}^{\prime}}=2b\lambda/\hat{\sigma}_{b}^{(i)}$,
$\boldsymbol{\Omega}^{(i)}$ is a diagonal matrix whose entries on diagonal are
$\sqrt{w_{l}^{(i)}}$ and $w_{l}^{(i)}$ is given by,

$\begin{split}w_{l}^{(i)}=\frac{\big{[}w_{\tau}^{(i)}\psi_{0}(\tilde{r}_{l}^{(i)})+\psi_{1}(\tilde{r}_{l}^{(i)})\big{]}}{\check{r}_{l}^{(i)}},\\\
w_{\tau}^{(i)}=\frac{\sum_{l=1}^{b}\big{[}2\rho_{1}(\tilde{r}_{l}^{(i)})-\psi_{1}(\tilde{r}_{l}^{(i)})\tilde{r}_{l}^{(i)}\big{]}}{\sum_{l=1}^{b}\psi_{0}(\tilde{r}_{l}^{(i)})\tilde{r}_{l}^{(i)}}.\end{split}$
(8)

where the notation $\check{r}_{l}^{(i)}$ is a shorthand for
$\check{r}_{l}^{(i)}(\boldsymbol{\beta})$ and
$\tilde{r}_{l}^{(i)}=\check{r}_{l}^{(i)}/\hat{\sigma}_{b}^{(i)}$. In the
spirit of iteratively reweighted least-squares (IRLS) [32]-[33], we use
iteratively reweighted Lasso (IR-LASSO) alternating between finding the weight
matrices $\boldsymbol{\Omega}^{(i)}$, refining
$\hat{\boldsymbol{\beta}}^{(i)}$ and updating $\hat{\sigma}_{b}^{(i)}$. The
M-scale estimates are calculated via fixed-point iterations at each step of
IR-LASSO.

####  III-A2 Fusing variable selections

Once parameter estimation and variable selection are performed at each node,
the chosen variables are communicated to the fusion center or cloud. In the
fusion center, a percentage-based voting rule is used to select the relevant
variables for the entire large-scale data set. The selection results of all
nodes are combined according to the following rule, that is, , if a parameter
is in the support within $100\times K$ percent of subsets, it is selected to
the support for the complete data set,
$\hat{\mathcal{S}}=\\{j:\sum_{i=1}^{s}\mathbf{1}(\hat{\beta}^{(i)}_{j}\neq
0)/s\geq K\\}$. The chosen variables $\hat{\mathcal{S}}$ are broadcast into
computation and storage nodes and deployed in the second stage of inference
which uses the selected variables only.

Consider we skip the fusion and directly pass down the variables selected at
each of $s$ nodes to the next stage. It is likely that models selected at some
subsets of data are either overfitting or undefitting, which results in
inaccurate inference results, thereby an appropriate rule of fusion is
recommended.

###  III-B Stage II: Robust Inference

####  III-B1 Construction of CIs (RSOB-$\tau$)

In this part, we perform the actual inference over the chosen variables from
the model selection stage by using the RSOB-$\tau$. The derived robust
$\tau$-estimation equations are used to compute bootstrap replicates instead
of MM-estimation equations in [22]. Scalability is achieved by conducting
inference in parallel on each distinct subset of data. Moreover, applying the
RSOB-$\tau$ on only the selected variables eliminates the bias introduced by
the regularization term.

In order to develop low-complexity inference procedures using bootstrap, it is
required for the underlying estimator to be expressed as a fixed-point
problem. It is well-known that $\tau$-estimator fulfills the above requirement
and can be represented as a solution of fixed-point equations as follows:

$\boldsymbol{\hat{\theta}}_{b}=\mathbf{f}(\boldsymbol{\hat{\theta}};\tilde{\mathbf{Y}}_{b})$
(9)

where
$\mathbf{f}:\mathbb{R}^{|\hat{\mathcal{S}}|+1}\rightarrow\mathbb{R}^{|\hat{\mathcal{S}}|+1}$
denotes a smooth function,
$\tilde{\mathbf{Y}}_{b}=\big{(}\check{\mathbf{y}},\tilde{\mathbf{X}}\big{)}$
with the explaining variables chosen in the model selection stage and
$\boldsymbol{\hat{\theta}}_{b}$ is a fixed-point of $\mathbf{f}$. The terms
$\boldsymbol{\hat{\theta}}_{b}$ and
$\mathbf{f}(\boldsymbol{\hat{\theta}}_{b};\tilde{\mathbf{Y}}_{b})$ are
defined, respectively, as follows:

$\begin{split}\boldsymbol{\hat{\theta}}_{b}&=\begin{bmatrix}\hat{\boldsymbol{\beta}}_{b}\\\
\\\ \hat{\sigma}_{b}\end{bmatrix}\\\
\mathbf{f}(\boldsymbol{\hat{\theta}}_{b};\tilde{\mathbf{Y}})&=\begin{bmatrix}\Big{(}\sum_{l=1}^{b}\hat{w}_{l}\tilde{\mathbf{x}}_{[l]}\tilde{\mathbf{x}}_{[l]}^{T}\Big{)}^{-1}\Big{(}\sum_{l=1}^{b}\hat{w}_{l}\check{y}_{l}\tilde{\mathbf{x}}_{[l]}\Big{)},\\\
\\\ \sum_{l=1}^{b}\hat{v}_{l}\hat{r}_{l},\end{bmatrix}\end{split}$ (10)

where

$\begin{split}\hat{v}_{l}&=\frac{1}{b\delta_{2}}\times\frac{\rho_{0}\big{(}\tilde{r}_{l}\big{)}}{\tilde{r}_{l}},\\\
\hat{w}_{l}&=\frac{\hat{w}_{\tau}\rho_{0}^{{}^{\prime}}\big{(}\tilde{r}_{l}\big{)}+\rho_{1}^{{}^{\prime}}\big{(}\tilde{r}_{l}\big{)}}{\hat{r}_{l}},\\\
\hat{w}_{\tau}&=\frac{\sum_{l=1}^{b}\Big{[}2\rho_{1}\big{(}\tilde{r}_{l}\big{)}-\rho_{1}^{{}^{\prime}}\big{(}\tilde{r}_{l}\big{)}\tilde{r}_{l}\Big{]}}{\sum_{l=1}^{b}\rho_{0}^{{}^{\prime}}\big{(}\tilde{r}_{l}\big{)}\tilde{r}_{l}},\\\
\hat{r}_{l}&=\check{y}_{l}-\tilde{\mathbf{x}}_{[l]}^{T}\hat{\boldsymbol{\beta}}_{b},\\\
\tilde{r}_{l}&=\frac{\hat{r}_{l}}{\hat{\sigma}_{b}}.\end{split}$ (11)

$\hat{w}_{l}$ down-weights the outlying observations. In order to construct
confidence intervals of the selected variables, the bootstrap replications of
$\boldsymbol{\hat{\theta}}_{b}$ are computed as follows:

$\boldsymbol{\hat{\theta}}_{n,b}^{\star}=\mathbf{f}(\boldsymbol{\hat{\theta}}_{n,b}^{\star};\tilde{\mathbf{Y}}^{\star}),$
(12)

where
$\tilde{\mathbf{Y}}^{\star}=\big{(}\tilde{\mathbf{Y}}_{b},\boldsymbol{\omega}^{\star}\big{)}$
denotes a bootstrap sample of size $n$ randomly drawn with replacement from
the given subset of data. The multiplicity of observations is determined by
the random weight vector $\boldsymbol{\omega}^{\star}\in\mathbb{R}^{b}$ drawn
from a multinomial distribution $\big{(}n,(1/b)\mathbf{1}_{b}\big{)}$. Instead
of computing a fully iterating bootstrap replicate
$\boldsymbol{\hat{\theta}}_{n,b}^{\star}$, we can approximate it via a one-
step iteration $\boldsymbol{\hat{\theta}}_{n,b}^{1\star}$ as follows:

$\boldsymbol{\hat{\theta}}_{n,b}^{1\star}=\mathbf{f}(\boldsymbol{\hat{\theta}}_{b};\tilde{\mathbf{Y}}^{\star}),$
(13)

However, the distribution of $\boldsymbol{\hat{\theta}}_{n,b}^{1\star}$ may
not exhibit the actual variability of the sampling distribution of
$\boldsymbol{\hat{\theta}}_{b}$, mainly because all bootstrap replicates
$\boldsymbol{\hat{\theta}}_{n,b}^{1\star}$ are calculated starting from the
same initial value $\boldsymbol{\hat{\theta}}_{b}$. It has been shown that
one-step iteration of bootstrap replications for many estimators could be
adjusted by a correction term to provide asymptotically true estimate of
bootstrap distribution [34]. In order to achieve asymptotically correct
bootstrap estimates, the one-step improvement of
$\boldsymbol{\hat{\theta}}_{n,b}^{1\star}$ with a linear correction term can
be written as [21], [35]:

$\boldsymbol{\hat{\theta}}_{n,b}^{R\star}=\boldsymbol{\hat{\theta}}_{b}+\Big{[}\mathbf{I}-\nabla\mathbf{f}(\boldsymbol{\hat{\theta}}_{b};\tilde{\mathbf{Y}}_{b})\Big{]}^{-1}\Big{(}\boldsymbol{\hat{\theta}}_{n,b}^{1\star}-\boldsymbol{\hat{\theta}}_{b}\Big{)},$
(14)

where
$\boldsymbol{\hat{\theta}}_{n,b}^{R\star}\in\mathbb{R}^{|\hat{\mathcal{S}}|+1}$
denotes the linearly corrected one-step bootstrap replication of
$\boldsymbol{\hat{\theta}}_{b}$ and $\nabla\mathbf{f}(\cdot)$ is a gradient
matrix with respect to $\boldsymbol{\theta}$. The linear correction term can
be obtained by inverting the matrix
$\Big{[}\mathbf{I}-\nabla\mathbf{f}(\boldsymbol{\hat{\theta}}_{b};\tilde{\mathbf{Y}}_{b})\Big{]}$
via the block matrix inversion lemma as described in Appendix.

Computational efficiency is attained because the correction factor and the
initial estimate $\boldsymbol{\hat{\theta}}_{b}$ are computed only once for
each subset of data. Furthermore, one-step bootstrap replications of
$\boldsymbol{\hat{\theta}}_{b}$ are computationally inexpensive.We show in
theorem 2 that $\boldsymbol{\hat{\theta}}_{n,b}^{R\star}$ would estimate the
same limiting distribution as the actual bootstrap distribution
$\boldsymbol{\hat{\theta}}_{n,b}^{\star}$ under certain regularity conditions.

The algorithm begins with generating $B$ bootstrap samples for each distinct
subset of data
$\tilde{\mathbf{Y}}^{(i)}_{b}=\big{(}\check{\mathbf{y}}^{(i)},\tilde{\mathbf{X}}^{(i)}\big{)}\in\mathbb{R}^{b\times(|\hat{\mathcal{S}}|+1)},i=1,\cdots,s$.
Across all computing and storage nodes, linearly corrected one-step bootstrap
replicates $\hat{\boldsymbol{\beta}}_{n,b}^{R\star,(ij)},j=1,\cdots,B$ are
computed. Then, confidence intervals associated with selected variables from
stage 1 are constructed by using bootstrap percentile method for each subset
of data. The estimated confidence intervals are communicated from each
computing node to the fusion center for performing inference on the complete
large-scale data. In the fusion center, the confidence intervals for the
complete large-scale data are produced by applying coordinate-wise averaging
over upper bounds and lower bounds of transmitted confidence intervals as
follows:

$\begin{split}\overline{\text{CI}}_{j}=\frac{1}{s}\sum_{i=1}^{s}\overline{\text{CI}}_{j}^{\star(i)}\\\
\underline{\text{CI}}_{j}=\frac{1}{s}\sum_{i=1}^{s}\underline{\text{CI}}_{j}^{\star(i)}\end{split}$
(15)

where $\overline{\text{CI}}_{j}$ and $\underline{\text{CI}}_{j}$ denote the
lower bound and upper bound of confidence interval associated with entry $j$
of the non-zero parameter.

Herein, we assumed that the proportion of outliers withing each subset of data
remains below $50\%$, implying the estimated confidence intervals are not
corrupted. Therefore, we may avoid the undesirable effects of outliers on the
veracity of confidence intervals. In case more than half of observations
within each subset of data are contaminated by outliers, we recommend using
the adaptive trimmed mean [36] or classical trimmed mean with manually set
trimming ratio [30].

Data: $s$ distinct subsamples

Output: $\mathbf{CI}$

1 for  _each subsample_ do

2 Generate $B$ bootstrap samples of size $n$ by randomly drawing with
replacement as follows:
$\tilde{\mathbf{Y}}^{\star}=\big{(}\tilde{\mathbf{Y}}_{b},\boldsymbol{\omega}^{\star}\big{)}$

3 Find the initial estimate $\boldsymbol{\hat{\theta}}_{b}$ by solving the
fixed-point problem given in equation (9)

4 Calculate the one-step linearly corrected bootstrap replication
corresponding to each bootstrap sample
$\hat{\boldsymbol{\beta}}_{n,b}^{R\star}$ from equation (33)

5 Compute confidence intervals $\mathbf{CI}^{\star(i)}$ for each subset of
data via bootstrap percentile method

6

7 end for

Combine the confidence intervals by coordinate-wise averaging

Algorithm 1 The RSOB-$\tau$ procedure

##  IV Statistical Robustness and Consistency Properties

In this section, the robustness properties of $\tau$-Lasso estimators and
linearly corrected one-step replications using RSOB-$\tau$ are characterized
by deriving their breakdown points. Furthermore, asymptotic normality of the
linearly corrected one-step bootstrap replications is established under
certain regularity conditions.

###  IV-A Robustness Properties of $\tau$-Lasso Estimators

The finite-sample breakdown point of $\tau$-Lasso is derived. It measures the
largest proportion of observations when arbitrarily replaced by outliers does
not cause unbounded maximum bias or equivalently does not break down. Given a
subset of data
$\check{\mathbf{Y}}=\big{(}\check{\mathbf{y}},\check{\mathbf{X}}\big{)}\in\mathbb{R}^{b\times(p+1)}$
randomly drawn from the original data
$\mathbf{Y}=\big{(}\mathbf{y},\mathbf{X}\big{)}\in\mathbb{R}^{n\times(p+1)}$,
the replacement finite-sample breakdown point (FBP)
$\epsilon^{*}\big{(}\hat{\boldsymbol{\beta}};\check{\mathbf{Y}}\big{)}$ of a
regression estimator $\hat{\boldsymbol{\beta}}\in\mathbb{R}^{p}$ is defined as

$\epsilon^{*}\big{(}\hat{\boldsymbol{\beta}};\check{\mathbf{Y}}\big{)}=\max\\{\frac{m}{b}:\sup_{\check{\mathbf{Y}}_{m}\in\mathcal{Y}_{m}}\|\hat{\boldsymbol{\beta}}\big{(}\check{\mathbf{Y}}_{m}\big{)}\|_{\ell_{2}}<\infty\\}$
(16)

where the set $\mathcal{Y}_{m}$ contains all datasets $\check{\mathbf{Y}}_{m}$
with $m$ $(0<m<b)$ out of the original $b$ observations replaced by arbitrary
values. The bounded supremum of $\ell_{2}$ term in the definition is
equivalent to having a bounded maximum bias. The following theorem proves the
results for finite-sample breakdown point of robust $\tau$-Lasso estimator by
extending the theoretical result for that of robust S-PENSE estimator
demonstrated in [17] to the class of $\tau$-Lasso estimators.

Theorem 1: Suppose $m(\delta)$ is the largest integer smaller than
$b\min(\delta,1-\delta)$ for a subset of data
$\check{\mathbf{Y}}=\big{(}\check{\mathbf{y}},\check{\mathbf{X}}\big{)}\in\mathbb{R}^{b\times(p+1)}$
and $\delta$ defined by equation (5). Then the finite-sample breakdown point
of the $\tau$-Lasso estimator is bounded from above and below as follows:

$\frac{m(\delta)}{b}\leq\epsilon^{*}(\hat{\boldsymbol{\beta}};\check{\mathbf{Y}})\leq\delta$
(17)

where $\hat{\boldsymbol{\beta}}$ denotes the $\tau$-Lasso estimator.

Proof: The complete proof is given in Supplemental Materials.

With minor modification of the proof given in the Supplemental Material, we
can extend the above theorem for data models including an intercept term.

###  IV-B Asymptotic Properties of TSRD-$\tau$

In this section, we show that the proposed inference procedure will enjoy
desirable asymptotic properties for fixed $p$ under certain regularity
conditions if the model selection procedure is consistent. The model selection
consistency in the first stage would imply that the relevant variables are
selected with probability converging to $1$ and the asymptotic distribution of
modified bootstrap $\tau$-estimates for non-zero coefficients would be the
same asymptotic normal distribution as if the true non-zero coefficients were
known in advance. It seems reasonable to conjecture that the model selection
procedure is consistent, extending the concept of Bolasso [27] to robust
$\ell_{1}$-penalized estimators. We explore this possibility in simulations.
The asymptotic properties of FRB for $\tau$-regression estimator were proven
in Theorem 1 of [35]. However, the inaccuracy in the linear correction factor
in the proof is rectified here. The asymptotic properties of BLFRB and FRB for
MM-regression estimator are established in [22] and [21]. Also note that in
this theorem, $\boldsymbol{\beta}_{0}\in\mathbb{R}^{|\hat{\mathcal{S}}|}$ is a
shorthand for $[\boldsymbol{\beta}_{0}]_{\hat{\mathcal{S}}}$ where parameter
vector $\boldsymbol{\beta}_{0}\in\mathbb{R}^{p}$ as given in equation (1) and
$\xrightarrow{P}$ denotes convergence in probability.

Theorem 2 : Let $\rho_{0}$ and $\rho_{1}$ be bounded $\rho$-functions
satisfying the properties of bounded $\rho$-function defined by Maronna et al.
[30] and have continuous third derivatives. Assume the model selection stage
produces consistent estimates, i.e.,
$\mathbb{P}(\hat{\mathcal{S}}=\mathcal{S})=o(1)$. Let
$\hat{\boldsymbol{\beta}}_{b}$ be the $\tau$-regression estimator for a subset
of data
$\big{(}\check{\mathbf{y}},\tilde{\mathbf{X}}\big{)}\in\mathbb{R}^{b\times(|\hat{\mathcal{S}}|+1)}$
randomly drawn from the original data
$\big{(}\mathbf{y},\underaccent{\bar}{\bX}\big{)}\in\mathbb{R}^{n\times(|\hat{\mathcal{S}}|+1)}$
and $\hat{\sigma}_{b}$ be the M-scale of residuals for the given subset of
data and assume that they are consistent estimators, that is,
$\hat{\boldsymbol{\beta}}_{b}\xrightarrow{P}\boldsymbol{\beta}_{0}$ and
$\hat{\sigma}_{b}\xrightarrow{P}\sigma_{0}$. Given the following regularity
conditions hold:

  1. 1.

The following vectors and matrices exist and are finite:

    1. 1.1

$\mathbb{E}\big{[}\big{(}\bar{w}_{\tau}\rho_{0}^{{}^{\prime}}(r)+\rho_{1}^{{}^{\prime}}(r)\big{)}/r\mathbf{x}\mathbf{x}^{T}\big{]}^{-1}$

    2. 1.2

$\mathbb{E}\big{[}\rho_{0}^{{}^{\prime}}(r)\mathbf{x}\big{]}$

    3. 1.3

$\mathbb{E}\big{[}\big{(}\bar{w}_{\tau}\rho_{0}^{{}^{\prime\prime}}(r)+\rho_{1}^{{}^{\prime\prime}}(r)\big{)}\mathbf{x}\mathbf{x}^{T}\big{]}$

    4. 1.4

$\mathbb{E}\big{[}\big{(}\bar{w}_{\tau}\rho_{0}^{{}^{\prime\prime}}(r)+\rho_{1}^{{}^{\prime\prime}}(r)\big{)}r\mathbf{x}\big{]}$

  2. 2.

$\mathbb{E}[\rho_{0}^{{}^{\prime}}(r)r]\neq 0$ and finite,

  3. 3.

$\rho_{0}^{{}^{\prime}}(u)/u$, $\rho_{1}^{{}^{\prime}}(u)/u$,
$\big{(}\rho_{0}^{{}^{\prime}}(u)-\rho_{0}^{{}^{\prime\prime}}(u)u\big{)}/u^{2}$
and
$\big{(}\rho_{1}^{{}^{\prime}}(u)-\rho_{1}^{{}^{\prime\prime}}(u)u\big{)}/u^{2}$
are continuous.

where
$\bar{w}_{\tau}=\big{(}2\mathbb{E}[\rho_{1}(r)]-\mathbb{E}[\rho_{1}^{{}^{\prime}}(r)r]\big{)}/\mathbb{E}[\rho_{0}^{{}^{\prime}}(r)r]$.
Then, the distribution of
$\sqrt{n}\big{(}\hat{\boldsymbol{\beta}}_{n,b}^{R\star}-\hat{\boldsymbol{\beta}}_{b}\big{)}$
converges weakly to the limiting distribution of
$\sqrt{b}\big{(}\hat{\boldsymbol{\beta}}_{b}-\boldsymbol{\beta}_{0}\big{)}$
and consequently to the limiting distribution of
$\sqrt{n}\big{(}\hat{\boldsymbol{\beta}}_{n}-\boldsymbol{\beta}_{0}\big{)}$ as
$n$ and $b$ approach infinity.

Proof: The complete proof is presented in Supplemental Materials.

###  IV-C Statistical Robustness of RSOB-$\tau$

In this section, we are interested in studying the robustness properties of
the proposed RSOB-$\tau$. The confidence intervals for regression parameters
were constructed by using RSOB-$\tau$ quantiles and thus, the breakdown point
of quantiles gives an insight into the robustness and reliability of the
inferences made using the proposed method. Here, we derive some theoretical
results about the breakdown point of quantile estimates of RSOB-$\tau$ and FRB
using $\tau$-estimators.

Before proceeding, we define few concepts essential in comprehending the
robustness results. Given $t\in\big{(}0,1\big{)}$, $\hat{q}_{t}$ is defined as
the $t_{th}$ upper quantile of a statistic $\hat{\boldsymbol{\beta}}_{b}$ such
that $P[\hat{\boldsymbol{\beta}}_{b}>\hat{q}_{t}]=t$. According to [37], the
upper breakdown point of a bootstrap estimate $\hat{q}_{t}^{\star}$ is defined
as the minimum proportion of arbitrarily large outliers in the subset of data
$\tilde{\mathbf{Y}}$ that can drive $\hat{q}_{t}^{\star}$ into infinity. In
what follows theoretical results on robustness properties of estimated
RSOB-$\tau$ and FRB quantiles using robust $\tau$-estimator are demonstrated.
The following theorems are proved by using similar techniques as in the proofs
given for BLFRB and FRB quantiles using robust MM-estimator [22] and [21].
Note that the theorems below also hold for subsets of data
$\tilde{\mathbf{Y}}\in\mathbb{R}^{b\times|\hat{\mathcal{S}}|}$ having
$|\hat{\mathcal{S}}|$ predictors selected from model selection stage.

Theorem 3: Suppose
$\mathbf{Y}=\big{(}\mathbf{y},\mathbf{X}\big{)}\in\mathbb{R}^{n\times(p+1)}$
is a large-scale data following the linear model given in equation (1).
Except, here no sparsity assumption is imposed on the parameter vector
$\boldsymbol{\beta}_{0}$. Assume that $\mathbf{Y}$ is in general position,
i.e., any subset of $p$ observations will result in a unique determination of
$\boldsymbol{\beta}_{0}$. Let $\hat{\boldsymbol{\beta}}_{n}$ be a robust
$\tau$-estimate of $\boldsymbol{\beta}_{0}$ based on the data $\mathbf{Y}$
whose breakdown point is $\epsilon$. Then, the breakdown point of the $t_{th}$
FRB quantile estimate of the regression parameters $[\beta_{0}]_{(l)}\mbox{,
}l=1,\cdots,p$ is determined by $\min(\epsilon_{t}^{\star},\epsilon)$ where
$\epsilon_{t}^{\star}$ is the smallest $\delta_{c}\in[0,1]$ that satisfies the
following inequality,

$P\big{[}Binomial(n,1-\delta_{c})<p\big{]}\geq t$ (18)

Note that $\delta_{c}$ is different from $\delta_{0}$ and $\delta_{1}$.

Proof: The details of the poof are given in Appendix.

Theorem 4: Suppose
$\check{\mathbf{Y}}=\big{(}\check{\mathbf{y}},\check{\mathbf{X}}\big{)}\in\mathbb{R}^{b\times(p+1)}$
is a subset of the large-scale data
$\mathbf{Y}=\big{(}\mathbf{y},\mathbf{X}\big{)}\in\mathbb{R}^{n\times(p+1)}$
formed by random resampling without replacement of the original full data set.
Except, here no sparsity assumption is imposed on the parameter vector
$\boldsymbol{\beta}_{0}$. Assume that $\check{\mathbf{Y}}$ is in general
position, i.e., any subset of $p$ observations will determine a bounded least-
square estimate. Let $\hat{\boldsymbol{\beta}}_{b}$ be a robust
$\tau$-estimator of $\boldsymbol{\beta}_{0}$ based on the bag of data
$\check{\mathbf{Y}}$,
$\epsilon^{*}(\hat{\boldsymbol{\beta}}_{b},\check{\mathbf{Y}})$ denotes the
finite-sample breakdown point of $\hat{\boldsymbol{\beta}}_{b}$ and
$\delta_{1}=0.5$. Then, all the RSOB-T bootstrap quantiles estimated using the
$\tau$-estimator that are constructed over $\check{\mathbf{Y}}$ will have
equal asymptotic breakdown point to $\hat{\boldsymbol{\beta}}_{b}$.

Proof: A detailed proof is given in Appendix.

##  V SIMULATIONS AND RESULTS

In this section, the performance of the proposed method is investigated in
simulations considering both variable selection and inference. In particular,
correct identification of sparse basis, statistical robustness of bootstrap
estimates, the quality of the parameter estimates and confidence intervals are
studied. The results in Theorem 2 is validated through computer simulations,
indicating the distribution of bootstrap replications by using RSOB-$\tau$
asymptotically converges to the sampling distribution of $\tau$-estimator. The
performance of the proposed method is assessed through computer simulations by
using different proportions of outliers, large-scale data in both low-
dimensional $(p<b\text{ or }n)$ versus high-dimensional $(p\approx n\text{ or
}\gg b)$ settings. It is assumed the large-scale data follow a linear
regression model where the parameter vector
$\boldsymbol{\beta}_{0}\in\mathbb{R}^{p}$ is sparse with $k_{s}$ non-zero
entries. The measurement noise vector, $\mathbf{v}$, is an additive white
Gaussian noise with a variance (AWGN)
$\sigma_{v}^{2}=\|\mathbf{X}\boldsymbol{\beta}_{0}\|_{\ell_{2}}^{2}10^{-\mbox{SNR}/10}/n$
(SNR in dB).

###  V-A Simulation Setting

Throughout simulations, the confidence intervals are reported in a nominal
level of $100\times(1-\alpha_{cl})\%=90\%$ with $\alpha_{cl}=0.1$. We create a
decreasing grid of $70$ lambda values with logarithmic spacing of $1.1$,
spanning $\lambda_{1}$ to $\lambda_{70}$ where $\lambda_{1}$ is set to
$\lambda_{\max}$. The maximum number of iteration in $\tau$ and MM-Lasso
estimators is fixed at $30$. The robustness of DPD-SIS procedure is adjusted
by a tuning parameter $\alpha$ set to $0.4$ (stable for a range of
contamination levels) to provide robustness without significant loss in
efficiency. The proposed algorithm was implemented in MATLAB except for the
estimation of initial S-Lasso which was done in R using PENSE [38]. It
provides a good initial estimate by constructing clean subsamples of data,
potentially removing outlying observations. This is achieved by using the
principal sensitivity components (PSCs) for EN estimator and removing the
observations with most extreme PSCs from the subsamples in an iterative
manner. A detailed description of initialization for PENSE can be found in
Supplementary Materials of [39]. In addition, we used the Dual Augmented
Lagrangian [40] implementation to solve the IR-LASSO, and MM- and S-estimators
as in [41].

The outliers in all simulations except for Section V-I are introduced by
randomly choosing the observations in y and replacing them with random values
chosen from a standard Gaussian distribution with $\sigma_{e}=250$.

###  V-B Data Standardization

Across all simulations for estimation problems using robust Lasso, it is
assumed the linear regression model has an intercept component and all columns
of the augmented regression matrix $[\mathbf{1}_{b\times
1}\check{\mathbf{X}}^{(i)}]$ except the first one are robustly standardized by
centring the columns using a bisquare location estimator and scaling them
using bisquare scale estimators [42]. The response vector
$\check{\mathbf{y}}^{(i)}$ is centred using the bisquare location estimator
[42].

###  V-C Calibration of Tuning Constants $c_{0}$ and $c_{1}$

In order to tune $c_{0}$ for the M-scale of residuals within all the robust
estimators discussed in this paper, we set $\delta_{i}=0.5$ for $i=0,1,2$ to
achieve the maximum robustness against outliers. $\delta_{i}$ controls the
breakdown point according to Theorem 4.1 in [17] and $\delta_{0}$ is
associated with initial S-Lasso estimator. It’s worth mentioning $0.5$ is the
largest value $\delta_{i}$ can take on. Note that $c_{0}$ is tuned for the
desired breakdown point with the assumption $\lambda=0$. In particular, it is
recommended to have maximum robustness because bootstrapping may exacerbate
the malicious behaviour of outliers in resampled data sets. The tuning
constant $c_{1}$ of $\tau$-Lasso and MM-Lasso is respectively adjusted to
$6.08$ and $4.68$ to attain $95\%$ efficiency under Gaussian errors when
$\lambda=0$. Likewise, $c_{1}=6.08$ for $\tau$-estimator and $c_{1}=4.68$ for
MM-estimator are adjusted to provide $95\%$ efficiency under Gaussian errors
for bootstrap replications.

###  V-D Choice of Fusion Parameter $K$

$K$ denotes the proportion of data subsets classified the given variable as
relevant. Herein, the fusion parameter $K$ for voting rule is set to $0.5$,
indicating the variables selected within at least $50\%$ of subsets are
regarded as relevant variables, i.e, majority voting scheme. The specific
value of $0.5$ is a compromise between reducing the false positive rate and
maximizing the selection of true relevant variables.

###  V-E Choice of Batch Size $b$

The choice of batch size $b$ is a design parameter, depending on computational
and storage resources available plus either model selection performance or
prediction performance. On one hand, randomly partitioning the entire large-
scale data of size $n$ into subsets of size $b$ and independently analyzing
smaller subsets of data leads to significant computational speed-up and data
storage benefits [20, 26]. On the other hand, the choice of batch size $b$
trades off between overfitting and underfitting in the variable selection
procedure in additional to statistical correctness of confidence intervals.
Setting the batch size $b$ to very large values results in rejecting relevant
variables. This limits the invertibility of sub-design matrices as the
eigenvalues are no longer strictly positive. However, one shall tune $b$ to be
large enough to discard the irrelevant variables. Moreover, very small batch
size $b$ equals to large number of data subsets results in wasted computation
as we are occupying unnecessarily many computation and storage nodes while
only a small proportion of them would be needed. The known information-
theoretic results [43] impose a condition on the batch size $b$ to ensure the
support recovery for any method, $b\geq c_{s}k_{s}\log(p)$ for a sufficiently
large constant $c_{s}$. On the other hand, BLB [20] requires the batch size to
some values $n^{0.6}<b<n^{0.9}$ such that the statistical correctness is
obtained. Combining these two constraints, one may choose the batch size in
the range of $\max(n^{0.6},c_{s}k_{s}\log(p))\leq b\leq n^{0.9}$. We then need
to take into account practical implications of computational constraints and
choose an appropriate batch size $b$.Given the batch size $b$, one can obtain
the number of subsamples $s$ and $\gamma$ in a straightforward manner.

We provide an example describing an appropriate batch size $b$. Let
$n=1000000$, $b=500$ and $k_{s}=50$, we then obtain that the minimum batch
size from support recovery point of view is $b\geq 311c_{s}$. Setting the
constant $c_{s}=30$, sufficiently large for sparse recovery, implies that
$b\geq 9400$. On the other hand, reliable statistical correctness is attained
by BLB for $b\geq(n^{0.6}\approx 4000)$. We shall occupy many resources over
$100$ computing cores. If we set $b=n^{0.75}\approx 32000$, we shall need $33$
computing cores and computations over $32000$ is manageable by single
computing cores, which gives a good compromise between computational burden
and statistical efficiency. Setting $b\approx n^{0.9}$ would mean almost $4$
subsets of size $250000$ which require higher computational power for each
node and not very desirable.

###  V-F Tuning of regularization parameter

In order to calibrate the robust penalized estimation problem,
$\lambda_{\max}$ is initially estimated by using the method introduced by Khan
et al. [44] and then improved upon via a binary search [16]. A set of
candidate lambdas in decreasing order starting from $\lambda_{max}$ is formed.
Selection of $\lambda$ is carried out through a robust version of Bayesian
Information Criterion (BIC) [45], defined as

$\mbox{RBIC}(\lambda)=b\log(\hat{\sigma}_{b}^{2}(\lambda))+C(b,p)\|\hat{\boldsymbol{\beta}}(\lambda)\|_{\ell_{0}}$
(19)

where $\hat{\sigma}_{b}$ denotes the robust M-scale of residuals and $C(b,p)$
is set to $\log(b)$ for settings where the dimensionality $p$ is much smaller
than the sample size $n$. $C(b,p)=\log(\log(b))\log(p)$ is chosen for
$\log(p)/b\rightarrow 0$ as $p\rightarrow\infty$ to attain better empirical
performance [46, 47, 48]. Here, standard BIC is modified by replacing the non-
robust estimate of scale with the robust M-estimate of scale to deal with
outliers. Finally, the optimal tuning parameter $\lambda$ minimizes the RBIC
over the pre-defined grid of lambdas

$\lambda^{*}=\operatorname*{argmin}_{\lambda\in\Lambda}\mbox{RBIC}(\lambda).$
(20)

###  V-G Scenarios

We consider the following scenarios for which simulation studies are carried
out.

  * •

Scenario 1: We set the simulation parameters as follows: $n=27000$, $p=30000$,
$b=900$ ($p/b=33.33$), $\gamma=0.6667$, $\mbox{SNR}=15\mbox{ dB}$, $B=1000$.
$\boldsymbol{\beta}_{0}\in\mathbb{R}^{p}$ is sparse with $k_{s}=40$ non-zero
entries. $[\boldsymbol{\beta}_{0}]_{\mathcal{S}}$ is set to
$3\times\mathbf{1}_{\mathcal{S}}$ and their positions are chosen randomly. The
covariate vectors $\mathbf{x}_{[i]},i=1,\cdots,n$ are drawn independently from
a multivariate Gaussian distribution
$\mathcal{N}\big{(}\mathbf{0},\boldsymbol{\Sigma}\big{)}$ with
$\Sigma_{ij}=\rho^{|i-j|}$ (Toeplitz covariance structure, $\rho=0.5$).

  * •

Scenario 2: The simulation parameters are set as follows: $n=2000000$, $p=80$,
$b=40000$ ($p/b=1/500$), $\gamma=0.730367$, $\mbox{SNR}=30\mbox{ dB}$,
$B=400$. $\boldsymbol{\beta}_{0}\in\mathbb{R}^{p}$ is sparse with $k_{s}=20$
non-zero entries. $[\boldsymbol{\beta}_{0}]_{\mathcal{S}}$ is set to
$3\times\mathbf{1}_{\mathcal{S}}$ and their positions are chosen randomly.
Explaining variables $\mathbf{x}_{[i]}$ in the regression matrix are i.i.d,
randomly drawn from a multivariate Gaussian distribution
$\mathcal{N}\big{(}\mathbf{0},\mathbf{I}_{p}\big{)}$.

  * •

Scenario 3: The simulation parameters are set as follows: $n=80000$, $p=100$,
$b=4000$ ($p/b=1/40$), $\gamma=0.73466$, $\mbox{SNR}=15\mbox{ dB}$,
$B\\!=\\!300$. The regression matrix $\mathbf{X}$ is randomly generated from
mutually independent observations drawn from
$\mathcal{N}(\mathbf{0},\boldsymbol{\Sigma})$ with $\Sigma_{ij}=\rho^{|i-j|}$
(Toeplitz covariance structure, $\rho\\!=\\!0.5$).
$\boldsymbol{\beta}_{0}\in\mathbb{R}^{p}$ is sparse with $k_{s}=10$ non-zero
entries. $[\boldsymbol{\beta}_{0}]_{\mathcal{S}}$ is set to
$3\times\mathbf{1}_{\mathcal{S}}$ and their positions are chosen randomly.

  * •

Scenario 4: We set the simulation parameters as follows: $n=4900$, $p=6000$,
$b=350$ ($p/b=17.14$), $\gamma=0.6895$, $\mbox{SNR}=15\mbox{ dB}$, $B=1000$.
$\boldsymbol{\beta}_{0}\in\mathbb{R}^{p}$ is sparse with $k_{s}=10$ non-zero
entries. $\boldsymbol{\beta}_{0}$ is set to

$\boldsymbol{\beta}_{0}=[2.5,2.5,2.5,2,3,3,3,3.5,3.5,3.5,\mathbf{0}_{p-k_{s}}^{T}]^{T}$

The covariate vectors $\mathbf{x}_{[i]},i=1,\cdots,n$ are drawn independently
from a multivariate Gaussian distribution
$\mathcal{N}\big{(}\mathbf{0},\boldsymbol{\Sigma}\big{)}$ with
$\Sigma_{ij}=\rho^{|i-j|}$ (Toeplitz covariance structure, $\rho=0.5$).

  * •

Scenario 5: We set the simulation parameters as follows: $n=20000$, $p=80$,
$\mbox{SNR}=15\mbox{ dB}$. $\boldsymbol{\beta}_{0}\in\mathbb{R}^{p}$ is sparse
with $k_{s}=15$ non-zero entries. $\boldsymbol{\beta}_{0}$ is set to

$\boldsymbol{\beta}_{0}=[3.5,3.5,3.5,5,5,5,2.5,2.5,2.5,1.5,2\times\mathbf{1}_{5}^{T},\mathbf{0}_{p-k_{s}}^{T}]^{T}$

The covariate vectors $\mathbf{x}_{[i]},i=1,\cdots,n$ are drawn independently
from a multivariate Gaussian distribution
$\mathcal{N}\big{(}\mathbf{0},\boldsymbol{\Sigma}\big{)}$ with
$\Sigma_{ij}=\rho^{|i-j|}$ (Toeplitz covariance structure, $\rho=0.5$).

###  V-H Variable Selection Performance with Different $\ell_{1}$-penalized
Estimators

In this subsection, we first substitute different $\ell_{1}$-penalized
estimators with $\tau$-Lasso estimators used in the model selection stage of
TSRD-$\tau$ introduced in Section III. We then carry out simulations to
compare the variable selection performance of the resulting procedures with
that of TSRD-$\tau$ and TSRD-MM. For this purpose, we use $\ell_{1}$-penalized
estimators such as RA-Lasso [19] and Sparse-LTS [18]. The former is
implemented in MATLAB by following exactly the computation algorithm described
in [19] whereas we use the well-known R package robustHD for Sparse LTS
regression [49]. We run the simulations on the synthetic data set described by
Scenario 5 except for setting $\text{SNR}=10\text{ dB}$ with normal errors. In
this experiment, we introduce contamination in the regression matrix
$\mathbf{X}$ and the response vector $\mathbf{y}$, simultaneously. For each
simulation, we consider four contamination schemes as follows:

  * •

Scheme 1: $10\%$ of observations in the response vector $\mathbf{y}$ are
replaced with random values drawn from standard Gaussian with
$\sigma_{e}=250$. We also replace the corresponding observations in
$\mathbf{X}$ with random values chosen from standard multivariate Gaussian
distribution with $\mathbf{\Sigma}=\sigma_{e}^{2}\mathbf{I}_{p}$. (large-
variance outliers)

  * •

Scheme 2: $10\%$ of observations in the response vector $\mathbf{y}$ are
replaced with random values drawn from a Gaussian distribution
$\mathcal{N}(250,1)$. We also replace the corresponding observations in
$\mathbf{X}$ with random values chosen from a multivariate Gaussian
distribution $\mathbf{N}(50\times\mathbf{1}_{p},\mathbf{I}_{p})$. (gross
outliers)

  * •

Scheme 3: The additive noise is heavy-tailed Student’s $t$-distribution with
one degree of freedom.

  * •

Scheme 4: The outlier contamination follows the same procedure as in Scheme 1
with the assumption of heavy-tailed additive errors, thereby combining Scheme
1 and 3.

Across all schemes, the same set of observations are replaced with outliers.
Herein, we perform a Monte-Carlo study of 20 trials where a random realization
of outlier in $\mathbf{y}$ is used at each trial. Table I shows the result of
variable selection method with different $\ell_{1}$-penalized estimators for
four contamination schemes, averaged over 20 trials. We observe that the
proposed variable selection algorithms using $\tau$-Lasso estimator and MM-
Lasso estimator perfectly recover the true relevant variables while keeping
the false positive rates low. In contrast, variable selection method using
Sparse-LTS results in significantly larger false positive rates for smaller
subsample sizes under contamination schemes 1 and 2. Except for contamination
scheme 2 (heavy-tailed errors), the variable selection with RA-Lasso results
in extremely overfitted models.

TABLE I:  Comparison of Variable Selection Performance with Different $\ell_{1}$\- penalized Estimators: The Proposed Algorithms using MM-Lasso and $\tau$-Lasso Exhibit a Reliable Performance in Recovering The True Sparse Basis (TPR=1) while Keeping False Positive Rates low. Variable Selection Method using Sparse-LTS fails to suppress the False Positives for Schemes 1 and 2 (Smaller subsample Size). Variable Selection with RA-Lasso Performs Reliably Only for Scheme 3, Heavy-Tailed Errors with No Outlier Contamination | scheme 1 |  | scheme 2 |  | scheme 3 |  | scheme 4
---|---|---|---|---|---|---|---
| TP | FP | CER |  | TP | FP | CER |  | TP | FP | CER |  | TP | FP | CER
$b=625$ |  |  |  |  |  |  |  |  |  |  |  |
TSRD-$\tau$ | 1 | 0.0085 | 0.0069 |  | 1 | 0 | 0 |  | 1 | 0 | 0 |  | 1 | 0 | 0
TSRD-MM | 1 | 0 | 0 |  | 1 | 0 | 0 |  | 1 | 0 | 0 |  | 1 | 0 | 0
Sparse-LTS | 1 | 0.14 | 0.1138 |  | 1 | 0.2008 | 0.1631 |  | 1 | 0.0008 | 0.0006 |  | 1 | 0.0023 | 0.0019
RA-Lasso | 1 | 1 | 0.8125 |  | 1 | 0.1461 | 0.1187 |  | 1 | 0 | 0 |  | 1 | 1 | 0.8125
$b=800$ |  |  |  |  |  |  |  |  |  |  |  |
TSRD-$\tau$ | 1 | 0 | 0 |  | 1 | 0 | 0 |  | 1 | 0 | 0 |  | 1 | 0 | 0
TSRD-MM | 1 | 0 | 0 |  | 1 | 0 | 0 |  | 1 | 0 | 0 |  | 1 | 0 | 0
Sparse-LTS | 1 | 0.0146 | 0.0119 |  | 1 | 0.05 | 0.0406 |  | 1 | 0 | 0 |  | 1 | 0 | 0
RA-Lasso | 1 | 1 | 0.8125 |  | 1 | 0.2477 | 0.2019 |  | 1 | 0 | 0 |  | 0.9966 | 0.9938 | 0.8081
$b=1000$ |  |  |  |  |  |  |  |  |  |  |  |
TSRD-$\tau$ | 1 | 0.0015 | 0.0013 |  | 1 | 0 | 0 |  | 1 | 0 | 0 |  | 1 | 0 | 0
TSRD-MM | 1 | 0 | 0 |  | 1 | 0 | 0 |  | 1 | 0 | 0 |  | 1 | 0 | 0
Sparse-LTS | 1 | 0.0023 | 0.0019 |  | 1 | 0.0077 | 0.0062 |  | 1 | 0 | 0 |  | 1 | 0 | 0
RA-Lasso | 0.92 | 0.9231 | 0.7650 |  | 1 | 0.4661 | 0.3787 |  | 1 | 0 | 0 |  | 0.3733 | 0.3569 | 0.4075

###  V-I Robustness of Bootstrap Replications

In this subsection, the robustness of RSOB-$\tau$ is quantified by an
uncertainty measure in comparison to BLB. In particular, we calculate the
standard deviation based on the bootstrap replications produced by both
methods and verify the results in Theorem 4 by assessing the relative error.
The bootstrap estimate of standard deviation is computed as follows:

$\widehat{\mbox{SD}}(\hat{\boldsymbol{\beta}}_{n})=\frac{1}{|\hat{\mathcal{S}}|}\sum_{l=1}^{|\hat{\mathcal{S}}|}\Bigg{(}\frac{1}{s}\sum_{i=1}^{s}\bigg{(}\sum_{j=1}^{B}\frac{\Big{(}[\hat{\boldsymbol{\beta}}_{n,b}^{\star(ij)}]_{l}-[\hat{\boldsymbol{\beta}}_{n,b}^{\star(i.)}]_{l}\Big{)}^{2}}{B-1}\bigg{)}^{1/2}\Bigg{)}$
(21)

where $[\hat{\boldsymbol{\beta}}_{n,b}^{\star(i.)}]_{l}$ is given by:

$[\hat{\boldsymbol{\beta}}_{n,b}^{\star(i.)}]_{l}=\frac{1}{B}\sum_{j=1}^{B}[\hat{\boldsymbol{\beta}}_{n,b}^{\star(ij)}]_{l}$
(22)

In order to measure how accurate the bootstrap estimate of standard deviation
approximates the average standard deviation of $\hat{\boldsymbol{\beta}}_{n}$,
we use the relative error criterion that is defined as follows:

$\varepsilon=\frac{\widehat{\mbox{SD}}(\hat{\boldsymbol{\beta}}_{n})-\overline{\mbox{SD}}(\hat{\boldsymbol{\beta}}_{n})}{\overline{\mbox{SD}}(\hat{\boldsymbol{\beta}}_{n})}$
(23)

where the average standard deviation
$\overline{\mbox{SD}}(\hat{\boldsymbol{\beta}}_{n})$ is defined as
$\sigma/\sqrt{n\mathcal{O}}$ based on asymptotic covariance of
$\tau$-estimator [29]. Here, $\mathcal{O}$ is set to $1$ for least square
estimator in BLB and $\mathcal{O}=0.95$ for $\tau$-estimator in RSOB-$\tau$
tuned to have $95\%$ Gaussian efficiency.

First, we show that even one outlying observation could drive
$\widehat{\mbox{SD}}(\hat{\boldsymbol{\beta}}_{n})$ based on bootstrap
estimates obtained by BLB into infinity whereas those based on bootstrap
replications obtained by the RSOB-$\tau$ remain resistant to outlier. Herein,
we run the simulations on the synthetic data set described by Scenario 2. In
this experiment, a data point within the original data set is randomly drawn
and its response vector is multiplied by an extreme value $\alpha_{o}$ powers
of $10$ to imitate the situation where outlier is introduced by misplacement
of decimal point. In regard to model selection, the robust and non-robust two-
stage methods exhibit reliable performance in selecting the true sparse basis
of the parameter vector with a $\text{TP}=1$ and $\text{CER}=0$. One might
have expected that the non-robust inference would fail in model selection.
However, only the bag of data containing the outlier yielded unreliable
estimates but the voting scheme in the fusion center reduces the adversarial
effect of outlier.

In the stage 2, the bootstrap estimates of
$\widehat{\mbox{SD}}(\hat{\boldsymbol{\beta}}_{n})$ are computed based on the
data set generated by the selected predictors from the model selection stage.
As it is observed in Fig. 3, both two-stage algorithms TSRD-$\tau$ and TSLL
perform remarkably well in terms of relative error when there are no outliers
present within the data set. However, the bootstrap estimates of standard
deviation obtained by BLB are severely influenced by the presence of even one
outlier. As the magnitude of $\alpha_{o}$ increases, the relative error gets
larger, implying that BLB is not robust to outliers. On contrary, the relative
error of standard deviation obtained by TSRD-$\tau$ is not influenced at all
by the presence of one outlier regardless of its magnitude.

Figure 3: The presence of only one extreme outlier in the data can drive the
bootstrap estimate of standard deviation obtained by BLB in two-stage Lasso-LS
into infinity. However, the bootstrap estimates of standard deviation obtained
by RSOB-$\tau$ in TSRD-$\tau$ remains almost unaltered to the presence of one
extreme outlier. The curves produced by TSRD-$\tau$ for all different values
of $\alpha$ overlap, implying that one outlier has almost no effect on
bootstrap estimates of standard deviation based on TSRD-$\tau$.

Theorem 4 states the upper breakdown point of RSOB-$\tau$ bootstrap quantile
estimates is $0.499$ for the simulation set-up described above. In order to
examine the robustness of bootstrap replications, we show that the TSRD-$\tau$
bootstrap replications are robust in the face of outliers even if the data is
contaminated severely by outliers. In this experiment, the outliers are
introduced by randomly choosing a percentage of the observations in
$\mathbf{y}$ and multiplying them with a random value $\alpha_{o}=100000$. The
proposed method TSRD-$\tau$ correctly recovers the true sparse basis with zero
false positive rate even in the presence of severe contamination. As shown in
Fig. 4 , the bootstrap estimate of standard deviation based on the RSOB-$\tau$
is only slightly influenced by the outliers at contamination levels as high as
$40\%$, hence verifying the results in Theorem 4. In other words, the impact
of outliers on the bootstrap estimates is bounded. Note that the curves in
Fig. 3 and Fig. 4 are obtained by averaging over 15 trials of Monte Carlo
simulations

Figure 4: The proposed inference method RSOB-$\tau$ at the stage 2 of
TSRD-$\tau$ exhibit strong resilience to outlier even when $40\%$ of
observations are contaminated by outliers and relative error of standard
deviation remains bounded.

###  V-J Statistical Convergence

In this subsection, the correctness of Theorem 2 is verified by computer
simulations. In other words, we show the distribution of
$\sqrt{n}\big{(}\hat{\boldsymbol{\beta}}_{n,b}^{R\star}-\hat{\boldsymbol{\beta}}_{b}\big{)}$
converges to the limiting distribution of
$\sqrt{n}\big{(}\hat{\boldsymbol{\beta}}_{n}-\boldsymbol{\beta}_{0}\big{)}$ as
$n$ and $b$ approach infinity. Here, we run the simulations on the synthetic
dataset described by Scenario 2, the same settings as in section V-I are used
for the sake of convenience to study the statistical convergence. We assume
that large-scale data is not contaminated by outliers. However, the number of
bootstrap samples within each subset of data is set to $1000$. Under the
condition $\tau$-estimator is tuned for $95\%$ normal efficiency, the limiting
distribution of
$\sqrt{n}\big{(}\hat{\boldsymbol{\beta}}_{n}-\boldsymbol{\beta}_{0}\big{)}$
obeys a multivariate Gaussian distribution
$\mathcal{N}\big{(}\mathbf{0},(\sigma^{2}/0.95)\mathbf{I}_{k_{s}}\big{)}$. The
distribution of
$\sqrt{n}\big{(}\hat{\boldsymbol{\beta}}_{n,b}^{R\star}-\hat{\boldsymbol{\beta}}_{b}\big{)}$
is formed by randomly drawing a subset
$\big{(}\check{\mathbf{y}},\tilde{\mathbf{X}}\big{)}$ from the original data
set
$\big{(}\mathbf{y},\underaccent{\bar}{\bX}\big{)}\in\mathbb{R}^{n\times|\hat{\mathcal{S}}|+1}$,
computing the initial $\tau$-estimate $\hat{\boldsymbol{\beta}}_{b}$ and
performing one-step linear correction of initial $\tau$-estimates for
bootstrap samples by using the derived equations. The plot on the right-hand
side of Fig. 5 shows the empirical distribution of
$\sqrt{n}\big{(}\hat{\boldsymbol{\beta}}_{n,b}^{R\star}-\hat{\boldsymbol{\beta}}_{b}\big{)}$
overlaps the true limiting distribution of
$\sqrt{n}\big{(}\hat{\boldsymbol{\beta}}_{n}-\boldsymbol{\beta}_{0}\big{)}$
for all elements of $\hat{\boldsymbol{\beta}}_{n,b}^{R\star}$. On contrary,
the plot on the left-hand side of Fig. 5 shows the empirical distribution of
$\sqrt{n}\big{(}\hat{\boldsymbol{\beta}}_{n,b}^{1\star}-\hat{\boldsymbol{\beta}}_{b}\big{)}$
underestimates the variability of the true limiting distribution of
$\sqrt{n}\big{(}\hat{\boldsymbol{\beta}}_{n}-\boldsymbol{\beta}_{0}\big{)}$
for all elements of $\hat{\boldsymbol{\beta}}_{n,b}^{1\star}$. Therefore we
can conclude the distribution of
$\sqrt{n}\big{(}\hat{\boldsymbol{\beta}}_{n,b}^{R\star}-\hat{\boldsymbol{\beta}}_{b}\big{)}$
provides a reliable approximation of the distribution of
$\sqrt{n}\big{(}\hat{\boldsymbol{\beta}}_{n}-\boldsymbol{\beta}_{0}\big{)}$.

Figure 5: The distribution of bootstrap replicates produced by linearly
corrected one-step $\tau$-estimator in RSOB-$\tau$ overlaps the true limiting
distribution of $\tau$-estimator whereas the distribution of bootstrap
replicates produced by one-step $\tau$-estimator underestimates the
variability of the true limiting distribution of $\tau$-estimator.

###  V-K Computational Complexity

In this subsection, we compare the computational complexity of the proposed
distributed inference method, RSOB-$\tau$ at stage 2 of TSRD-$\tau$ to a
robust realization of BLB method employing $\tau$-estimator for computing
bootstrap replicates. We run the simulations on the synthetic data illustrated
by Scenario 3 where the proportion of outliers is set to $10\%$. The
experiment was conducted in parallel on a single node of a high-performance
computing cluster (Triton) where $22$ computing cores and $14$ GB of memory
were requested and a Dell PowerEdge C$4130$ node was granted. The cumulative
processing time is recorded after each iteration where new set of bootstrap
samples are successively added to the bags. The cumulative processing time of
RSOB-$\tau$ versus robustified BLB is demonstrated in Fig. 6. As the number of
bootstrap samples increases, the proposed RSOB-$\tau$ requires significantly
less processing time in comparison to robustified BLB. This implies the
RSOB-$\tau$ achieves remarkably higher computational efficiency.

Figure 6: The RSOB-T method employing linearly corrected one-step
$\tau$-estimator is significantly faster than robustified BLB.

###  V-L Overfitting versus Underfitting and their implications on inference

In order to make valid inferences, we shall ensure the model selection
procedure performs reliably. Herein, we examine how confidence intervals are
influenced by overfitting and underfitting in the model selection. To do so,
we run the simulations on the synthetic data described by Scenario 5 where the
proportion of outliers is set to $10\%$. We set the number of data partition
$s=40$. We study the validity of the statistical inference procedure
TSRD-$\tau$ when model selection can not perfectly recover the true support,
resulting in either overfitted models or underfitted models. We first consider
a scenario where the model selection procedure fails to reject all irrelevant
variables associated with zero coefficients of the parameter vector, resulting
in an extremely overfitted model. We further consider the case that the model
selection procedure fails to select all relevant variables and $10$ out of
$15$ relevant variables are classified as irrelevant variables, resulting in
an underfitted model.

We observe that confidence intervals for non-zero coefficients of underfitted
model are much larger than that of overfitted model, as indicated by Fig. 7.
In particular, confidence interval for one coefficient of underfitted model is
very biased. Therefore, we do not only lose information about the relevant
variables not chosen by the model selection in underfitted model. Also, the
inference results may not reliably reflect confidence intervals for the non-
zero coefficients of selected model. Moreover, we observe from Fig. 8 that the
confidence intervals for the given coefficients of parameter vector associated
with false positives cover $0$. We can identify these coefficients as zero and
variables corresponding to them as irrelevant. In this experiment, the
confidence intervals associated with only 5 coefficients did not contain zero.
Hence, the variable selection can still be improved even when one faces
overfitting.

Figure 7:  Confidence intervals for non-zero coefficients of underfitted model
are much larger than that of overfitted model. Besides, confidence interval
for one non-zero coefficient of the underfitted model is extremely biased.
(the dash-dot lines indicate the true value of parameter vector for the
corresponding entry)  Figure 8: In the overfitted model, we observe that
confidence intervals for false positives associated with zero coefficients of
the parameter vector cover zero for the given entries. This implies that we
can identify the variables associated with confidence intervals covering zero
as irrelevant. ( the dash-dot lines indicate the true value of parameter
vector for the corresponding entry)

###  V-M Fusing the Variables Selected by S-Lasso Estimates

We now study how the variable selection algorithm performs when instead of the
variables selected by $\tau$-Lasso estimates, the variables selected by the
initial S-Lasso estimates across nodes are aggregated via the majority voting
scheme. To do so, we carry out a series of Monte-Carlo simulations of 20
trials where a random realization of outlier is used at each trial. We run the
simulations on the synthetic data set described by Scenario 3 for low-
dimensional regime and the synthetic data set described by Scenario 4 for
high-dimensional regime.

Table II shows the result of simulations. We observe that if we aggregate the
variables selected by the initial S-Lasso estimates instead of those selected
by $\tau$-Lasso estimates, the performance will remain almost the same. This
could be explained by the fact that both $\tau$-Lasso and S-Lasso estimators
promote sparsity and are robust to gross outliers. We suspect that low
Gaussian efficiency has little effect on the performance of sparse recovery.

TABLE II:  Fusing the Variables Selected by S-Lasso Estimates: If we use the variables selected by S-Lasso estimates instead of those selected by $\tau$-Lasso, the model selection will remain almost unchanged. |  | TSRD-$\tau$ |  | Fusion / S-Lasso
---|---|---|---|---
| $\xi_{o}$ | TP | FP | CER |  | TP | FP | CER
LD | $0.1$ | 1 | 0 | 0 |  | 1 | 0.0011 | 0.0010
$0.2$ | 1 | 0 | 0 |  | 1 | 0.0017 | 0.0015
HD | $0.1$ | 1 | 0 | 0 |  | 1 | 0 | 0
$0.2$ | 1 | 0 | 0 |  | 1 | 0 | 0

##  VI Conclusion

This paper introduced robust and distributed inference procedures for large
scale data where data exhibits an underlying low-dimensional structure and is
contaminated by outliers. We propose two-stage inference procedures called
TSRD-$\tau$ and TSRD-MM. The former employs the class of robust
$\tau$-estimators whereas the latter employs that of MM-estimators. In the
first stage, active explaining variables are selected by local variable
selection employing robust Lasso estimators. The selections from each node are
combined by applying a fusion rule at the fusion center or cloud. The
selection is broadcast to the computational nodes. In the second stage, actual
inferences on the selected variables are performed by using the robust and
computationally efficient bootstrap procedures. Confidence intervals are
constructed, parameter estimates are found, and standard deviations are
quantified. The favorable statistical properties including consistency and
robustness of the proposed method were established using analytical methods
and verified in simulations. Moreover, the quantitative robustness properties
of robust $\tau$-Lasso were established, in particular its finite-sample
breakdown point.

Future directions of research include extending the classical asymptotic
analysis to high-dimensional asymptotic analysis where $p$ grows with $n$ to
infinity. It is also an open question how one can establish asymptotic results
for local optima. Finally, it would be interesting to devise a rigorous proof
for the model consistency of proposed TSRD-$\tau$ and TSRD-MM procedures.

##  Appendix A Derivation of the Linear Correction Term

The linear correction may be derived by inverting a block-matrix as follows:

$\Big{[}\mathbf{I}-\nabla\mathbf{f}(\boldsymbol{\hat{\theta}}_{b};\tilde{\mathbf{Y}})\Big{]}^{-1}=\begin{bmatrix}\mathcal{A}&\boldsymbol{\eta}\\\
\boldsymbol{\zeta}&a\end{bmatrix}^{-1},$ (24)

where $\mathcal{A}$, $\boldsymbol{\eta}$, $\boldsymbol{\zeta}$ and $a$ are
given by

$\begin{split}\mathcal{A}&=\big{(}\hat{\mathbf{A}}_{b}\big{)}^{-1}\Bigg{[}\mathcal{A}_{2}-\mathcal{A}_{1}\Bigg{]},\\\
\boldsymbol{\eta}&=\big{(}\hat{\mathbf{A}}_{b}\big{)}^{-1}\Bigg{[}\mathcal{\eta}_{2}-\mathcal{\eta}_{1}\Bigg{]},\\\
\boldsymbol{\zeta}&=\frac{1}{b\delta_{2}}\sum_{l=1}^{b}\rho_{0}^{{}^{\prime}}(\tilde{r}_{l})\tilde{\mathbf{x}}_{[l]}^{T},\\\
a&=\frac{1}{b\delta_{2}}\sum_{l=1}^{b}\rho_{0}^{{}^{\prime}}(\tilde{r}_{l})\tilde{r}_{l}.\end{split}$
(25)

$\mathcal{A}_{1}$, $\mathcal{A}_{2}$, $\boldsymbol{\eta}_{1}$ and
$\boldsymbol{\eta}_{2}$ are calculated as follows:

$\begin{split}\mathcal{A}_{1}&=\frac{1}{b}\sum_{l=1}^{b}\tilde{\mathbf{x}}_{[l]}\nabla_{\boldsymbol{\beta}}w_{\tau}\rho_{0}^{{}^{\prime}}(\tilde{r}_{l}),\\\
\mathcal{A}_{2}&=\frac{1}{b\hat{\sigma}_{b}}\sum_{l=1}^{b}\Big{[}\hat{w}_{\tau}\rho_{0}^{{}^{\prime\prime}}(\tilde{r}_{l})+\rho_{1}^{{}^{\prime\prime}}(\tilde{r}_{l})\Big{]}\tilde{\mathbf{x}}_{[l]}\tilde{\mathbf{x}}_{[l]}^{T},\\\
\mathcal{\eta}_{1}&=\frac{1}{b}\sum_{l=1}^{b}\nabla_{\sigma}w_{\tau}\rho_{0}^{{}^{\prime}}(\tilde{r}_{l})\tilde{\mathbf{x}}_{[l]},\\\
\mathcal{\eta}_{2}&=\frac{1}{b\hat{\sigma}_{b}}\sum_{l=1}^{b}\big{[}\hat{w}_{\tau}\rho_{0}^{{}^{\prime\prime}}(\tilde{r}_{l})+\rho_{1}^{{}^{\prime\prime}}(\tilde{r}_{l})\Big{]}\tilde{\mathbf{x}}_{[l]}\tilde{r}_{l},\end{split}$
(26)

where $\nabla_{\boldsymbol{\beta}}w_{\tau}$ and $\nabla_{\sigma}w_{\tau}$ are
shorthands for $\partial
w_{\tau}(\boldsymbol{\hat{\theta}}_{b})/\partial\boldsymbol{\beta}$ and
$\partial w_{\tau}(\boldsymbol{\hat{\theta}}_{b})/\partial\sigma$,
respectively. $\hat{\mathbf{A}}_{b}$, $\nabla_{\boldsymbol{\beta}}w_{\tau}$
and $\nabla_{\sigma}w_{\tau}$ are given by

$\hat{\mathbf{A}}_{b}=\frac{1}{b}\sum_{l=1}^{b}\hat{w}_{l}^{(i)}\tilde{\mathbf{x}}_{[l]}\tilde{\mathbf{x}}_{[l]}^{T},$
(27)

$\begin{split}\nabla_{\boldsymbol{\beta}}w_{\tau}=\frac{\sum_{l=1}^{b}\Big{[}\rho_{1}^{{}^{\prime\prime}}(\tilde{r}_{l})\tilde{r}_{l}-\rho_{1}^{{}^{\prime}}(\tilde{r}_{l})\Big{]}\tilde{\mathbf{x}}_{[l]}^{T}/\hat{\sigma}_{b}}{\sum_{l=1}^{b}\rho_{0}^{{}^{\prime}}(\tilde{r}_{l})\tilde{r}_{l}},\\\
+\frac{\sum_{l=1}^{b}\Big{[}\rho_{0}^{{}^{\prime\prime}}(\tilde{r}_{l})\tilde{r}_{l}+\rho_{0}^{{}^{\prime}}(\tilde{r}_{l})\Big{]}\tilde{\mathbf{x}}_{[l]}^{T}/\hat{\sigma}_{b}}{\sum_{l=1}^{b}\rho_{0}^{{}^{\prime}}(\tilde{r}_{l})\tilde{r}_{l}}\hat{w}_{\tau}.\end{split}$
(28)

$\begin{split}\nabla_{\sigma}w_{\tau}=\frac{\sum_{l=1}^{b}\Big{[}\rho_{1}^{{}^{\prime\prime}}(\tilde{r}_{l})\tilde{r}_{l}-\rho_{1}^{{}^{\prime}}(\tilde{r}_{l})\Big{]}\tilde{r}_{l}/\hat{\sigma}_{b}}{\sum_{l=1}^{b}\rho_{0}^{{}^{\prime}}(\tilde{r}_{l})\tilde{r}_{l}}\\\
+\frac{\sum_{l=1}^{b}\Big{[}\rho_{0}^{{}^{\prime\prime}}(\tilde{r}_{l})\tilde{r}_{l}+\rho_{0}^{{}^{\prime}}(\tilde{r}_{l})\Big{]}\tilde{r}_{l}/\hat{\sigma}_{b}}{\sum_{l=1}^{b}\rho_{0}^{{}^{\prime}}(\tilde{r}_{l})\tilde{r}_{l}}\hat{w}_{\tau}\end{split}$
(29)

On the other hand,

$\begin{bmatrix}\mathcal{A}&\boldsymbol{\eta}\\\
\boldsymbol{\zeta}&a\end{bmatrix}^{-1}=\begin{bmatrix}\mathbf{M}_{b}&\mathbf{d}_{b}\\\
\mathbf{N}_{b}&\mathbf{q}_{b}\end{bmatrix}.$ (30)

##  Appendix B Derivation of the One-step Bootstrap Replicates

The one-step bootstrap replicates $\hat{\boldsymbol{\beta}}_{n,b}^{1\star}$
and $\hat{\sigma}_{n,b}^{1\star}$ are calculated as follows:

$\begin{split}\hat{\boldsymbol{\beta}}_{n,b}^{1\star}&=\Big{(}\sum_{l=1}^{b}\omega_{l}^{\star}w_{l}^{\star}\tilde{\mathbf{x}}_{[l]}\tilde{\mathbf{x}}_{[l]}^{T}\Big{)}^{-1}\sum_{l=1}^{b}\omega_{l}^{\star}w_{l}^{\star}\check{y}_{l}\tilde{\mathbf{x}}_{[l]},\\\
\hat{\sigma}_{n,b}^{1\star}&=\sum_{l=1}^{b}\omega_{l}^{\star}\check{v}_{l}^{\star}\big{(}\check{y}_{l}-\tilde{\mathbf{x}}_{[l]}^{T}\hat{\boldsymbol{\beta}}_{b}\big{)},\end{split}$
(31)

where $v_{l}^{\star}$, $w_{l}^{\star}$ and $w_{\tau}^{\star}$ are computed as
follows:

$\begin{split}v_{l}^{\star}&=\frac{b}{n}\hat{v}_{l},\\\
w_{l}^{\star}&=\frac{w_{\tau}^{\star}\rho_{0}^{{}^{\prime}}(\tilde{r}_{l})+\rho_{1}^{{}^{\prime}}(\tilde{r}_{l})}{\hat{r}_{l}},\\\
w_{\tau}^{\star}&=\frac{\sum_{l=1}^{b}\omega_{l}^{\star}\Big{[}2\rho_{1}(\tilde{r}_{l})-\rho_{1}^{{}^{\prime}}(\tilde{r}_{l})\tilde{r}_{l}\Big{]}}{\sum_{l=1}^{b}\omega_{l}^{\star}\rho_{0}^{{}^{\prime}}(\tilde{r}_{l})\tilde{r}_{l}}.\end{split}$
(32)

Therefore, the linearly corrected one-step bootstrap replications using
$\tau$-estimators are calculated as follows:

$\begin{split}\hat{\boldsymbol{\beta}}_{n,b}^{R\star}&=\hat{\boldsymbol{\beta}}_{b}+\mathbf{M}_{b}\Big{(}\hat{\boldsymbol{\beta}}_{n,b}^{1\star}-\hat{\boldsymbol{\beta}}_{b}\Big{)}+\mathbf{d}_{b}\Big{(}\hat{\sigma}_{n,b}^{1\star}-\hat{\sigma}_{b}\Big{)},\\\
\mathbf{M}_{b}&=\Big{(}\mathcal{A}-\boldsymbol{\eta}a^{-1}\boldsymbol{\zeta}\Big{)}^{-1},\\\
\mathbf{d}_{b}&=-\mathcal{A}^{-1}\boldsymbol{\eta}\Big{(}a-\boldsymbol{\zeta}\mathcal{A}^{-1}\boldsymbol{\eta}\Big{)}^{-1}.\end{split}$
(33)

## References

  * [1] K. Slavakis, G. B. Giannakis, and G. Mateos, “Modeling and optimization for big data analytics:(statistical) learning tools for our era of data deluge,” _IEEE Signal Processing Magazine_ , vol. 31, no. 5, pp. 18–31, 2014.
  * [2] R. Tibshirani, “Regression shrinkage and selection via the lasso,” _Journal of the Royal Statistical Society: Series B (Methodological)_ , vol. 58, no. 1, pp. 267–288, 1996.
  * [3] A. Chatterjee and S. N. Lahiri, “Bootstrapping lasso estimators,” _Journal of the American Statistical Association_ , vol. 106, no. 494, pp. 608–625, 2011.
  * [4] A. Javanmard and A. Montanari, “Confidence intervals and hypothesis testing for high-dimensional regression,” _The Journal of Machine Learning Research_ , vol. 15, no. 1, pp. 2869–2909, 2014.
  * [5] H. Liu, B. Yu _et al._ , “Asymptotic properties of lasso+ mls and lasso+ ridge in sparse high-dimensional linear regression,” _Electronic Journal of Statistics_ , vol. 7, pp. 3124–3169, 2013.
  * [6] A. Chatterjee, S. N. Lahiri _et al._ , “Rates of convergence of the adaptive lasso estimators to the oracle distribution and higher order refinements by the bootstrap,” _The Annals of Statistics_ , vol. 41, no. 3, pp. 1232–1259, 2013.
  * [7] L. Wasserman and K. Roeder, “High dimensional variable selection,” _Annals of statistics_ , vol. 37, no. 5A, p. 2178, 2009.
  * [8] H. Liu, X. Xu, and J. J. Li, “A bootstrap lasso+ partial ridge method to construct confidence intervals for parameters in high-dimensional sparse linear models,” _arXiv preprint arXiv:1706.02150_ , 2017.
  * [9] J. D. Lee, D. L. Sun, Y. Sun, J. E. Taylor _et al._ , “Exact post-selection inference, with application to the lasso,” _Annals of Statistics_ , vol. 44, no. 3, pp. 907–927, 2016.
  * [10] R. Berk, L. Brown, A. Buja, K. Zhang, L. Zhao _et al._ , “Valid post-selection inference,” _The Annals of Statistics_ , vol. 41, no. 2, pp. 802–837, 2013.
  * [11] C.-H. Zhang and S. S. Zhang, “Confidence intervals for low dimensional parameters in high dimensional linear models,” _Journal of the Royal Statistical Society: Series B: Statistical Methodology_ , pp. 217–242, 2014.
  * [12] R. Lockhart, J. Taylor, R. J. Tibshirani, and R. Tibshirani, “A significance test for the lasso,” _Annals of statistics_ , vol. 42, no. 2, p. 413, 2014\.
  * [13] R. F. Barber, E. J. Candès _et al._ , “Controlling the false discovery rate via knockoffs,” _Annals of Statistics_ , vol. 43, no. 5, pp. 2055–2085, 2015.
  * [14] P. Bühlmann _et al._ , “Statistical significance in high-dimensional linear models,” _Bernoulli_ , vol. 19, no. 4, pp. 1212–1242, 2013.
  * [15] M. Martinez-Camara, M. Muma, B. Bejar, A. M. Zoubir, and M. Vetterli, “The regularized tau estimator: A robust and efficient solution to ill-posed linear inverse problems,” _arXiv preprint arXiv:1606.00812_ , 2016.
  * [16] E. Smucler and V. J. Yohai, “Robust and sparse estimators for linear regression models,” _Computational Statistics & Data Analysis_, vol. 111, pp. 116–130, 2017.
  * [17] G. V. C. Freue, D. Kepplinger, M. Salibián-Barrera, and E. Smucler, “Pense: A penalized elastic net s-estimator,” 2017.
  * [18] A. Alfons, C. Croux, and S. Gelper, “Sparse least trimmed squares regression for analyzing high-dimensional large data sets,” _The Annals of Applied Statistics_ , pp. 226–248, 2013.
  * [19] Q. Zheng, C. Gallagher, and K. Kulasekera, “Robust adaptive lasso for variable selection,” _Communications in Statistics-Theory and Methods_ , vol. 46, no. 9, pp. 4642–4659, 2017.
  * [20] A. Kleiner, A. Talwalkar, P. Sarkar, and M. I. Jordan, “A scalable bootstrap for massive data,” _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , vol. 76, no. 4, pp. 795–816, 2014.
  * [21] M. Salibian-Barrera, R. H. Zamar _et al._ , “Bootrapping robust estimates of regression,” _The Annals of Statistics_ , vol. 30, no. 2, pp. 556–582, 2002.
  * [22] S. Basiri, E. Ollila, and V. Koivunen, “Robust, scalable, and fast bootstrap method for analyzing large scale data,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 4, pp. 1007–1017, 2015.
  * [23] E. Mozafari-Majd and V. Koivunen, “Robust, sparse and scalable inference using bootstrap and variable selection fusion,” in _2019 IEEE 8th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP)_. IEEE, 2019, pp. 271–275.
  * [24] ——, “Robust variable selection and distributed inference using $\tau$-based estimators for large-scale data,” in _2020 28th European Signal Processing Conference (EUSIPCO)_. IEEE, 2021, pp. 2453–2457.
  * [25] A. Ghosh and M. Thoresen, “A robust variable screening procedure for ultra-high dimensional data,” _arXiv preprint arXiv:2004.14851_ , 2020.
  * [26] X. Chen and M.-g. Xie, “A split-and-conquer approach for analysis of extraordinarily large data,” _Statistica Sinica_ , pp. 1655–1684, 2014.
  * [27] F. R. Bach, “Bolasso: model consistent lasso estimation through the bootstrap,” in _Proceedings of the 25th international conference on Machine learning_ , 2008, pp. 33–40.
  * [28] M. Martinez-Camara, M. Muma, A. M. Zoubir, and M. Vetterli, “A new robust and efficient estimator for ill-conditioned linear inverse problems with outliers,” in _2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2015, pp. 3422–3426.
  * [29] V. J. Yohai and R. H. Zamar, “High breakdown-point estimates of regression by means of the minimization of an efficient scale,” _Journal of the American statistical association_ , vol. 83, no. 402, pp. 406–413, 1988.
  * [30] R. A. Maronna, R. D. Martin, V. J. Yohai, and M. Salibián-Barrera, _Robust statistics: theory and methods (with R)_. John Wiley & Sons, 2019.
  * [31] F. H. Clarke, _Optimization and nonsmooth analysis_. Siam, 1990, vol. 5.
  * [32] C. S. Burrus, J. Barreto, and I. W. Selesnick, “Iterative reweighted least-squares design of fir filters,” _IEEE Transactions on Signal Processing_ , vol. 42, no. 11, pp. 2926–2936, 1994.
  * [33] P. W. Holland and R. E. Welsch, “Robust regression using iteratively reweighted least-squares,” _Communications in Statistics-theory and Methods_ , vol. 6, no. 9, pp. 813–827, 1977.
  * [34] W. R. Schucany and S. Wang, “One-step bootstrapping for smooth iterative procedures,” _Journal of the Royal Statistical Society: Series B (Methodological)_ , vol. 53, no. 3, pp. 587–596, 1991.
  * [35] M. Salibian-Barrera, S. Van Aelst, and V. J. Yohai, “Robust tests for linear regression models based on $\tau$-estimates,” _Computational Statistics & Data Analysis_, vol. 93, pp. 436–455, 2016.
  * [36] R. V. Hogg, “Some observations on robust estimation,” _Journal of the American Statistical Association_ , vol. 62, no. 320, pp. 1179–1186, 1967.
  * [37] K. Singh _et al._ , “Breakdown theory for bootstrap quantiles,” _The annals of Statistics_ , vol. 26, no. 5, pp. 1719–1732, 1998.
  * [38] D. Kepplinger, M. Salibián-Barrera, and G. Cohen Freue, “https://cran.r-project.org/package=pense,” 2021.
  * [39] G. V. C. Freue, D. Kepplinger, M. Salibián-Barrera, and E. Smucler, “Robust elastic net estimators for variable selection and identification of proteomic biomarkers,” _The Annals of Applied Statistics_ , vol. 13, no. 4, pp. 2065–2090, 2019.
  * [40] R. Tomioka, T. Suzuki, and M. Sugiyama, “Super-linear convergence of dual augmented lagrangian algorithm for sparsity regularized estimation,” _Journal of Machine Learning Research_ , vol. 12, no. May, pp. 1537–1586, 2011.
  * [41] C. Croux, G. Dhaene, and D. Hoorelbeke, “Robust standard errors for robust estimators,” _CES-Discussion paper series (DPS) 03.16_ , pp. 1–20, 2004\.
  * [42] R. A. Maronna, “Robust ridge regression for high-dimensional data,” _Technometrics_ , vol. 53, no. 1, pp. 44–53, 2011.
  * [43] M. J. Wainwright, “Information-theoretic limits on sparsity recovery in the high-dimensional and noisy setting,” _IEEE transactions on information theory_ , vol. 55, no. 12, pp. 5728–5741, 2009.
  * [44] J. A. Khan, S. Van Aelst, and R. H. Zamar, “Robust linear model selection based on least angle regression,” _Journal of the American Statistical Association_ , vol. 102, no. 480, pp. 1289–1299, 2007.
  * [45] G. Schwarz _et al._ , “Estimating the dimension of a model,” _Annals of statistics_ , vol. 6, no. 2, pp. 461–464, 1978.
  * [46] Y. Fan and C. Y. Tang, “Tuning parameter selection in high dimensional penalized likelihood,” _Journal of the Royal Statistical Society: SERIES B: Statistical Methodology_ , pp. 531–552, 2013.
  * [47] Y. Kim, S. Kwon, and H. Choi, “Consistent model selection criteria on high dimensions,” _The Journal of Machine Learning Research_ , vol. 13, no. 1, pp. 1037–1057, 2012.
  * [48] A. Ghosh and S. Majumdar, “Ultrahigh-dimensional robust and efficient sparse regression using non-concave penalized density power divergence,” _IEEE Transactions on Information Theory_ , vol. 66, no. 12, pp. 7812–7827, 2020.
  * [49] A. Alfons, “robusthd: An r package for robust regression with high-dimensional data,” _Journal of Open Source Software_ , vol. 6, no. 67, p. 3786, 2021\.
  * [50] Q. Li, Z. Zhu, and G. Tang, “The non-convex geometry of low-rank matrix optimization,” _Information and Inference: A Journal of the IMA_ , vol. 8, no. 1, pp. 51–96, 2019.
  * [51] Y. Chi, Y. M. Lu, and Y. Chen, “Nonconvex optimization meets low-rank matrix factorization: An overview,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 20, pp. 5239–5269, 2019.
  * [52] V. J. Yohai and R. H. Zamar, “High breakdown-point estimates of regression by means of the minimization of an efficient scale,” _University of Washington, Dept. of Statistics. Technical Report 84_ , pp. 1–41, 1986.
  * [53] G. V. C. Freue, D. Kepplinger, M. Salibián-Barrera, and E. Smucler, “Proteomic biomarker study using novel robust penalized elastic net estimators,” _the Annals of Applied Statistics, submitted_.
  * [54] J. Shao, _Mathematical statistics: exercises and solutions_. Springer Science & Business Media, 2006.
  * [55] P. J. Bickel and D. A. Freedman, “Some asymptotic theory for the bootstrap,” _The annals of statistics_ , pp. 1196–1217, 1981.
  * [56] P. Rousseeuw and V. Yohai, “Robust regression by means of s-estimators,” in _Robust and nonlinear time series analysis_. Springer, 1984, pp. 256–272.

|  Emadaldin Mozafari-Majd
---|---
|  Visa Koivunen
---|---

This supplemental material contains a section on initial variable screening
(Preprocessing) for very high-dimensional settings along with additional
simulation results and technical proofs of the theorem discussed in the paper.

## Initial Variable Screening (Preprocessing): Data with Very High-Dimensional
Subsets

In order to reduce the computational burden in settings with very high-
dimensional subsets, an initial variable screening procedure called Density
Power Divergence-SIS (DPD-SIS) [25] is employed to further reduce the model
complexity, the number of variables to an order of sample size $n$. In
distributed and parallel architecture, we use the variable screening only in
very high-dimensional settings where $p$ is much larger than $b$. Basically,
the DPD-SIS extends the Sure Independence Screening (SIS) to address
robustness in the presence of outliers. The robust screening procedure assigns
a certain score to each predictor and then predictors are ranked in descending
order based on the calculated score. In order to compute the score, the
marginal estimate of each regression coefficient is obtained via a minimum DPD
estimator and then its absolute value determines the score associated with
each predictor. At each node, the DPD-SIS variable screening procedure is
utilized to discard a certain number of irrelevant predictors and then
distinct subsets with reduced dimensionality is passed down to next step for
further processing, that is, $\check{\mathbf{X}}\in\mathbb{R}^{b\times
p}\rightarrow\bar{\mathbf{X}}\in\mathbb{R}^{b\times q}$ with $(q\ll p)$ an
order of sample size. A reasonable choice for $q$ is the number of
observations within subsets of data, that is, $q=b$. Note that variable
screening procedure is dispensable in low-dimensional models, i.e. $q=p$.
Subdividing a high-dimensional data into smaller subsets may be allowed as
long as the batch size $b$ satisfies the requirements for sparse recovery and
statistical correctness of bootstrap computations as specified in Section V.E,
Choice of Batch Size, within the main body of the paper. The algorithmic
details of DPD-SIS can be found in [25].

The DPD-SIS is initially applied to the distinct subsets of data and top $b$
predictors, the exact order of sample size, within each subset of data are
kept based on their score and the remaining predictors are discarded. We then
form a set of predictors appearing within at least half of data subsets. We
now place these set of predictors on top of the ranking within each of $s$
data subsets, keep the top $b$ predictors and reject the remaining ones. The
model selection proceeds with the data set of reduced dimensionality.

## Additional Simulation Results

###  B-A High-dimensional: Model Selection and Inference

In this part, we study the performance of the proposed methods, TSRD-$\tau$
and TSRD-MM, for a large-scale high-dimensional data set ($p>b$) in terms of
model selection and robustness of confidence intervals to outliers. We run the
simulations on the synthtetic dataset described by Scenario 1. In the current
high-dimensional setting, an initial variable screening procedure is employed
to reduce the dimensionality prior to model selection. The DPD-SIS is
initially applied to the distinct subsets of data and top $b$ covariates, the
exact order of sample size, within each subset of data are kept based on their
score and the remaining covariates are discarded. The model selection proceeds
with the data set of reduced dimensionality.

TABLE III:  Model Selection in High-Dimensional: The Proposed Two-Stage Robust
Inference Methods Achieve a Perfect Recovery of True Sparse Basis (TP=1) with
Small Number of False Positives under Zero Contamination to Moderate
Contamination. In Contrast, the Two-Stage Lasso-LS Completely Fails at
Recovering the Sparse Basis ($\text{TP}=0$) for all Scenarios Except for
Outlier-Free.

| TSRD-$\tau$ |  | TSRD-MM |  | TSLL
---|---|---|---|---|---
$\xi_{o}$ | TP | FP | CER |  | TP | FP | CER |  | TP | FP | CER
$0$ | 1 | 0.0004 | 0.0004 |  | 1 | 0.0013 | 0.0013 |  | 1 | 0 | 0
$0.1$ | 1 | 0.0003 | 0.0003 |  | 1 | 0.0009 | 0.0009 |  | 0 | 0 | 0.0013
$0.2$ | 1 | 0.0002 | 0.0002 |  | 1 | 0.0002 | 0.0002 |  | 0 | 0 | 0.0013

  * 1

Although no truly active variables are selected by TSLL in the model selection
stage, one might be mislead by low CER. This can be attributed to fact that
the number of truly active variables $k_{s}$ are insignificant compared to
$p$.

The proposed two-stage robust inference algorithms are compared to their two-
stage non-robust counterpart employing Lasso in the first stage and BLB based
on least square estimator in the second stage. For the sake of brevity, the
non-robust two-stage inference method is regarded as two-stage Lasso-LS
(TSLL). Before using Lasso, variable screening is carried out by the DPD-SIS
to ensure the data with the reduced dimensionality contains relevant
variables. In regard to model selection, the performance is quantified by
using confusion matrix and CER as represented in Table III. The model
selection algorithms in the stage 1 of TSRD-$\tau$ and TSRD-MM could perfectly
identify all true non-zero parameters for different proportions of outliers at
the highly underdetermined setting, $p/b=33.33$. In contrast, the non-robust
TSLL algorithm fails completely at identifying the sparse basis of the
parameter vector except for outlier-free scenario.

Figure 9:  When there are no outlying observations, the confidence intervals
based on bootstraps methods in TSRD-$\tau$ and TSRD-MM provide reliable
estimates of the CIs based on bootstrap percentiles of least-square estimator.
The confidence intervals produced by TSRD-T and TSRD-MM methods exhibit
robustness to outliers and their lengths are slightly affected by an increase
in the proportion of outliers. Thus, reliable parameter estimates and
confidence intervals are obtained even in the presence of outliers (the dash-
dot lines indicate the true value of non-zero entries of parameter vector).

The selected variables from the first stage are used to construct confidence
intervals based on the bootstrap methods RSOB-$\tau$ and BLFRB. The confidence
intervals constructed for the first $15$ selected variables are shown in Fig.
9. The CIs formed by robust inference methods at the stage 2 of TSRD-$\tau$
and TSRD-MM remain resistant to contamination and length of CIs are slightly
inflated with an increase in the proportion of outliers. In regard to outlier-
free scenario, the robust bootstrap methods provide reliable estimates of the
CIs constructed by bootstrap percentiles of least-square estimator. Hence, it
can be concluded that the proposed two-stage inference methods, TSRD-$\tau$
and TSRD-MM, can be used to perform robust statistical inference for large-
scale high-dimensional data sets. Note that no comparison is made to two-stage
Lasso-LS in the presence of outlying observations due to zero true positive
rate, i.e., no element of sparse basis was identified.

###  B-B Effect of Substituting $\tau$-Lasso with non-regularized $\tau$\- and
MM-estimators on Variable Selection

We study the effect of substituting the $\tau$-Lasso estimator with
$\tau$-estimator and MM-estimator under the settings described by Scenario 3
and $10\%$ outlier ratio. We conduct the above experiment with $20$ trials for
which a random realization of outlier is used at each trial. The $\tau$\- and
MM-estimators were initialized with S-Lasso estimates of the parameter vector
$\boldsymbol{\beta}$. We observe in our simulations that both $\tau$\- and MM-
estimators fail in recovering the correct support and all irrelevant variables
are selected within the estimated support In contrast, TSRD-$\tau$ employing
the $\tau$-Lasso estimator succeeds in recovering the correct support,
indicated by the following contingency table.

TABLE IV:  Model Selection: The Proposed TSRD-$\tau$ Method Achieves a Perfect Recovery of True Sparse Basis (TP=1) with no False Positives ($\text{FP}=0$). In Contrast, all irrelevant variables are selected within the estimated support when one replaces the $\tau$-Lasso estimator with non-regularized estimators | TSRD-$\tau$ |  | $\tau$ |  | MM
---|---|---|---|---|---
$\xi_{o}$ | TP | FP | CER |  | TP | FP | CER |  | TP | FP | CER
$0.1$ | 1 | 0 | 0 |  | 1 | 1 | 0.9 |  | 1 | 1 | 0.9

Note that we obtained the results demonstrated in the above table by averaging
over 20 trials.

###  B-C Effect of Initial Estimate on Variable Selection

In order to explore the influence of initial estimates on the variable
selection procedure, we perform a series of simulations with high-dimensional
and low-dimensional data. We run the $\tau$-Lasso estimator when randomly
initialized and then when initialized with non-regularized $\tau$-estimates
(applicable only to low-dimensional regime) and compare the results of model
selection with recommended procedure where the $\tau$-Lasso estimator is using
S-Lasso estimates as the initial point. In both high-dimensional and low-
dimensional data, we carry out a Monte-Carlo study of 20 trials where a random
realization of outlier is used at each trial. In case of $\tau$-Lasso
estimator with random initialization, we initiate the algorithm with a
randomly distributed multivariate Gaussian
$\mathcal{N}(1000\times\mathbf{1}_{p+1},(250)^{2}\times\textbf{I}_{p+1})$ for
each trial and batch of data, chosen to be far from the true coefficient
$\boldsymbol{\beta}_{0}$. We run the simulations on the synthetic data set
described by Scenario 3 for low-dimensional regime and the synthetic data set
described by Scenario 4 for high-dimensional regime.

As shown in Tables V-VI, the model selection algorithm using $\tau$-Lasso
achieves the exact support recovery regardless of initialization across all
trials for the low-dimensional and high-dimensional regimes. Although the
model selection procedure succeeds in perfectly recovering the true support
when $\tau$-Lasso estimator initialized randomly. In fact, the $\tau$-Lasso
optimization problem is solved via alternating minimization where the sub-
problems are non-convex themselves. Recent results show that many well-known
nonconvex optimization problems possess a well behaved landscape where all
second-order stationary points are global minima [50] and [51]. We conjecture
that the variable succeeds in recovering the true support when the
$\tau$-Lasso optimization problem is initialized randomly due to potential
nice landscape of optimization problem in conjunction with collaborative
nature of fusion procedure. We suspect this result may not entirely generalize
to all scenarios and this topic requires further study from the optimization
perspective.

TABLE V:  Influence of Initial Estimate on Model Selection in Low-Dimensional Data: Regardless of How the $\tau$-Lasso Estimator is Initialized, the Model Selection Method Succeeds in Perfectly Recovering the True Sparse Basis (TP=1) with no False Positives ($\text{FP}=0$) under Contamination. | TSRD-$\tau$ |  | random init. |  | init. by $\tau$-estimates
---|---|---|---|---|---
$\xi_{o}$ | TP | FP | CER |  | TP | FP | CER |  | TP | FP | CER
$0.1$ | 1 | 0 | 0 |  | 1 | 0 | 0 |  | 1 | 0 | 0
$0.2$ | 1 | 0 | 0 |  | 1 | 0 | 0 |  | 1 | 0 | 0
TABLE VI:  Influence of Initial Estimate on Model Selection in High-Dimensional Data: Regardless of How the $\tau$-Lasso Estimator is Initialized, the Model Selection Method Succeeds in Perfectly Recovering the True Sparse Basis (TP=1) with no False Positives ($\text{FP}=0$) under Contamination. | TSRD-$\tau$ |  | random init.
---|---|---|---
$\xi_{o}$ | TP | FP | CER |  | TP | FP | CER
$0.1$ | 1 | 0 | 0 |  | 1 | 0 | 0
$0.2$ | 1 | 0 | 0 |  | 1 | 0 | 0

###  B-D The Effect of Sample Size, Dimensionality , and Number of Subsamples
on Computational Complexity of RSOB-$\tau$

We conduct a number of experiments to examine how the computational complexity
of RSOB-$\tau$ scales with sample size, dimensionality and the number of
subsamples. To do so, we plot the processing time against one of the above
parameters while keeping the remaining parameters at their defaults. This
process continues until all parameters have been allowed to vary for a range
of values. We run the simulations on the data set described by Scenario 3
where $10\%$ of observations are contaminated by outliers. We then compare
them to the result of simulations obtained using robust realization of BLB
method based on $\tau$-estimator and report the corresponding results in Fig.
10-12. We conduct the experiment on a single node of a high-performance
computing cluster (Triton) where 22 computing cores and $25$ GB of memory were
used. Note that we report the processing time of RSOB-$\tau$ versus that of
robustified BLB where a distinct subsample is allocated to each computing
core. The number of bootstrap samples is fixed at $B=300$.

Fig. 10 shows that as we increase the sample size $n$, the processing time
associated with RSOB-$\tau$ grows at much slower rate than the processing time
associated with the BLB method based on $\tau$-estimator, thereby achieving
higher computational efficiency for larger sample size. In Fig. 11, we obtain
similar results when plotting the processing time against dimensionality $p$.
That is, RSOB-$\tau$ achieves higher computational efficiency for larger
dimensions $p$.

Figure 10: The processing time of RSOB-$\tau$ method employing linearly
corrected one-step $\tau$-estimator scales much better with sample size than
the robustified BLB method. Figure 11: The processing time of RSOB-$\tau$
method employing linearly corrected one-step $\tau$-estimator scales much
better with dimensionality than the robustified BLB method.

Fig. 12 shows that as we increase the number of subsets (data partitions), the
processing time associated with both RSOB-$\tau$ and the robustified BLB
decreases and simultaneously the gap between the two curves shrinks.
Therefore, there would little benefit in excessively increasing the number of
subsets when considering the processing time. We assume that each subsample is
processed by allocating one computing core to each subsample. Setting $s$ to
very large values leads to wasted computation as we are occupying so many
resources while a small portion of each resource is needed. One could use
$\lfloor n/n^{0.9}\rfloor<s<\lfloor{n/\max(c_{s}k_{s}\log(p),n^{0.6})}\rfloor$
as a crude approximation and choose $b=n/s$ so that it uses the computational
and storage capabilities of the computational nodes efficiently. As it would
be difficult to allocate computing cores for large numbers, we record the
processing time associated with each subsample and their maximum is considered
to be the processing time for better interpretability.

Figure 12: The above plot indicates if one excessively increases the number of
data partitions, small computational gains are made by using RSOB-$\tau$ than
using the robustified BLB.

## Proof of Theorem 1

The finite-sample breakdown point of $\tau$-Lasso estimator can be derived in
two separate stages: first stage is dedicated to show the boundedness of
finite-sample breakdown point from above and second stage is dedicated to show
the boundedness of finite-sample breakdown point from below. Once the
boundedness from above and below is established, the above theorem is proved
immediately. Before proceeding any further, we extend the lemma 5.1 [30] for
M-scale estimators to $\tau$-scale estimators as follows:

Lemma 1: Consider any sequence of samples
$\Big{(}\hat{\mathbf{Y}}^{(k)}=\big{(}\hat{\mathbf{y}}^{(k)},\hat{\mathbf{X}}^{(k)}\big{)}\in\mathbb{R}^{b\times(p+1)}\Big{)}_{k\in\mathbb{N}}=\Big{(}\hat{\mathbf{Y}}^{(1)},\hat{\mathbf{Y}}^{(2)},\hat{\mathbf{Y}}^{(3)},\cdots\Big{)}$
and corresponding residual vector
$\mathbf{r}^{(k)}=\hat{\mathbf{y}}^{(k)}-\hat{\mathbf{X}}^{(k)}\hat{\boldsymbol{\beta}}^{(k)}$
for $\hat{\mathbf{Y}}^{(k)}$. Suppose $r^{(k)}_{l}$ denotes the residual for
$\big{(}\hat{y}^{(k)}_{l},(\hat{\mathbf{x}}_{[l]}^{(k)})^{T}\big{)}$ a row in
$\hat{\mathbf{Y}}^{(k)}$, Then

  1. 1.

Let $C=\\{l:|r_{l}^{(k)}|\rightarrow\infty\\}$, if $\\#(C)>b\delta$, then
$\hat{\sigma}_{\tau}(\mathbf{r}^{(k)})\rightarrow\infty$ as
$k\rightarrow\infty$.

  2. 2.

Let $D=\\{l:|r_{l}^{(k)}|\mbox{ is bounded }\\}$, if $\\#(D)>b-b\delta$, then
$\hat{\sigma}_{\tau}(\mathbf{r}^{(k)})$ is bounded.

where $\rightarrow\infty$ denotes the left-hand side of arrow tends to
infinity. The part $1$ of the lemma implies that if the number of unbounded
entries of $\mathbf{r}^{(k)}$ exceeds $b\delta$, then $\tau$-scale estimate of
$\mathbf{r}^{(k)}$ goes to infinity. On the other hand, if the number of
bounded entries of $\mathbf{r}^{(k)}$ exceeds $b-b\delta$, then $\tau$-scale
estimate of $\mathbf{r}^{(k)}$ remains bounded. To prove the above lemma, we
use the results from lemma 3.2 in [52], the $\tau$-scale is bounded from above
and below as follows:

$\bar{c}_{i}\hat{\sigma}_{\text{M}}(\mathbf{u})\leq\hat{\sigma}_{\tau}(\mathbf{u})\leq\sqrt{\sup_{t\in\mathbb{R}}\rho_{1}(t)}\hat{\sigma}_{\text{M}}(\mathbf{u})\quad\forall\mathbf{u}\in\mathbb{R}^{b},$
(34)

where $\bar{c}_{i}$ is a positive constant, $\hat{\sigma}_{\text{M}}$ denotes
the M-scale estimate and $\hat{\sigma}_{\tau}$ denotes the $\tau$-scale
estimate. Now, we can take advantage of the inequality associating M-scale
estimate with $\tau$-scale estimate. We extend the lemma 5.1 in p.184 of [30]
for M-scale estimators and derive similar results for $\tau$-scale estimators.
According to part 1 of lemma 5.1 in [30],

$\begin{split}\mbox{if
}\\#(C)>b\delta\Rightarrow\hat{\sigma}_{\text{M}}(\mathbf{r}^{(k)})\rightarrow\infty,\\\
\bar{c}_{i}\times\infty\leq\hat{\sigma}_{\tau}(\mathbf{r}^{(k)})\leq\sqrt{\sup_{t\in\mathbb{R}}\rho_{1}(t)}\times\infty,\end{split}$
(35)

which implies both upper and lower bound of $\tau$-scale estimate goes to
infinity. This proves part 1 of the lemma, if $\\#(C)>b\delta$ then
$\hat{\sigma}_{\tau}(\mathbf{r}^{(k)})\rightarrow\infty$. To prove the second
part of lemma, we know from part 2 of lemma 5.1 [30],

$\begin{split}\\#(D)>b-b\delta\Rightarrow\hat{\sigma}_{\text{M}}(\mathbf{r}^{(k)})\mbox{
is bounded},\\\
\bar{c}_{i}\times\mbox{bounded}\leq\hat{\sigma}_{\tau}(\mathbf{r}^{(k)})\leq\sqrt{\sup_{t\in\mathbb{R}}\rho_{1}(t)}\times\mbox{bounded},\end{split}$
(36)

which implies both lower and upper bounds of $\tau$-scale estimates are
bounded. This proves the second part of the lemma, if $\\#(D)>b-b\delta$, then
$\hat{\sigma}_{\tau}(\mathbf{r}^{(k)})$ is bounded. Now, we need to prove in
two separate stages that the finite-sample breakdown point of the $\tau$-Lasso
estimator is bounded from the above and the below.

###  B-E Bounded From Below

In order to establish the boundedness of $\hat{\boldsymbol{\beta}}$ from
below, we need to show that the sequence of $\tau$-scale estimates
$\big{(}\hat{\boldsymbol{\beta}}^{(k)}\in\mathbb{R}^{p}\big{)}_{k\in\mathbb{N}}$
is bounded for any arbitrary sequence of contaminated samples
$\big{(}\hat{\mathbf{Y}}_{m}^{(k)}\big{)}_{k\in\mathbb{N}}$ with $m\leq
m(\delta)$. The method of proof by contradiction can be used to establish the
boundedness from below. That is, it is assumed the sequence
$\big{(}\hat{\boldsymbol{\beta}}^{(k)}\big{)}_{k\in\mathbb{N}}$ is unbounded
and then shown that $\hat{\boldsymbol{\beta}}^{(k)}$ violates the optimality
condition where it is assumed to attain the minimum of the $\tau$-Lasso
objective function. Thus, $\hat{\boldsymbol{\beta}}^{(k)}$ should be bounded
to be a minimizer of the objective function.

Suppose $\tilde{\boldsymbol{\beta}}\in\mathbb{R}^{p}$ has a bounded
$\ell_{1}$-norm such that
$\|\tilde{\boldsymbol{\beta}}\|_{\ell_{1}}=K_{b}<\infty$. For the
uncontaminated observations $(\check{y}_{l},\check{\mathbf{x}}_{[l]}^{T})$
within the contaminated sample $\hat{\mathbf{Y}}_{m}^{(k)}$, the corresponding
residuals
$|r_{l}^{(k)}\big{(}\tilde{\boldsymbol{\beta}},\hat{\mathbf{Y}}_{m}^{(k)}\big{)}|=|\check{y}_{l}-\check{\mathbf{x}}_{[l]}^{T}\tilde{\boldsymbol{\beta}}|<\infty$
are bounded based on the triangle inequality. Without loss of generality, it
is assumed $b\delta=b\min(\delta,1-\delta)$. Since the inequality $m\leq
m(\delta)\leq b\delta$ holds based on the theorem assumption, it can be shown
the number of bounded residuals $\\#(D)\geq b-m\geq b-b\delta$. Therefore, we
can conclude
$\hat{\sigma}_{\tau}\Big{(}\mathbf{r}^{(k)}\big{(}\tilde{\boldsymbol{\beta}},\hat{\mathbf{Y}}_{m}^{(k)}\big{)}\Big{)}$
can be large but still bounded based on the above lemma,

$\sup_{k\in\mathbb{N}}\hat{\sigma}_{\tau}\Big{(}\mathbf{r}^{(k)}\big{(}\tilde{\boldsymbol{\beta}},\hat{\mathbf{Y}}_{m}^{(k)}\big{)}\Big{)}<\infty.$
(37)

Now, let the sequence
$\big{(}\|\hat{\boldsymbol{\beta}}^{(k)}\|_{\ell_{1}}\big{)}_{k\in\mathbb{N}}$
be unbounded. An unbounded sequence does not converge, i.e. sequences that
contain arbitrarily large numbers. Hence, there exists a sequence index
$k_{0}$ such that
$\|\hat{\boldsymbol{\beta}}^{(k_{0})}\|_{\ell_{1}}>K_{b}+\frac{1}{\lambda}\sup_{k\in\mathbb{N}}\hat{\sigma}_{\tau}^{2}\Big{(}\mathbf{r}^{(k)}\big{(}\tilde{\boldsymbol{\beta}},\hat{\mathbf{Y}}_{m}^{(k)}\big{)}\Big{)}$
and as a result we can say for every $k^{{}^{\prime}}\geq k_{0}$,

$\begin{split}\overbrace{\hat{\sigma}_{\tau}^{2}\Big{(}\mathbf{r}^{(k^{{}^{\prime}})}\big{(}\hat{\boldsymbol{\beta}}^{(k^{{}^{\prime}})},\hat{\mathbf{Y}}_{m}^{(k^{{}^{\prime}})}\big{)}\Big{)}+\lambda\|\hat{\boldsymbol{\beta}}^{(k^{{}^{\prime}})}\|_{\ell_{1}}}^{\mathcal{L}\big{(}\hat{\boldsymbol{\beta}}^{(k^{{}^{\prime}})},\hat{\mathbf{Y}}_{m}^{(k^{{}^{\prime}})}\big{)}}\\\
>\hat{\sigma}_{\tau}^{2}\Big{(}\mathbf{r}^{(k^{{}^{\prime}})}\big{(}\hat{\boldsymbol{\beta}}^{(k^{{}^{\prime}})},\hat{\mathbf{Y}}_{m}^{(k^{{}^{\prime}})}\big{)}\Big{)}\\\
+\lambda\Bigg{(}K_{b}+\frac{1}{\lambda}\sup_{k\in\mathbb{N}}\hat{\sigma}_{\tau}^{2}\Big{(}\mathbf{r}^{(k)}\big{(}\tilde{\boldsymbol{\beta}},\hat{\mathbf{Y}}_{m}^{(k)}\big{)}\Big{)}\Bigg{)}\\\
\geq\hat{\sigma}_{\tau}^{2}\Big{(}\mathbf{r}^{(k^{{}^{\prime}})}\big{(}\tilde{\boldsymbol{\beta}},\hat{\mathbf{Y}}_{m}^{(k^{{}^{\prime}})}\big{)}\Big{)}+\lambda
K_{b}=\mathcal{L}\big{(}\tilde{\boldsymbol{\beta}},\hat{\mathbf{Y}}_{m}^{(k^{{}^{\prime}})}\big{)}\Rightarrow\\\
\mathcal{L}\big{(}\hat{\boldsymbol{\beta}}^{(k^{{}^{\prime}})},\hat{\mathbf{Y}}_{m}^{(k^{{}^{\prime}})}\big{)}\geq\mathcal{L}\big{(}\tilde{\boldsymbol{\beta}},\hat{\mathbf{Y}}_{m}^{(k^{{}^{\prime}})}\big{)}\end{split}$
(38)

where $\lambda$ determines the amount of regularization imposed by the
$\ell_{1}$-norm and $\mathcal{L}(\cdot)$ is $\tau$-Lasso objective function as
given in equation (38). The above result contradicts the fact that
$\hat{\boldsymbol{\beta}}^{(k^{{}^{\prime}})}$ is the minimizer of the
objective function, as the loss function is larger for
$\hat{\boldsymbol{\beta}}^{(k^{{}^{\prime}})}$. Thus,
$\hat{\boldsymbol{\beta}}^{(k)}$ should be bounded for $m\leq m(\delta)$ and
the boundedness from below is proved.

###  B-F Bounded From Above

The boundedness from above can be established by showing that the estimator
breaks down for $m>b\delta$. In order to prove such a property, proof by
contradiction is used in this work. That is, it is assumed that the sequence
of estimates $\big{(}\hat{\boldsymbol{\beta}}^{(k)}\big{)}_{k\in\mathbb{N}}$
minimizing the $\tau$-Lasso objective function over the contaminated sample
$\hat{\mathbf{Y}}_{m}^{(k)}$ is bounded. and then shown the $\tau$-Lasso
objective function evaluated at $\hat{\boldsymbol{\beta}}^{(k)}$ achieves
larger value than the $\tau$-Lasso objective function evaluated at a given
unbounded $\tilde{\boldsymbol{\beta}}^{(k)}$. This contradicts the assumption
that $\hat{\boldsymbol{\beta}}^{(k)}$ is the minimum of $\tau$-Lasso objective
function. Thus, it can be concluded $\hat{\boldsymbol{\beta}}^{(k)}$ must be
unbounded for $m>b\delta$ establishing the boundedness from above.

As the $\tau$-Lasso objective function is comprised of three primary
components, its evaluation is carried out in three separate steps as follows:

  1. 1.

M-scale of residuals, component 1

  2. 2.

$\frac{1}{b}\sum_{l=1}^{b}\rho_{1}\Bigg{(}\frac{\hat{r}_{l}^{(k)}}{\hat{\sigma}_{\text{M}}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}\Bigg{)}$,
component 2

  3. 3.

$\lambda\|\hat{\boldsymbol{\beta}}^{(k)}\|_{\ell_{1}}$, component 3

Once these components are evaluated, the $\tau$-Lasso objective function is
formed by a simple addition and multiplication operation over these
components. In order to proceed with the proof, the $\tau$-Lasso objective
function is evaluated for the bounded sequences
$\hat{\boldsymbol{\beta}}^{(k)}$ and unbounded sequences
$\tilde{\boldsymbol{\beta}}^{(k)}$ and then compared with each other to arrive
at a contradiction.

####  B-F1 Evaluation of $\tau$-Lasso Objective Function for
$\hat{\boldsymbol{\beta}}^{(k)}$

##### Component 1

Suppose the set $C\subset\\{1,\cdots,b\\}$ denotes the indices of the
observations within the original contaminated-free sample $\check{\mathbf{Y}}$
replaced by outliers to construct the contaminated sample
$\hat{\mathbf{Y}}_{m}^{(k)}$ with $m=\\#(C)$. To simplify the proof without
loss of the generality, we choose an arbitrary
$\mathbf{x}_{0}\in\mathbb{R}^{p}$ with unit $\ell_{2}$-norm,
$\|\mathbf{x}_{0}\|_{\ell_{2}}=1$. Now, we construct a contaminated sequence
of samples with $m$ outliers given by

$\big{(}\hat{y}_{l}^{(k)},\hat{\mathbf{x}}_{[l]}^{(k)}\big{)}=\begin{cases}\big{(}k^{\nu+1},\mathbf{x}_{0}k\big{)}&l\in
C\\\ \big{(}\check{y}_{l},\check{\mathbf{x}}_{[l]}\big{)}&l\notin
C\end{cases},$ (39)

where $0<\nu\leq 1$ and $\big{(}k^{\nu+1},\mathbf{x}_{0}k\big{)}$ are chosen
to account for outlying observations. The sequence for outlying observations
diverges as $k$ goes to infinity.

First, we assume that $\hat{\boldsymbol{\beta}}^{(k)}$ is bounded in norm and
consequently, we have
$|r_{l}^{(k)}\big{(}\hat{\boldsymbol{\beta}}^{k},\hat{\mathbf{Y}}_{m}^{(k)}\big{)}|=|\hat{y}_{l}-(\hat{\mathbf{x}}_{[l]}^{(k)})^{T}\hat{\boldsymbol{\beta}}^{(k)}|<\infty$
is bounded for $l\notin C$ and $k\in\mathbb{N}$. $\hat{r}_{l}^{(k)}$ is a
shorthad for
$r_{l}^{(k)}\big{(}\hat{\boldsymbol{\beta}}^{k},\hat{\mathbf{Y}}_{m}^{(k)}\big{)}$.
On the other hand, the residuals corresponding to contaminated observations,
$l\in C$ are lower bounded by

$\begin{split}|\hat{r}_{l}^{(k)}|=|k^{\nu+1}-k\mathbf{x}_{0}^{T}\hat{\boldsymbol{\beta}}^{(k)}|=k|k^{\nu}-\mathbf{x}_{0}^{T}\hat{\boldsymbol{\beta}}^{(k)}|\mbox{
and }\mathbf{x}_{0}^{T}\hat{\boldsymbol{\beta}}^{(k)}\\\
\leq\|\mathbf{x}_{0}\|_{\ell_{1}}\|\hat{\boldsymbol{\beta}}^{(k)}\|_{\ell_{1}}\Rightarrow|k^{\nu}-\mathbf{x}_{0}^{T}\hat{\boldsymbol{\beta}}^{(k)}|\geq|k^{\nu}-\|\mathbf{x}_{0}\|_{\ell_{1}}\|\hat{\boldsymbol{\beta}}^{(k)}\|_{\ell_{1}}|\\\
\Rightarrow|\hat{r}_{l}^{(k)}|\geq
k|k^{\nu}-\|\mathbf{x}_{0}\|_{\ell_{1}}\|\hat{\boldsymbol{\beta}}^{(k)}\|_{\ell_{1}}|\end{split}$
(40)

In addition, the right-hand side of inequality goes to infinity as $k$
approaches infinity which implies the residuals $\hat{r}_{l}^{(k)}$ go to
infinity for $l\in C$ as well. Based on the above lemma and lemma 5.1 in p.184
of [30] , we conclude that both
$\hat{\sigma}_{\tau}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}$ and
$\hat{\sigma}_{\text{M}}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}$ go to infinity
for $\\#(C)>b\delta$. As a result, we can decompose the M-estimation of scale
equation as follows:

$\sum_{l\notin
C}\rho_{0}\Bigg{(}\frac{\hat{r}_{l}^{(k)}}{\hat{\sigma}_{\text{M}}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}\Bigg{)}+\sum_{l\in
C}\rho_{0}\Bigg{(}\frac{\hat{r}_{l}^{(k)}}{\hat{\sigma}_{\text{M}}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}\Bigg{)}=b\delta$
(41)

Recalling the proof of Theorem 4.1 in [53], it follows

$\rho_{0}\big{(}\frac{1}{\gamma}\big{)}=\frac{b\delta}{m}$ (42)

where

$\gamma=\lim_{k\rightarrow\infty}\frac{\hat{\sigma}_{\text{M}}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}{k^{\nu+1}}$
(43)

##### Component 2

On the other hand, we have

$\begin{split}\frac{1}{b}\sum_{l=1}^{b}\rho_{1}\Bigg{(}\frac{\hat{r}_{l}^{(k)}}{\hat{\sigma}_{\text{M}}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}\Bigg{)}=\frac{m}{b}\rho_{1}\Bigg{(}\frac{1-\mathbf{x}_{0}^{T}\hat{\boldsymbol{\beta}}^{(k)}/k^{\nu}}{\hat{\sigma}_{\text{M}}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}/k^{\nu+1}}\Bigg{)}\Rightarrow\\\
\lim_{k\rightarrow\infty}\frac{m}{b}\rho_{1}\Bigg{(}\frac{1-\mathbf{x}_{0}^{T}\hat{\boldsymbol{\beta}}^{(k)}/k^{\nu}}{\hat{\sigma}_{\text{M}}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}/k^{\nu+1}}\Bigg{)}=\frac{m}{b}\rho_{1}\Big{(}\frac{1}{\gamma}\Big{)}\end{split}$
(44)

##### Component 3

$\lim_{k\rightarrow\infty}\|\hat{\boldsymbol{\beta}}^{(k)}\|_{\ell_{1}}$
remains bounded according to the assumption of the proof.

##### Deriving $\tau$-Lasso Objective Function

Now, we can evaluate the $\tau$-Lasso loss function according to equation (38)
as follows:

$\begin{split}\lim_{k\rightarrow\infty}\frac{\hat{\sigma}_{\tau}^{2}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}{k^{2\nu+2}}=\lim_{k\rightarrow\infty}\frac{\hat{\sigma}_{\text{M}}^{2}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}{k^{2\nu+2}}\times\frac{m}{b}\rho_{1}\Bigg{(}\frac{1-\mathbf{x}_{0}^{T}\hat{\boldsymbol{\beta}}^{(k)}/k^{\nu}}{\hat{\sigma}_{\text{M}}^{2}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}/k^{\nu+1}}\Bigg{)}\\\
\Rightarrow\lim_{k\rightarrow\infty}\frac{\hat{\sigma}_{\tau}^{2}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}{k^{2\nu+2}}=\frac{m\gamma^{2}}{b}\rho_{1}\Big{(}\frac{1}{\gamma}\Big{)}\Rightarrow\lim_{k\rightarrow\infty}\frac{\mathcal{L}\big{(}\hat{\boldsymbol{\beta}}^{(k)},\hat{\mathbf{Y}}_{m}^{(k)}\big{)}}{k^{2\nu+2}}\\\
=\lim_{k\rightarrow\infty}\Bigg{(}\frac{\hat{\sigma}_{\tau}^{2}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}{k^{2\nu+2}}+\lambda\frac{\|\hat{\boldsymbol{\beta}}^{(k)}\|_{\ell_{1}}}{k^{2\nu+2}}\Bigg{)}=\frac{m\gamma^{2}}{b}\rho_{1}\Big{(}\frac{1}{\gamma}\Big{)}\end{split}$
(45)

where
$\lim_{k\rightarrow\infty}\|\hat{\boldsymbol{\beta}}^{(k)}\|_{\ell_{1}}/k^{2\nu+2}=0$
due to the bounded norm of $\hat{\boldsymbol{\beta}}^{(k)}$.

####  B-F2 Evaluation of $\tau$-Lasso Objective Function for
$\tilde{\boldsymbol{\beta}}^{(k)}$

Now, we evaluate the $\tau$-Lasso objective function for an unbounded sequence
$\tilde{\boldsymbol{\beta}}^{(k)}$ in a three-step procedure as follows:

##### Component 1

The unbounded sequence is assumed to be
$\tilde{\boldsymbol{\beta}}^{(k)}=\frac{k^{\nu}}{2}\mathbf{x}_{0}$. The
residuals for this sequence become

$\hat{r}_{l}^{(k)}=\begin{cases}k^{\nu+1}-\frac{k^{\nu+1}}{2}\mathbf{x}_{0}^{T}\mathbf{x}_{0}&l\in
C\\\
\check{y}_{l}-\frac{k^{\nu}}{2}\mathbf{x}_{0}^{T}\check{\mathbf{x}}_{[l]}&l\notin
C\end{cases}=\begin{cases}\frac{k^{\nu+1}}{2}&l\in C\\\
\check{y}_{l}-\frac{k^{\nu}}{2}\mathbf{x}_{0}^{T}\check{\mathbf{x}}_{[l]}&l\notin
C\end{cases}$ (46)

where $\mathbf{x}_{0}^{T}\mathbf{x}_{0}=\|\mathbf{x}_{0}\|_{\ell_{2}}^{2}=1$
and all residuals go to infinity as $k\rightarrow\infty$. Hence, both
$\hat{\sigma}_{\tau}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}$ and
$\hat{\sigma}_{\text{M}}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}$ tend to
infinity. The decomposition of M-estimation of scale equation yields,

$\sum_{l\notin
C}\rho_{0}\Bigg{(}\frac{\check{y}_{l}-\frac{k^{\nu}}{2}\mathbf{x}_{0}^{T}\check{\mathbf{x}}_{[l]}}{\hat{\sigma}_{\text{M}}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}\Bigg{)}+\sum_{l\in
C}\rho_{0}\Bigg{(}\frac{k^{\nu+1}/2}{\hat{\sigma}_{\text{M}}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}\Bigg{)}=b\delta$
(47)

Using the proof of Theorem 4.1 in [53] , it can be inferred

$\rho_{0}\Bigg{(}\frac{1}{\lim_{k\rightarrow\infty}\frac{\hat{\sigma}_{\text{M}}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}{k^{\nu+1}/2}}\Bigg{)}=\frac{b\delta}{m}$
(48)

where

$\lim_{k\rightarrow\infty}\frac{\hat{\sigma}_{\text{M}}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}{k^{\nu+1}/2}=\gamma$
(49)

##### Component 2

On the other hand, we have

$\begin{split}\frac{1}{b}\sum_{l=1}^{b}\rho_{1}\Bigg{(}\frac{\hat{r}^{(k)}_{l}}{\hat{\sigma}_{\text{M}}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}\Bigg{)}=\frac{m}{b}\rho_{1}\Bigg{(}\frac{1}{\frac{\hat{\sigma}_{\text{M}}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}{k^{\nu+1}/2}}\Bigg{)}\Rightarrow\\\
\lim_{k\rightarrow\infty}\frac{m}{b}\rho_{1}\Bigg{(}\frac{1}{\frac{\hat{\sigma}_{\text{M}}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}{k^{\nu+1}/2}}\Bigg{)}=\frac{m}{b}\rho_{1}\Big{(}\frac{1}{\gamma}\Big{)}\end{split}$
(50)

##### Component 3

$\lim_{k\rightarrow\infty}\|\tilde{\boldsymbol{\beta}}^{(k)}\|_{\ell_{1}}$
diverges as $k$ goes to infinity.

##### Deriving $\tau$-Lasso objective function

Now, we can evaluate the $\tau$-Lasso loss function according to equation (38)
as follows:

$\begin{split}\lim_{k\rightarrow\infty}\frac{\hat{\sigma}_{\tau}^{2}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}{k^{2\nu+2}}=\lim_{k\rightarrow\infty}\frac{\hat{\sigma}_{\text{M}}^{2}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}{k^{2\nu+2}}\times\frac{m}{b}\rho_{1}\Bigg{(}\frac{1}{\frac{\hat{\sigma}_{\text{M}}\big{(}\hat{\mathbf{r}}^{(k)}\big{)}}{k^{\nu+1}/2}}\Bigg{)}\\\