On the benefits of output sparsity for multi-label classification
Evgenii Chzhen∗1,2,3 , Christophe Denis†1 , Mohamed Hebiri‡1 , and Joseph Salmon§2

arXiv:1703.04697v1 [math.ST] 14 Mar 2017

1 LAMA,

Université Paris-Est Marne-la-Vallée, Université Paris-Est, 5 Boulevard Descartes,
77420, Champs-sur-Marne, France
2 LTCI, Télécom ParisTech, Université Paris-Saclay, 46 rue Barrault, 75013, Paris, France,
3 Institute for Information Transmission Problems of the Russian Academy of Sciences
(Kharkevich Institute), Bolshoy Karetny per. 19, build.1, 127051, Moscow, Russia

Abstract
The multi-label classification framework, where each observation can be associated with a set of
labels, has generated a tremendous amount of attention over recent years. The modern multi-label
problems are typically large-scale in terms of number of observations, features and labels, and the
amount of labels can even be comparable with the amount of observations. In this context, different
remedies have been proposed to overcome the curse of dimensionality. In this work, we aim at
exploiting the output sparsity by introducing a new loss, called the sparse weighted Hamming loss.
This proposed loss can be seen as a weighted version of classical ones, where active and inactive
labels are weighted separately. Leveraging the influence of sparsity in the loss function, we provide
improved generalization bounds for the empirical risk minimizer, a suitable property for large-scale
problems. For this new loss, we derive rates of convergence linear in the underlying output-sparsity
rather than linear in the number of labels. In practice, minimizing the associated risk can be performed
efficiently by using convex surrogates and modern convex optimization algorithms. We provide
experiments on various real-world datasets demonstrating the pertinence of our approach when
compared to non-weighted techniques.

1

Introduction

The field of Extreme Multi-Label Classification (EMLC), is concerned with classification problems
where the number of classes is large and possibly comparable with the number of observations, with one
specific characteristic: a data point is tagged with several labels. It has recently attracted a vast amount
of contributions due to the variety of problems that EMLC can model: text categorization [Gao et al.,
2004], functional genomics [Barutcuoglu et al., 2006], image classification [Li et al., 2014] to name a few.
The objective in multi-label classification is to predict a binary vector Y ∈ {0, 1}L for a given observation
X ∈ X , where L is the number of labels available and X = RD is the feature space.
A pioneering approach, called Binary Relevance (BR), requires to build L binary classifiers [Tsoumakas et al., 2010] in a one-versus-all fashion. Yet, when the number of classes is large,
possibly of the same order as the available observations, BR becomes too expensive. Moreover, BR
degrades when the labels are highly correlated and poor performance were established [Dembczyński
∗ evgenii.chzhen@univ-paris-est.fr

† christophe.denis@u-pem.fr
‡ mohamed.hebiri@u-pem.fr

§ joseph.salmon@telecom-paristech.fr

1

et al., 2012] in this context. Such limitations led to different workarounds such as probabilistic classifier
chains [Cheng et al., 2010], graph based algorithms [Tan et al., 2015] and correlated logistic regression [Li
et al., 2016].
Our motivation in this work starts from the observation that for a given X ∈ X the output label vector
Y is often sparse: for each data point only a few labels are generally active in real-world scenarios (i.e., Y
has few non-zero coordinates, say K, with K  L). To the best of our knowledge, exploiting the output
sparsity for improved learning appeared first in [Hsu et al., 2009] for EMLC, where the authors relied
on compressed sensing [Donoho, 2006, Candès et al., 2006] to approximate the output by K log(D)
measurements. In terms of theory Hsu et al. [2009] defined the notion of sparsity for E[Y |X]. Hence
the observed outcomes might not necessary be sparse. Meanwhile, in several real-world applications of
EMCL the output vectors are indeed sparse, see Table 3. Of course this can be due to errors in the labeling
process: some meaningful labels are expected to be missed by human annotators when the number of
labels L is huge. Consequences of such errors were investigated in [Jain et al., 2016], where the sparsity
was defined in terms of relevant labels for a each data point.
Here, we aim at providing, in a principle way, a statistical framework for EMLC by exploiting ideas
of sparsity and some heuristics proposed in [Hsu et al., 2009, Jain et al., 2016]. We introduce a novel
sparsity assumption, which reflects the nature of multi-label classification problems. Our assumption,
unlike the one proposed in [Hsu et al., 2009], restricts the possibility to observe a large amount of ones
in the output labels for any observed data point. We also propose a new loss, tailored to handle output
sparsity. The proposed loss function aims at producing a sparse prediction with the possibility to make a
small amount of mistakes on the zero labels of the underlying signal. Such a feature is also addressed
by Jain et al. [2016] who proposed to omit zero labels to build their classifier. This allows to handle
scenarios where the annotator is not malicious, and does not add wrong tags: this is typically the case for
a dataset like Wikipedia where the community would automatically erase wrong labels, while missing
ones might be common due to the millions of possible tags. Another benefit of our approach is that the
naive classifier that outputs only ones, is not achieving the best performance, as it is the case in [Jain
et al., 2016].
The organization of the paper is as follows. In Section 2, we introduce our framework, our output
sparsity model and our new loss function for EMLC. We also discuss the advantages of the proposed
loss over standard ones. In Section 3, we adapt classical generalization bounds for this loss. These
bounds are advantageous due to their linear dependence on the output sparsity level K. As a consequence
of our analysis we introduce a novel notion of Rademacher complexity. It can be seen as the classical
Rademacher complexity separated into two parts, associated with either active labels (the ones) or inactive
labels (the zeros). Our analysis is somehow similar to the one introduced in [Yu et al., 2014], though
our assumption and the loss function of Section 2 allows us to develop faster rates of convergence.
The final contribution of this section is an adaptation of Zhang’s lemma [Zhang, 2004] to this context.
In Section 4.1, we perform a numerical study in order to show the advantage of our loss function over
classical ones, for both synthetic and real data.

2

Framework and Notation

Throughout the paper, the following notation is used: for any N , L, D ∈ N, [D] B {1, . . . , D}, we write
N ∨ L B max{N , L} and N ∧ L B min{N , L}. We use N for the number of observations, L for the
number of labels and D for the dimension of the feature space X . For a matrix Q, we denote by Q>
its transpose. The Euclidean norm is written k · k2 and for any a, b ∈ RM we denote by ha, bi = a> b
their inner product. For any A, B ∈ RN ×L we denote by hA, Bi = Tr(A> B) their Hilbert-Schmidt
inner product. For a given matrix A ∈ RN ×L and for any p ≥ 1, its Schatten p-norm [Bhatia, 1997] is
 L∧N
1/p
P p
written kAkσ ,p B
σi (A)
, with σi (A) being the singular values of A in decreasing order. We
i=1

2

use kAkσ ,∞ B σ1 (A) to denote the spectral norm of the matrix A. We write X ∈ X = RD for a feature
vector and Y ∈ {0, 1}L for its associated label vector (Y l being the l th component of Y ), and P be the
underlying probability distribution of X × {0, 1}L .
Definition 1 (Sparse vectors). We say that a binary vector Y ∈ {0, 1}L is K-sparse for K ∈ [L] if
kY k0 B #{l ∈ [L] : Y l = 1} = K ,

(1)

where #A refers to the cardinality of a finite set A.

Sparsity is an important property that often arises in statistics and signal processing to bypass the
curse of dimensionality, especially in high dimensional regression or inverse problems [Elad, 2010,
Bühlmann and van de Geer, 2011]. In this work, though, we are rather interested in sparsity of the output
variable Y , a somehow unconventional context.
Assumption 1 (Label sparsity). There exists a positive integer K  L, s.t. ,
L
X

P
1{Y l =1} > K|X = 0, a.s. .
l=1

We emphasize that K  L and restrict our attention only to this extreme case. Note that Assumption 1
is different from the one presented in [Hsu et al., 2009] and is in fact stronger. Indeed, under AssumpP
P
tion 1, E[ Ll=1 1{Y l =1} |X] is bounded by K. Meanwhile, assuming that E[ Ll=1 1{Y l =1} |X] ≤ K does not
necessary mean that every realization of Y would be sparse.
Definition 2 (Active labels). We say that a label l ∈ [L] is active for a binary vector Y ∈ {0, 1}L if Y l = 1.
If Y l = 0 we say that l is inactive. The set of all active (inactive) labels is called active (inactive) set of
the vector Y .
The choice of the loss function is crucial for EMLC. One of the most popular one is the Hamming
loss, however this loss does not take into account label sparsity, since it treats all labels equally. In [Jain
et al., 2016], the authors proposed propensity loss functions, that exploit the observation that the label
vector Y for a given data point X, often does not consists of all possible active labels. Being more
precise, the authors introduced another label vector Y ∗ , which is a noisy version of the true label vector
Y . Assuming, that the true label vector Y is deterministic, the authors defined a label to be relevant, if it
is active in the underlying deterministic label vector Y .
Here we consider the following observation framework. Let DN = {(X1 , Y1 ), . . . , (XN , YN )} be
i.i.d. copies of (X, Y ). Let PN be the empirical probability distribution. Let F be a class of hypothesis,
i.e., each element f ∈ F is a measurable function from X to RL . With weights p0 , p1 ≥ 0 such that
p0 + p1 = 1 we can consider the following asymmetric risk
L 
X
h
i
h
i
R(f ) =
p0 E L0 (f l (X))1{Y l =0} + p1 E L1 (f l (X))1{Y l =1} .

(2)

l=1

The empirical counterpart of the risk in Eq. (2) is given by
R̂N (f ) =

N L 

1 XX
p0 L0 (f l (Xi ))1{Y l =0} + p1 L1 (f l (Xi ))1{Y l =1} .
i
i
N

(3)

i=1 l=1

Note that we also write R̂N (f , DN ) for R̂N (f ) when the dependency on DN is needed. In this work we
are interested in the empirical risk minimizer (ERM)
fˆ ∈ arg min R̂N (f ) ,
f ∈F

3

(4)

Proposed loss
Hamming loss: p0 = p1 = 0.5
[Jain et al., 2016]: p0 = 0, p1 = 1
2K
Our choice: p0 = 2K
L , p1 = 1 − L

Ŷ0 ≡ 0
p1 K

Ŷ1 ≡ 1
p0 (L − K)

Ŷ2K
p0 K

Ŷw ≡ 1 − Y
p1 K + p0 (L − K)

K
2
K − 2KL

0
2
2K − 2KL

2K 2
L

0

K
2
3K − 4KL

K
2

L−K
2

K
2

L
2

Table 1: Loss examples costs for several classifiers with underlying true label being K-sparse, with
K  L: Ŷ0 ≡ 0: output no label, Ŷ1 ≡ 1 output all labels, Ŷ2K : output correct active set plus K mistakes
on inactive set, Ŷw ≡ 1 − Y : always wrong
and to compare fˆ with f ∗ , the best predictor1 in F
f ∗ ∈ arg min R(f ) .

(5)

f ∈F

2.1

Motivation and examples

In this section we provide some motivational examples for the risk in Equation (2). First, we notice
that this risk is essentially asymmetric, and is motivated by the fact that under Assumption 1 the output
label vector consists of at most K active labels. Treating Y l = 1 and Y l = 0 in the same way might force
classifiers to overfit on inactive labels only.
We handle this phenomenon by adding weights to our loss function, and such weights only rely on
the (global) sparsity of the labels. This approach is different from the one presented in [Jain et al., 2016],
where the authors proposed to add weights that are (local) label dependent. Intuitively p0 should be a
small weight to reflect the presence of inactive labels and p1 should be almost 1 to give priority to active
labels over inactive ones.
l
A particular case of our loss includes Lj (Ŷ l ) = L0/1
j (Ŷ ) = 1{Y l ,j} for j ∈ {0, 1}, where Ŷ is a binary
prediction vector. This leads to introducing a weighted Hamming loss written as
Lw
0/1 (Y , Ŷ ) =

L 

X
p0 1{Ŷ l =1} 1{Y l =0} + p1 1{Ŷ l =0} 1{Y l =1} .

(6)

l=1

We remind that the Hamming loss can be written as
LH (Y , Ŷ ) =

L
X
l=1

1{Y l ,Ŷ l } .

(7)

Let us now study the advantage of the weighted Hamming loss with specific weights.
One of the drawbacks of the Hamming loss is that it treats all labels equally and therefore, does not
take into account label sparsity. For more insight, let us consider the scenario where Y is exactly K-sparse
and let us analyze the following classifiers (a synthesis is also given in Table 1):
•
•
•
•

Ŷ0 ≡ 0: predicts all labels inactive,
Ŷ1 ≡ 1: predicts all labels active,
Ŷw ≡ 1 − Y : misspredicts all labels,
Ŷ2K : correctly predicts the active set of Y and makes exactly K mistakes on its inactive set.

Intuitively, one would like to build a loss which is able to differentiate between the first three predictions
and Ŷ2K . Indeed, in large-scale problems the predictions similar to Ŷ2K provide more valuable insights
compared to the first three ones.
1 One does not have access to f ∗ since P is unknown.

4

With our choice of weights
p0 =

2K
,
L

p1 = 1 −

2K
,
L

(8)

the introduced loss function treats Ŷ0 , Ŷ1 and Ŷw almost equally and it promotes predictions with small
amount of mistakes on inactive sets and correct predictions on active sets. Meanwhile, the Hamming loss
does not make any difference between Ŷ2K and Ŷ0 , and the loss considered in [Jain et al., 2016] gives a
high promotion to naive classifiers like Ŷ1 .
Remark 1. In practice, the weights in Eq. (8) rely on the unknown sparsity constant K in Assumption 1.
Since this quantity is unknown to the practitioner, a simple strategy consists in performing a rough
estimation based on the observed labels. Hence, we consider
p̂0 =

2K̂
,
L

p̂1 = 1 − p̂0 ,

(9)

where we estimate the output sparsity level by the maximal sparsity on the observations:
K̂ = max
i∈[N ]

L
X

Yil .

l=1

3 Generalization Bounds
In this section we study generalization properties of empirical risk minimization (ERM) under sparsity
Assumption 1. The generalization bounds usually have the following form: with high probability, the
error of the estimator is bounded by an empirical estimate of the error plus a residual term depending on
the complexity of the hypothesis class investigated [Bartlett et al., 2005].
There are various notions of complexities, one of the most popular being the Vapnik-Chervonenkis
(VC) dimension [Vapnik, 1998]. However, VC-dimension is distribution-free, hence model independent,
so it is not well tailored to our context. An important alternative technique to estimate complexity of
the hypothesis class relies on Rademacher averages, which provide data driven bounds. In the context
of multilabel classification it was first used by [Yu et al., 2014] that additionally accounted for possible
missing labels.
In Section 3.1 we focus on adapting standard analysis based on Rademacher complexities [Koltchinskii, 2006, Bartlett et al., 2002, Bartlett and Mendelson, 2002] to provide generalization bounds linear
in K. In Section 3.1.1 we adapt Talagrand’s concentration inequality in the form of [Bousquet, 2002].
This result is a first step towards bounds based on local Rademacher complexities [Bartlett et al., 2005],
but further development on this road is outside the scope of this paper. In Section 3.2 we prove a
generalization of Zhang’s Lemma [Zhang, 2004] for the case of weighted loss functions.

3.1

Rademacher complexities

The standard technique of comparing the ERM with the best classifier in the class F [Bartlett and
Mendelson, 2002] relies on the following inequality
R(fˆ) ≤ R(f ∗ ) + 2 sup R(f ) − R̂N (f ) .

(10)

f ∈F

From now on, w.l.o.g. we assume L0 , L1 are bounded by 1. Hence, our next step is to notice that the
function supf ∈F R(f ) − R̂N (f , DN ) has bounded differences and that McDiarmid’s inequality [McDiarmid, 1989], reminded in the supplementary material, can be applied to control the supremum in terms
of its expectation. With Assumption 1, we can state
5

Proposition 1. Consider the following function
ϕ(DN ) = sup R(f ) − R̂N (f , DN ) ,
f ∈F

0
0
and let DN
be a second set of observations such that DN
and DN differ only in one observation (Xk , Yk )
0
0
0
0
and (Xk , Yk ), where (Xk , Yk ) is a copy of (X, Y ) independent from DN . Then,
0
ϕ(DN ) − ϕ(DN
) ≤

L
2K
p0 +
p , a.s. .
N
N 1

(11)

Proof. See supplementary material for the proof.
Notice that choosing p0 = 0, p1 = 1 leads to the best upper bound in terms of K. However, this
choice of weights is useless in practice: Ŷ1 the naive classifier that always predicts a label as active, is
a minimizer of the empirical risk. To avoid this degeneracy, we recommend p0 > 0. In particular, the
proposed choice in Equation (8) allows to bound the differences of ϕ almost surely by a factor 4K/N .
Without Assumption 1 the differences would be bounded by a constant of order L/N , a possibly not
informative bound in large-scale scenarios with N and L being large.
We also emphasize that McDiarmid’s inequality consists in bounding differences of the supremum
over the whole domain. However in our case, we can get a better bound almost surely, following [Combes,
2015]. In the light of the above discussion we have the following lemma
Lemma 1. Setting weights as in Eq. (8) and δ ∈]0, 1[, then with probability at least 1 − δ the following
inequality holds
r
K
1
sup R(f ) − R̂N (f ) ≤ EN sup R(f ) − R̂N (f ) + 4 √
log .
δ
N
f ∈F
f ∈F
Using a standard symmetrization inequality to control the expected deviation of the empirical mean
from the true one in terms of Rademacher averages, we can get the following proposition
Proposition 2. For any class of hypothesis F of measurable functions
EN sup R(f ) − R̂N (f ) ≤ 2p0 R0 (L0 ◦ F ) + 2p1 R1 (L1 ◦ F ) ,
f ∈F

where for j = 0 and j = 1
Rj (F ) =

N X
L
X
1
EN Eε sup
εil f l (Xi )1{Y l =j} ,
i
N
f ∈F

(12)

i=1 l=1

are two Rademacher complexities of the class F and εil are i.i.d. Rademacher variables.
Assuming that L0 and L1 are C-Lipschitz and using contraction property [Ledoux and Talagrand,
1991, Theorem 4.12.] we can write
EN sup R(f ) − R̂N (f ) ≤ 2Cp0 R0 (F ) + 2Cp1 R1 (F ) .
f ∈F

Combining Lemma 1, Proposition 2 and the contraction property [Ledoux and Talagrand, 1991, Theorem
4.12.] we get the following theorem

6

Theorem 1. Setting weights as in Eq. (8) and δ ∈]0, 1[, and assuming that L0 and L1 are C-Lipschitz,
then with probability at least 1 − δ, the following bound holds
r
K
1
∗
ˆ
R(f ) ≤ R(f ) + 4Cp0 R0 (F ) + 4Cp1 R1 (F ) + 8 √
log .
δ
N
Notice that, the previous theorem provides rates with linear dependence on K, hence our next goal is
to control R1 (F ) and R0 (F ), with possibly the following rate w.r.t. to K, L and N
L
R0 (F ) . √ ,
N

K
R1 (F ) . √ .
N

Such bounds√are needed in order to have linear dependency
on K. Indeed this would lead to p1 R1 (F ) '
√
R1 (F ) . K/ N and p0 R0 (F ) ' KR0 (F )/L . K/ N .
3.1.1 Refined bound: Talagrand’s inequality
Similar but somehow sharper bounds on the risk of the empirical minimizer fˆ can be obtained if we
manage to bound the largest variance of the class members f ∈ F [Bartlett et al., 2005]. This technique
relies on Talagrand’s concentration inequality [Bousquet, 2002] for empirical processes. The only
difference from the previous analysis holds in a refined version of Lemma 1.
Lemma 2. Setting weights as in Eq. (8), δ ∈]0, 1[ and letting r > 0 be such that
∀j ∈ {0, 1},

∀f ∈ F ,

E[max L2j (f l (X))] ≤ r ,
l∈[L]

then with probability at least 1 − δ the following inequality holds
r
sup R(f ) − R̂N (f ) ≤ 2E sup R(f ) − R̂N (f ) + K
f ∈F

f ∈F

(13)

32r
1 10K
1
log +
log .
N
δ 3N
δ

Proof. See supplementary material for details.
Applying Proposition 2 and the contraction principle [Ledoux and Talagrand, 1991, Theorem 4.12.]
we get the following theorem
Theorem 2. Setting weights as in Eq. (8), δ ∈]0, 1[ and letting r > 0 be such that the Conditions (13) are
satisfied, we have with probability at least 1 − δ
r
128r
1 40K
1
∗
ˆ
R(f ) ≤ R(f )+8Cp0 R0 (F ) + 8Cp1 R1 (F ) + K
log +
log .
N
δ 3N
δ
It is important to point out that this bound does not give an improvement in K, however it can be made
tighter if the variance bound r decreases w.r.t. N . As mentioned in [Bartlett et al., 2005], Theorem 2
is not useful when applied to the whole class of functions F . Meanwhile, by applying Theorem 2 to
a subset F 0 ⊂ F or to a modified version of F , a better result could be obtained [Bartlett et al., 2005,
Theorem 3.3.]. In multi-label settings, this type of bounds was first considered in [Xu et al., 2016],
where the authors developed a specific regularization for ERM-type algorithms. Though, they used local
Rademacher complexities and did not investigate the influence of K on their bound.

7

3.2

Margin-based performance bound

Previous sections have focused on establishing bounds on the risk (2) of some estimator fˆ that lies in
some class of functions F . Here, we are more concerned by controlling the (weighted) excess risk of
classifiers, a more natural goal in multi-label setting.
A classifier is a measurable function g : X → {0, 1}L , and here we restrict our attention to classifiers
of the form



f l (x) ≥ 0
1,
gfl (x) = 
,
(14)

0,
otherwise
where2 f : X 7→ RL . A natural measure of performance in the context of multi-label classification with
sparse labels is the Weighted Hamming risk discussed in Section 2.1
Rw
0/1 (f ) =

L 
X
l=1

h
i
h
o
p0 E 1{f l (X)≥0} 1{Y l =0} + p1 E 1{f l (X)<0} 1{Y l =1} ,

(15)

0/1 l
l
with L0/1
0 (f ) = 1{f l (X)≥0} , L1 (f ) = 1{f l (X)<0} and weights p0 , p1 as in Eq. (8). Let η(X) =
(η 1 (X), . . . , η L (X))> be the regression function, where each component is given by η l (X) = P{Y l = 1|X}.
The Bayes classifier is given by

Proposition 3. The minimizer of Eq. (15) is given for all l ∈ [L] by


l

1, η (x) ≥ p0
(gη∗ )l (x) = 
.

0, η l (x) < p0
w
l
We denote by Rw∗
0/1 = R0/1 (η (x) − p0 ) the minimum of Eq. (15) and define the Weighted Hamming
excess risk of a classifier f as
w
w∗
E0/1
(f ) = Rw
0/1 (f ) − R0/1 .

(16)

l
Since the cost functions L0/1
j (f ) used in the risk (15) are not C-Lipschitz, the analysis provided in the
previous sections is not applicable. The purpose of the present section is to develop a different analysis to
bound the above excess risk.
Given the set of observations DN , an estimator of the Bayes classifier can be defined as the minimizer
of the empirical counterpart of the above risk

R̂w
N ,0/1 (f ) =

N L 

1 XX
p0 1{f l (Xi )≥0} 1{Y l =0} + p1 1{f l (Xi )<0} 1{Y l =1} .
N

(17)

i=1 l=1

Minimization the former is a computationally difficult problem, since indicator functions are not convex.
l
Therefore some convex surrogate of L0/1
j (f ) can be introduced. Establishing a generalization of Zhang’s
Lemma [Zhang, 2004] for the case of weighted loss functions, we are able to show good performance in
terms of the Weighted Hamming excess risk Eq. (16) while dealing with its convexification [Zhang, 2004,
Bartlett et al., 2006]. Due to major advances in convex optimization theory [Boyd and Vandenberghe,
2004] it has become possible to solve large-scale convex problems efficiently. Convex optimization in machine learning and more specifically in the classification literature has proved its efficiency. It is reflected
by the popularity of such methods as Boosting [Freund and Schapire, 1997], logistic regression [Friedman
et al., 2000] and support vector machine [Vapnik, 1998].
2 for instance, f can be defined as in the previous sections

8

Let φ : R 7→ R+ be a non-negative convex surrogate. We consider a convexified form of Equation (15)
Rw
φ (f ) =

L 
X
h
i
h
i
p0 E φ(−f l (X))1{Y l =0} + p1 E φ(f l (X))1{Y l =1} ,

(18)

l=1

where f : X 7→ RL (one may keep in mind to Equation (14) to deduce a classifier g from f ). Now,
the aim is to show that optimizing an empirical version of a convexified risk can be beneficial in terms
l l
of weighted Hamming loss. For convenience, we also introduce the following notation Rw
φ (η , f ) =
p1 η l φ(f l ) + p0 (1 − η l )φ(−f l ), so that

Rw
φ (f ) =

L
X
l=1

h
i
l
l
EX Rw
φ (η (X), f (X)) .

We write fφ∗ (η) : X 7→ RL for the minimizer
fφ∗ (η) ∈ arg min Rw
φ (η, f ) ,

(19)

f

l l
and ∆Rw
φ (η , f ) for the label-wise excess risk
l l
w l l
w l
∗ l
∆Rw
φ (η , f ) = Rφ (η , f ) − Rφ (η , (fφ ) ) .

Now, we are ready to state the main result of this section
Lemma 3 (Weighted Zhang’s Lemma). Assume that for all l ∈ [L] the function φ is classificationcalibrated, i.e., (fφ∗ )l (η l ) > 0 when η l > p0 and (fφ∗ )l (η l ) < 0 when η l < p0 . Assume there exist c > 0 and
s ≥ 1 such that for all l ∈ [L] and for all η l ∈ [0, 1],
s

l
p0 − η l ≤ cs ∆Rw
φ (η , 0) .

(20)

Then for any measurable vector function f ,
Rw (f ) − Rw∗ ≤ cL

s−1 
s Rw (f
φ

∗
) − Rw
φ (fφ )

1
s

.

(21)

Proof. See supplementary material for details.
Notice that in the context of multi-label classification, convexification deteriorates the rate of convergence in terms of its dependence on the number of labels. While in previous sections the rate was linear in
K, we now pay an additional L(s−1)/s term (replacing the total number of labels L by the sparsity constant
w ∗
K is a point we plan to address in future works). The term on the right hand side Rw
φ (f ) − Rφ (fφ ) can be
decomposed into estimation and approximation errors, both depending on the complexity of the class F .
However, the approximation error could also be bounded at the rate K/N up to a log factor (depending
on the complexity of F ) by adapting more refined techniques, see for instance in [Bartlett et al., 2006,
Theorem 4]. Since our loss function is asymmetric we assume (fφ∗ )l (η l ) < 0 when η l < p0 in Lemma 3.
In the original paper [Zhang, 2004] the author assumed that the surrogate function φ was invariant with
respect to reflection of coordinates. Note that standard surrogate functions satisfy all the assumptions
of Lemma 3:
Example 1. Here we provide some examples of convex functions φ which satisfy the assumptions of
Lemma 3:
9

Settings
N = 2D = 2L = 200
K =2
K =6
K = 10
K = 14
K = 18

Median output sparsity
Weighted
Classical
2.47(0.02)
0.04(0.0)
6.83(0.01) 0.43(0.01)
9.85(0.03) 1.81(0.04)
10.90(0.14) 4.11(0.15)
10.98(0.21) 6.61(0.2)

Recall (micro)
Weighted Classical
1.0(0.0) 0.02(0.0)
1.0(0.0) 0.07(0.0)
0.90(0.0) 0.18(0.0)
0.72(0.0) 0.29(0.0)
0.58(0.0) 0.36(0.0)

Precision (micro)
Weighted Classical
0.80(0.0) 1.0(0.03)
0.88(0.0) 1.0(0.0)
0.91(0.0) 1.0(0.0)
0.93(0.0) 0.99(0.0)
0.95(0.0) 0.99(0.0)

Table 2: Numerical experiments with synthetic data. Evaluation of each quantity is averages over 100
simulations and the variance is reported in brackets. The regularization parameter is fixed λ = 30.
√
• Square loss : φ(v) = (1 − v)2 with s = 2, c = 2 and
(fφ∗ )l (η l ) =

η l − p0
,
η l (1 − 2p0 ) + p0

√
• Boosting loss: φ(v) = e−v with s = 2, c = 2 and
(fφ∗ )l (η l ) =

(1 − p0 )η l
1
log
.
2
(1 − η l )p0

Proof. See supplementary material for details.

4

Linear prediction approach and experiments

We now provide explicit bounds on Rademacher complexities for linear approximation of Y . This setting
has been recently studied by Yu et al. [2014], Xu et al. [2016], Chen and Lin [2012], though the choice of
loss function considered by the authors does not handle output sparsity. We start with a generalization
bound for the linear prediction case.
Theorem 3. Let X = RD and f (Xi ) = W > Xi , where W ∈ RD×L , W = {W ∈ RD×L : kW kσ ,1 ≤ λ}.
Let W ∗ be the minimum of the risk Eq. (2) over W , and let Ŵ be the minimum of its empirical
counterpart Eq. (3). Assume that EX [kXk2 ] ≤ 1, and weights as in Eq. (8). Then, with probability at
least 1 − δ we have
r
r

 K
K
K 
1
∗

R(Ŵ ) ≤ R(W )+4λC 
+√
log .
 + 8 √
N
δ
NL
N
Proof. See supplementary material for details.
Theorem 3 is an important case of the analysis developed in Section 3.1, illustrating √
that with convex
surrogates, we√can still get rate of convergence linear in K. Also notice that the term K/ N L is small as
compared to K/N , due to our sparse scenario K  L.
√
This result also provides a theoretical way to choose the parameter λ ' K controlling the size of W .
By increasing λ and, as a consequence the set W , we obtain a better approximation of the Bayes rule
in Proposition 3. At the same time the increase of λ leads to a worse estimation error in Theorem 3.

10

Dataset: bibtex

Dataset: delicious
0.25

Mean test F1-micro

Mean test F1-micro

0.20

0.15

0.10

0.05

Weighted square loss: threshold=0.5
Classical square loss: threshold=0.5

0.00
0

25

50

75

100

125

lambda

150

175

0.20
0.15
0.10
0.05

Weighted square loss: threshold=0.5
Classical square loss: threshold=0.5

0.00

200

0

Dataset: stackex cs

0.35

25

50

75

100

125

lambda

150

175

200

Dataset: stackex philosophy
0.30

Mean test F1-micro

Mean test F1-micro

0.30
0.25
0.20
0.15
0.10
0.05

Weighted square loss: threshold=0.5
Classical square loss: threshold=0.5

0.00
0

25

50

75

100

125

lambda

150

175

0.25
0.20
0.15
0.10
0.05

Weighted square loss: threshold=0.5
Classical square loss: threshold=0.5

0.00

200

0

25

50

75

100

125

lambda

150

175

200

Figure 1: Mean test F1 -micro score for various datasets, with θ fixed to 0.5 for predicting with
Equation (23).

4.1

Experiments

Let us introduce some notation to compute the estimator analyzed in Theorem 3. First, let W =
[W 1 , . . . , W L ] ∈ RD×L , W = {W ∈ RD×L : kW kσ ,1 ≤ λ} and Xi ∈ RD for each i ∈ [N ]. We denote by
X = [X1 , . . . , XN ]> ∈ RN ×D the design matrix and by Y = (Yil )i=1,...,N ,l=1,...,L ∈ RN ×L the label matrix.
The estimator we considered is the solution of the convex problem
Ŵ ∈ arg min
W ∈W

where we denote by A
defined by

1
P̂
N

(Y − XW )

2
σ ,2

,

(22)

B the Hadamard product of two matrices, and by P̂ ∈ RN ×L the weight matrix
Pˆil

p


 p̂0
p
=

 p̂1

if Yil = 0
,
if Yil = 1

Notice that by choosing p̂0 = p̂1 = 1/2, Ŵ becomes a usual Frobenius norm minimizer under nuclear
norm constraint [Xu et al., 2016, Yu et al., 2014, Chen and Lin, 2012]. However our choice of the
weight matrix P̂ allows to handle sparse label space. The minimization problem in (22) is solved with
a Frank-Wolf algorithm [Jaggi and Sulovský, 2010, Jaggi, 2013], reminded in supplementary material,
since it is a simple and efficient algorithm, which requires to compute only the largest singular value of the
gradient at each iteration. Once the matrix Ŵ is obtained, we form the final prediction by thresholding
each entry of the matrix X Ŵ , i.e., for a threshold θ
Ŷil = 1{hXi ,Ŵ l i≥θ} .
11

(23)

Receiver operating characteristic: bibtex

Receiver operating characteristic: delicious
1.0

True Positive Rate

True Positive Rate

1.0

0.8

0.6

0.4

0.2

0.0

Weighted: (area = 0.80)
Classical: (area = 0.79)
0.0

0.2

0.4

0.6

False Positive Rate

0.8

Receiver operating characteristic: cs

0.4

0.2

Weighted: (area = 0.68)
Classical: (area = 0.74)
0.0

0.2

0.4

0.6

False Positive Rate

0.8

1.0

Receiver operating characteristic: philosophy
1.0

True Positive Rate

True Positive Rate

0.6

0.0

1.0

1.0

0.8

0.6

0.4

0.2

0.0

0.8

Weighted: (area = 0.83)
Classical: (area = 0.91)
0.0

0.2

0.4

0.6

False Positive Rate

0.8

1.0

0.8

0.6

0.4

0.2

0.0

Weighted: (area = 0.90)
Classical: (area = 0.88)
0.0

0.2

0.4

0.6

False Positive Rate

0.8

1.0

Figure 2: Micro averaged ROC curve for the best classifier obtained via cross-validation over λ
in Equation (22).
To get a relevant regularization parameter, we have performed 5-fold cross-validation on each realworld dataset, and for each scoring function with regularizers selected among the following values:
λ ∈ {0.1, 0.5, 1, 10, 20, 40, 60, 80, 100, 120, 140, 160, 180, 200}. To evaluate the influence of the loss
functions on the method, we have considered the following score functions, see for instance [Natarajan
and Jain, 2016]:
•
•
•
•

F1 -micro score,
AUC (micro averaged),
Recall (micro averaged),
Precision (micro averaged).

The goal of our experiments is to show the advantage of the weighted scheme proposed in Equation (22)
over non-weighted techniques, in terms of different performance measures. The benefit of our introduced
loss is reflected in

4.2

Data

For our empirical study we consider both synthetic and real-world datasets. By controlling the ration K/L
on synthetic dataset (see supplementary material for details), we can show empirical evidence to support
our framework. We also considered the following real-world datasets, bibtex [Katakis et al., 2008],
delicious [Tsoumakas et al., 2008], stackex cs, philosophy [Charte et al., 2015]. The real-world datasets
are obtained via mldr [Charte and Charte, 2015] R package, and normalized to satisfy the assumption
of the Theorem 3. Table 3 presents a description of real-world datasets with a focus on the estimated
sparsity.

12

Dataset
bibtex
delicious
stackex philosphy
stackex cs

N
7395
16105
3971
9270

D
1836
500
842
635

L
159
983
233
274

K̂
28
25
5
5

p̂0
0.352
0.051
0.043
0.036

Table 3: Datasets summary: Number of observations (N ), features dimensions (D), number of labels (L),
maximal observed sparsity (K̂) and p̂0 evaluated as Eq. (9).

4.3

Results

The performance of the proposed method on synthetic data is reported in Table 2. We observe a decreasing
performance of our method as K grows, but note that for the large sparsity constants (e.g., K ≥ 14), our
scenario K  L no longer holds.
The performance of out algorithm on the real-world datasets is reported on Figure 1 and 2. The
proposed algorithm achieves better performance over classical (non-weighted) in terms of F1 -micro score.
It is important to point out that our introduced loss is admitting a prediction to be good if it has few
amount of mistakes on the inactive set of the underlying signal Y . Since, ROC curves consider only True
Positive Rates and False Positive Rates, the method we proposed may show worse performance in terms
of AUC. Indeed, by allowing some amount of mistakes on the inactive set of the vector Y , we increase
the False Positive Rate and hence, decrease the performance in terms of the AUC score. This effect is
illustrated on Figure 2. To plot the ROC curves, we consider the classifier with λ fixed achieving the best
performance in terms of F1 -micro score, and each point on the curve corresponds to different threshold
levels θ in Eq. (9).

5

Conclusion

We have introduced a new framework for multi-label classification with large number of classes. We
have defined a natural output sparsity assumption and provided an associated weighted loss (and risk)
leveraging this sparsity property. Interestingly, we could prove generalization bounds linear in the sparsity
of the labels. We could also control the performance of a minimizer of a convexified version of our
proposed risk. This leads to good practical performance, not only for the introduced risk but also for
standard scores in multi-label setting such as recall or F1 -scores.

6 Acknowledgement
We would like to thank Jean-Baptiste Schiratti, Mónika Csikós and Guillaume Papa for their useful
comments and remarks. This work was partially supported by ”laboratoire d’excellence Bézout of
Université Paris Est” and by ”Chair Machine Learning for Big Data at Télécom ParisTech”.

References
P. L. Bartlett and S. Mendelson. Rademacher and Gaussian complexities: risk bounds and structural
results. J. Mach. Learn. Res., 3:463–482, 2002. 5
P. L. Bartlett, S. Boucheron, and G. Lugosi. Model selection and error estimation. Mach. Learn., 48(1-3):
85–113, 2002. 5

13

P. L. Bartlett, O. Bousquet, and S. Mendelson. Local Rademacher complexities. Ann. Statist., 33(4):
1497–1537, 2005. 5, 7
P. L. Bartlett, M. I. Jordan, and J. D. Mcauliffe. Convexity, classification, and risk bounds. J. Am. Statist.
Assoc., 101(473):138–156, 2006. 8, 9
Z. Barutcuoglu, R. E. Schapire., and O. G. Troyanskaya. Hierarchical multi-label prediction of gene
function. Bioinformatics, 22(7):830–836, 2006. 1
R. Bhatia. Matrix analysis, volume 169 of Graduate Texts in Mathematics. Springer-Verlag, New York,
1997. 2
O. Bousquet. A Bennett concentration inequality and its application to suprema of empirical processes.
C. R. Acad. Sci. Paris, Ser. I, 334:495–500, 2002. 5, 7, 16
S. Boyd and L. Vandenberghe. Convex optimization. Cambridge University Press, Cambridge, 2004. 8
P. Bühlmann and S. van de Geer. Statistics for high-dimensional data. Springer Series in Statistics.
Springer, Heidelberg, 2011. Methods, theory and applications. 3
E. J. Candès, J. Romberg, and T. Tao. Robust uncertainty principles: Exact signal reconstruction from
highly incomplete frequency information. IEEE Trans. Inf. Theory, 52(2):489–509, 2006. 2
F. Charte and D. Charte. Working with multilabel datasets in R: The mldr package. The R Journal, 7(2):
149–162, dec 2015. 12
F. Charte, A. J. Rivera, M. J. del Jesus, and F. Herrera. QUINTA: a question tagging assistant to improve
the answering ratio in electronic forums. In EUROCON, pages 1–6, 2015. 12
Y.-N. Chen and H.-T. Lin. Feature-aware label space dimension reduction for multi-label classification.
In NIPS, pages 1529–1537, 2012. 10, 11
W. Cheng, E. Hüllermeier, and K. Dembczynski. Bayes optimal multilabel classification via probabilistic
classifier chains. In ICML, pages 279–286, 2010. 2
R. Combes. An extension of McDiarmid’s inequality. arXiv preprint arXiv:1511.05240, 2015. 6
K. Dembczyński, W. Waegeman, W. Cheng, and E. Hüllermeier. On label dependence and loss minimization in multi-label classification. Mach. Learn., 88(1-2):5–45, 2012. 1
D. L. Donoho. Compressed sensing. IEEE Trans. Inf. Theory, 52(4):1289–1306, 2006. 2
M. Elad. Sparse and redundant representations. From theory to applications in signal and image
processing. New York, NY: Springer, 2010. 3
M. Frank and P. Wolfe. An algorithm for quadratic programming. Naval research logistics quarterly, 3
(1-2):95–110, 1956. 21
Y. Freund and R. E. Schapire. A decision-theoretic generalization of on-line learning and an application
to boosting. Journal of computer and system sciences, 55(1):119–139, 1997. 8
J. Friedman, T. Hastie, and R. R. Tibshirani. Additive logistic regression: a statistical view of boosting.
Ann. Statist., 28(2):337–407, 2000. 8
S. Gao, W. Wu, C. H. Lee, and T. S. Chua. A MFoM learning approach to robust multiclass multi-label
text categorization. In ICML, pages 329–336, 2004. 1
14

D.J. Hsu, M. Sham Kakade, J. Langford, and T. Zhang. Multi-label prediction via compressed sensing.
In NIPS, pages 772–780, 2009. 2, 3
M. Jaggi. Revisiting Frank-Wolfe: Projection-free sparse convex optimization. In ICML, pages 427–435,
2013. 11, 21, 22
M. Jaggi and M. Sulovský. A Simple Algorithm for Nuclear Norm Regularized Problems. ICML, 2010.
11, 21
H. Jain, Y. Prabhu, and M. Varma. Extreme multi-label loss functions for recommendation, tagging,
ranking & other missing label applications. In KDD, pages 935–944, 2016. 2, 3, 4, 5
I. Katakis, G. Tsoumakas, and I. Vlahavas. Multilabel text classification for automated tag suggestion. In
Proc. ECML/PKDD 2008 Discovery Challenge, pages 75–83, 2008. 12
V. Koltchinskii. Rademacher penalties and structural risk minimization. IEEE Trans. Inf. Theory, 47(5):
1902–1914, September 2006. 5
M. Ledoux and M. Talagrand. Probability in Banach spaces. Isoperimetry and processes. Springer-Verlag,
Berlin, 1991. 6, 7
Q. Li, B. Xie, J. You, W. Bian, and D. Tao. Correlated logistic model with elastic net regularization for
multilabel image classification. IEEE Trans. Image Process., 25(8):3801–3813, Aug 2016. 2
X. Li, F. Zhao, and Y. Guo. Multi-label image classification with a probabilistic label enhancement model.
In UAI, pages 430–439, 2014. 1
C. McDiarmid. On the method of bounded differences. In Surveys in Combinatorics, number 141 in
London Mathematical Society Lecture Note Series, pages 148–188. Cambridge University Press, 1989.
5, 16
N. Natarajan and P. Jain. Regret bounds for non-decomposable metrics with missing labels. In NIPS,
pages 2874–2882. 2016. 12
M. Tan, Q. Shi, A. van den Hengel, C. Shen, J. Gao, F. Hu, and Z. Zhang. Learning graph structure for
multi-label image classification via clique generation. In CVPR, pages 4100–4109, 2015. 2
G. Tsoumakas, I. Katakis, and I. Vlahavas. Effective and efficient multilabel classification in domains
with large number of labels. In Proc. ECML/PKDD 2008 Workshop on Mining Multidimensional Data,
2008. 12
G. Tsoumakas, I. Katakis, and I. Vlahavas. Mining Multi-label Data, pages 667–685. Springer US,
Boston, MA, 2010. 1
V. N. Vapnik. Statistical learning theory. Wiley, Wiley New York, 1998. 5, 8
C. Xu, T. Liu, D. Tao, and C. Xu. Local Rademacher complexity for multi-label learning. IEEE Trans.
Image Process., 25(3):1495–1507, 2016. 7, 10, 11
H.-F. Yu, P. Jain, P. Kar, and I. S. Dhillon. Large-scale multi-label learning with missing labels. In ICML,
pages 593–601, 2014. 2, 5, 10, 11, 20
T. Zhang. Statistical behavior and consistency of classification methods based on convex risk minimization.
Ann. Statist., 32(1):56–85, 02 2004. 2, 5, 8, 9

15

A

Appendix

A.1

Generalization bound

Theorem 4 (McDiarmid’s inequality [McDiarmid, 1989]). Let Z1 , . . . , ZN be independent random
variables taking values in Z. Further, let ϕ : Z N 7→ R be a function of Z1 , . . . , ZN that satisfies
∀i ∈ [N ] ∀z1 , . . . , zN , zi0 ∈ Z,
Then for all t > 0,

ϕ(z1 , . . . , zi , . . . , zN ) − ϕ(z1 , . . . , zi0 , . . . , zN ) ≤ ci .


2t 2
P{φ ≤ Eφ + t} ≤ exp  − PN



 .
2

i=1 ci

Proof of Proposition 1. Let ϕ(DN ) = supf ∈F R(f ) − R̂N (f , DN ) and denote by f¯ ∈ F the classifier3
such that ϕ(DN ) = R(f¯) − R̂N (f¯, DN ) . Let D0 , DN be two sets of observations such that D0 and DN
N

N

differ only in one observation (Xk , Yk ) and (Xk0 , Yk0 ), where (Xk0 , Yk0 ) is a copy of (X, Y ) independent from
DN , hence
0
ϕ(DN ) − ϕ(DN
) ≤

L 

p1 X
L1 (f¯l (Xk ))1{Y l =1} − L1 (f¯l (Xk0 ))1{(Y 0 )l =1}
k
k
N
l=1

L 

p0 X
+
L0 (f¯l (Xk ))1{Y l =0} − L0 (f¯l (Xk0 ))1{(Y 0 )l =0}
k
k
N
l=1

L 

p X
L
L1 (f l (Xi ))1{Y l =1} − L1 (f¯l (Xk0 ))1{(Y 0 )l =1} + p0
≤ 1
k
k
N
N
l=1

2K
L
≤
p1 + p0 ,
N
N
PL
where in the last inequality we used that l=1 (L1 (f¯l (Xk ))1{Y l =1} − L1 (f¯l (Xk0 ))1{(Y 0 )l =1} ) consists of at
k
k
most 2K non-zero terms almost surely under Assumption 1.
A.1.1 Refined bound
The following theorem [Bousquet, 2002] is used to prove Lemma 2.
Theorem 5 ([Bousquet, 2002]). Let X1 , . . . , XN be a sequence of independent random variables with
values in some polish space X and distributed according to P. Let F be a P-measurable function from
X N to R. Let Z = F(X1 , . . . , XN ), let Ak be the sigma algebra generated by (X1 , . . . , Xk ) for all k ∈ [N ]
and let Akn be the sigma field generated by (X1 , . . . , Xk−1 , Xk+1 , . . . , XN ) for all k ∈ [N ]. We denote by
0
Ekn [·] the expectation taken conditionally to Akn . Let (Z, Z10 , . . . , ZN
), be a sequence of A-measurable
random variables and let (Zk )k=1,...,N be a sequence of random variables respectively Akn -measurable.
Assume that the following inequalities are satisfied
Zk0 ≤ Z − Zk ≤ 1 a.s., Ekn [Zk0 ] ≥ 0 and Zk0 ≤ 1 a.s.
(24)
PN
Let r > 0 be a real satisfying k=1 Ekn [(Zk0 )2 ] ≤ N r almost surely. If the following condition holds
N
X
k=1

(Z − Zk ) ≤ Z a.s.,

(25)

3 We suppose that the supremum is reached in the definition of ϕ(D ). When this is not the case an  argument would be
N

needed that we leave to the reader.

16

hence for all δ > 0 we obtain


p
log(1/δ)
≤ 1−δ .
P Z ≤ E[Z] + (4E[Z] + N r) log(1/δ) +
3

(26)

The following result is the consequence of Theorem 5.
Theorem 6. Let r > 0 be such that
∀j ∈ {0, 1},

∀f ∈ F ,

E[max L2j (f l (X))] ≤ r ,

(27)

l∈[L]

and
Z = sup R(f ) − R̂N (f ) .
f ∈F

With p0 , p1 as in Eq. (8) hence under Assumption (1) with probability at least 1 − δ we have
r
8K
32K 2 r
4K
Z ≤ E[Z]+ E[Z]
log(1/δ) +
log(1/δ) +
log(1/δ) .
N
N
3N
Proof of Theorem 6. To prove Theorem 6 we need to justify the conditions of Theorem 5, with appropriate
choice of random variables. We define the following variables
Z=

N
sup R(f ) − R̂N (f , DN ) ,
4K f ∈F

(28)

and
Zk =

N
−k
sup R(f ) − R̂N (f , DN
),
4K f ∈F

Zk0 =

N
R(fk ) − R̂N (fk , DN ) − Zk ,
4K

(29)

−k
−k
N
) . We denote by f¯ the
R(fk ) − R̂N (fk , DN
where DN
= DN \ {(Xk , Yk )} and fk is such that Zk = 4K
function for which the supremum is reached in Z. Hence we get

N
N
−k
R(f¯) − R̂N (f¯, DN ) −
R(f¯) − R̂N (f¯, DN
)
4K
4K
L
1 X
p0 L0 (f¯l (Xk ))1{Y l =0} + p1 L1 (f¯l (Xk ))1{Y l =1}
≤
k
k
4K

Zk0 ≤ Z − Zk ≤

l=1

≤ 1 a.s. ,
where in the last inequality we used Assumption 1 and boundedness of L0 and L1 . Moreover we have
Ekn [Zk0 ] ≥

h
i
N
Ekn R(fk ) − R̂N (fk , DN ) − Zk = 0 ,
4K

We also have
(N − 1)Z =

N
N
N
X
N X
N X
−k
−k
(R(f¯) − R̂N (f¯, DN
)) ≤
R(f¯) − R̂N (f¯, DN
) ≤
Zk ,
4K
4K
k=1

k=1

hence
N
X
k=1

(Z − Zk ) ≤ Z ,
17

k=1

and finally we have
N
X
k=1

Ekn [(Zk0 )2 ] ≤

N
X

Ekn

≤

1
8K 2

L
2 
1 X
p0 L0 (fkl (Xk ))1{Y l =0} + p1 L1 (fkl (Xk ))1{Y l =1}
i
k
4K
l=1

k=1

1
≤
16K 2
≤



N
X

k=1
N
X
k=1

Ekn


2 
l
l
p0 L max L0 (fk (Xk )) + Kp1 max L1 (fk (Xk ))
l∈[L]

l∈[L]





Ekn p02 L2 max L20 (fkl (Xk )) + Ekn K 2 p12 max L21 (fkl (Xk ))
l∈[L]

l∈[L]

!

 Nr
Nr  2 2
L2
p02 2 + p12 ≤ N r
p0 L + K 2 p12 =
2
8
8K
K

a.s. ,

and we conclude.
Proof of Lemma 2. It is sufficient to apply the following elementary result to the Theorem 6, to obtain
the proof of Lemma 2.
Lemma 4. For every u, v ≥ 0,

√

√
√
u + v ≤ u + v,
√
2 uv ≤ u + v.

(30)
(31)

Using previous inequality to bound product we can rewrite Equation (28) as follows
r
32r
10
Z ≤ 2E[Z]+K
log(1/δ) + K
log(1/δ) .
N
3N

A.2

(32)

On convexification

Proof of Lemma 3.
Rw (f ) − Rw∗ = EX

L
X

p0 (1 − η l (X))(1{f l (X)≥0} − 1{η l (X)≥p0 } ) + p1 η l (X)(1{f l (X)<0} − 1{η l (X)<p0 } )
|
{z
}
l=1
=1{η l (X)≥p } −1{f l (X)≥0}
0

=

L
X
l=1

≤
≤

L
X

EX (p0 − η l (X))(1{f l (X)≥0} − 1{η l (X)≥p0 } )
EX [1{(η l (X)−p0 )f l (X)≤0} η l (X) − p0 ]

l=1
L 
X

s

EX [1{(η l (X)−p0 )f l (X)≤0} η l (X) − p0 ]

l=1
L 
X

≤c

l=1

1/s

1/s
l
EX [1{(η l (X)−p0 )f l (X)≤0} ∆Rw
.
φ (η , 0)]

(by Jensen’s inequality)
(thanks to Eq. (20))

l
Now it is sufficient to show that for all l ∈ [L] having (η l (X) − p0 )f l (X) ≤ 0 implies ∆Rw
φ (η , 0) ≤

w l
w l l
l l
l
l
∆Rw
φ (η , f ). We only need to prove that (η − p0 )p ≤ 0 implies Rφ (η , 0) ≤ Rφ (η , p ). To see this, we
consider the following cases:

18

• η > p0 : By assumption, we have (fφ∗ )l (η l ) > 0. Moreover, (η − p0 )pl ≤ 0 implies pl ≤ 0. Hence
w l
l
0 ∈ [pl , (fφ∗ )l (η l )] and by convexity of the real function t 7→ Rw
φ (η , t) we have Rφ (η , 0) ≤

w l
w l l
∗ l l
l l
max{Rw
φ (η , p ), Rφ (η , (fφ ) (η ))} = Rφ (η , p ).

• η < p0 : Note that we require (fφ∗ )l (η l ) < 0. Therefore the same analysis as above brings the desired
result.
• η l = p0 : Since (fφ∗ )l (η l ) = 0 we have the desired result.
Hence we have proved the following inequality
Rw (f ) − Rw∗ ≤ c

L 
X

1/s
l l
EX ∆Rw
(η
,
f
)
.
φ

l=1

Then, by concavity of x1/s we conclude that
w

w∗

R (f ) − R

1−1/s

≤ cL

L
X

l l
EX ∆Rw
φ (η , f )

1/s

1−1/s

= cL

l=1


1/s
w
w ∗
Rφ (f ) − Rφ (fφ )
.

Proof for Examples 1. We first notice, that the minimization in Equation (19) is separable, i.e., fφ∗ (η) =

((fφ∗ )1 (η 1 ), . . . , (fφ∗ )L (η L ))> , where for all l ∈ [L]

l l
(fφ∗ )l (η l ) ∈ arg min Rw
φ (η , f ) .

(33)

f

l l
l
• Let φ(v) = (1 − v)2 and fix l ∈ [L], minimizing Rw
φ (η , f ) over f we obtain

(fφ∗ )l (η l ) =

p0 − η l
l
η (1 − 2p0 ) + p0

,

l
substituting to ∆Rw
φ (η , 0) we get
l
∆Rw
φ (η , 0) =

p0 − η l

2

η l (1 − 2p0 ) + p0

,

and finally, since η l (1 − 2p0 ) + p0 ≤ 2 we conclude by
p0 − η

l 2

≤2

p0 − η l

2

η l (1 − 2p0 ) + p0

.

l l
l
• Let φ(v) = e−v and fix l ∈ [L], minimizing Rw
φ (η , f ) over f we obtain

(fφ∗ )l (η l ) =

 l

η (1 − p0 )
1
log
,
2
p0 (1 − η l )

l
substituting to ∆Rw
φ (η , 0) we get
l
∆Rw
φ (η , 0) =

q
q
2
(1 − p0 )η l − p0 (1 − η l ) ,
19

and finally, since

p
p
2
(1 − p0 )η l + p0 (1 − η l ) ≤ 2 we conclude by

p0 − η l

2

=

q
q
q
2 q
2
(1 − p0 )η l − p0 (1 − η l )
(1 − p0 )η l + p0 (1 − η l )

≤2

A.3

q
q
2
(1 − p0 )η l − p0 (1 − η l ) .

Linear prediction approach

Proof of Theorem 3. Our proof follows the analysis presented in [Yu et al., 2014] with slight modification
and use of Assumption 1. It is sufficient to provide bound for
R1 (W ) =

N X
L
X

1
εil f l (Xi )1{Y l =1} ,
EN Eε sup
i
N
f ∈F
i=1 l=1

R0 (W ) =

1
E E sup
N N ε f ∈F

to this end let V1 = [V11 , . . . , V1L ], where V1l =

N
P
i=1

N X
L
X
i=1 l=1

εil f l (Xi )1{Y l =0}



i

,

εil Xi 1{Y l =1} , hence
i

1
λ
EN Eε sup hW , V1 i ≤ EN Eε kV1 kσ ,∞
N
N
W ∈W
v
u
t
q
L
X
λ
λ
2
EN Eε kV1 k2σ ,2 =
EN Eε
V1l 2 ,
≤
N
N

R1 (W ) =

l=1

using expression for V1l we can write
EN Eε

L X
N
X
l=1

i=1

2

εil Xi 1{Y l =1}
i

= EN Eε

L X
N
X
l=1 i=1

2

= EN

N
X
i=1

kXi k22

kXi k22 1{Y l =1} + EN Eε

L X
X

i

L
X
l=1

1{Y l =1} ≤ KEN
i

N
X
i=1

l=1 i,i 0

εil εil 0 1{Y l =1} 1{Y l0 =1} hXi , Xi 0 i
i

i

kXi k22 ≤ KN .

where in the second equality we used the definition of Rademacher variables and in the second last
P
inequality we used the fact that Ll=1 1{Y l =1} ≤ K almost surely under Assumption 1, therefore we have
i

r

R1 (W ) ≤ λ

K
.
N

For the second Rademacher complexity we can write
1
λ
E E sup hW , V0 i ≤ EN Eε kV0 kσ ,∞
N N ε W ∈W
N
v
u
t
q
L
X
λ
λ
2
2
≤
EN Eε kV0 kσ ,2 =
EN Eε
V0l 2 ,
N
N

R0 (W ) =

l=1

20

Algorithm 1 Frank-Wolfe algorithm for minimizing F(X) with a nuclear norm constraint.
Result: Ŵ
k=0
W0 = 0
while Stopping criteria is not achieved do
[uk , sk , vk ] = −svds(−∇F(Wk ), 1)
2
αk = k+2
Wk+1 = (1 − αk ) · Wk + αk · λ · uk vk>
k = k+1
where V0 = [V01 , . . . , V0L ], where V0l =

EN Eε

N
L X
X
l=1

i=1

2

εil Xi

= EN Eε

N
P
i=1

εil Xi 1{Y l =0} , hence
i

L X
N
X
l=1 i=1

2

= EN

N
X
i=1

kXi k22

kXi k22 1{Y l =0} + EN Eε
i

L
X
l=1

1{Y l =0} ≤ LEN
i

N
X
i=1

X
i,i 0

εil εil 0 1{Y l =0} 1{Y l0 =0} hXi , Xi 0 i
i

i

kXi k22 ≤ LN ,

Hence we have
r

K
,
R1 (W ) ≤ λ
N
r
L
R0 (W ) ≤ λ
.
N

(34)
(35)

We conclude by substituting p0 , p1 and Equation (34) into Theorem 1 and we get with probability at least
1−δ
r
n
o
K
K p
K
) + 4√
log 1/δ .
(36)
sup R(W ) − R̂N (W ) ≤ 2Cλ(
+√
N
NL
N
W ∈W

A.4

Frank-Wolfe algorithm

Let F : RD×L 7→ R be convex function, hence we would like to solve the following minimization problem
min

W ∈RD×L

s.t.

F(W )
kW kσ ,1 ≤ λ.

One of the possible approach is to consider the Frank-Wolfe4 algorithm [Frank and Wolfe, 1956] (see
also [Jaggi, 2013, Jaggi and Sulovský, 2010] for applications to similar matrix optimization problems).
This algorithm is attractive due to its cheap iteration cost, i.e., on each iteration we need to evaluate
only the top singular value of the gradient of F (in contrast to forward-backward variants that would
require a full SVD evaluation at each iteration). The Frank-Wolfe algorithm applied to Problem (A.4) is
summarized in Algorithm 1. The function svds(−∇F(Wk ), 1) computes the largest singular value and the
corresponding singular vectors.
4 also referred to as conditional gradient

21

As a stopping criterion we follow the one proposed in [Jaggi, 2013] that can be specified in our
context as:
|hWk , ∇F(Wk )i + λ · sk )| ≤ ,

(37)

for some  > 0 This duality gap criterion ensures that F(Wk ) − F(Ŵ ) ≤  for the output of the algorithm.

A.5

Synthetic data

First, we remind that the sigmoid function is given by:
S(t) =

1
.
1 + e−t

(38)

We consider the following framework to simulate synthetic datasets with controlled sparsity constant K.
Each entry of the design matrix X ∈ RN ×D is simulated according to a uniform distribution on [−1, 1].
We form a matrix of coefficients B ∈ RL×D , in which every entry is either 2 or −2 with probability 0.5.
We transform the matrix Z = XB > into probabilities by applying the sigmoid function Eq. (38) to each
entry. Having the matrix Z of probabilities we pick the K/2 highest values in each row of Z and set
them to 1. Moreover, we set the L − K lowest values in each row of Z to 0. Now, we have a matrix with
K/2 ones, K/2 value between zero and one (probabilities) and L − K zeros in each row. Each element of
the label matrix Y is obtained as Bernoulli random variable with probability from the described matrix,
see Equation (39).


l

1, with probability Zi
l
Yi = 
.
(39)

0, with probability 1 − Zil
This framework, allows to obtain Y , in which rows are at least K/2 sparse and at most K sparse.

22

