Multi-Task Feature Selection on Multiple Networks via Maximum Flows∗†
Mahito Sugiyama‡§

Chloé-Agathe Azencott‡Yoshinobu Kawahara∥

Dominik Grimm‡¶

Karsten M. Borgwardt‡¶

Abstract
We propose a new formulation of multi-task feature selection coupled with multiple network regularizers, and showthat the problem can be exactly and eﬃciently solved by
maximum flow algorithms. This method contributes to oneof the central topics in data mining: How to exploit structural information in multivariate data analysis, which has
numerous applications, such as gene regulatory and social
network analysis. On simulated data, we show that the proposed method leads to higher accuracy in discovering causal
features by solving multiple tasks simultaneously using networks over features. Moreover, we apply the method tomulti-locus association mapping with Arabidopsis thaliana
genotypes and flowering time phenotypes, and demonstrate
its ability to recover more known phenotype-related genes
than other state-of-the-art methods.
Keywords: Multi-task learning; Feature selection; Network
regularizer; Multi-locus association mapping

1 Introduction
Knowledge discovery from structured data is one of thecentral topics in data mining. In particular, networks,
or graphs, have attracted considerable attention in thecommunity, as they may represent molecular, biological,
social, or other types of systems whose functionality andmechanisms are far from being completely understood.
Large amounts of data are now available as networks
across a wide range of domains, from biological pathways in KEGG to chemical compounds in PubChemand social networks on the web.
A crucial concern when studying such systems is
to determine which part of the network is responsible
for performing a particular function. For instance,∗ Supplementary information is available at

http://webdav.tuebingen.mpg.de/u/karsten/
Veroeﬀentlichungen/papers/sdm14supp.pdf† M.S. is funded by the Alexander von Humboldt Foundation.
‡ Machine Learning and Computational Biology ResearchGroup, Max Planck Institutes Tübingen, Germany§ Contact: mahito.sugiyama@tuebingen.mpg.de¶ Zentrum für Bioinformatik, Eberhard Karls UniversitätTübingen, Germany∥ The Institute of Scientific and Industrial Research (ISIR),
Osaka University, Japan

medical researchers are interested in determining thesubset of proteins in an interaction network that governs
the response to a particular treatment. Structural
biologists seek to find out which part of the 3D structureof a protein, which can be modeled by a graph, iscorrelated with a particular function. Neuroscientists
search for subgraphs in brain connectivity networks
from functional MRI screens that correlate with certain
types of behavior or cognitive tasks. Geneticists areinterested in sets of mutations in interacting genes that
may be associated with heritable diseases.
Hence the general problem of feature selection on
networks is of broad interest across disciplines. In
these networks, features coincide with vertices (nodes)
and the network topology can be viewed as a priori
knowledge about relationships between features.
The common approach to this problem is to use
Lasso-based regression [28] with an ℓ1 -regularizer of theweight vector and additional structured regularizers that
represent relationships between features. Examples include supervised [3, 14, 15] and unsupervised [26] learning methods and applications in computer vision [9] andstatistical genetics [6, 18, 19]. Many of those methods
have also been proposed for multi-task learning, where
the commonalities between related tasks are leveraged
to improve the quality of models on individual tasks.
In spite of their success, we see a number of drawbacks to regression-based approaches in this context.
First, they do not scale to millions or even hundredsof thousands of features, although such a setting is common, for instance, in genetics. Second, regression-based
approaches concentrate on optimizing a prediction loss,
while the problem to solve is often formulated in termsof finding features that are relevant for, correlated to or
associated with a property of interest.
These two issues have been addressed by our recent
work in statistical genetics, which proposes a new
formulation of network-constrained feature selection
called SConES [2]. This method directly maximizes ascore of association rather than minimizing a prediction
error. Its optimization scheme is exact and eﬃcient,
thanks to a minimum-cut reformulation, and it has been
empirically shown to recover more causal features than

199

Copyright © SIAM.
Unauthorized reproduction of this article is prohibited.

as multiple tasks to increase the accuracy in retrieving
its regression-based counterparts.
An additional issue arises in multi-task settings: phenotype-related SNPs. We confirm in this paper the
Most current methods [16, 32, 33] assume that the same eﬃcacy of our approach compared to the state-of-thefeatures should be selected across all tasks. While this is art methods in feature selection on Arabidopsis thaliana
reasonable for some application domains, one can think SNPs and flowering time phenotypes.
This paper is organized as follows: We present ourof numerous examples where this assumption is violated.
For instance, lung diseases such as asthma and chronic approach to multi-task feature selection across several
obstructive pulmonary disease may be linked to a set of networks in Section 2, discuss related work in Section 3,
common mutations, but there is no indication that the evaluate our method on synthetic and real data in
Section 4, and summarize our contribution in Section 5.
exact same mutations are causal in both diseases.
Moreover, to the best of our knowledge, none of themulti-task approaches incorporating structured regular- 2 Feature Selection on Networks
izers make it possible to consider diﬀerent structural We directly perform feature selection on networks. Forconstraints for diﬀerent tasks. However, we may want mally, we only require, for each task, a network over theto consider diﬀerent biological pathways for diﬀerent dis- features, that is, a weighted graph G = (V, E) with aeases, or to highlight diﬀerent parts of brain connectivity set of vertices (features) V and edges E, and a function
networks for diﬀerent correlated behaviors.
q : V → R that assigns to each vertex v a quantity q(v),To address these two issues, we propose a new measuring its relevance for the problem at hand.
formulation of SConES for multi-task feature selection
This is diﬀerent from the typical setting of multivaricoupled with multiple network regularizers to improve ate data analysis, where a design matrix X ∈ RN ×|V |feature selection in each task by combining and solving and a response vector y ∈ RN about N individuals are
multiple tasks simultaneously. The key to extending given. In such a framework, feature selection is usuSConES to multiple tasks, which is the main technical ally solved as a regularized linear regression problem,
contribution in this paper, is the unification of multiple where one tries to determine a subset of features of X
networks into a single network. This strategy enables us which minimizes a prediction error of y. In this conto solve the multiple tasks as a single task, and hence we text a relevance score q(v) can be easily obtained bystill obtain the exact solution through a minimum cut measuring the association between y and each featurereformulation. Thanks to the eﬃciency of the maximum of X. Many techniques are available for that purpose:
flow algorithm that solves this problem, our approach Pearson’s correlation coeﬃcient or the cosine similarity
is still tractable for large-scale networks. To the best for linear associations, and the Hilbert-Schmidt indepenof our knowledge, this is the first non-regression-based dence criterion (HSIC) [11], sequence kernel association
formulation of multi-task feature selection on multiple test (SKAT) [29], or maximum information coeﬃcientfeature space networks.
(MIC) [25] for non-linear associations.
The motivating application underlying this paper
is to give an eﬃcient method for multi-locus associa- 2.1 Single Task. We first introduce the single-tasktion mapping in genetics. This problem, which aims at formulation of SConES [2] for feature selection. Let
explaining the genetic basis of diseases and other ob- f : 2V → R be additive in the sense that it is defined asserved traits, is receiving growing interest in the con∑f (S) := v∈S q(v).
text of genome-wide association studies (GWAS) [21]. (2.1)
The goal of GWAS is to find single-nucleotide polyThis function measures the goodness of a subset S ofmorphisms (SNPs), single positions that diﬀer in thefeatures via q(v) for each feature v ∈ S.
genomes of diﬀerent individuals, which are significantly
Our goal is to find a subset S ⊂ V which maximizes
associated with variance in phenotype (diseases or otherf (S) under the constraints that the cardinality of S
observed traits). Since SNPs can be viewed as features,is small and its elements tend to be connected to one
feature selection techniques have been widely developed
another. As conducting an exhaustive search over all
and applied to this problem (e.g. [1]). The proposedconnected subnetworks is not feasible, we formulate themethod holds the following advantages over the typicalproblem as follows by focusing on local connectivity:
GWAS setting: First, the method can exploit a priori biological knowledge about network structures over (2.2)
argmax f (S) − g(S),S⊂V
SNPs, derived for example from protein-protein inter∑g(S):= λ e∈B w(e) + η|S|,
action networks. Second, we often have not only one
but several related phenotypes for each set of SNPs,
where B = { {v, u} ∈ E | v ∈ V \ S, u ∈ S } is the set
and hence our method can treat them simultaneouslyof edges located at the boundary of S and w : E → R+

200

Copyright © SIAM.
Unauthorized reproduction of this article is prohibited.

is a weighting function. The first term in the penalty
function g enforces the connectivity of S, as it penalizes
selecting a vertex without selecting all of its neighbors.
The second term enforces its sparsity and the cardinalityof S is penalized. The two real-valued parameters λ andη control these constraints.
A major advantage of this formulation is that
Equation (2.2) can be exactly solved by maximum
flow algorithms by adding source and sink nodes to
the given network G (Supplementary Note A or [2]).
The smallest known time complexity of these algorithms is O( |V ||E| log(|V |2 /|E|) ) [10] and the BoykovKolmogorov algorithm [4] is more eﬃcient in practice.
Regularization Path. An interesting property of
the regularization parameter η, which was not analyzed
in [2], is its anti-monotonicity with respect to thenumber of selected features. Specifically, if we denotethe selected features for each η by S(η), we have S(η) ⊂S(η ′ ) if and only if η > η ′ . Moreover, we can easily
check that our formulation satisfies all assumptions
to apply the parametric maximum flow algorithm [8](Supplementary Note B). With this algorithm, we canobtain the entire regularization path [13] along with thechanges in η without increasing the time complexity.
In practice, this property of η is particularly interesting when we are given cardinality constraints a prioriover the size of the set of selected features. Then we candirectly pick from the regularization path the solutions
that fulfill these constraints.
2.2 Multiple Tasks. Our main contribution isa new, generalized formulation of SConES (Equation (2.2)) to achieve feature selection for multiple tasks
simultaneously. In what follows, we assume that the setof vertices (features) V is shared all over K tasks, andfor each task i we have a network Gi = (V, Ei ) associated with a respective scoring function qi . Given such a
set of K networks G = {G1 , G2 , . . . , GK }, the multi-task
feature selection is formulated as

to select the same features across tasks.
Here we show that this problem can be reduced
to a single-task feature selection similar to that ofEquation (2.2) and thus can also benefit from maximum
flow algorithms. We show an example for K = 2
in Figure 1a. First we replicate the vertices of each
network Gi so that all sets of vertices are disjoint, that
is, G′i = (Vi′ , Ei′ ) such that Vi′ ∩ Vj′ = ∅ for every
i, j ∈ {1, . . . , K} with i ̸= j. All edges are copied on thereplicated set Vi′ and assume that vertices are indexedfrom 1 to n in each network Gi , where vertices have thesame index if they are identical in the original set V .
The m-th vertex of a network Gi is denoted by vim . We
then construct a unified network U (G) = (Ṽ , Ẽ) fromthe set of K networks G = {G1 , . . . , GK } by connecting
each pair of replicated vertices in the following manner:∪K∪n∪K
Ṽ := i=1 Vi′ , Ẽ := i=1 Ei′ ∪ m=1 Am , where{ m m}
Am := {vi , vj } i, j ∈ {1, . . . , K}, i ̸= j .
The weight w̃ of edges is given as w̃(e) = wi (e) if e ∈ Ei′
and w̃(e) = µ/λ otherwise. Thus U (G) has |Ṽ | = Kn∑K
vertices and |Ẽ| = i=1 |Ei | + nK(K − 1)/2 edges.
Theorem 2.1. Given a set of K networks G ={G1 , . . . , GK }. For every subset S ⊂ Ṽ in the unified
network U (G), we havef (S) =K∑
i=1

fi (Si ), g(S) =

K
∑
i=1

gi (Si ) +∑h(Si , Sj ),i<jwhere f, g are defined over U (G) as in Eq. (2.1), (2.2).

Proof. Suppose that K = 2 for simplicity. We can easily
generalize the following proof for the case of K > 2 byconsidering sums over pairs.
The first equation f (S) = f1 (S1 ) + f2 (S2 ) directly
follows from the definition of S1 and S2 . Thus we
focus on the second equation g(S) = g1 (S1 ) + g2 (S2 ) +h(S1 , S2 ). Since V1′ ∩ V2′ = ∅ and therefore
|S| =
∑
|S∑
|
+
|S
|,
all
we
have
to
prove
is
λ
w(e)
=
1
2
e∈B
∑
λ e∈B1 w1 (e) + λ e∈B2 w2 (e) + µ|S1 △ S2 |. From
K
∑
) ∑
(
the definition of the unified network, we have B =(2.3) argmaxh(Si , Sj ),fi (Si ) − gi (Si ) −
′
′
′
′
B
1 ∪ B2 ∪ BA , where Bi = { {v, u} ∈ Ei | v ∈S1 ,...,SK ⊂V i=1i<j′
Vi \ S, u ∈ S }, which is a set of edges on the∑∑qi (v), gi (Si ) := λwi (e) + η|Si |.
fi (Si ) :=unified network corresponding to Bi , and BA is given
v∈Si
e∈Bi
as B = { {vim , vjm } | vim ∈ S and vjm ̸∈ S, i ̸= j } ⊂∪n A
m=1 Am . We show an example of B in Figure 1b.
We introduce a new penalty function h : 2V × 2V → We
see that the cardinality of BA is the same as that
R defined as h(S, S ′ ) := µ| S △ S ′ |, where µ is a of the set S △ S . Thus we have λ ∑12e∈B w(e) =∑real-valued regularization parameter and S △ S ′ is the λ( ∑
w
(e)
+
w
(e)
+
(µ/λ)|B□
1
2
A | ).
e∈B1′
e∈B2′
symmetric diﬀerence between them, that is, S △ S ′ =(S ∪S ′ )\(S ∩S ′ ). The penalty function h represents our The multiple task problem in (2.3) is therefore exactly
belief that similar networks should be associated with equivalent to the single task feature selection probrelated features, and the larger µ, the more we enforce lem (2.2) over the unified network U (G), and can therethis belief. A large µ is thus better when it is desirable fore be solved directly by the maximum flow algorithm.

201

Copyright © SIAM.
Unauthorized reproduction of this article is prohibited.

a

b

3

S

B
G

BA
G
B

Figure 1: (a) Example of two networks (left) whichshare vertices and have diﬀerent edges (solid and dottedlines), and the unified network (right), where vertices
are duplicated and new edges (dotted lines) are added.
(b) Example of two duplicated networks G′1 and G′2 .
Note that this formulation applies even if somefeatures are missing for some task. In that case Vi′
contains only vertices corresponding to the features
available for task i and A contains only edges {vim , vjm },where the feature m is available for both tasks i and j.
2.3 Spectral Analysis. We describe our formulation from a graph spectral point of view. We denoteby f ∈ {0, 1}|V | the indicator vector of a subset S ⊂ V :fv is set to 1 if v ∈ S and 0 otherwise, and denote byc ∈ R|V | the vector composed of values q(v).
The function f in Equation (2.1) can be rewritten
as f (S) = cT f . Let L be the Laplacian matrix ofa network G. Then g(S) = λf T Lf + η∥f ∥0 holds
for the penalty
term in Equation (2.2) since we have∑f T Lf = {v,u}∈E (fv − fu )2 if all edge weights are 1,
where (fv − fu )2 is the XOR of fv and fu and coincides
with the cardinality of B in Equation (2.2). The
same relationship holds if the edges are weighted. The
Laplacian graph regularizer f T Lf is often used in theliterature for penalizing disconnected features (e.g. [18]).
Thus we can state the formulation of single-task feature
selection as: argmaxf ∈{0,1}|V | cT f − λf T Lf − η∥f ∥0 .
In the multi-task framework, the formulation isnaturally extended to)∑K (Targmax i=1 cT
i fi − λfi Li fi − η∥fi ∥0
f1 ,...,fK−∑2
i<j µ∥fi − fj ∥2 ,∑where ∥f − f ′ ∥22 = v∈V (fv − fv′ )2 . Suppose that f˜and c̃ are concatenations of the vectors f∑
1 . . . , fK andK
c1 . . . , cK , respectively. We directly have i=1 ( cTi fi −T ˜˜η∥fi ∥0 ) = c̃ f − η∥f ∥0 . Moreover, we can prove that
∑
∑K
T
2
˜T ˜
i=1 λfi Li fi +i<j µ∥fi − fj ∥2 = λf L̃f , whereL̃ is the Laplacian matrix of the unified network U (G).
Thus this becomes single-task feature selection:argmaxf˜∈{0,1}K|V | c̃T f˜ − λf˜T L̃f˜ − η∥f˜∥0 .

Related Work

Current work on feature selection with network information typically concentrates on regularized linear regression (an excellent overview is given by [30]).
In a single-task context, the Lasso regression
model [28] minimizes the prediction error together withthe ℓ1 -norm of the regression parameter vector, which
encourages sparse solutions. This fact popularized theuse of Lasso for feature selection. Note however thatthe ℓ1 -norm is a convex relaxation of the cardinality
constraint one would really want to enforce.
Group Lasso [31] partitions features into groups
and encourages to select entire such groups via anℓ1 /ℓ2 penalty. If a graph over the features is available,groups can be defined as pairs of connected features [15].
However, the number of groups becomes prohibitively
massive on large-scale networks, which grows exponentially once one starts to consider connected subsets ofhigher cardinality. A popular instance of Group Lasso
is Elastic Net [34], in which all features belong to asingle group; Elastic Net is particularly suited when thenumber of features is larger than that of samples andwhen several correlated features should be selected.
Grace (Graph-constrained estimation) [18, 19] adds
a Laplacian graph regularizer analogous to that of
SConES to the objective. While our method is anassociation-based approach, Grace aims at minimizinga prediction error in a Lasso framework. In practice, it
is at least one order of magnitude slower than SConES
on a network of the same size [2]. aGrace (adaptive
Grace) [19] and GOSCAR [30] employ refined types of
network regularizers which allow connected features tohave eﬀects of opposite directions.
Several multi-task versions of Lasso, in which related tasks are coupled with each other, have been
proposed as well: Multi-Task Lasso [23] uses an ℓ2 norm on each weight across all tasks to reward solutions, where the same features are selected for all tasks.
Graph-Guided Fused Lasso [16] extends this idea bycoupling the weight vectors of correlated tasks: the morecorrelated two tasks are, the more solutions in which
they have similar weight vectors are rewarded. Recent
developments combine multi-task learning with structured regularization on the input features (e.g. [6, 17]).
Finally, [7] incorporates network-structured feature selection to multi-task learning. This Lasso approach has
the advantage to integrate task covariance, which weleave to future work. However, it uses a single network over the features and its extension to multiple taskdependent networks is not straightforward.
Note that “network feature selection” [12] can also
refer to the quite diﬀerent class of problems where theobjects, not their features, are connected over a graph.

202

Copyright © SIAM.
Unauthorized reproduction of this article is prohibited.

easily be shown that this gives exactly the same solution
We evaluate the proposed method, which we refer to as Grace. As the incidence matrix can be constructedas Multi-SConES, on both synthetic and real data. in linear time in the number of vertices and edges, thisThroughout all experiments, q(v) is set to the absolute is a much faster implementation.
Parameter selection. In feature selection, there
value of Pearson’s correlation coeﬃcient between aisgenerallyno ground truth to validate selected featuresfeature v and the response y. This setting is slightly
intraining.
Thus one must use a proxy to evaluate
diﬀerent from that in [2], where linear SKAT was used.
therelative
qualityof the solutions given by diﬀerent
Environment. We used Mac OS X version 10.7.4
parameter
values.
One
such proxy, used in [2], is thewith 2 × 3 GHz Quad-Core Intel Xeon CPU andstabilityoftheselection.
Another possibility, which we
16 GB of memory. SConES and Multi-SConES were
used
here,
is
to
considerpredictivity, which is still an
implemented in R, version 2.15.1. All experiments wereindicatorofthequalityoftheselected features although
performed in the R environment.
we
do
not
wish
to
optimizeitdirectly.
Comparison partners. We evaluate our method
For
every
method,we
performed
10-fold crossin both single-task and multi-task contexts. We systemvalidationandselected
optimalparametersthat yieldatically use as a baseline the top-k features ranked by qthelowestmean
squarederror(MSE).
alone, where k is set to be the same number of features
Evaluation criteria. Our main goal is to recoveras selected by SConES or its multi-task version. This
truly
causal features, or, in other words, to accurately
allows us to evaluate the impact of the structured reguclassify
the features into causal and non-causal. As thislarizer. We refer to this method as Correlation Ranking.
binaryclassification
problem is imbalanced, we evaluIn single-task feature selection, we also compare
ate
performance
using
Matthews correlation coeﬃcient
SConES with two standard feature selection algorithms,(MCC[22]).
MCC
rangesfrom −1 to 1, 1 being best.
Lasso [28] and Elastic Net [34], as well as three state-ofIn
experimentsonsynthetic
data, we also evaluate
the-art structured regularizer methods, group Lasso [27]thepredictivityofthe
selectedfeaturesand report MSE(with groups formed by edges as suggested by [15]),(ranging
from
0
to
1,
0
being
best)
on
a test set using
Grace, and aGrace [18, 19]. Grace and aGrace, whichridge
regression
ontheselectedfeatures.
use a Laplacian graph regularizer, can be considered
4

Experiments

Lasso equivalents of SConES. In Lasso and Elastic Net,the network information was just ignored.
In multi-task settings, we compare Multi-SConES
to multi-task Lasso [23] and multi-task Grace. Since
Grace is for single-task feature selection, we construct
an artificial dataset including a given network using thereformulation in Lemma 1 of [18], followed by applying
multi-task Lasso to the dataset (Supplementary Note C).
As we cannot determine sign changes in the objectivefunction of aGrace [19, Section 2.2] over multiple tasks,
we do not use a multi-task version of aGrace. Since thesemethods cannot treat diﬀerent networks for diﬀerent
tasks, we are limited, for our experiments in this paper,to cases where all tasks share the same network.
In addition, as group Lasso clearly underperforms
Grace in our single-task experiments, we keep Graceas the sole structured regularized Lasso comparison
partner in the multi-task setting.
We used the glmnet package in R for Lasso, Elastic
Net, and multi-task Lasso, and the SGL package in R
for group Lasso. We implemented Grace and aGracein R, based on the reformulation in Lemma 1 of [18].
As this method requires to compute the single value
decomposition (SVD) of the network’s Laplacian, it
does not scale to large networks. We therefore adopt
a new reformulation, replacing the matrix obtained by
SVD with the incidence matrix of the network. It can

4.1 Evaluation on Synthetic Data. We first evaluate (Multi-)SConES on synthetic data. We simulate
four types of gene regulatory networks, models 1 to 4,
which are exactly the same as described in [18] (Supplementary Note D) and use their combination in themulti-task setting. In each network there are 2, 200 features and the first 44 features are causal to the response.
Models 1 and 3 (resp. 2 and 4) are positively (resp. negatively) correlated networks. Model 3 (resp. model 4) is
identical to model 1 (resp. model 2), but the connection
in models 3 and 4 is weaker than in models 1 and 2.
For each model, we generate training and test
datasets of 100 samples each, and report the average
MCC as well as the test MSE over 50 repetitions.
Eﬃciency. We first analyze the runtime of MultiSConES with respect to the number of tasks. We create
multi-task problems with varying number of tasks from
1 to 100 by repeatedly combining model 1 with itself,
and report the runtimes of multi-task Lasso, MultiGrace and Multi-SConES in Figure 2a.
Empirically, the runtime of Multi-SConES increases
cubically with the number of tasks. While this is suboptimal, in particular compared with Multi-Grace, wemust remember that Multi-Grace cannot use diﬀerent
networks for diﬀerent tasks. Moreover, Multi-SConES
is still eﬃcient enough to make it possible to analyze

203

Copyright © SIAM.
Unauthorized reproduction of this article is prohibited.

2

5
10
20
Number of tasks

50

100

MCC

MCC

0.9
0.8
0.7
0.6
0.5
0.4
0.3

400

400

400

300

300

300

300

200

100
0.15 0.20 0.25 0.30

η

0
0.01

MSE

MSE

200

200

100
1.00

λ

100

sharing all features

200

100

0
0.01 0.10 1.00 10.00

μ

MSE

400

300

MSE

400

0

1

0.9
0.8
0.7
0.6
0.5
0.4
0.3

300

100

0.01

0.9
0.8
0.7
0.6
0.5
0.4
0.3

400

200

0.10

0.9
0.8
0.7
0.6
0.5
0.4
0.3

MSE

1.00

0.9
0.8
0.7
0.6
0.5
0.4
0.3

MCC

MCC

10.00

MSE

Running time (s)100.00

0.8
0.7
0.6
0.5
0.4
0.3

MCC

b 0.9
Lasso
Grace
SConES

MCC

a

200

100

100

0
0.01 0.10 1.00 10.00

00.01 0.10 1.00 10.00

00.01 0.10 1.00 10.00

sharing 3/4 features

sharing 1/2 features

sharing no feature

μ

μ

μ

Figure 2: (a) Running time with respect to changes in number of tasks under fixed regularization parameters. Data
are mean ± SEM. (b) Feature selection performance with respect to changes in regularization parameters η, λ, and µ.
Note that x-axes for λ and µ have logarithmic scales. The eﬀects of changes in µ are reported for various feature-sharing
scenarios. Two parameters η and λ behave identically independently of the amount of true causal features shared by thetasks and corresponding plots are therefore not reported. Data are means ± SEM.
hundreds of thousands of features for the order of a these diﬀerent scenarios and is therefore not reported.
In contrast, MCC (resp. MSE) shows a concave
dozen of tasks, which matches the statistical genetics(resp. convex) response to η, which fits to our theoretisetting that motivates our study.
In the particular case where we want to select the cal analysis (Section 2.1, Regularization path), althoughsame features for all tasks and have a single network here we evaluate solutions in terms of generalization erstructure over the features, Multi-SConES reduces to a ror on independently generated test datasets.
In practice, this means that in cases well-behaved
single-task problem over a network of same size: q(v) is∑enough,we do not need to carefully tune the regularizasimply replaced with K
q
(v)
and
µ
set
to
zero.
In
i
i=1
that case Multi-SConES is as eﬃcient as the single-task tion parameters λ and µ. As the entire regularization
SConES, which is much more eﬃcient than Grace over path with respect to η can be obtained without increasing the time complexity, finding optimal parameters forthe same network (as corroborated by [2]).
Parameter sensitivity. Next, we analyze the Multi-SConES becomes attractively inexpensive.
Feature selection performance. We then evalubehavior of Multi-SConES with respect to changes in
the regularization parameters λ, η, and µ. For that ate the feature selection performance of Multi-SConES
purpose, we fix two of those parameters, and run Multi- in both single-task and multi-task settings (Figure 3).
In the single-task setting, SConES shows much betSConES for two-task feature selection over models 1
terperformance
than the baseline Correlation Rankand 2. To understand parameter sensitivity with respect
ing.
This
means
that the penalty function g in Equato the amount of causal features shared across tasks,
tion
(2.2)
works
well
to select connected features. Morewe perform experiments for four cases: models 1 and 2
over,
SConES
outperforms
all the other methods inshare all, 3/4, half, or none of their features. We use
terms
of
MCC,
showingthat
it is better at recovering
λ = 1, µ = 1, and η = 0.2 when they are fixed.
truecausalfeatures.
Only
Lasso
(and, in one case, ElasResults are shown in Figure 2b. Multi-SConES isticNet)
outperforms
SConESinterms
of predictivity ofsensitive to η, while more robust to µ and robust to λ iftheselectedfeatures.
Howeverthefeaturesit selects are
it is set large enough. The robustness of Multi-SConES
too
sparseanddisconnected,
resultinginnotably
worsewith respect to λ can be understood as follows: once λMCCscoresanddiﬃculties
ininterpretability.
Theseis large enough to cause the true causal features to be
selected, if they form a subnetwork disconnected from results are consistent with the behavior reported in [2].
To evaluate Multi-SConES, we create multi-task
the rest of the network, the corresponding penalty term
problems
by combining the models. More precisely,
becomes 0, and increasing λ will not aﬀect the objective.
formodels1 and 2, two-task problems are created by
Similarly, if the true causal features are identical
combiningbothmodels, and three-task problems byacross all tasks, the penalty term controlled by µ is also
combiningthemwithmodel 3; for models 3 and 4, two0, and varying µ will not aﬀect the objective. However,task
problems
arecreatingby combining both models,
if the causal features are not shared across all tasks,
andthree-task
problems
bycombiningthem with model
setting µ too large enforces the selection of too many
1.
The
four-task
problemcombinesallfour models.
identical features and leads to poor solutions. The
Multi-SConES
outperforms
SConESin all cases.
behaviors of λ and η, however, remain unchanged across204

Copyright © SIAM.
Unauthorized reproduction of this article is prohibited.

400

0.8

300

0.6

MSE

MCC

Model 1
1.0

0.4

100

0.2
0
Model 2
1.0

CR LA EN GL GR AG SC

CR LA GR SC

CR LA GR SC

CR LA GR SC

Single task

Two tasks

Three tasks

Four tasks

0

MSE

MCC

0.4

CR LA EN GL GR AG SC

CR LA GR SC

CR LA GR SC

CR LA GR SC

Single task

Two tasks

Three tasks

Four tasks

MSE

MCC

0.6
0.4
0.2
CR LA GR SC

CR LA GR SC

CR LA GR SC

Single task

Two tasks

Three tasks

Four tasks

Four tasks

CR LA EN GL GR AG SC

CR LA GR SC

CR LA GR SC

CR LA GR SC

Single task

Two tasks

Three tasks

Four tasks

CR LA EN GL GR AG SC

CR LA GR SC

CR LA GR SC

CR LA GR SC

Single task

Two tasks

Three tasks

Four tasks

CR LA EN GL GR AG SC

CR LA GR SC

CR LA GR SC

CR LA GR SC

Single task

Two tasks

Three tasks

Four tasks

100
80

0.6

60

MSE

MCC

CR LA EN GL GR AG SC

140
120
100
80
60
40
20
0

0.8

0.4

40
20

0.2
0

CR LA GR SC

Three tasks

100

0

0.8

Model 4
1.0

CR LA GR SC

Two tasks

50

0.2

0

CR LA GR SC

Single task
150

0.6

Model 3
1.0

CR LA EN GL GR AG SC

200

0.8

0

200

CR LA EN GL GR AG SC

CR LA GR SC

CR LA GR SC

CR LA GR SC

Single task

Two tasks

Three tasks

Four tasks

0

Figure 3: Feature selection performance for synthetic data. MCC (left column) should be maximized and MSE(right column) should be minimized. CR: the ranking of correlations (baseline), LA: Lasso, EN: the elastic net,
GL: group Lasso, GR: Grace, AG: aGrace, and SC: SConES. Data are means ± SEM.
Moreover, performance (MCC and MSE) improves with
the number of tasks. This confirms that our multi-task
formulation on feature networks is eﬀective comparedto solving each task independently. Furthermore, MultiSConES achieves significantly better MCC than all of its
comparison partners and is also now superior in terms
of predictivity. Our method is therefore eﬀective formulti-task feature selection on networks.
We also examine Multi-SConES when causal features are not exactly shared. For each data model, we
perform two-task feature selection (models 1 + 2 and3 + 4) assuming that a fraction σ (σ = 3/4, 1/2 or 0) of
the causal features are shared between both tasks. Results are shown in Figure 4. Once again, Multi-SConES
clearly outperforms all other methods in terms of MCC.
The features it selects are also more predictive, except
for those selected by multi-task Lasso when half of thefeatures are shared between tasks. The more featuresare shared, the better Multi-SConES is at recovering
causal, explanatory features. This holds for all multitask methods and is typical in multi-task learning.

4.2 Multi-Locus Association Mapping. As a realworld application, we performed large-scale multi-locus
mapping of Arabidopsis thaliana flowering time phenotypes. Our goal here is to uncover which SNPs are associated with flowering time, using a network over SNPs
derived from biological properties.
Data preparation. We used the Arabidopsis
thaliana GWAS data1 collected by [1]. This dataset contains 216,130 SNPs (features) for 199 individuals, withmissing values, and there are 23 flowering phenotypes
in total. It is well known that A. thaliana is susceptible
to population structure confounding, that is, the existence of subpopulations (clusters) of individuals due to,for instance, diﬀerent demographic histories or diverse
environmental influences, induces fake correlations between genotypes and phenotypes. We corrected it using
principal component analysis as in [2, 24].
To evaluate the quality of the selected SNPs (features), we used the 282 candidate (causal) genes for flow1 https://cynin.gmi.oeaw.ac.at/home/resources/atpolydb/
genomic-polymorphism-data-in-arabidopsis-thaliana

205

Copyright © SIAM.
Unauthorized reproduction of this article is prohibited.

400

0.6

300

MSE

500

0.8
MCC

Model 1
1.0

0.4
0.2

100

0

CR LA GR SC

Model 2
1.0

CR LA GR SC

0CR LA GR SC

sharing 3/4 features sharing 1/2 features sharing no feature

MCC

MSE
CR LA GR SC

CR LA GR SC

0CR LA GR SC

Model 3 sharing 3/4 features sharing 1/2 features sharing no feature
1.0

MSE

0.6
0.4
0.2
CR LA GR SC

CR LA GR SC

CR LA GR SC

Model 4 sharing 3/4 features sharing 1/2 features sharing no feature
1.0

140
120
100
80
60
40
20
0

CR LA GR SC

CR LA GR SC

CR LA GR SC

CR LA GR SC

CR LA GR SC

sharing 3/4 features sharing 1/2 features sharing no feature
100
80

0.6

60

MSE

0.8

0.4

CR LA GR SC

sharing 3/4 features sharing 1/2 features sharing no feature

0.8
MCC

100
50

0.2

MCC

CR LA GR SC

150

0.4

0.2
0

CR LA GR SC

200

0.6

0

CR LA GR SC

sharing 3/4 features sharing 1/2 features sharing no feature

0.8

0

200

40
20

CR LA GR SC

CR LA GR SC

CR LA GR SC

sharing 3/4 features sharing 1/2 features sharing no feature0

CR LA GR SC

CR LA GR SC

CR LA GR SC

sharing 3/4 features sharing 1/2 features sharing no feature

Figure 4: Feature selection performance in two tasks forsynthetic data. A fraction σ (= 3/4, 1/2 or 0) of thecausal features are shared. Data are means ± SEM.

Once again, Multi-SConES shows much better performance in terms of MCC than its competitors. In
addition, the proportion of SNPs near candidate genes
among the selected SNPs is higher for Multi-SConES
than those for the other methods. Finally, combining several phenotypes helps recovering more candidate
causal genes. Altogether, this means that our multitask feature selection strategy can also be eﬀectively
employed for the important real-world problem of multilocus association mapping in A. thaliana.
5 Conclusion
In this paper we have proposed Multi-SConES, a new
formulation of multi-task feature selection with multiple network regularizers. We directly optimize feature
relevance scores and exactly solve the formulation by
maximum flow algorithms. Compared to the typicalstructured Lasso approaches, Multi-SConES shows improved ability to discover causal features in simulated
and real-world experiments.
Unlike existing structured sparsity multi-task feature selection methods, Multi-SConES can use diﬀerent
networks for diﬀerent tasks, and yields a clear, binaryclassification of features. Another attractive property of
our approach is the possibility to incorporate cardinality
constraints on the size of the solution set.
Currently, we model the relationship between tasks
with a single parameter µ, which controls how coupled
the solutions are. Some multi-task Lasso models can
include more detailed structures of correlation between
tasks (e.g. [16]). Others also consider relationships
between tasks using a task covariance matrix [7]. In
future work, we will study how to integrate these typesof more complex task relationships into Multi-SConES.
Another interesting direction is to incorporate proximal methods with structured norms, which is also shown
to be related to the maximum flow problem [20].

ering time listed in [5, (Table 12 in Dataset S1)], originally provided in [1], as an approximation of the gold
standard. For each selected SNP, we checked whether
or not it is located within 20 kb of one of the 282 candidate genes as in [2, 5]. If a SNP belongs to more than
two genes, we assigned it to the closest gene.
We derived a network over SNPs from theA. thaliana protein-protein interaction network between
genes from TAIR2 (The Arabidopsis Information Resource). SNPs were connected with a weight of 1 ifthey belong to the same gene or connected genes. In
addition, we connected each pair of SNPs adjacent on
the genomic sequence with a small weight of 0.01.
Results. For each of two phenotypes (2W and References
LDV), which have high correlations to other phenotypes[1] Atwell, S., et al.: Genome-wide association study ofin average, we picked two additional phenotypes they
107 phenotypes in Arabidopsis thaliana inbred lines.
are highly correlated with (4W and FTGH for 2W, 0W
Nature 465(7298), 627–631 (2010)
and FT10 for LDV) and checked whether or not these[2] Azencott, C.A., Grimm, D., Sugiyama, M., Kawahara,
additional phenotypes improve performance of MultiY., Borgwardt, K.M.: Eﬃcient network-guided multiSConES and its competitors.
locus association mapping with graph cuts. Bioinformatics 29(13), i171–i179 (2013)For each method, we determined the optimal pa[3] Bach, F.: Structured sparsity-inducing norms throughrameters by 10-fold cross-validation and run it on thesubmodular functions. In: NIPS. pp. 118–126 (2010)
full data to get a final set of selected SNPs. We report MCC as well as the ratios of candidate SNPs (resp. [4] Boykov, Y., Kolmogorov, V.: An experimental comparison of min-cut/max-flow algorithms for energy minigenes) retrieved with respect to the number of selectedmization in vision. TPAMI 26(9), 1124–1137 (2004)
SNPs (resp. genes) in Table 1.
[5] Brachi, B., et al.: Linkage and association mappingof Arabidopsis thaliana flowering time in nature. PLoS
Genetics 6(5), e1000940 (2010)

2 http://www.arabidopsis.org/

206

Copyright © SIAM.
Unauthorized reproduction of this article is prohibited.

Table 1: Results of multi-locus mapping for Arabidopsis thaliana genotypes and flowering time phenotypes.
Phenotype
Lasso

MCC
Grace

SConES

Hit ratio of SNPs
Lasso
Grace
SConES

Hit ratio of genes
Lasso Grace SConES

2W
2W + 4W
2W + FT GH
2W + 4W + FT GH

0.001−0.001
0.001
0.005−0.001−0.003
0.000
0.002

0.014
0.016
0.024
0.027

7/126
7/175
9/173
15/183

4/98
6/198
7/146
16/26542/33881/802
106/818
101/679

2/112
2/163
9/162
6/174

1/91
2/191
7/135
3/256

7/124
11/24013/25013/208LDV
LDV + 0W
LDV + FT10
LDV + 0W + FT10

0.001
0.005
0.001
0.003

0.000
0.007
0.001
0.002

0.016
0.020
0.021
0.023

6/116
16/196
12/214
18/283

7/144
19/206
10/191
19/32373/66786/702
92/762
81/482

2/107
2/183
1/199
2/265

2/131
2/187
1/181
1/307

9/202
10/20910/221
10/153[6] Chen, X., et al.: A two-graph guided multi-task lasso
approach for eQTL mapping. In: AISTATS. pp. 208–
217 (2012)[7] Fei, H., Huan, J.: Structured feature selection and task
relationship inference for multi-task learning. KAIS
35(2), 345–364 (2013)[8] Gallo, G., Grigoriadis, M.D., Tarjan, R.E.: A fastparametric maximum flow algorithm and applications.
SICOMP 18(1), 30–55 (1989)[9] Gao, S., Tsang, I., Chia, L.: Laplacian sparse coding,
hypergraph Laplacian sparse coding, and applications.
TPAMI 35(1), 92–104 (2013)[10] Goldberg, A.V., Tarjan, R.E.: A new approach to the
maximum-flow problem. JACM 35(4), 921–940 (1988)[11] Gretton, A., Bousquet, O., Smola, A., Schölkopf, B.:
Measuring statistical dependence with Hilbert-Schmidt
norms. In: ALT. pp. 63–77. Springer (2005)[12] Gu, Q., Han, J.: Towards feature selection in network.
In: CIKM. pp. 1175–1184 (2011)[13] Hastie, T., Rosset, S., Tibshirani, R., Zhu, J.: The entire regularization path for the support vector machine.
JMLR 5, 1391–1415 (2005)[14] Huang, J., Zhang, T., Metaxas, D.: Learning withstructured sparsity. JMLR 12, 3371–3412 (2011)[15] Jacob, L., Obozinski, G., Vert, J.P.: Group lasso withoverlap and graph lasso. In: ICML. pp. 433–440 (2009)[16] Kim, S., Sohn, K.A., Xing, E.: A multivariate regression approach to association analysis of a quantitative
trait network. Bioinformatics 25(12), i204–i212 (2009)[17] Lee, S., Xing, E.P.: Leveraging input and output
structures for joint mapping of epistatic and marginal
eQTLs. Bioinformatics 28(12), i137–i146 (2012)[18] Li, C., Li, H.: Network-constrained regularizationand variable selection for analysis of genomic data.
Bioinformatics 24(9), 1175–1182 (2008)[19] Li, C., Li, H.: Variable selection and regression analysis
for graph-structured covariates with an application to
genomics. Ann Appl Stat 4(3), 1498–1516 (2010)[20] Mairal, J., Jenatton, R., Obozinski, G., Bach, F.:
Convex and network flow optimization for structured
sparsity. JMLR 12, 2681–2720 (2011)[21] Marchini, J., Donnelly, P., Cardon, L.R.: Genome-

wide strategies for detecting multiple loci that influence
complex diseases. Nat Genet 37(4), 413–417 (2005)[22] Matthews, B.W.: Comparison of the predicted andobserved secondary structure of T4 phage lysozyme.
Biochimica et Biophysica Acta 405(2), 442–451 (1975)[23] Obozinski, G., Wainwright, M.J., Jordan, M.I.: Highdimensional union support recovery in multivariate
regression. In: NIPS. pp. 1217–1224 (2008)[24] Price, A.L., Patterson, N.J., Plenge, R.M., Weinblatt,M.E., Shadick, N.A., Reich, D.: Principal components
analysis corrects for stratification in genome-wide association studies. Nat Genet 38(8), 904–909 (2006)[25] Reshef, D.N., et al.: Detecting novel associations inlarge data sets. Science 334(6062), 1518–1524 (2011)[26] Saha, B., Pham, D.S., Phung, D., Venkatesh, S.:
Sparse subspace clustering via group sparse coding. In:
SDM. pp. 130–138 (2013)[27] Simon, N., Friedman, J., Hastie, T., Tibshirani, R.: Asparse-group lasso. J Comput Graph Stat 22(2), 231–
245 (2013)[28] Tibshirani, R.: Regression shrinkage and selection via
the lasso. J Roy Stat Soc B 58(1), 267–288 (1996)[29] Wu, M.C., Lee, S., Cai, T., Li, Y., Boehnke, M.,
Lin, X.: Rare-variant association testing for sequencing
data with the sequence kernel association test. Am J
Hum Genet 89(1), 82–93 (2011)[30] Yang, S., Yuan, L., Lai, Y.C., Shen, X., Wonka, P., Ye,J.: Feature grouping and selection over an undirected
graph. In: KDD. pp. 922–930 (2012)[31] Yuan, M., Lin, Y.: Model selection and estimation inregression with grouped variables. J Roy Stat Soc B 68,
49–67 (2006)[32] Zhang, Y., Yeung, D.Y., Xu, Q.: Probabilistic multitask feature selection. In: NIPS. pp. 2559–2567 (2010)[33] Zhou, Y., Jin, R., Hoi, S.C.H.: Exclusive lasso for
multi-task feature selection. In: AISTATS. vol. 9, pp.
988–995 (2010)[34] Zou, H., Hastie, T.: Regularization and variable selection via the elastic net. J Roy Stat Soc B 67(2), 301–320(2005)

207

Copyright © SIAM.
Unauthorized reproduction of this article is prohibited.