Sélection automatique de blocs de variables
Pierre Alquier
Laboratoire de Probabilités et Modèles AléatoiresUniversité Paris 7, site Chevaleret, case 7012
175, rue du Chevaleret75205 Paris Cedex 13 - FRANCE&Laboratoire de Statistique du CREST
Timbre J350
3, Avenue Pierre Larousse
92240 Malakoff - FRANCE
Résumé: En génomique apparaissent certains modèles de classification et de régression
avec un paramètre β parcimonieux (c’est-à-dire que la plupart des coordonnées βi sontnulles1 ) et constant par morceaux, voir par exemple les articles de Rapaport, Barillot et Vert (2008) et Huang, Salim, Lei, O’Sullivan et Pawitan (2009). Des méthodesd’estimation ont été proposées dans ce contexte, par exemple le fused LASSO de Tibshirani, Saunders, Rosset, Zhu et Knight (2005). L’idée de cet estimateur est d’adapter lapénalité dans l’estimateur LASSO de Tibshirani (1996) en forçant βi à être égal à βi+1
dans beaucoup de cas, ce qui permet de sélectionner des groupes de variables sans lesconnaı̂tre a priori. Dans l’article Alquier (2010), on a proposé une méthode dont l’idéeest d’utiliser une adaptation de la forme duale du LASSO, donnée dans Osborne, Presnell
et Turlach (2000). On présente ici brièvement cette méthode, puis on donne quelquesrésultats d’application à des données CGH.
Abstract: In genomics applications, it makes sense to consider classification or regression
estimation with a parameter β “sparse” (most of its coordinates are 0) and “blocky” (βiand βi+1 are likely to be equal). We refer the reader for example to the paper of Rapaport,
Barillot and Vert (2008) or Huang, Salim, Lei, O’Sullivan and Pawitan (2009). Someestimation methods taking this information into account are available, such as the fused
LASSO, by Tibshirani, Saunders, Rosset, Zhu, and Knight (2005). This method relies onan adaptation of the penalty in the well known LASSO, defined by Tibshirani (1996), toselection of unknown blocks of features. In the paper Alquier (2010), we proposed anotherestimation method, that relies on an adapatation of the dual form of the LASSO providedby Osborne, Presnell and Turlach (2000). We describe shortly this method, and provide
some applications to array-CGH data.
Mots-clés: sélection de variables, modèles parcimonieux, sparsité, régression linéaire,groupes de variables, tableaux CGH.
1

Il est devenu courant d’utiliser l’anglicisme sparse dans ce cas, et sparsité pour parcimonie.

11Introduction

On se place dans le modèle de régression linéaire gaussienneY ∼ N Xβ, σ 2 In

(1)où pour la simplicité de la présentation, X est une matrice déterministe n × p, et β =(β1 , ..., βp )0 avec possiblement p > n (on utilise la convention que si u ∈ Rp , u est unvecteur colonne et u0 sa transposée). On note X1 , ..., Xp les colonnes de X, et P la loi
de Y donnée par (1). On suppose que les données sont renormalisées de façon à avoir
Xj0 Xj /n = 1 pour tout j. Soit le critère d’erreurL(b) = kX(b − β)k22 .
Sans hypothèses supplémentaires, il est impossible d’estimer correctement β. Une
hypothèse courante est de supposer que le nombre de coordonnées non-nulles dans β est
petit (hypothèse de sparsité). Dans ce cas, l’estimateur LASSO"#pX
β̂LASSO (λ) = arg min kY − Xbk22 + λ|bj |bj=1

pour un λ > 0, défini par Tibshirani (1996), ou le Dantzig selector, de Candès et Tao(2007), ont de bonnes performances en pratique et sont bien compris en théorie - voir
par exemple Candès et Tao (2007) et Bickel, Ritov et Tsybakov (2009) pour les résultatsthéoriques. La sparsité a de plus souvent un sens pour le praticien, en terme de sélection
de variables.
En génomique, dans certains modèles, cette hypothèse est naturelle si elle est combinéeavec l’existence de “blocs de coefficients” dans β, c’est-à-dire que le nombre de j tels que
βj 6= βj+1 est aussi petit. Par exemple,
β = (0, 0, 0, 0, 0, 0, 0, 0, 5, 5, 5, 5, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 2, 2, 0, 0, 0, 0)0 .
Ceci est par exemple le cas dans les données de type CGH, où les index j sont des positions
de “sondes” le long de chromosomes. On observe la quantité d’ARN correspondant àchaque sonde de façon très bruitée: Yj = βj + εj (dans ce cas, X = In ). Dans des cellules
cancéreuses, des portions entières de chromosomes peuvent subir des amplifications ou
des délétions, ce qui se traduira respectivement par une augmentation où une diminution
de toute une série de βj consécutifs. Pour un chromosome d’une cellule saine, βj = 0 pour
tout j.
Pour traiter ce cas, adapter la pénalité du LASSO de la façon suivante semble naturel:"#ppXXβ̂F −LASSO (λ, µ) = arg min kY − Xbk22 + λ|bj | + µ|bj − bj−1 |
b

j=1

2

j=2

pour λ, µ > 0. Cet estimateur a été introduit par Tibshirani, Saunders, Rosset, Zhu, etKnight (2005) sous le nom de fused LASSO. On pourra cependant remarquer la ressemblance avec les méthodes de débruitage d’images utilisant la variation totale, cf. Mumford
et Shah (1989) ou Rudin, Osher et Fatemi (1992). Dans des articles de bioinformatiqueet de biologie, de telles pénalisations ont effectivement été utilisées, cf. par exemple lesarticles de Rapaport, Barillot et Vert (2008) et Huang, Salim, Lei, O’Sullivan et Pawitan (2009).
On peut comparer par exemple le fused LASSO au group LASSO introduit par Yuan etLin (2006): contrairement au fused LASSO, le group LASSO ne peut fonctionner si on neconnait pas a priori la position des groupes de variables. La même remarque s’appliqueá la méthode de Zhao, Rocha et Yu (2009). En contrepartie, les garanties théoriques
sur le group LASSO sont plus générales, voir Chesnau et Hebiri (2008). Dans le cas
du fused LASSO, les seuls résultats disponibles sont dans le cas X = In , prouvés par
Rinaldo (2009).
En pratique, le fused LASSO peut être calculé rapidement par une optimization coordonnée par coordonnée, cf. Friedman, Hastie, Höfling et Tibshirani (2007). Plusrécemment, des algorithmes capables de calculer rapidement toutes les solutions (pourtoutes les valeurs de λ et µ) sont apparus. Toutes reposent sur une généralisation del’algorithme LARS de Efron, Hastie, Johnstone, et Tibshirani (2004): Hoefling (2009) ou
Tibshirani et Taylor (2009). Toujours en pratique, on pourra aussi considérer le smooth
LASSO de Hebiri (2008), étudié en détail par Hebiri et Van de Geer (2010), comme unealternative au fused LASSO, qui, bien que conçu pour un objectif légèrement différent,est souvent plus stable numériquement:"#ppXXβ̂F −LASSO (λ, µ) = arg min kY − Xbk22 + λ|bj | + µ(bj − bj−1 )2 .
bj=1

j=2

Cet estimateur peut aussi être vu comme un cas particulier de l’estimateur strictured
elastic net de Slawski, zu Castell et Tutz (2010).
Dans l’article Alquier (2010), on propose un estimateur alternatif dans ce contexte.
Celui-ci ne repose pas sur une modification de la pénalisation du LASSO, mais de laforme duale du LASSO donnée par Osborne, Presnell et Turlach (2000). On présente ici
cet algorithme (Section 2), ensuite, on rappele un résultats théoriques sur cet algorithme(Section 3). On donne enfin une illustration sur des données réelles.

2

L’algorithme

L’algorithme est itératif: il part de β̂ (0) = (0, . . . , 0) et à chaque étape, il met à jour ungroupe de coordonnées. On décrit maintenant cette mise à jour. On suppose donc queβ̂ (m) est donné et on va définir β̂ (m+1) .
3On commence par quelques notations. Soit la fonction de seuillage “doux” ∀x ∈R, ∀u ≥ 0, γ(x, u) = sign(x)(|x| − u)+ . Pour tout j ∈ {1, ..., p} et k ∈ {1, ..., p − j + 1} ondéfinit 1j,k ∈ Rp par (1j,k )i = 1 si i ∈ {j, ..., j + k − 1}, et (1j,k )i = 0 sinon, autrement
dit si i appartient au groupe de coordonnées de longueur k, commençant à la position j.
Choisissons un s > 0 (on discute ce choix dans la section suivante).
Pour tout groupe défini par la position initiale j et la longueur k, on définit le vecteur(m+1),j,kβ̃
par:
Pj+k−1(Y −X β̂ (m) )0 h=jXh(m+1),j,k(m)s, √10 X 0 X1β̃j= β̂j + γ
10j,k X 0 X1j,k
j,kj,k
..
(2).



P

j+k−1

(Y −X β̂ (m) )0 h=j Xh(m+1),j,k(m), √10 Xs 0 X1
 β̃j+k−1 = β̂j+k−1 + γ
10j,k X 0 X1j,kj,kj,k(m+1),j,k(m)et β̃`= β̂` pour tout ` ∈/ {j, ..., j + k − 1}.
On doit ensuite choisir un couple (j, k), on pose simplement:β̂ m = arg max kX(β̃ (m+1),j,k − β̂ (m) )k22 .
j,k3

Un premier résultat théorique

Le résultat suivant, prouvé dans Alquier (2010), constitue un premier pas dans l’analysedes propriétés statistiques de cet algorithme. Il garantit que l’on peut itérer autant defois que l’on veut l’étape de mise à jour sans risquer le surapprentissage.
Théorème: On a, pour tout M ∈ N,"#s2

P L(β̂ (M ) ) ≤ L(β̂ (M −1) ) ≤ · · · ≤ L(β̂ (0) ) ≥ 1 − p2 e− 2σ2 .

Ceci signifie que si l’on choisit s = 2σ[log(p/ε)]1/2 , alors avec probabilité au moins
1 − ε, à chaque étape le risque de l’estimateur diminue.
La preuve repose sur les considérations géométriques sur le LASSO donnés dans Osborne, Presnell et Turlach (2000) et Alquier (2008a,2008b).

4

Test sur des données simuléesPour des raisons de place, on ne présente pas ici en détail les simulations effectuées dansAlquier (2010). Notons que si cette étude expérimentale est assez brève, il en ressort
que dans certains cas, les performances de l’algorithme de sélection itétarive de blocs devariables présenté ici sont une réelle amélioration par rapport à celles du fused LASSO.
45Applications sur des données CGH

La Figure 1 montre un exemple d’itération sur des données réelles: le chromosome 17 d’unpatient atteint d’un type particulier de cancer. On est ici dans le cas particulier présenté
en introduction Yi = βi + εi , 1 ≤ i ≤ n = 7728.
Figure 1: Un exemple d’itération de l’algorithme.

Tous les calculs ont été effectués avec le logiciel R, cf. R Development Core Team(2008). Les codes sont disponibles auprès de l’auteur. Les données ont été fournies parle Dr. Jean-Paul Feugeas (INSERM U944, Institut Universitaire d’Hématologie) que jetiens a remercier particulièrement pour avoir motivé ce travail.

Bibliographie[1] Alquier, P. (2008a), Iterative feature selection in least square regression estimation,
Annales de l’IHP B (Probas. Stats.), 44(1), 47–88.
[2] Alquier, P. (2008b), LASSO, iterative feature selection and the correlation selector:
oracle inequalities and numerical perfomances, Elec. Journ. of Statistics, 2, 1129–1152.
[3] Alquier, P. (2010), An algorithm for iterative selection of blocks of features, Proceedingsof ALT’10, M. Hutter, F. Stephan, V. Vovk & T. Zeugmann Edts., LNAI, Springer, 35–49.
[4] Bickel, P., Ritov, Y. and Tsybakov, A. (2009), Simultaneous analysis of LASSO andDantzig selector, The Annals of Statistics, 37(4), 1705–1732.
[5] Candès, E. and Tao, T. (2007), The Dantzig selector: statistical estimation when p ismuch larger than n, The Annals of Statistics, 35(6), 2313–2351.
5

[6] Chesnau, C. and Hebiri, M. (2008), Some theoretical results on the grouped variables
LASSO, Mathematical Methods in Statistics, 17(4), 317–326.
[7] Efron, B., Hastie, T., Johnstone, I. and Tibshirani, R. (2004), Least angle regression,
The Annals of Statistics, 32(2), 407–499.
[8] Friedman, J., Hastie, T., Höfling, H. and Tibshirani, R. (2007), Pathwise coordinate
optimization, The Annals of Applied Statistics, 1(2), 302–332.
[9] Hebiri, M. (2008), Regularization with the smooth LASSO procedure, Preprint, arXiv:0803.0668.
[10] Hebiri, M. and Van de Geer, S. (2010), The smooth LASSO and other `1 +`2 -penalizedmethods, Preprint, arXiv: 1003.4885.
[11] Hoefling, H. (2009), A path algorithm for the fused LASSO signal approximator,Preprint, arXiv: 0910.0526.
[12] Huang, J., Salim, A., Lei, K., O’Sullivan, K. and Pawitan, Y. (2009), Classification ofarray CGH data using smoothed logistic regression model, Statistics in Medicine, 8(30),3798–3810.
[13] Mumford, D. and Shah, J. (1989), Optimal approximations by piecewise smooth functions and associated variational problems, Communications on Pure and Applied Mathematics, 42(5), 57–685.
[14] Osborne, M., Presnell, B. and Turlach, B. (2000), On the LASSO and its dual, Journalof Computational and Graphical Statistics, 9(2), 319–337.
[15] R Development Core Team (2008), R: A Language and Environment for StatisticalComputing, R Foundation for Statistical Computing, Vienna, http://www.R-project.org[16] Rapaport, F., Barillot, E. and Vert, J.-P. (2008), Classification of array-CGH datausing fused SVM, Bioinformatics, 24(13), 1375–1382.
[17] Rinaldo, A. (2009), Properties and Refinements of the fused LASSO, The Annals ofStatistics, 37(5B), 2922–2952.
[18] Rudin, L., Osher, S. and Fatemi, E. (1992), Nonlinear total variation based noise
removal algorithm, Physica D (Nonlinear Phen.), 60, 259–268.
[19] Slawski, M., zu Castell, W. and Tutz, G. (2010), Feature selection guided by structuralinformation, The Annals of Applied Statistics, 4(2), 1056–1080.
[20] Tibshirani, R. (1996), Regression shrinkage and selection via the LASSO, JRSS B(Method.), 58(1), 267–288.
[21] Tibshirani, R., Saunders, M., Rosset, S., Zhu, J. and Knight, K. (2005), Sparsity andsmoothness via the fused LASSO, JRSS B (Method.), 67(1), 91–108.
[22] Tibshirani, R. J. and Taylor, J. (2009), Regularization path for least squares problems
with generalized `1 penalties, Preprint.
[23] Yuan, M. and Lin, Y. (2006), Model selection and estimation in regression withgrouped variables, JRSS B (Method.), 68(1), 49–67.
[24] Zhao, P. and Rocha, G. and Yu, B. (2009), The composite absolute penalties forgrouped and hierarchical variable selection, The Annals of Statistics, 37(6A), 3468–3497.

6