1

On a class of optimization-based robust estimators

arXiv:1705.02837v1 [] 8 May 2017

Laurent Bako

Abstract—We consider in this paper the problem of estimating
a parameter matrix from observations which are affected by two
types of noise components: (i) a sparse noise sequence which,
whenever nonzero can have arbitrarily large amplitude (ii) and a
dense and bounded noise sequence of "moderate" amount. This is
termed a robust regression problem. To tackle it, a quite general
optimization-based framework is proposed and analyzed. When
only the sparse noise is present, a sufficient bound is derived on
the number of nonzero elements in the sparse noise sequence that
can be accommodated by the estimator while still returning the
true parameter matrix. While almost all the restricted isometrybased bounds from the literature are not verifiable, our bound
can be easily computed through solving a convex optimization
problem. Moreover, empirical evidence tends to suggest that it is
generally tight. If in addition to the sparse noise sequence, the
training data are affected by a bounded dense noise, we derive
an upper bound on the estimation error.

I. I NTRODUCTION
In many engineering fields such as control system design,
signal processing, machine learning or statistics, one is frequently confronted with the problem of empirically uncovering
a mathematical relationship between a number of signals of
interest. The usual method to achieve this goal is to run an
experiment during which one measures (a finite number of)
samples of the relevant signals and proceed with fitting a
certain model structure to the experimental data samples. This
process is known as system identification [11], [19]. A issue of
critical importance during this process is that the experimental
data samples might be contaminated by a measurement noise
of relatively high level due for example to intermittent sensor
failures or various communication disruptions. To cope with
the troublesome effects of the noise, the model estimation must
be designed with care.
In this paper we consider the situation where the data are
corrupted by two types of noise: a sparse noise sequence which
shows up only intermittently in time but can take on arbitrarily
large values whenever it is nonzero; and a more standard dense
noise component of moderate amount.
II. T HE

ROBUST REGRESSION PROBLEM

Consider a system described by an equation of the form
yt = Ao xt + ft + et

(1)

where yt ∈ Rm and xt ∈ Rn are respectively the output and
the regressor vector at time t; Ao ∈ Rm×n is an unknown
parameter matrix; ft and et are some noise terms which are
unobserved.
L. Bako is with Laboratoire Ampère – Ecole Centrale de Lyon – Université
de Lyon – 36 avenue Guy de Collongue, 69134 Ecully, France – E-mail:
laurent.bako@ec-lyon.fr

N

Problem. Given a finite collection {xt , yt }t=1 of measurements obeying the relation (1), the robust regression problem
of interest here is the one of finding an estimate of the
parameter matrix Ao under the assumptions that {et } and {ft }
are unknown but enjoy the following (informal) properties:
• {et } is a dense noise sequence with bounded elements
accounting for moderate model mismatches or measurement noise.
• {ft } is such that the majority of its elements are equal
to zero while the remaining nonzero elements can be of
arbitrarily large magnitude. The nonzero elements of that
sequence are usually termed gross errors or outliers. They
can account for possible intermittent sensor faults. We
will refer to {ft } as the sequence of sparse noise.
For the time being, these are just informal descriptions of the
characteristics of the sequences {ft } and {et }. They will be
made more precise whenever necessary in the sequel for the
need of stating more formal results.
Let Y ∈ Rm×N and X ∈ Rn×N be data matrices formed respectively with N output measurements and regressor vectors.
Then it follows from (1) that
Y = Ao X + E + F,
m×N

(2)

m×N

where E ∈ R
and F ∈ R
are unknown noise
components. The matrices Y and X can be structured or
not, depending on whether the system (1) is dynamic or not.
For example, when the model (1) is of MIMO FIR type, Y
contains a finite collection of output measurements while X
is a Hankel matrix containing lagged inputs of the system. In
this case Y and X take the form


Y = y1 y2 · · · yN ,


u1
u2
···
uN
 u0
u1
· · · uN −1 


X= .
.
..  .
..
 ..
···
. 
u1−nf

u2−nf

···

uN −nf

where {ut } and {yt } stand respectively for the input and
output of the system and the maximum lag nf is called the
order of the model. In the sequel, the notations of the type yt
and xt with subindex t ∈ I , {1, . . . , N } refer to the columns
of Y and X respectively.
Relevant prior works. The so formulated regression problem
is called a robust regression problem in connection with the
fact that the error matrix F assume columns of (possibly)
arbitrarily large amplitude. It has applications in e.g., the identification of switched linear systems [1], [15], [14], subspace
clustering [2], etc. Existing approaches for solving the robust
regression problem can be roughly divided into two groups:
methods from the field of robust statistics [17], [12], [9] which
have been developed since the early 60s and a class of more

2

recent methods inspired by the compressed sensing paradigm
[3], [4], [18], [21], [13]. The first group comprises methods
such as the least absolute deviation (LAD) estimator [8], the
least median of squares [16], the least trimmed squares [17],
the family of M-estimators [9]. The latter group can be viewed
essentially as a refreshed look at the so-called least absolute
deviation method. There has been however a fundamental
shift of philosophy in the analysis. While in the framework
of robust statistics, robustness of an estimator is measured
in terms of the breakdown point (the asymptotic minimum
proportion of points which cause the estimation error induced
by an estimator to be unbounded if they were to be arbitrarily
corrupted by gross errors), in the compressed-sensing-inspired
category of robust methods, the analysis aims generally at
characterizing properties of the data that favor exact recovery
of the true parameter matrix Ao . In this latter group, the LAD
estimator is sometimes regarded as a convex relaxation of a
combinatorial sparse optimization problem.
To the best of our knowledge, only the papers [18] provides
an explicit bound on the estimation error induced by the
LAD estimator. However that bound does not fully apply to
the current setting since the estimators although similar are
of different natures. Indeed, the LAD estimator stands only
as a special case of the current framework. Moreover the
bound in [18] is not easily computable while ours is. The
references [4] and [13] provide some bounds for a noise-aware
version of the LAD estimator which are based respectively
on the Restricted Isometry Property (RIP) and a measure of
subspace angles. Unfortunately numerical evaluation of those
bounds is a process of exponential complexity, a price that is
unaffordable in practice.
A related but different problem from the regression problem
considered here is that of sparse signal recovery studied in the
field of compressed sensing [5], [7]. This is about finding the
sparsest solution to an underdetermined set of linear equations.
Various analysis approaches have been devised which rely on
the RIP constant, the mutual coherence, the nullspace property,
to name but a few. Again, these analysis results either cannot
be extended efficiently to the robust regression problem or lead
to bounds that are NP-hard to compute [20], [10], [6].
Contributions. In this paper we propose and analyze a class
of optimization-based robust estimators. It is shown that the
robust properties of the estimators are essentially inherited
from a key property of the to-be-optimized performance function (or loss function) called column-wise summability. The
proposed framework admits the LAD estimator and its usual
variants as special cases. Moreover it applies to both SISO and
MIMO systems. When the dense noise component E in (2)
is identically equal to zero, we derive bounds on the number
of gross errors (nonzero columns of F ) that the estimator is
able to accommodate while still returning the true parameter
matrix Ao . In comparison with the existing literature, the
proposed bounds have the important advantage that they are
numerically computable through convex optimization. When
both E and F are active, exact recovery of the true parameter
matrix is no longer possible. In this scenario, we derive upper
bounds on the parametric estimation error in function of the

amplitude of E and the number of nonzero columns of F .
Again, computable but (possibly) looser versions of those
bounds are obtainable.
The current paper can be viewed as a generalization of our
previous work reported in [3]. While [3] provides an analysis
of mostly a single estimator (namely the LAD estimator)
relying on nonsmooth optimization theory, we focus here
on a much larger class of optimization-based robust estimators by highlighting some key robustness-inducing properties.
Moreover, we provide, for the considered class of estimators,
stability results which permit the estimation of parametric error
bounds.
Outline. The rest of the paper is organized as follows. Section
III defines the optimization-based approach to the robust
regression problem. Section IV discusses the properties of the
proposed estimation framework. Section V provides further
comments. Section VI reports some numerical experiments.
Lastly, Section VII contains some concluding remarks.
Notations. R is the set of real numbers; R≥0 (respectively
R>0 ) is the set of nonnegative (respectively positive) real
numbers ; RN is the space of N -tuples (vectors) of real
numbers. For any vector x = [x1 · · · xN ]⊤ ∈ RN , the
p-norm of x with p ∈ {1, . . . , ∞} is defined by kxkp =
PN
p 1/p
. A special case is the limit case p = ∞
i=1 |xi |
in which kxk∞ = maxi=1,...,N |xi |. For any matrix A =
[a1 · · · aN ] with ai ∈ Rm , the induced p-norm of A is
defined by kAkp = supx∈RN ,kxkp =1 kAxkp .
Cardinality of a finite set. Throughout the paper, whenever S
is a finite set, the notation |S| will refer to the cardinality of
S. However, for a real number x, |x| will denote the absolute
value of x.
Submatrices and subvectors. Let X ∈ Rn×N and I =
{1, . . . , N } be the index set for the columns of X. If I ⊂ I,
the notation XI denotes a matrix in Rn×|I| formed with the
columns of X indexed by I. We will use the convention that
XI = 0 ∈ Rn when the index set I is empty.
III. A

CLASS OF ROBUST ESTIMATORS

Let DN be the set of N data points generated by system
(1) for any possible values of the noise sequences, i.e.,
n
DN = (Y, X) ∈ Rm×N × Rn×N :
o
f
e
∃(E, F ) ∈ GN
× GN
, (2) holds ,
f
e
with GN
⊂ Rm×N and GN
⊂ Rm×N denoting the set of dense
and sparse noise matrices respectively. The estimation problem
aims at determining the unknown parameter matrix Ao given
a point (Y, X) in DN . Of course, this quest would not make
much sense if the noises E and F were completely arbitrary
since in this case, we would have DN = Rm×N × Rn×N
hence losing any informativity concerning the data-generating
system. Therefore some minimum constraints need to be put
on E and F as informally described above.
With respect to the estimation problem just stated, an estimator
is a set-valued map Ψ : DN → P(Rm×n ), (Y, X) 7→
Ψ(Y, X) which is defined from the data space DN to the power

3

set P(Rm×n ) of the parameter space. For (Y, X) generated
by a system of the form (1), one would like to design an
estimator achieving, whenever possible, the ideal property that
Ψ(Y, X) = {Ao }. In default of that ideal situation, a more
pragmatic goal is to search for a Ψ so that Ao ∈ Ψ(Y, X) and
Ψ(Y, X) is of small size in some sense despite the troublesome
effects of the unknown noise components E and F . The design
of an optimal estimator requires specifying a performance
index (usually called a loss function) which is to be minimized.
In this paper, we study the properties of the estimator of
the parameter matrix Ao in (2) defined by
Ψ(Y, X) = arg min ϕ(Y − AX)

(3)

A∈Rm×n

where ϕ : M (R) → R≥0 is a convex function defined on the
set M (R) of (all) real matrices. It is assumed that ϕ has the
following properties:
P1. For all B, C ∈ M (R) of compatible dimensions,
ϕ([B

C]) = ϕ(B) + ϕ(C)

(5)

P3. There exists a constant real number ε ≥ 0 such that for
all B ∈ M (R) with n rows and N columns,
ℓ(B) − |Iεc (B)| ε ≤ ϕ(B) ≤ ℓ(B)

IV. P ROPERTIES

(6)

where

We first study the conditions under which the true parameter
matrix Ao in (1) can be exactly recovered. Theorem 1 and
Theorem 2 stated next show that if the number of nonzero
columns in the matrix V , E + F is less than a certain
threshold, then Ψ(Y, X) = {Ao }.
Theorem 1 (A necessary and sufficient condition). Let ϕ
be a function satisfying (4)-(6) with ε = 0 and Ψ be
defined as in (3). Let d be an integer and assume that
rank(X) = n. For any A ∈ Rm×n and Y ∈ Rm×N , let
Ic (Y − AX) = {t ∈ I : yt − Axt 6= 0}. Then the following
statements are equivalent.
(i)


Iεc (B) = i ∈ {1, . . . , N } : ℓ(bi ) > ε

∀A ∈ Rm×n , ∀Y ∈ Rm×N , |Ic (Y − AX)| ≤ d

and |Iεc (B)| is the cardinality of Iεc (B) and bi ∈ Rn is
the ith column of the (n, N )-matrix B.
The property (4) will be called column-wise summability.
Since ϕ is a function defined over the space of real matrices of any dimensions, it is also defined for n-dimensional
vectors of real numbers. Hence according to property (4), if
B = [b1 · · · bN ] with column vectors bi ∈ Rn , then
ϕ(B) =

N
X

ϕ(bi ).

i=1

The so-defined function ϕ is not necessarily a norm. For any
εo ≥ 0 and any vector norm ℓo , it can be verified that the
function ϕ defined by
ϕ(B) =

N
X

max(0, ℓo (bi ) − εo )

OF THE ROBUST ESTIMATORS

A. Exact recoverability
(4)

with [B C] denoting the matrix formed by concatenating column-wise B and C.
P2. There exists a matrix norm ℓ : M (R) → R≥0 such that
for all B, C ∈ M (R), conformable for addition,
ϕ(B) ≤ ϕ(B − C) + ℓ(C)

enjoys some impressive robustness properties with respect to
the sparse noise matrix F . The term sparse is used here to
mean that a relatively large proportion of the column vectors
of F are equal to zero. And saying that Ψ is robust with
respect to F means that Ψ(Y, X) does not depend on (or
is insensitive to) the magnitudes of the nonzero columns of
F under the sparsity condition. Therefore those few columns
which are nonzero can have arbitrarily large magnitude. As
will be shown in the sequel, the robustness properties of Ψ are
inherited from the properties P1-P3 of the objective function
ϕ. In the special case where ϕ is a norm, the properties P2P3 are automatically satisfied so that P1 becomes the only
key property required. As to the convexity of ϕ, it is intended
just for computational reasons as it eases the solving of the
optimization problem in (3).

(7)

i=1

is positive and convex and satisfies properties (4)-(6) but it
is not a norm for εo > 0 since in this case, ϕ(B) = 0 does
not imply that B = 0. But if εo = 0 in (7), then ϕ = ℓ
by (6) soP
that ϕ corresponds to the matrix norm defined by
N
ϕ(B) = i=1 ℓo (bi ). We note in this latter case that (6) is
trivial while (5) reduces to the triangle inequality.
We will show in the sequel that the estimator Ψ in (3)

⇒
(ii)
max
c

max

I ⊂I: Λ∈Rm×n
|I c |=d
Λ6=0




Ψ(Y, X) = A
(8)


1
ϕ(ΛXI c )
<
ϕ(ΛX)
2

(9)

Here and in the following, the notation I , {1, . . . , N } is used
to denote the index set for the columns of the data matrices.
Proof: We first note that the rank assumption on X is
intended to insure that (9) is well-defined since then, with ϕ
being a norm, ϕ(ΛX) 6= 0 whenever Λ 6= 0.
(i) ⇒ (ii): Assume that (i) holds.
Consider an arbitrary subset I c of I such that |I c | = d.
Let Λ be any matrix in Rm×n satisfying Λ 6= 0. Finally,
consider a matrix Y ∈ Rm×N defined by YI c = 0 and
YI 0 = ΛXI 0 where I 0 = I \ I c . Then Ic (Y − ΛX) ⊂ I c and
so |Ic (Y − ΛX)| ≤ d. Hence by (i) {Λ} = arg minH ϕ(Y −
HX) which means that ϕ(Y − ΛX) < ϕ(Y − HX) for any
H ∈ Rm×n , H 6= Λ. In particular, by taking H = 0 we get
ϕ(Y − ΛX) < ϕ(Y ). It follows from the property (4) that
ϕ(YI c − ΛXI c ) + ϕ(YI 0 − ΛXI 0 ) < ϕ(YI c ) + ϕ(YI 0 ).
Using now the relations YI c = 0 and YI 0 = ΛXI 0
yields ϕ(ΛXI c ) < ϕ(ΛXI 0 ) or, equivalently, ϕ(ΛXI c ) <

4

1/2ϕ(ΛX). Eq. (9) then follows from the fact that I c and
Λ are arbitrary.
(ii) ⇒ (i): To begin with, note that if Eq. (9) holds for some
d, then it holds also for any d0 ≤ d. As a result, the equality
|I c | = d in (9) can be changed to |I c | ≤ d. Assuming
(ii), let A ∈ Rm×n and Y ∈ Rm×N be matrices satisfying
|Ic (Y − AX)| ≤ d. Set I c = Ic (Y − AX) and I 0 = I \ I c .
Then for all Λ ∈ Rm×n such that Λ 6= 0,
2ϕ(ΛXI c ) < ϕ(ΛX) = ϕ(ΛXI c ) + ϕ(ΛXI 0 ),
where the equality is obtained by the property (4) of ϕ. It
follows that
ϕ(ΛXI c ) < ϕ(YI 0 − (A + Λ)XI 0 ).

(10)

On the other hand, we know by (5) that
ϕ(YI c − AXI c )−ϕ(YI c − (A + Λ)XI c ) ≤ ϕ(ΛXI c ).
Combining with the inequality (10) yields
ϕ(Y − AX) < ϕ(Y − (A + Λ)X).
Since Λ is an arbitrary nonzero matrix, this inequality says
that A is the unique minimizer of V (H) = ϕ(Y − HX).
Consider a data pair (Y, X) generated by (1). By letting

πϕc (X) = max d : Eq. (9) holds ,
(11)

and assuming that πϕc (X) > 0 we can see that whenever
|Ic (Y − Ao X)| ≤ πϕc (X), Ao can be exactly recovered by
computing Ψ(Y, X). Of course this is likely to hold only
if the dense noise component E does not exist. So in the
situation where E = 0, the theorem says that Ao can be
uniquely obtained by convex optimization provided that the
number of outliers (nonzero columns of F ) is less than or
equal to πϕc (X). For the condition of exact recoverability to
be checkable we must be able to compute πϕc (X). The bad
news are that evaluating numerically such a number is likely
to be NP-hard in most cases.
In the sequel, we investigate sufficient conditions of exact recovery which are more tractable from a numerical standpoint.
For this purpose let us introduce some definitions.
Definition 1. A matrix X = [x1 · · · xN ] ∈ Rn×N is said
to be self-decomposable if rank(X) = n and for all k ∈ I,
xk ∈ im(X6=k ) where X6=k , XI\{k} is the matrix obtained
from X by removing its k-th column and im(·) refers to range
space.
For a matrix to be self-decomposable it is enough that X6=k
be full row rank for any k ∈ I. Achieving this condition in
practice seems easy provided that the number N of measurements is large enough compared to the dimension n of X.
Definition 2 (self-decomposability amplitude). Let X ∈
Rn×N be a self-decomposable matrix. We call selfdecomposability amplitude of X, the number ξ(X) defined
by
o
n
(12)
kγk k∞ : xk = X6=k γk .
ξ(X) = max min
k∈I γk ∈RN −1

The so-defined ξ(X) constitutes a quantitative measure of

richness (or genericity) of the regressor matrix X. By richness
it is meant here how much, in a global sense, the columns of
X are linearly independent. ξ(X) is expected to be small if
the columns of X are somehow strongly linearly independent.
Remark 1. If for some k the norm of xk was to be considerably large in comparison to the norm of the other columns of
X, then ξ(X) would get large hence reducing recoverability
capacity of the considered class of estimators (see also Eq.
(9)). Such situations can be alleviated by normalizing each
column of X, i.e., for example by replacing (yk , xk ) by
(ỹk , x̃k ) , (yk / kxk k , xk / kxk k) under the assumption that
xk 6= 0 for all k ∈ I.
With the help of the device of self-decomposability amplitude (12), we can state a condition for exact recovery of the
parameter matrix Ao by solving the optimization problem in
(3). A similar result was proven in [3] for the Least Absolute
Deviation (LAD) estimator.
Theorem 2 (A sufficient condition for exact recovery). Let ϕ
be a function satisfying (4)-(6) with ε = 0 and Ψ be defined as
in (3). Assume that X is self-decomposable. Then the following
statement is true:
∀A ∈ Rm×n , ∀Y ∈ Rm×N ,

|Ic (Y − AX)| < T ξ(X)
⇒


(13)
Ψ(Y, X) = A .

where T : R>0 → R>0 is the function defined by T (α) =
1
1
1+
.
2
α
Proof: The proof is completely parallel to that of Theorem
11 in [3]. From the assumptions, each xk , k ∈ I, can be
written as a linear combination of the columns of X6=k . Let
γk ∈ RN −1 be any vector satisfying xk = X6=k γk . It follows
that for any Λ ∈ Rm×n ,
X

γk,t Λxt
ϕ(Λxk ) = ϕ
t∈I\{k}

with γk,t denoting the entry of γk ∈ RN −1 indexed by t.
Under the assumptions of the theorem, ϕ is a norm. So, it
is positive and satisfies the triangle inequality property. As a
result we can write
X
ϕ(Λxk ) ≤
|γk,t | ϕ(Λxt ) ≤ kγk k∞ (ϕ(ΛX) − ϕ(Λxk ))
t6=k

where the rightmost term follows from the property (4) of ϕ.
Since this holds for any γk such that xk = X6=k γk , it holds
also for
n
o
γk⋆ = arg min kγk∞ : xk = X6=k γ .
γ∈RN −1

Hence,
ϕ(Λxk ) ≤ ξ(X) (ϕ(ΛX) − ϕ(Λxk )) ∀k ∈ I, ∀Λ ∈ Rm×n .
(14)
or equivalently,
ϕ(Λxk ) ≤

ξ(X)
ϕ(ΛX) ∀k ∈ I, ∀Λ ∈ Rm×n .
1 + ξ(X)

Let I c be any subset of I and pose |I c | = d. Summing the

5

previous inequality over the set I c yields
max
Λ6=0

ϕ(ΛXI c )
1
 |I c |
≤
ϕ(ΛX)
2T ξ(X)

relation
(15)

Note that the term on the right hand side is well-defined since
by the self-decomposability assumption, rank(X) = n which
implies that ϕ(ΛX)
 6= 0 whenever Λ 6= 0. Therefore (9) holds
if |I c | < T ξ(X) and the conclusion follows from Theorem
1.
It is worth noting that the threshold T (ξ(X)) on the number
of correctable outliers does not depend on ϕ. Hence this
threshold is valid when the estimator is defined from any
matrix norm obeying (4).
Remark 2. The statement of Theorem 2 still holds true if we
replace ξ(X) with the ϕ-dependent number δϕ (X) defined by
ϕ(Λxk )
δϕ (X) = max sup
k∈I Λ6=0 ϕ(ΛX6=k )

(16)

when it is assumed that ϕ is a norm and rank(X6=k ) = n
for all k. Doing so will give a less conservative condition
for exact recovery. However δϕ (X) seems much harder to
evaluate numerically than ξ(X).
Remark 3 (A few useful properties of ξ(X)).
•

•

For any nonsingular matrix R ∈ Rn×n , ξ(RX) = ξ(X).
It follows that the number ξ(X) depends only on the
subspace spanned by the rows of the regressor matrix X.
For any self-decomposable X ∈ Rn×N , ξ(X) is lowerbounded in the following sense
1
,
N −1
This follows from the more general observation that
ξ(X) ≥

ξ(X) ≥ max P
k∈I

kxk k
t6=k kxt k

for any vector norm k·k. As a result, T (ξ(X) is upperbounded as follows
N
.
2
Theorem 2 provides a sufficient condition for exact recovery
in the situation where the function ϕ is a norm. Next, another
condition is stated which holds in the general case.
T (ξ(X)) ≤

ϕ(YI c − AXI c ) − ϕ(YI c − (A + Λ)XI c ) < ϕ(ΛXI 0 ).
By (5), we can note that ϕ(YI c − AXI c ) − ϕ(YI c − (A +
Λ)XI c ) ≤ ℓ(ΛXI c ). It then follows that
ℓ(ΛXI c ) < ϕ(ΛXI 0 )
is a sufficient condition for Ψ(Y, X) = {A}. Finally, invoking
(6) allows us to observe that ℓ(ΛXI 0 ) − |Iεc (ΛXI 0 )| ε ≤
ϕ(ΛXI 0 ) which implies that ℓ(ΛXI c ) < ℓ(ΛXI 0 ) −
|Iεc (ΛXI 0 )| ε is a sufficient condition for Ψ(Y, X) = {A}.
We have hence proved the proposition.
B. Uncertainty set induced by dense noise
When both E and F are nonzero in the data-generating
system (1), Ψ(Y, X) is likely to be a non-singleton subset of
P(Rm×n ) especially if we consider all possible realizations of
the unknown components E and F . In this case the desirable
properties of the estimator are in default of better (i) that it
contains Ao and (ii) that its size with respect to some metric
is as small as possible. In this section we are interested in
estimating the size of Ψ(Y, X) when both dense noise E and
sparse noise F are active in the data-generating system (1).
A notion of estimator gain. Similarly to the concept of
system gain in control [22], one could define the gain of an
estimator, that is, a quantitative measure of the sensitivity
of the estimator with respect to the perturbations affecting
the measurements. Consider a data pair (Y, X) generated
by a system of the form (1) with Ao being the parameter
matrix sought for. Let us fix the sparse noise matrix F or
view it somehow as part of the data-generating system. This
consideration proceeds from the fact that Ψ can be insensitive
to F (when acting alone) under, for example, the condition
derived in Theorem 2. Let E be bounded in the sense that
ℓ(E) is finite with ℓ being the norm appearing in (6). Then
we can define a gain of the estimator with respect to the dense
noise component E. More specifically, an (ℓ, q)-gain of the
estimator Ψ with respect to the dense noise E may be defined
by
kA⋆ − Ao kq
.
(18)
gℓ,q (Y, X) =
sup
ℓ(E)
A⋆ ∈Ψ(Y,X)
0<ℓ(E)<∞
F sparse

Proposition 1. Consider a triplet (ϕ, ℓ, ε) satisfying (4)-(6).
For A ∈ Rm×n and Y ∈ Rm×N , pose I c = Ic (Y − AX),
I 0 = I \ I c = {t ∈ I : yt − Axt = 0} and Iεc (ΛXI 0 ) =

t ∈ I 0 : ℓ(Λxt ) > ε . Then Ψ(Y, X) = {A} if

Here k·kq denotes matrix q-norm. The so-defined number
gℓ,q (Y, X) provides an upper bound on the distance from the
set Ψ(Y, X) to Ao in function of the amount of dense noise.
The following theorem and its corollaries show that if the
number of nonzero columns in F is no larger than a certain
threshold, then gℓ,q (Y, X) exists and is finite.

∀Λ ∈ Rm×n , Λ 6= 0.

Theorem 3. Let (Y, X) be the data generated by system (1)
subject to the noise components E and F . Consider a triplet
(ϕ, ℓ, ε) satisfying (4)-(6). Let S 0 ⊂ I be a set such that FS 0 =
0 and let S c = I \ S 0 . Assume that the matrix X and the
partition (S 0 , S c ) are such that there exists α > 0 such that

|Iεc (ΛXI 0 )| ε < ℓ(ΛXI 0 ) − ℓ(ΛXI c )

(17)

Proof: Ψ(Y, X) = {A} is equivalent to
ϕ(Y − AX) < ϕ(Y − (A + Λ)X)
for any Λ ∈ Rm×n , Λ 6= 0. Using the definitions of the sets
I 0 and I c and applying property (4) of ϕ yields the equivalent

ℓ(ΛXS 0 ) − ℓ(ΛXS c ) ≥ α kΛkq ∀Λ ∈ Rm×n ,

(19)

6

with k·kq denoting some matrix q-norm.
Then for any A⋆ ∈ Ψ(Y, X), it holds that


1
kA⋆ − Ao kq ≤
2ℓ(ES 0 ) + |Iεc | ε
(20)
c
γℓ,q (X, S )

with1 Iεc = Iεc (YS 0 − A⋆ XS 0 ) = t ∈ S 0 : ℓ(yt − A⋆ xt ) > ε
and
ℓ(ΛXS 0 ) − ℓ(ΛXS c )
(21)
γℓ,q (X, S c ) = inf
m×n
kΛkq
Λ∈R
Λ6=0

if the dense noise matrix E is such that ℓo (et ) ≤ εo for
all t ∈ I, then by taking ε = εo the set Iεc defined in the
statement of Theorem 3 corresponds to the empty set so that
(22) holds as well in this case. In connection with the concept
of estimator gain discussed earlier, one can interpret the factor
2/γℓ,q (X, S c ) as an estimate of the gain (of the estimator Ψ)
with respect to dense noise.
Lastly, it is interesting to see that when ϕ is a norm, if E =
0 then the result of Theorem 3 implies that Ψ(Y, X) = {Ao }
provided (19) is true.

where k·kq refers to matrix q-norm.

V. D ISCUSSIONS

Proof: By definition of Ψ(Y, X) in (3),

For the purpose of illustrating the extent of the results above,
let us discuss further the situation where ϕ reduces to a norm.

ϕ(Y − A⋆ X) ≤ ϕ(Y − AX) ∀A ∈ Rm×n
By letting Λ = A − Ao , Λ⋆ = A⋆ − Ao and applying (2), the
last inequality takes the form
⋆

ϕ(F + E − Λ X) ≤ ϕ(F + E − ΛX) ∀Λ ∈ R

m×n

.

In particular, for Λ = 0, we get ϕ(F + E − Λ⋆ X) ≤ ϕ(F + E)
which, thanks to property (4) of ϕ, takes the form
ϕ(FS c + ES c − Λ⋆ XS c )+ϕ(ES 0 − Λ⋆ XS 0 )
≤ ϕ(FS c + ES c ) + ϕ(ES 0 ).
Now applying property (5) to the first member of the left hand
side and rearranging yields
⋆

A. Scenario when the loss function is a norm
Corollary 1. Let (Y, X) be the data generated by system (1)
subject to the noise components E and F . Let S 0 and S c be
defined as in the statement of Theorem 3. Assume that ϕ is a
norm i.e., it satisfies (4)-(6) with ε = 0.

If X is self-decomposable and |S c | < T ξ(X) , then for any
A⋆ ∈ Ψ(Y, X),
kA⋆ − Ao kq ≤ Bϕ,q (|S 0 |, X)ϕ(ES 0 )
2
,
h
N −r i
σϕ,q (X) 1 −
T (ξ(X))
ϕ(ΛX)
σϕ,q (X) = inf
Λ6=0 kΛkq

Bϕ,q (r, X) =

Using (6) gives
ℓ(ES 0 − Λ⋆ XS 0 ) − |Iεc | ε − ℓ(Λ⋆ XS c ) ≤ ϕ(ES 0 ) ≤ ℓ(ES 0 ).
Here we used the fact that Iεc (ES 0 − Λ⋆ XS 0 ) is equal to the
set Iεc defined in the statement of the theorem.
Applying the triangle inequality property of ℓ, it can be seen
that ℓ(Λ⋆ XS 0 ) − ℓ(ES 0 ) ≤ ℓ(ES 0 − Λ⋆ XS 0 ). Combining with
the previous inequality yields
ℓ(Λ⋆ XS 0 ) − ℓ(Λ⋆ XS c ) ≤ 2ℓ(ES 0 ) + |Iεc | ε.
Finally, it follows from the definition of γℓ,q (X, S c ) in (21)
that


γℓ,q (X, S c ) kΛ⋆ kq ≤ 2ℓ(ES 0 ) + |Iεc | ε .

The condition (19) guarantees that γℓ,q (X, S c ) is well-defined
and is positive. Hence the statement of the theorem is established.
Theorem 3 constitutes an interesting stability result in that
it provides a finite upper bound on the distance from Ao to
the set Ψ(Y, X) as a function of the amplitude of the dense
noise matrix E. It applies to any estimator Ψ defined as in (3)
with ϕ a function obeying (4)-(6). In particular, in the situation
where ϕ is a norm (in which case ε can be taken equal to zero
in (6)), the inequality in (20) simplifies to
2
ℓ(ES 0 ).
γℓ,q (X, S c )

(22)

If ϕ is defined as in (7) (which, recall, is not a norm) and
1 The

notation Iεc is used for simplicity reasons.

(23)

where

⋆

ϕ(ES 0 − Λ XS 0 ) − ℓ(Λ XS c ) ≤ ϕ(ES 0 ).

kA⋆ − Ao kq ≤

ON SOME SPECIAL CASES

(24)

(25)

Proof: The principle of the proof is to show that
γℓ,q (X, S c ) is well-defined and then find a positive underestimate of it. Using the property (4) of ϕ and the fact that
ϕ = ℓ, we can write


2ϕ(ΛX) 1 ϕ(ΛXS c )
ℓ(ΛXS 0 ) − ℓ(ΛXS c )
.
=
−
kΛkq
kΛkq
2
ϕ(ΛX)
On the other hand we know from the proof of Theorem 2 (see
Eq. (15)) that
1
ϕ(ΛXS c )
≤
|S c |
ϕ(ΛX)
2T (ξ(X))
so that


1−


|S c |
ϕ(ΛX)
ℓ(ΛXS 0 ) − ℓ(ΛXS c )
≤
T (ξ(X)) kΛkq
kΛkq

Taking now the infimum on both sides of the inequality symbol
over all nonzero matrices Λ ∈ Rm×n yields


|S c |
≤ γℓ,q (X, S c ).
σϕ,q (X) 1 −
T (ξ(X))
It follows from the rank condition imposed on X (by the selfdecomposability assumption) that σϕ,q (X) > 0. This shows
that γℓ,q (X, S c ) is well defined and is strictly positive. Finally,
since ϕ = ℓ, invoking (22) gives the result.
Two important comments can be made at this stage.

7

First it is interesting to note that the bound Bϕ,q (r, X)
is an increasing function of ξ(X). Therefore it is all the
smaller as ξ(X) is small. That is, the error bound will be
small if the data matrix X is rich enough.
• Second, Bϕ,q (r, X) is a decreasing function of r. This
means that the upper bound on the estimation error
decreases when the number of gross error columns in F
decreases. In the extreme case where S 0 = N (no gross
error), Bϕ,q (|S 0 |, X) in (23) reduces to 2/σϕ,q (X).
Beyond these observations it should be noted that
 a key
assumption of Corollary 1 is that |S c | < T ξ(X) with S c
being the index set of the nonzero columns in F . Realizing
this condition requires on the one hand that the number of
nonzero columns in the sparse noise matrix F be small and
on the other hand that ξ(X) be small2 (which means that the
data must be generic). Indeed this condition is not necessarily
as strong as it might appear to be at first sight. For example,
it can be relaxed as follows. Observe that the sum E + F
is not uniquely defined from model (2). Taking advantage of
this, one can always absorb in E all nonzero columns of F
whose magnitude does not exceed a certain level. To see this,
let I = {t ∈ S c : ℓ(et + ft ) ≤ εo } where εo = maxt∈I ℓ(et ).
Then we can define Ẽ and F̃ such that E + F = Ẽ + F̃ and
F̃S 0 ∪I = 0 that is, we set ẽt = ft +et and f˜t = 0 for any t ∈ I
and (ẽt , f˜t ) = (et , ft ) otherwise. As a consequence, E and F
in Corollary 1 can be replaced by Ẽ and F̃ respectively so
that |S| and |S c | are replaced by |S| + |I| and |S c | − |I|. The

condition of the corollary then becomes |S c |−|I| < T ξ(X) ,
which is potentially easier to fulfill.
•

Remark 4 (sum of p-norms). Evaluating numerically the
bound Bϕ (r, X) might prove to be a hard problem due to
the potential difficulty in computing the term σϕ,q (X) in (25).
A particular case of interest is when ϕ consists of a sum
of p-norms
the column vectors, i.e. when it is defined by
Pof
N
bN ]. In this case
ϕ(B) =
i=1 kbi kp for B = [b1 · · ·
if we take q = 2 in (23) and (25), it is easy to see that
1/2
1/2
λmin (XX ⊤) ≤ σϕ,2 (X) with λmin (·) denoting the square
root of the minimum eigenvalue. Replacing σϕ,2 (X) with
1/2
λmin (XX ⊤) in (24) yields an overestimate of Bϕ (r, X) which
is computable.
Remark 5. Corollary 1 still holds true if one replaces
T (ξ(X)) with πϕc (X) defined in (11). As shown in [18], the
number πϕc (X) in (11) is computable although at the price
of a combinatorial complexity. However if the n-dimension of
X is small enough the complexity of the algorithm proposed
there can be affordable. Then by using our formula (24) and
Remark 4 above, it is possible therefore to obtain a smaller
bound on the estimation error.
B. Single output case: ℓ1 norm
In this section, we discuss for an illustrative purpose, the
applicability of Theorem 3 to the case of single-output systems. This is an interesting case to highlight since it represents
2 Recall that T is a decreasing function hence implying that T (ξ(X)) is
large when ξ(X) is small.

the most classical situation. Consider the single-output system
defined by
yt = (θo )⊤ xt + ft + et
(26)
where yt , et , ft are scalars and xt and θo are n-dimensional
vectors. By letting Y = [y1 · · · yN ] ∈ R1×N and defining
E and F similarly, we obtain
Y = (θo )⊤ X + F + E.

(27)

This last equation corresponds indeed to (2) where the matrix
Ao reduces
to the row vector (θo )⊤ . In this case, if we let
PN
n
ϕ(B) =
i=1 kbi k2 then for any θ ∈ R , the columns of
(the row vector) Y − AX are scalars so that
ϕ(Y − θ⊤ X) =

N
X

yt − θ⊤ xt

t=1

2

=

N
X

yt − θ⊤ xt . (28)

t=1

As a result, Ψ coincides in this case with the Least Absolute
Deviation (LAD) estimator. The following corollary specializes the result of Theorem 3 to the LAD estimator.
1×N
Corollary 2. Let (Y, X)
× Rn×N be generated by
 ∈R
c
c
model (26). Let S = t ∈ I : ft 6= 0 , S 0 = I \ S
 . Assume
c
that X is self-decomposable and |S | < T ξ(X) . Then for
any θ⋆ ∈ arg min Y − θ⊤ X 1 ,
θ∈Rn

where


kθ⋆ − θo k2 ≤ B1,2 |S 0 |, X kES 0 k1
2
,
h
N −r i
σ1,2 (X) 1 −
T (ξ(X))
X ⊤η 1
.
σ1,2 (X) = inf
η6=0
kηk2

B1,2 (r, X) =

Again here the bound B1,2 (r, X) can be numerically overestimated by following the idea of Remark 4.
VI. N UMERICAL

ILLUSTRATIONS

The performance of the estimator Ψ has been extensively
tested in some existing papers in the special case of the LAD
(see e.g., [3]) . We therefore concentrate here on evaluating
numerically an estimate of the gain of the estimator based on
Corollary 1 and Remark 4. The estimation is carried out for
the case where ϕ consists in the sum of 2-norms and q = 2.
Four different cases are studied:
(a) Static data: X ∈ R2×200 is sampled from a Gaussian distribution N (0, I2 ) with zero-mean and identitycovariance.
(b) Dynamic data generated by a switched linear system:
X ∈ R2×200 is formed with the regressors (yt−1 , ut−1 )
generated by a switched linear system composed of 3
subsystems of order 1. This is a switched ARX system defined by yt = aσ(t) yt−1 + bσ(t) ut−1 with the
switching signal σ(t) ∈ {1, 2, 3} generated from a
uniform distribution and input ut being a white noise
with Gaussian distribution; (a1 , b1 ) = (−0.40, −0.15),
(a2 , b2 ) = (1.55, −2.10) and (a3 , b3 ) = (1, −0.65).

8

3

3
T (ξ(X))
c (X)
πϕ

T (ξ(X))
c (X)
πϕ

Bound

Bound

2

2

1
1

0
0

2

4

6

8

10

12

14

16

18

20

22

24

26

28

0

2

Percentage of nonzeros [%]
(a) static system: ξ(X) = 0.0083

4

6

8

10

12

14

3

18

20

3
T (ξ(X))
c (X)
πϕ

T (ξ(X))
c (X)
πϕ

2

2

Bound

Bound

16

Percentage of nonzeros [%]
(b) switched system: ξ(X) = 0.0127

1

1

0

2

4

6

8

10

12

Percentage of nonzeros [%]
(c) linear system: ξ(X) = 0.0188

0
0

2

4

6

8

10

12

14

16

18

20

22

Percentage of nonzeros [%]
(d) nonlinear system: ξ(X) = 0.0107

Fig. 1: An overestimate of Bϕ using respectively πϕc (X) and T (ξ(X)) for a data matrix X ∈ R2×200 : (a) static data sampled
from a Gaussian distribution; (b) data generated by a switched system; (c) data generated by a linear dynamic system ; (d) data
generated by a dynamic nonlinear system. In each case, the x-axis is limited to the range of nonzero gross errors proportions
which statisfy the stability condition |S c | /N < T ξ(X) /N (see e.g., Corollary 1).
then the so-defined estimator will inherit robustness properties.
Considering a data set generated by a linear model subject to
both sparse and dense noises, we showed that the estimator is
insensitive to the sparse noise when this latter is acting alone
and provided that the number of its nonzero components is no
larger than a certain (computable) threshold. Conditions are
proposed for the exact recovery of the true parameter matrix
when only the sparse noise is active. When both types of noises
affect the measurements we propose computable bounds on
the parametric estimation error. By assuming stochasticity of
the dense noise sequence, the obtained bounds are probably
improvable by exploiting appropriately the statistics of the
dense noise. This is a matter than can be investigated in future
research.

(c) Dynamic data generated by a linear ARX system defined
by yt = a1 yt−1 + b1 ut−1 with the (a1 , b1 ) defined above
in case (b).
(d) Dynamic data generated by a nonlinear NARX system
2
defined by yt = (yt−1 + 2.5)/(1 + yt−1
) + ut−1 .
Following Remark 1, the columns of all data matrices X have
been normalized to unit 2-norm before being processed.
Figure 1 plots the obtained estimate of the estimator gain
against the proportion of correctable outliers. As remarked in
Section V, the gain estimate increases as the proportion of
outliers gets larger. But the growth rate of the gain estimate
depends on the genericity of the data matrix X. The more
generic the columns of X are, the smaller the growth rate
of the estimation error is when regarded as a function of
the proportion of outliers. The experiment confirms also the
intuition according to which static data tend to be more generic
than data generated by a dynamic system. Among the three
cases of dynamic systems, the linear system appears to be the
one generating the least generic data.

The author is grateful to the Associate Editor and the
anonymous reviewers for constructive feedback.

VII. C ONCLUSIONS

R EFERENCES

In this paper we have discussed a somewhat general framework for designing a robust estimator. Given the training data,
the estimator is defined as the minimizing set of a certain
performance index applying to the data. We have shown that
if the performance function possesses some key properties,

[1] L. Bako. Identification of switched linear systems via sparse optimization. Automatica, 47:668–677, 2011.
[2] L. Bako. Subspace clustering through parametric representation and
sparse optimization. IEEE Signal Processing Letters, 21:356–360, 2014.
[3] L. Bako and H. Ohlsson. Analysis of a nonsmooth optimization approach
to robust estimation. Automatica, 66:132–145, 2016.

ACKNOWLEDGEMENT

9

[4] E. Candès and P. A. Randall. Highly robust error correction by convex
programming. IEEE Transactions on Information Theory, 54:2829–
2840, 2006.
[5] E. J. Candès and M. B. Wakin. An introduction to compressive sampling.
IEEE Signal Processing Society, 25:21–30, 2008.
[6] D. L. Donoho, M. Elad, and V. Temlyakov. Stable recovery of
sparse overcomplete representations in the presence of noise. IEEE
Transactions on Information Theory, 52:6–18, 2006.
[7] S. Foucart and H. Rauhut. A mathematical introduction to compressive
sensing. Birkhäuser, 2013.
[8] P. J. Huber. The place of l1 -norm in robust estimation. Computational
Statistics and Data Analysis, 5:255–262, 1987.
[9] P. J. Huber and E. M. Ronchetti. Robust Statistics. A. John Wiley &
Sons, Inc. Publication (2nd Ed), 2009.
[10] A. Juditsky and A. Nemirovski. On verifiable sufficient conditions for
sparse signal recovery via ℓ1 minimization. Mathematical Programming,
127:57–88, 2011.
[11] L. Ljung. System Identification: Theory for the user (2nd Ed.). PTR
Prentice Hall., Upper Saddle River, USA, 1999.
[12] R. A. Maronna, R. D. Martin, and V. J. Yohai. Robust Statistics: Theory
and Methods. John Wiley & Sons, Inc., 2006.
[13] K. Mitra, A. Veeraraghavan, and R. Chellappa. Analysis of sparse
regularization based robust regression approaches. IEEE Transactions
on Signal Processing, 61:1249–1257, 2013.

[14] N. Ozay and M. Sznaier. Hybrid system identification with faulty measurements and its application to activity analysis. In IEEE Conference
on Decision and Control and European Control Conference, Orlando,
FL, USA, 2011.
[15] N. Ozay, M. Sznaier, C. Lagoa, and O. Camps. A sparsification approach
to set membership identification of a class of affine hybrid systems. IEEE
Transactions on Automatic Control, 57:634–648, 2012.
[16] P. J. Rousseeuw. Least median of squares regression. Journal of the
American Statistical Association, 79:871–880, 1984.
[17] P. J. Rousseeuw and A. M. Leroy. Robust Regression and Outlier
Detection. John Wiley & Sons, Inc., 2005.
[18] Y. Sharon, J. Wright, and Y. Ma. Minimum sum of distances estimator:
Robustness and stability. In American Control Conference, St Louis,
MO, USA, 2009.
[19] T. Soderstrom and P. Stoica. System identification. Prentice Hall, Upper
Saddle River, USA, 1989.
[20] A. M. Tillmann and M. E. Pfetsch. The computational complexity
of the restricted isometry property, the nullspace property, and related
concepts in compressed sensing. IEEE Transactions on Information
Theory, 60:1248–1259, 2014.
[21] W. Xu, E.-W. Bai, and M. Cho. System identification in the presence
of outliers and random noises: A compressed sensing approach. Automatica, 50:2905–2911, 2014.
[22] G. Zames. Feedback and optimal sensitivity: Model reference transformations, multiplicative seminorms, and approximate inverses. IEEE
Transactions on Automatic Control, 26:301–320, 1981.

