IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, XXX XXXX

1

A New Sparse and Robust Adaptive Lasso
Estimator for the Independent Contamination Model

arXiv:1705.02162v1 [] 5 May 2017

Jasin Machkour, Michael Muma, Member, IEEE, Bastian Alt, Abdelhak M. Zoubir, Fellow, IEEE

Abstract—Many problems in signal processing require finding
sparse solutions to under-determined, or ill-conditioned, linear
systems of equations. When dealing with real-world data, the
presence of outliers and impulsive noise must also be accounted
for. In past decades, the vast majority of robust linear regression
estimators has focused on robustness against rowwise contamination. Even so called ‘high breakdown’ estimators rely on the
assumption that a majority of rows of the regression matrix is
not affected by outliers. Only very recently, the first cellwise
robust regression estimation methods have been developed. In
this paper, we define robust oracle properties, which an estimator
must have in order to perform robust model selection for
under-determined, or ill-conditioned linear regression models
that are contaminated by cellwise outliers in the regression
matrix. We propose and analyze a robustly weighted and adaptive
Lasso type regularization term which takes into account cellwise
outliers for model selection. The proposed regularization term is
integrated into the objective function of the MM-estimator, which
yields the proposed MM-Robust Weighted Adaptive Lasso (MMRWAL), for which we prove that at least the weak robust oracle
properties hold. A performance comparison to existing robust
Lasso estimators is provided using Monte Carlo experiments.
Further, the MM-RWAL is applied to determine the temporal
releases of the European Tracer Experiment (ETEX) at the source
location. This ill-conditioned linear inverse problem contains
cellwise and rowwise outliers and is sparse both in the regression
matrix and the parameter vector. The proposed RWAL penalty
is not limited to the MM-estimator but can easily be integrated
into the objective function of other robust estimators.
Index Terms—Sparse and Robust Estimation, Outlier, Lasso,
Independent Contamination Model, Robust Oracle Properties,
Atmospheric Emissions.

I. I NTRODUCTION
Many of today’s signal processing problems can be formulated as a linear regression
y = Xβ + u,

(1)

where we assume that the regressor matrix X ∈ Rn×p , the
errors u ∈ Rn×1 and observations y ∈ Rn×1 are independent
and identically distributed (iid) random variables, β ∈ Rp×1
are the unknown parameters of interest, and X and u are
mutually independent.
The presence of outliers and impulsive noise has been
reported in applications as diverse as wireless communication,
ultrasonic systems, computer vision, electric power systems,
automated detections of defects, biomedical signal analysis,
genomics and the estimation of the temporal releases of
J. Machkour, M. Muma and A.M. Zoubir are with the Signal Processing
Group Technische Universität Darmstadt. Merckstr. 25, 64283 Darmstadt,
Germany. Email: {machkour,muma,zoubir}@spg.tu-darmstadt.de. B. Alt is
with the Bioinspired Communication Systems, Rundeturmstr. 12, 64283
Darmstadt, Germany. Email: bastian.alt@bcs.tu-darmstadt.de

a pollutant to the atmosphere. See [1], [2], [3], [4], [5],
[6], [7], [8], [9] and references therein. Violations of the
Gaussian assumption cause a drastic performance drop for the
commonly used least-squares estimator (LSE) [10], [11], [12]
β̂LSE = arg min ky − xβk22 .

(2)

β

For decades, the vast majority of robust linear regression
estimators has focused on robustness against ’rowwise’ contamination. Under this so-called Tukey-Huber contamination
model (THCM) [11], a small fraction of rows of X may be
contaminated. Even ’high-breakdown’ regression estimators,
such as the LTS-, S-, MM-, and τ -estimators [1], [11] rely
on the THCM. In [13], Rousseuw and Van den Bossche
state that recently researchers have come to realize that the
outlying rows paradigm is no longer sufficient for modern
high-dimensional data sets. It often happens that most data
cells (entries) in a row are regular and just a few of them are
anomalous.
The case that independent cells of X are outliers is referred
to as the independent contamination model (ICM) [14], [15],
[16]. Only very recently, the first ‘cellwise robust’ regression
estimation methods have been developed [15], [16]. The extension of existing estimators to other contamination models,
such as the ICM, and even the development of completely
new robust estimators is necessary to solve many real-world
problems. For example, the estimation of the spatio-temporal
emissions of a pollutant, given noisy observations, can be
formulated as a linear inverse problem with the help of an
atmospheric dispersion model [6]. The data of the European
Tracer Experiment (ETEX) which was conducted in Monterfil,
Brittany in 1994, where Perfluorocarbon (PFC) tracers were
released into the atmosphere, for instance, contains both
cellwise and rowwise outliers.
Additionally to the robustness considerations, atmospheric
inverse problems, like many other problems in signal processing, require finding sparse solutions to under-determined, or
ill-conditioned, linear systems of equations. For example, handling large datasets in terms of model interpretation, including
the case where the number of explanatory variables p is larger
than the sample size n, requires penalized estimators, such as
the classical least absolute shrinkage and selection operator
(Lasso) [17]
β̂Lasso = arg min ky − xβk22 + λkβk1 .

(3)

β

with λ ∈ R+ .
Many other regularizations have been proposed [18], [19],
[20], [21], [22]. In this paper the focus lies on Lasso estimation, to select a robust and interpretable model in high

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, XXX XXXX

dimensional settings. Zou [23] showed that the Lasso variable
selection can be inconsistent, so that the oracle properties do
not hold and proposed the adaptive Lasso
ad

β̂ Lasso = arg min ky − xβk22 + λ
β

p
X
j=1

ŵj |βj |,

(4)

where ŵj = 1/|β̂j |γ (γ > 0) are non-negative weights
√
depending on β̂, which is a n-consistent estimator of β.
Just like the LSE, the Lasso and the adaptive Lasso rely on
the Gaussian noise assumption and are sensitive to outliers.
In recent years, some robust and regularized approaches have
been proposed that replace the penalized square objective
function by a penalized bounded objective function [7], [9],
[24], [25], [26]. These methods, however, again, rely on the
THCM, and to date, no penalized robust regression method
exists that can handle cellwise and rowwise outliers.
Original Contributions: First, we give a weak and strong
definition of what we call the ’robust oracle properties’. These
are properties that estimators aiming at performing robust variable selection need to have. Next, we propose and analyze a
robustly weighted and adaptive Lasso-type regularization term,
which takes into account cellwise outliers for model selection.
The proposed regularization term is integrated into the objective function of the MM-estimator, which yields the proposed
MM-Robust Weighted Adaptive Lasso (MM-RWAL), for which
we prove that at least the weak robust oracle properties hold.
We would like to highlight, that the proposed Robust Weighted
Adaptive Lasso penalty can easily be integrated into the
objective function of other robust estimators. A performance
comparison to existing robust Lasso estimators is provided
using Monte Carlo experiments. Further, a challenging realdata application of estimating the sparse non-negative spatiotemporal emissions of a pollutant is considered, given noisy
observations y and an imprecisely estimated ill-conditioned
and sparse dispersion model X. This example contains both
cellwise and rowwise outliers.
Notation: Scalars are denoted by lowercase letters, e.g., x,
column vectors by bold-faced lowercase letters, e.g. x, matrices by bold-faced uppercase letters, e.g. X, sets are denoted
by calligraphic letters, e.g. X with associated cardinality |X |.
The jth column of a matrix X is denoted by xj while (x)i:j
denotes the vector that contains the entries i to j of vector x.
The ith element of vector x is denoted by xi , Ip is the
p-dimensional identity-matrix, 0p is the p-dimensional allzeros vector and diag(x) forms a matrix that contains the
entries of x as its diagonal. β̂ refers to the estimator (or
estimate) of the parameter vector β, (·)⊤ is the transpose
operator. The derivative of a function f with respect to its
argument is abbreviated by f ′ . P (X) is the probability of
event X. Bin(1, ǫ) denotes the binomial distribution with one
trial and a success probability of ǫ. Convergence to the normal
distribution with mean vector µ and covariance matrix Σ is
d
denoted by → N (0, Σ).
Organization: Section II discusses the Tukey-Huber and Independent Contamination models and motivates the use of
cellwise robust methods. Section III defines the robust oracle
properties, introduces the proposed estimator and provides

2

algorithms to compute the estimates. Section IV provides
numerical experiments, while Section V contains a challenging
real-data application of source estimation for an atmospheric
inverse problem. Finally, Section VI concludes the paper with
a brief outlook on future work.

II. T UKEY-H UBER AND I NDEPENDENT C ONTAMINATION
M ODEL
A. Tukey-Huber Contamination Model (THCM)
The THCM is based on the assumption that a majority (1 − ǫ) of the data points is not contaminated, while a
minority ǫ is contaminated. The univariate THCM is formulated as
x = (1 − b)y + bz,

(5)

where b, y, z are mutually independent scalar random variables. Let F be the true distribution of the data and let G be an
unspecified contaminating distribution, from which the outliers
are generated. Then, with y ∼ F , z ∼ G and b ∼ Bin(1, ǫ), the
distribution of the observed variable x becomes
H = (1 − ǫ)F + ǫG.

(6)

The multivariate THCM is defined by
x = (Ip − bIp )y + bIp z,

(7)

where y, x, z are p-dimensional random variables.
A highly valuable property of the THCM is that the percentage of contaminated rows in the data-matrix stays unchanged
under affine transformations, that is, if the random vector x
follows the THCM, then the affine transformed vector
x̃ = Ax + b

(8)

also follows the THCM. Thus, estimators designed for the
THCM can be affine equivariant. Additionally, many important robustness concepts such as the influence function and
the breakdown point are based on the THCM [10], [11],
[12], [1]. However, the disadvantages of this contamination
model, which occur especially in high dimensional settings,
are alerting. First of all, the assumption that a major fraction
(1 − ǫ) of the data points is outlier free, very unlikely holds
in higher dimensions. A mathematical legitimation of this
criticism will be given in Eq. (14) and the associated Fig. 2.
Secondly, as illustrated in Fig. 1a we see that a few cellwise
outliers in the THCM lead to flagging the whole corresponding
row of X as outliers. Although much information is lost,
such contamination may be handled by high breakdown point
estimators in low dimensions. However, Fig. 1b illustrates that
in cases where the number of predictors p exceeds the number
of rows n, it becomes more and more likely that a few highly
contaminated predictors force the THCM-based estimators to
flag all data points as outliers, which makes it impossible to
draw any inferences from the data.

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, XXX XXXX

3

of a row in a p-dimensional dataset being contaminated by the
formula
Pcont,row = 1 − (1 − ǫ)p .
(11)
So, for any high breakdown point estimator, tuned to have
the highest possible breakdown point of 50%, we obtain the
inequality
Pcont,row < 0.5
⇔ 1 − (1 − ǫ)p < 0.5
⇔

This means that the tolerable fraction of contamination ǫ in
every predictor - for simplicity, we assume that ǫ is equal for
every predictor - is bounded and depends on the dimension of
the dataset. From this point of view, the probability of having
only THCM-outliers in the data depends on the dimension of
the predictor-matrix and the number of rows as follows:

(a) p < n.

P (”THCM
ICM”)
 
 


−1
n
n
n
= 1·
·
· ...·
⌈nǫ2 ⌉
⌈nǫ3 ⌉
⌈nǫp ⌉

n!
n!
·
=
(n − nǫ2 )! · (nǫ2 )! (n − nǫ3 )! · (nǫ2 )!
−1
n!
... ·
(n − nǫ2 )! · (nǫp )!

(b) p > n.
Fig. 1. Illustration of how the THCM downweights outliers in the case of
(a) a few cellwise outliers for p < n and (b) a few highly contaminated
predictors for p < n.

B. Independent Contamination Model (ICM)
The ICM is defined by
x = (Ip − B)y + Bz,

(9)

where B = diag(b1 , b2 , . . . , bp ) and b1 , b2 , . . . , bp are independent Bernoulli random variables with success probability
ǫj , j = 1, . . . , p. Loosely speaking, each cell xij corresponding to the predictor j in every single row has a probability ǫj
of being contaminated. Notice that for P (b1 = b2 = . . . =
bp ) = 1, the ICM reduces to the THCM.
The main issue with the ICM is that it is not equivariant
under affine transformations. Let x be a random vector and
A an invertible quadratic p-dimensional matrix. Then, for an
affine transformation of x
x̃ = Ax + b = A(Ip − B)y + ABz + b

6= (Ip − B)Ay + BAz + b.

(12)

ǫ < 1 − 0.5p .

(10)

Therefore, if AB 6= BA, x̃ does not follow the ICM. The
lack of affine equivariance has a far reaching consequence
for the ICM, which is referred to as ’outlier propagation’.
Outlier propagation means that an outlying cell in a predictor
may spread over other components of the corresponding data
point, e.g. by linearly combining the predictors in a regression
model. From these considerations, we calculate the probability

(13)

Since the THCM assumes that ǫ1 = . . . = ǫp = ǫ, it follows
that
P (”THCM
ICM”)
−(p−1)

n!
=
(n − nǫ)! · (nǫ)!

−(p−1)
n
=
.
⌈nǫ⌉

(14)

Fig. 2 illustrates how rapidly the probability of having only
THCM-outliers decreases for fixed n and ǫ. This, again
supports the statement made in [13] that the outlying rows
paradigm is no longer suffcient for modern high-dimensional
data sets.
III. P ROPOSED M ETHODS
A. Proposed Definitions of Robust Oracle Properties
Extending the ideas of [23] to the ICM, we next introduce
what we call the ’robust oracle properties’. We propose a
strong and a weak version of the robust oracle properties.
Definition III.1. (Weak Robust Oracle Properties)
Let {x1 , x2 , . . . , xp } be the set of predictors,
A := {j : βj 6= 0 ∧ ǫj < 0.5}

(15)

the set of indices corresponding to the set of robust active
predictors, and A∗n := {j : β̂j∗n 6= 0} the set of indices
corresponding to the set of predictors chosen by a Lasso type
estimator to be active. Then, a Lasso type estimator needs to
have the following properties to be a robust oracle estimator:

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, XXX XXXX

1

n = 25,
n = 25,
n = 25,
n = 50,
n = 50,
n = 50,
n = 300,
n = 300,
n = 300,

10 -10
10 -20

= 0.1
= 0.2
= 0.3
= 0.1
= 0.2
= 0.3
= 0.1
= 0.2
= 0.3

10 -30

4

Intuitively speaking, the strong robust oracle properties hold
for any Lasso type estimator, whose adaptive ℓ1 -penalty term
ensures the penalization of the predictors in an ascending order
given by the predictor contamination order statistics
(ǫ)1:p ≤ (ǫ)2:p ≤ . . . ≤ (ǫ)p:p ,

(17)

while choosing the tuning parameter λ in the ℓ1 -penalty such
that predictors enter the model until the breakdown point of
the estimator is reached.

10 -40
B. MM-Robust Weighted Adaptive Lasso (MM-RWAL)

10 -50
1

5

10

15 20

Fig. 2. Probability of having only THCM-outliers for different values of n
and ǫ as a function of p.

1) Consistency in variable selection:
lim P (A∗n = A) = 1.
n→∞

2) Asymptotic normality:
∗n
√
d
n(β̂ A − β A ) → N (0, Σ∗ ), where Σ∗ is the covariance matrix, knowing the true subset model.
Loosely speaking, a Lasso type estimator is said to be a
robust oracle estimator if it selects only the active variables,
while especially leaving out highly outlier-contaminated ones.
We call a predictor highly contaminated if at least 50% of
its entries are outliers. The weak robust oracle properties
are necessary for any robust estimator, because a single
highly contaminated variable might introduce outliers to most
observations and lead to a breakdown of the estimator. Note
that we differentiate between the active set and the robust
active set of variables, since active variables should not be
chosen, when they are highly contaminated.

In this section, we introduce a new method called the MMRobust Weighted Adaptive Lasso (MM-RWAL). Recently,
the MM-Lasso and adaptive MM-Lasso were introduced to
robustify against outliers [26]. The objective function of the
MM estimator is
!
n
X
ri (β)
.
(18)
ρ
β̂ MM = arg min
sn (r(β̂ 1 ))
β
i=1
Here, ρ(·) is a robustifying function (see, e.g. [10], [11], [12]),
r(β̂ 1 ) = y − xβ̂ 1 is the residual of an S-estimator whose
estimates β̂ 1 have the property of minimizing a robust Mscale sn (r(β)) that satisfies


n
1X
ri (β)
ρ
= b,
n i=1
sn

where b is usually chosen such that consistency under the
Gaussian distribution is obtained. For the MM-(adaptive)
Lasso, (18) is extended by the penalty terms of (3) and (4),
respectively [26].
The proposed MM-RWAL estimator minimizes an MM
objective function to which a robust adaptive ℓ1 -penalty term
is added:
!
p
n
X
X
RWAL
ri (β)
ŵj |βj |.
+ λn
ρ
β̂MM = arg min
sn (r(β̂ 1 ))
β
j=1
i=1
(19)
Here, the MM-Lasso estimator [26] is used to calculate the
Lasso
weights according to ŵj = 1/|zj · β̂jMM
|, with zj defined in
Definition III.2. (Strong Robust Oracle Properties)
(23).
Let {x1 , x2 , . . . , xp } be the set of predictors,
To robustify the variable selection of the adaptive MM
A := {j : βj 6= 0 ∧ 1−Πkj=1 (1−(ǫ)j:p ) < 0.5, k ∈ {1, . . . , p}} Lasso, we propose to incorporate a measure of outlyingness for
(16) each predictor. We use the Stahel Donoho Outlyingness (SDO)
the set of indices corresponding to the set of robust active [27] and adjust it in a similar vein to the Adjusted Stahel
p
predictors, and A∗n := {j : β̂j∗n 6= 0} the set of indices Donoho Outlyingness of [28]. Let B = {b1 , b2 , . . . , bn } ⊂ R
corresponding to the set of predictors chosen by a Lasso type be a set of n observations. Then, the robust Stahel Donoho
estimator to be active. Then, a Lasso type estimator needs to outlyingness is given by
have the following properties to be a robust oracle estimator:
|a⊤ bi − med(a⊤ x)|
r(bi , x) = sup
, i = 1, . . . , n, (20)
mad(a⊤ x)
a∈Sp
1) Consistency in variable selection:
where Sp = {a ∈ Rp : kak2 = 1} and med(·) and mad(·)
lim P (A∗n = A) = 1.
n→∞
denote the median and the median absolute deviation (mad).
Since we assume that most rows flagged by the SDO as
2) Asymptotic normality:
√
∗n
outliers are not outlying in all of their components, the SDO
d
n(β̂ A − β A ) → N (0, Σ∗ ), where Σ∗ is the covariis extended by also taking into account the outlyingness of
ance matrix, knowing the true subset model.
the predictors. The idea that has been introduced in [28] in a

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, XXX XXXX

similar vein, is to adjust the SDO of every observation using
the outlyingness of every single predictor. This gives us the
Predictor Outlyingness (PO)
n
X
|xij − med(xj )|
,
cj =
mad(xj )
i=1

j = 1, . . . , p.

(21)

i = 1, . . . , n,

j = 1, . . . , p. (22)

We chose the tuning parameter α to be 0.5 throughout this
paper, to equally weigh the SDO and the PO, in order to
perform well in both contamination models, THCM and ICM.
Finally, by applying a weight function w(·), summing up the
rows of the outlyingness matrix and dividing by its cell sum,
we obtain the weights
n

zj =

X
p
w(rij ),
p
n
P P
i=1
w(rij )

j = 1, . . . , p,

(23)

i=1 j=1

where

designed outlyingness-matrix. Then, the weight of the jth
predictor is given by
n

zj =

Combining both, the SDO and the PO, we introduce an
outlyingness-matrix, whose (i, j)-th element is
rij = αri + (1 − α)cj ,

5

(25)

i=1 j=1

 2 

P
c
p · n(1 − ǫj ) +
r
ij
i∈Cj
=
 2 .
P
P P
c
pn − n
ǫl +
r
il
l∈E
l∈E i∈Cl

(26)

Now, we assume that the magnitude of outlyingness is fixed,
that is rij = r for all i, j. This yields

 2 

c
p · 1 − ǫj 1 −
r
.
zj =

 2 
P
c
ǫl
p− 1−
r
l∈E

(27)

Eq. (27) provides additional insights that help us understand
the behavior of our proposed weights. It shows us: the larger
ǫj the smaller zj . This leads to
Proposition III.3.
Let ǫj ∈ E, rij = r for all i, j and ǫ1 ≥ ǫ2 ≥ . . . ≥ ǫp . Then,
it follows that

Pp

j=1 zj = p.

In order to downweight cells, whose overall outlyingness
exceeds a certain threshold, we choose w(·) to be the Huber
weight function,
w(r) = 1(r≤c) + (c/r)2 1(r≤c) ,

X
p
w(rij )
p
n P
P
w(rij ) i=1

(24)

√
with c = min( χ2p (0.5), 4) as proposed in [28].

C. Analysis of the Proposed Predictor Weights
In this section, we analyze the behavior of the proposed
weights zj . For this purpose, we define two sets:
1. Cj = {i : w(rij ) < 1, i = 1, . . . , n} with cardinality
|Cj | = ⌈nǫj ⌉, is the set of indices corresponding to the
contaminated cells in predictor j.
2. E = {l : ǫl > 0, l = 1, . . . , p} with cardinality |E| =: γ
is the set of indices corresponding to the contaminated
predictors, so γ is the number of contaminated predictors.

z1 ≤ z2 ≤ . . . ≤ zp .
In this step, we assume that all contaminated predictors
contain the same fraction of contamination, which results in


 2 
c
p· 1−ǫ 1−ǫ
r
zj =
.
(28)

 2 
c
p− 1−
· γǫ
r
To carry out a plausibility analysis for the derived formulas,
we assume now that all p-predictors are contaminated (γ = p)
with the same fraction of outliers ǫ. Intuitively, we expect that
all predictors get the same weight and no one preferred over
an other one. Since our weights are designed to sum up to p,
we expect each predictor j ∈ {1, . . . , p} to receive the weight
zj = 1. Applying these assumptions to Eq. (28), we obtain

 2 
c
p − pǫ 1 −
r
zj =

 2  = 1,
c
p − pǫ 1 −
r

(29)

which confirms our expectations.

We will start with the most general formula that describes
the behavior of our weights and introduce step by step
assumptions, which simplify the equation.

The proof of the robust oracle properties in the next
section requires the weights zj to be smaller than one. Here,
we will prove this property only for Eq. (28).

Let ǫj be the fraction of contamination of the jth predictor,
and rij the magnitude of outlyingness of each cell in our

Proposition III.4.
If ǫj > 0, then zj ≤ 1 for all j.

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, XXX XXXX

Proof.
Let γ < p, then we obtain for fixed r and ǫ1 = ǫ2 = . . . =
ǫp = ǫ
 2 

c
p · 1 − ǫj
r
zj =

 2 
c
· γǫ
p− 1−
r

 2 
c
p − pǫ 1 −
r
=

 2 
c
p − γǫ 1 −
r

 2 
c
p − γǫ 1 −
r
<

 2  = 1.
c
p − γǫ 1 −
r

D. Proof of the Robust Oracle Properties for the Proposed
MM-RWAL
To prove that the MM-RWAL possesses at least the weak
robust oracle properties, we will need some assumptions. Let
G denote the distribution of the rows (xi1 , xi2 , . . . , xip ) in X
and let F be the distribution of the errors ui , which results in
the distribution H of the data points (xi1 , xi2 , . . . , xip , yi ) to
become
H(x, y) = G(x)F (y − x⊤ β).
(30)

6

we know that, with arbitrarily high probability, there exists a
constant M1 > 0 with

M1
kβ̂R − β0 k < √ .
n

(31)

Now, let
√
√ 

n
X
ri (β 0,I + u1 / n, β0,II + u2 / n)
ρ1
Gn (u1 , u2 ) =
sn (r(β̂ 1 ))
i=1
√
s
X
|β0,j + u1,j / n|
+ λn
|β̂2,j |
j=1
p
X |u2,j−s /√n|
.
+ λn
j=s+1 |zj · β̂2,j |
We obtain (β̂ R,I , β̂ R,II ) by minimizing Gn (u1 , u2 ), subject
to ku1 k + ku2 k2 ≤ M12 that we get from:
M1
kβ̂R − β 0 k ≤ √
n
⇔

u2
M1
u1
β 0,I + √ + β0,II + √ − β 0 ≤ √
n
n
n

A1. All occurring ρ(·)-functions are twice continuously
differentiable and Ψ = ρ′ .

⇔

u1 + u2
√
n

A2. The density f of the error terms u is an even and
monotonically decreasing function of |u|.

⇔ ku1 k2 + ku2 k2 ≤ M12 .

We now formulate our assumptions as follows:
M1
≤√
n

A3. The second moments of G exist.
√
RWAL
n-consistent
is a
A4. The estimator β̂ R := β̂MM
estimator of β R := βRWAL
MM , where β R complies with
Definition III.1.
A5. β̂ 2 is

√

n-consistent.

Note, that we will explicitely prove only the consistency in
variable selection here, because the weights zj do not affect
the asymptotic normality.

We next have to show hat Gn (u1 , u2 ) − Gn (u1 , 0p−s ) > 0
holds under the given condition and when kuk > 0. Let
Gn (u1 , u2 ) − Gn (u1 , 0p−s ) = D + E
with D and E being defined in Eq. (32).
With the Mean Value Theorem, we obtain

Proof. (Consistency in Variable Selection)
In this proof, we stick to the notation and structure of
Theorem 2 in [29] and Theorem 5 in [26].
With (A4.) and β0 = (β0,I , β0,II )⊤ being the true parameter
vector corresponding to Definition III.1, where I = A are
the s indices belonging to the robust active predictors and
II = {1, . . . , p} \ A = AC is the complementary set of A,

⊤

−1

D = (0s , u2 ) √
nsn (r(β̂ 1 ))

n
X
i=1

Ψ1



ri (θ∗n )
sn (r(β̂1 ))


xi ,

√
√
where θ∗n = (β 0,I +u1 / n, (1−αn )u2 / n) and αn ∈ [0, 1].

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, XXX XXXX

7

Gn (u1 , u2 ) − Gn (u1 , 0p−s ) =
√
√
√ 


p
n  
X
ri (β 0,I + u1 / n, u2 / n)
ri (β0,I + u1 / n, 0p−s )
λn X |u2,j−s |
ρ1
− ρ1
=
+√
n j=s+1 |zj · β̂2,j |
sn (r(β̂ 1 ))
sn (r(β̂ 1 ))
i=1
|
{z
} |
{z
}
=:D

=:E

Applying the Mean Value Theorem a second time yields


n
X
ri (θ∗n )
−1
Ψ1
xi
(0s , u2 )⊤ √
nsn (r(β̂ 1 )) i=1
sn (r(β̂1 ))
−1

⊤

= √
(0s , u2 )
nsn (r(β̂ 1 ))

+√

1
ns2n (r(β̂ 1 ))

n
X

Ψ1

i=1



ri (β 0 )
sn (r(β̂ 1 ))

≥ λn ·

(0s , u2 )⊤

Ψ′1

i=1



⊤

= √
(0s , u2 )
nsn (r(β̂ 1 ))

+

1
ns2n (r(β̂1 ))

n
X

Ψ1

i=1

(0s , u2 )⊤

n
X
i=1




xi


ri (θ∗∗
n )
xi x⊤
i
sn (r(β̂ 1 ))

ri (β 0 )
sn (r(β̂ 1 ))

Ψ′1




xi

ri (θ∗∗
n )
sn (r(β̂ 1 ))

E. Computation of the Weights for MM-RWAL

kθ∗n −β0 k.



|u2,j |
√
|z · nβ̂2,j |
{j:βj =0∧ǫj =0} j

|u2,j |
= λn ·
√
|z · nβ̂2,j |
{j:βj =0∧ǫj >0} j
X
+
X

|u2,j |
√
| nβ̂2,j |
{j:βj =0∧ǫj =0}

The main problem of calculating the weights zj is to compute the supremum in the SDO, because the cardinal number
of Sp is infinite and the objective function is non-convex.
Therefore, we need to apply a random search algorithm to
obtain an approximation of the supremum. We chose to take
a subsample from Sp by sampling from a (p − 1)-dimensional
unit-hypershere, as in [31]. We use the following algorithm
[32] to obtain Sp in (20).
Algorithm III.5. (Uniform Sampling from a p-Dimensional
Unit Hypersphere)

|u2,j |
√
|z
·
nβ̂2,j |
{j:βj =0∧ǫj >0} j
X
+
X

The above inequality follows from Proposition III.4 and
the last equation follows with assumption (A5.). Now, let
M2 , M3 > 0 be some real numbers with M3 λn > M2 , then
we have with arbitrarily high probability

and the proposition follows for sufficiently large n.

p
p
X
λn X |u2,j−s |
|u2,j−s |
E= √
= λn
√
n j=s+1 |zj · β̂2,j |
nβ̂2,j |
j=s+1 |zj ·



j=s+1

|u2,j |
= λn Ωp (ku2 k).
√
| nβ̂2,j |

= ku2 k · (−M2 + M3 λn ) > 0


xi x⊤
i

where
≤
The last equation follows from
the assumptions (A1.)-(A3.), Lemma 1 in [26] and Lemma 5
in the Technical Report associated with [30].
Additionally, we have that E is stochastically bounded from
below:

= λn ·

p
X



Gn (u1 , u2 ) − Gn (u1 , 0p−s ) > −M2 ku2 k + M3 λn ku2 k

· (u1 , (1 − αn )u2 ) = OP (ku2 k),
kθ∗∗
n −β0 k

|u2,j |
√
| nβ̂2,j |
{j:βj =0∧ǫj >0}
X

|u2,j |
+
√
| nβ̂2,j |
{j:βj =0∧ǫj =0}

√
√
· (u1 / n, (1 − αn )u2 / n)
−1



X

= λn
n
X

(32)



1) Generate p vectors with k entries
xj = (x1j , x2j , . . . , xkj )⊤ ,

j = 1, . . . , p,

(33)

where xij ∼ N (0, 1).
2) Calculate k p-dimensional vectors


ai =

p
X
j=1

xij
q
· ej ,
x2i1 + x2i2 + . . . + x2ip

where ej is the jth unit vector.
3) Set Sp := {ai ∈ Rp : i ∈ {1, . . . k}}.

i = 1, . . . , k,
(34)

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, XXX XXXX

IV. N UMERICAL E XPERIMENTS
A. Simulation Setup
Two different Monte Carlo studies are conducted, to assess
the performance of the proposed RWAD MM-Lasso.
Scenario 1: p > n, correlated predictors, cellwise outliers
A setup with p = 50 predictors and n = 30 observations
is considered. The regression parameters are defined by
βj = j/5 j ∈ {1, . . . , 5}, while βj = 0 j ∈ {6, . . . , 50}.
Correlated predictors xj , j ∈ {1, . . . , p} are generated by
sampling from a multivariate zero mean Gaussian distribution
with covariance matrix Σij = 0.5|i−j| , ∀i, j ∈ {1, . . . , p}.
The errors ui are zero mean i.i.d. Gaussian distributed with
variance σ 2 = 0.52 . The responses yi follow the linear model
yi = x⊤
i β + ui ,

i = 1, . . . , n.

To create cellwise outliers in X, ǫ = {0%, 10%, 20%, 30%}
of the predictors xj are contaminated. For the contaminated
predictors, 30% of the entries are independently and additively
contaminated by samples drawn from the distribution
xcont ∼ N (0, 1002 ).
Scenario 2: p > n, correlated predictors, cellwise outliers
and additive outliers
The setup is identical to Scenario 1 but additionally, 5% of
the responses are additively contaminated by samples drawn
from the distribution ycont ∼ N (0, 1002 ).
B. Performance Measures
To assess the performance in terms of parameter estimation
and model selection, we display the average mean squared
error (MSE)
MSE(β̂) =

p
R
1 X1X
(r)
(βj − β̂j )2 ,
R r=1 p j=1

(35)

the average false positive rate (FPR)
(r)

R
1 X |{j ∈ {1, . . . , p} : βj = 0 ∧ β̂j 6= 0}|
FPR(β̂) =
R r=1
|j ∈ {1, . . . , p} : βj = 0|
(36)
and the average false negative rate (FNR)
(r)
R
1 X |{j ∈ {1, . . . , p} : βj 6= 0 ∧ β̂j = 0}|
,
R r=1
|j ∈ {1, . . . , p} : βj 6= 0|
(37)
(r)
where β̂j refers to the parameter estimate of the rth Monte
Carlo experiment. All results represent averages over R = 100
Monte Carlo simulations.

FNR(β̂) =

C. Benchmark Methods and Choice of Parameters
The performance of the MM-RWAL is compared to:
• OLS Lasso [17]
• OCD Lasso [32]
• M-Lasso and adaptive M-Lasso [9]
• MM-Lasso and adaptive MM-Lasso [26]
• Sparse LTS [25]

8

For all methods, we use a grid of Nλ = 1000 candidate
regularization parameters, which are equally spaced on the
interval (0, λmax ]. Here, λmax is the value that results in a
Lasso estimate for which all regression parameters are equal
to zero, i.e., β̂ = 0. λmax is calculated by
λmax =

2
max y⊤ xj ,
n j∈{1,...,p}

(38)

using a robust correlation, as described in [25]. For all
methods, we choose the penalty parameter λ that provides
the lowest MSE so that the different methods are comparable, independent of the method that is used to select λ.
For the
 OCD Lasso
 estimator [32] we use a threshold of
chuber 2
tj = 10mad(x
. For all methods that include the SDO
j)
5
measure, 10 samples are used to approximate the supremum.
For the M-Lasso [9], Hubers ρ-function with a clipping point
of chuber = 1.215 is used, while for the adaptive M-Lasso we
use the bisquare ρ-function with a clipping of cbisquare = 3.44.
The initial estimate is the M-Lasso solution, as proposed in
[9]. For the sparse LTS, we use a subsample proportion of
75%. Further, 500 subsamples are used for the first two CSteps and the best 10 subsample sets are kept to carry out
the C-Steps until convergence [25]. For the MM-Lasso, the
S-Ridge estimator serves as initial estimate, as described in
[24]. For the MM-Lasso [26], we use a bisquare ρ-function
with clipping constant cbisquare = 3.44. The adaptive MMLasso and MM-RWAL use the same bisquare function and
are initialized with the the MM-Lasso estimate.
D. Simulation Results
Tables I and II display the estimation results for Scenarios 1
and Scenarios 2, respectively. Table III documents the average
computation time for one Monte Carlo Run of Scenario 2
TM
for the different methods, using an Intel R Core i7-4510U
with 8 GB RAM. In both setups, the OLS Lasso breaks
down even for low contamination in either the regressors or
the responses and FNR(β̂ OLS ) > 0.95since the OLS Lasso
selects a regularization parameter for which all coefficients
are equal to zero. The MM-RWAL performs best in terms
of robust model selection. In all experiments, it maximizes
the probability of correctly finding the indices of the non-zero
parameters, i.e. 1−(F P R+F N R), and in most cases provides
the best parameter estimation accuracy in terms of the MSE.
A drawback is that it is computationally heavy, which can be
contributed mainly to the MM-Lasso algorithm. Future work
will investigate combining the RWAL with computationally
more efficient estimators.
V. A R EAL DATA E XAMPLE OF S OURCE E STIMATION FOR
AN ATMOSPHERIC I NVERSE P ROBLEM
Quantifying the emissions of a pollutant into the atmosphere is essential, for example, in the case of nuclear power
plant accidents, volcano eruptions, or to track the releases of
greenhouse gases. In this paper, we apply penalized robust
estimation to determine the temporal releases of the particles
of the European Tracer Experiment (ETEX) at the source location. During the ETEX experiment tracers (perfluorocarbons)

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, XXX XXXX

ǫ = 0%

Estimator
OLS Lasso[17]

9

ǫ = 10%

ǫ = 20%

ǫ = 30%

n·MSE

FPR

FNR

n·MSE

FPR

FNR

n·MSE

FPR

FNR

n·MSE

FPR

FNR

0.14

0.11

0.74

1.88

0.019

0.95

2.8

0.019

0.97

2.76

0.018

0.97

OCD Lasso[32]

0.15

0.11

0.078

1.09

0.034

0.71

1.43

0.06

0.40

0.92

0.09

0.40

M-Lasso [9]

1.22

0.43

0.064

1.25

0.47

0.11

1.25

0.48

0.17

1.26

0.48

0.22

ad. M-Lasso [9]

0.15

0.066

0.12

0.36

0.07

0.24

0.41

0.07

0.30

0.53

0.05

0.41

MM-Lasso [26]

0.16

0.19

0.09

0.39

0.27

0.11

0.48

0.34

0.098

0.74

0.43

0.16

ad. MM-Lasso [26]

0.13

0.052

0.16

0.34

0.063

0.28

0.40

0.073

0.34

0.61

0.084

0.46

MM-RWAL

0.13

0.05

0.09

0.33

0.046

0.11

0.40

0.056

0.098

0.62

0.042

0.16

sparse LTS [25]

0.24

0.12

0.16

0.39

0.11

0.28

0.43

0.10

0.31

0.58

0.092

0.38

TABLE I
n· MSE, FPR AND FNR OF THE ESTIMATORS FOR S CENARIO 1, WITH ǫ CONTAMINATED PREDICTORS . B EST RESULTS AND PROPOSED ESTIMATOR ARE
HIGHLIGHTED WITH BOLD FONT.

ǫ = 0%

Estimator

ǫ = 10%

ǫ = 20%

ǫ = 30%

n·MSE

FPR

FNR

n·MSE

FPR

FNR

n·MSE

FPR

FNR

n·MSE

FPR

FNR

OLS Lasso [17]

2.24

0.013

0.98

1.88

0.019

0.95

2.8

0.019

0.97

2.76

0.018

0.97

OCD Lasso[32]

1.34

0.016

0.98

1.83

0.022

0.37

2.8

0.016

0.97

2.69

0.029

0.95

M-Lasso [9]

1.23

0.35

0.16

1.20

0.35

0.16

1.24

0.37

0.21

1.24

0.36

0.26

ad. M-Lasso [9]

0.46

0.031

0.16

0.52

0.031

0.40

0.59

0.036

0.43

0.75

0.023

0.56

MM-Lasso [26]

0.29

0.17

0.15

0.44

0.27

0.14

0.59

0.35

0.15

0.83

0.46

0.16

ad. MM-Lasso [26]

0.21

0.038

0.22

0.29

0.052

0.28

0.45

0.056

0.38

0.65

0.067

0.48

MM-RWAL

0.21

0.038

0.15

0.29

0.042

0.14

0.44

0.035

0.15

0.63

0.036

0.16

sparse LTS [25]

0.23

0.11

0.17

0.32

0.10

0.22

0.45

0.097

0.32

0.56

0.097

0.41

n·MSE, FPR AND FNR OF

THE ESTIMATORS FOR

TABLE II
S CENARIO 2, WITH ǫ CONTAMINATED PREDICTORS . B EST RESULTS AND PROPOSED ESTIMATOR ARE
HIGHLIGHTED WITH BOLD FONT.

ACT [s]
OLS Lasso [17]

0.02

OCD Lasso [32]

0.97

M-Lasso [9]

95.00

ad. M-Lasso [9]

97.51

MM-Lasso [26]

177.63

ad. MM-Lasso [26]

188.97

MM-RWAL

191.54

sparse LTS [25]

266.79

TABLE III
AVERAGE C OMPUTATION T IME (ACT) FOR DIFFERENT ESTIMATION
METHODS FOR ONE M ONTE -C ARLO RUN OF S CENARIO 2. B EST RESULTS
AND PROPOSED ESTIMATOR ARE HIGHLIGHTED WITH BOLD FONT.

were released into the atmosphere in Monterfil, Brittany in
1994. Hourly measurements were taken at 168 ground-level
sampling stations in 17 European countries as illustrated in
Fig. 3.
Atmospheric dispersion models, such as the Lagrangian
Particle Dispersion Model (LPDM) allow to formulate the
source estimation problem as a linear inverse problem according to Eq. (1) as follows. The regression matrix X is

Fig. 3. Locations of the 168 base stations used in the European Tracer
Experiment (ETEX).

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, XXX XXXX

2) Normalize the data robustly using the mad and the
median,
y − med(y)
y←
mad(y)

300
frequency=2013

250

10

200

xj ←

150
100
50
0
-4

-3

-2

-1

0

1

2

3

4

10 -13
Fig. 4. Histogram of the residuals.

30
25

xj − med(xj )
.
mad(xj )

3) When applying the median or the mad to the predictors
xj , only use samples which are greater than zero1 .
4) Apply a robust PCA [33] and reconstruct X using only
Np principal components, such that the mad of the Np
principal components corresponds to 90% of the total
mad2 .
5) Further, since we know that the number of particles
omitted by the source can only be positive βj ≥ 0, we
impose a non-negativity constraint on the parameters βj ,
j ∈ {1, . . . , p}.
This leads to the positive Lasso.

20

B. The Positive MM-RWAL Estimator

15

Definition V.1. (The Positive Lasso Estimator)
10

β̂ pos Lasso = arg min ||y − Xβ||22 ,

5
0

β

0

1

2

3

4

5

6

7

8

9

10 -14
Fig. 5. Histogram of the 68th predictor x68 (zero bin omitted).

estimated as in [6] by using the Flexible Particle Dispersion
Model and is formed by X = (X̃1 , . . . , X̃k , . . . , X̃K )⊤ ,
where the kth matrix describes the distribution of the particles
from the source to the kth sensor. Every regression parameter
βj , j ∈ {1, . . . , p} corresponds to the amount of PFC that is
released by the source at time instant j. The sampling interval
is one hour and 120 measurements are taken at each sampling
station, resulting in p = 120 unknown regression variables.
⊤ ⊤
The responses y = (y1 , . . . , yn )⊤ = (ỹ1⊤ , . . . , ỹk⊤ , . . . , ỹK
)
are a set of stacked observations of the K = 168 sensors.
The ETEX data is sparse in β, since only 12 of the
regression parameters, i.e. 10%, are unequal to zero. The
residuals r = y − Xβ, given the ground truth values of β,
are non-Gaussian, as displayed by the histogram in Fig. 4.
Additionally, the regression matrix X contains outliers, as
exemplified for the 68th predictor via the histogram in Fig. 5.
Furthermore, X is sparse, since most of the time, the PFC
particles do not reach a sensor, resulting in a regression matrix,
which is dominated by zero-valued cells.

A. Pre-Processing of the ETEX data
The following preprocessing steps are applied to the data:
1) Remove data points (xj , yj )⊤ , where all entries of the
predictor xj are equal to zero.

subject to ||β||1 ≤ t and βj ≥ 0 j ∈ {1, . . . , p}.

(39)

The positive Lasso can be calculated using a modified cyclic
coordinate descent algorithm as follows.
Algorithm V.2. (The Positive Lasso Estimator Using Cyclic
Coordinate Descent)
P
1) Standardize
the regressors so that i xij /n = 0 and
P 2
i xij = 1.
2) Initialize the regression parameters with an arbitrary
value, e.g.,
βj = 0 j ∈ {1, . . . , p}.
3) Calculate the update of the jth regression parameter by
keeping the parameters k 6= j fixed. Start with j = 1 and
end with j = p. If a regression parameter is negative
replace it by zero
 
 
n


X
X
β̂j = max S 
xij (yi −
xik βˆk ), λ , 0 .


i=1

k6=j

4) Repeat step 3 until convergence.

Here, S(x, λ) = sign(x)(|x| − λ)+ is the soft thresholding
function with (·)+ = max{·, 0}.
1 The predictors contain an overwhelming number of components which
are zero since X is highly sparse. Robust estimates like the median or the
mad will result in a value of zero if more than 50% of the entries of xj are
zero. Obviously, taking only the positive components into account leads to
estimates that are based on the distribution of the data that we are actually
interested in.
2 The reason for applying a robust PCA is to provide for all algorithms a
matrix that is not as badly conditioned as the original regression matrix. The
value 90% has been empirically determined from a range of possible values
between 85% − 99%.

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, XXX XXXX

MSE

FPR

FNR

1-(FPR+FNR)

OLS Lasso [17]

70.66

0.046

0.36

0.594

OCD Lasso [32]

3.75 · 104

0.029

0.438

0.533

Ad. M-Lasso [9]

68.23

0

0.25

0.75

MM-Lasso [26]

59.84

0.039

0.5

0.461

Ad. MM-Lasso [26]

61.88

0.053

0.143

0.804

MM-RWAL

62.1

0.053

0

0.947

TABLE IV
MSE, FPR AND FNR OF THE ETEX DATA USING DIFFERENT ESTIMATION
METHODS . B EST RESULTS AND PROPOSED ESTIMATOR ARE HIGHLIGHTED
WITH BOLD FONT.

The LARS algorithm can be modified in a similar way,
as proposed in [34]. Therewith, Algorithm V.2 is modified to
compute the positive MM-Lasso estimates,


p
n
X
X
yi − xi β
|βj |,
ρ
+λ
β̂ MM = arg min
σ̂
β
j=1
i=1
subject to βj ≥ 0 j ∈ {1, . . . , p}

and the extension to the positive MM-RWAL is straightforward.
C. Parameter Selection and Performance Metrics
All parameters are set as described in Section IV and the
optimal regularization parameter λ∗ is found by using the
parameter with the lowest MSE compared to the ground truth:
p

1X
(βj − β̂(λ))2
λ = arg min
p j=1
λ∈Λ
∗

The performance metrics are the mean squared error (MSE)
p

1X
MSE(β̂) =
(βj − β̂j )2 ,
p j=1

VI. C ONCLUSION
The problem of finding sparse solutions to underdetermined, or ill-conditioned, linear regression problems that
are contaminated by cellwise and rowwise outliers was investigated. We defined ’robust oracle properties’ that are required
to perform robust variable selection for such models. We
introduced and analyzed a robustly weighted and adaptive
Lasso type regularization term and integrated it into the objective function of the MM-estimator, resulting in the proposed
MM-Robust Weighted Adaptive Lasso (MM-RWAL) for which
we showed that at least the weak robust oracle properties
hold. An algorithm to compute the weights was proposed
and analyzed. The MM-RWAL outperformed existing robust
sparse estimators in numerical experiments and proved its usefulness in a real-data application of estimating the sparse nonnegative spatio-temporal emissions of a pollutant, given noisy
observations and an imprecisely estimated ill-conditioned and
sparse dispersion model containing cellwise and rowwise
outliers. In future work, the proposed RWAL penalty can easily
be integrated into the objective function of other rowwise
robust estimators to extend them to the cellwise contamination
framework.
ACKNOWLEDGMENT
The authors would like to thank Marta Martinez-Camara
and Martin Vetterli for making us aware of the ETEX experiment and for many interesting discussions on this application.

|{j ∈ {1, . . . , p} : βj = 0 ∧ β̂j 6= 0}|
|j ∈ {1, . . . , p} : βj = 0|

(41)

R EFERENCES

|{j ∈ {1, . . . , p} : βj 6= 0 ∧ β̂j = 0}|
.
|j ∈ {1, . . . , p} : βj 6= 0|

(42)

and the false negative rate (FNR)
FNR(β̂) =

Lasso estimate still has poor model selection properties. The
OCD Lasso estimates are dominated by outlying predictors
rendering this estimator useless for the above real-data application. The MM-RWAL is the only estimator, which correctly
detects all coefficients that are equal to zero, i.e. FNR=0.
The MM-RWAL also by far outperforms its competitors in
correctly finding the indices of the non-zero parameters, i.e.
1 − (FPR + FNR) = 0.947.

(40)

the false positive rate (FPR)
FPR(β̂) =

11

The false positive rate measures how many time instances
are flagged as containing a source emission, while there was
actually none. The false negative rate measures how many time
instances containing a source emission are falsely left out.
D. Results for the ETEX Experimental Data
Table IV displays the MSE, FPR and FNR for the ETEX
data set, while Fig. 6 shows the estimated source emissions
for the Lasso, Ad. M-Lasso, Ad. MM-Lasso and MM-RWAL.
In general, the difference in the MSE is not very large
between all estimators, except for the OCD-Lasso. For this
example, even the OLS Lasso estimate results in a reasonable
estimation, mainly thanks to robust PCA. However, the OLS

[1] A.M. Zoubir, V. Koivunen, Y. Chakhchoukh, and M. Muma, “Robust
Estimation in Signal Processing: A Tutorial-Style Treatment of Fundamental Concepts,” IEEE Signal Process. Mag, vol. 29, no. 4, pp. 61–80,
Jul 2012.
[2] E. Elhamifar and R. Vidal, “Sparse Subspace Clustering: Algorithm,
Theory, and Applications,” IEEE Trans. Pattern Anal. Mach. Intell, vol.
35, no. 11, pp. 2765–2781, Nov 2013.
[3] G. Liu, Z. Lin, S. Yan, J. Sun, Y. Yu, and Y. Ma, “Robust Recovery of
Subspace Structures by Low-Rank Representation,” IEEE Trans. Pattern
Anal. Mach. Intell, vol. 35, no. 1, pp. 171–184, Jan 2013.
[4] R. Couillet, F. Pascal, and J.W. Silverstein, “Robust Estimates of
Covariance Matrices in The Large Dimensional Regime,” IEEE Trans.
Inf. Theory, vol. 60, no. 11, pp. 7269–7278, Nov 2014.
[5] F. Pascal, Y. Chitour, and Y. Quek, “Generalized Robust Shrinkage
Estimator and Its Application to STAP Detection Problem,” IEEE Trans.
Signal Process, vol. 62, no. 21, pp. 5640–5651, Nov 2014.
[6] M. Martinez-Camara, B. Béjar Haro, A. Stohl, and M. Vetterli, “A Robust Method for Inverse Transport Modeling of Atmospheric Emissions
Using Blind Outlier Detection,” Geoscientific Model Develop, vol. 7,
no. 5, pp. 2303–2311, 2014.
[7] 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 Proc. IEEE Int. Conf. Acoust. Speech Signal
Process, Apr 2015, pp. 3422–3426.

90
80
70
60
50
40
30
20
10
0

12

30
β̂Lasso
Ground Truth

β̂M
Ground Truth

25
20

βj

βj

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. XX, NO. XX, XXX XXXX

15
10
5

0

20

40

60

80

100

0

120

0

20

40

j

100

120

80

100

120

30
Ad
β̂MM
Ground Truth

25

RWAL
β̂MM
Ground Truth

25
20

βj

20

βj

80

j

30

15

15

10

10

5

5

0

60

0

20

40

60

80

100

120

0

0

20

j

40

60

j

Fig. 6. Estimated source emissions for the Lasso, Ad. M-Lasso, Ad. MM-Lasso and MM-RWAL.

[8] W.-J. Zeng, H.-C. So, and Jiang X., “Outlier-Robust Greedy Pursuit
Algorithms in ℓp -Space for Sparse Approximation,” IEEE Trans. Signal
Process, vol. 64, no. 1, pp. 60–75, Jan 2016.
[9] E. Ollila, “Adaptive Lasso Based On Joint M-estimation of Regression
and Scale,” in Proc. Eur. Signal Process. Conf. IEEE, Aug 2016, pp.
2191–2195.
[10] F. R. Hampel, E. M. Ronchetti, P. J. Rousseeuw, and W. A. Stahel,
Robust Statistics: The Approach Based on the Influence Function, John
Wiley & Sons, 2005.
[11] R. A. Maronna, R. D. Martin, and V. J. Yohai, Robust Statistics, John
Wiley & Sons, 2006.
[12] P. J. Huber and E. M. Ronchetti, Robust Statistics, John Wiley & Sons,
2009.
[13] P. J. Rousseeuw and W. Van den Bossche, “Detecting Deviating Data
Cells,” arXiv preprint arXiv:1601.07251, 2016.
[14] F. Alqallaf, S. Van Aelst, V.J. Yohai, and R.H. Zamar, “Propagation of
Outliers in Multivariate Data,” Ann. Stat, pp. 311–331, Feb 2009.
[15] V. Öllerer, A. Alfons, and C. Croux, “The Shooting S-estimator for
Robust Regression,” Computation. Stat, vol. 31, no. 3, pp. 829–844,
2016.
[16] A. Leung, H. Zhang, and R. Zamar, “Robust Regression Estimation and
Inference in the Presence of Cellwise and Casewise Contamination,”
Comput. Stat. Data An, vol. 99, pp. 1–11, 2016.
[17] R. Tibshirani, “Regression Shrinkage and Selection via The Lasso,” J.
Roy. Stat. Soc. B. Met, pp. 267–288, 1996.
[18] J. Friedman, T. Hastie, and R. Tibshirani, “Regularization Paths for
Generalized Linear Models Via Coordinate Descent,” J. Stat. Softw,
vol. 33, no. 1, pp. 1, 2010.
[19] J. Huang, P. Breheny, and S. Ma, “A selective review of group selection
in high-dimensional models,” Statistical Science, vol. 27, no. 4, 2012.
[20] P.-L. Loh and M. J. Wainwright, “Regularized M-estimators with
nonconvexity: statistical and algorithmic theory for local optima,” in Adv.
Neural Inf. Process. Syst. 26, C. J. C. Burges, L. Bottou, M. Welling,
Z. Ghahramani, and K. Q. Weinberger, Eds., pp. 476–484. Curran
Associates, Inc., 2013.
[21] C.-H. Zhang and T. Zhang, “A general theory of concave regularization
for high-dimensional sparse estimation problems,” Statistical Science,
pp. 576–593, 2012.
[22] H. Lian and Y. Kim, “Nonconvex penalized reduced rank regression
and its oracle properties in high dimensions,” J. Multivariate Anal., vol.
143, pp. 383–393, 2016.

[23] H. Zou, “The Adaptive Lasso and Its Oracle Properties,” J. Am. Stat.
Assoc, vol. 101, no. 476, pp. 1418–1429, 2006.
[24] R. A. Maronna, “Robust Ridge Regression for High-Dimensional Data,”
Technometrics, vol. 53, no. 1, pp. 44–53, 2011.
[25] A. Alfons, C. Croux, and S. Gelper, “Sparse Least Trimmed Squares
Regression for Analyzing High-Dimensional Large Data Sets,” Ann.
Appl. Stat, vol. 7, no. 1, pp. 226–248, Mar 2013.
[26] E. Smucler and V. J. Yohai, “Robust and sparse estimators for linear
regression models,” Comp. Statist. Data Anal., vol. 111, pp. 116–130,
2017.
[27] D. Donoho, “Breakdown Properties of Multivariate Location Estimators,” Tech. Rep., Ph.D. qualifying paper, Harvard University, 1982.
[28] S. Van Aelst, E. Vandervieren, and G. Willems, “Stahel-Donoho
Estimators with Cellwise Weights,” J. Stat. Comput. Sim, vol. 81, no.
1, pp. 1–27, 2011.
[29] J. Huang, S. Ma, and C. Zhang, “Adaptive Lasso for Sparse HighDimensional Regression Models,” Stat. Sinica, pp. 1603–1618, 2008.
[30] V. J. Yohai, “High Breakdown-Point and High Efficiency Robust
Estimates for Regression,” Ann. Stat, pp. 642–656, 1987.
[31] W. P. Petersen and A. Bernasconi, “Uniform Sampling from an
n-Sphere,” Tech. Rep., Swiss Center for Scientific Computing, ETH,
Zürich.
[32] J. Machkour, B. Alt, M. Muma, and A. M. Zoubir, “The OutlierCorrected-Data Adaptive Lasso: A New Robust Estimator for the
Independent Contamination Model,” in Eur. Signal Process. Conf, to
be published, Aug 2017.
[33] C. Croux, P. Filzmoser, and M.R. Oliveira, “Algorithms for ProjectionPursuit Robust Principal Component Analysis,” Chemometr. Intell. Lab,
vol. 87, no. 2, pp. 218 – 225, 2007.
[34] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani, “Least Angle
Regression,” Ann. Stat, vol. 32, no. 2, pp. 407–499, Apr 2004.

