Efficient Multitask Feature and Relationship Learning

Han Zhao and Otilia Stretcu
Machine Learning Department
Carnegie Mellon University{han.zhao, ostretcu}@cs.cmu.edu

Alexander J. Smola
Amazon Web Services
alex@smola.org

Abstract
We consider a multitask learning problem, inwhich several predictors are learned jointly.
Prior research has shown that learning the relations between tasks, and between the input
features, together with the predictor, can lead tobetter generalization and interpretability, which
proved to be useful for applications in many domains. In this paper, we consider a formulationof multitask learning that learns the relationships both between tasks and between features,
represented through a task covariance and afeature covariance matrix, respectively. First,
we demonstrate that existing methods proposed
for this problem present an issue that may leadto ill-posed optimization. We then propose analternative formulation, as well as an efficientalgorithm to optimize it. Using ideas from optimization and graph theory, we propose an efficient coordinate-wise minimization algorithm
that has a closed form solution for each block
subproblem. Our experiments show that the
proposed optimization method is orders of magnitude faster than its competitors. We also provide a nonlinear extension that is able to achieve
better generalization than existing methods.

1

INTRODUCTION

In machine learning the goal is often to train predictivemodels for one or more tasks of interest. Making accurate
predictions relies heavily on the existence of labeled data
for the desired tasks. However, in real-world problems
data is often hard to acquire (e.g., medical domains) or
expensive to label (e.g., image segmentation). For manytasks, this makes it impractical or impossible to collectlarge volumes of labeled data. Multitask learning is a subproblem of the general transfer learning paradigm that

Geoffrey J. Gordon
Microsoft Research Montreal
Carnegie Mellon University
geoff.gordon@microsoft.com

aims to improve generalization performance in a learning task, by learning models for multiple related tasks
simultaneously. It has received considerable interest in
the past decades (Adel et al., 2017, Argyriou et al., 2007,2008, Caruana, 1997, Chen et al., 2011, 2012, Evgeniou
and Pontil, 2004, Jacob et al., 2009, Jawanpuria et al.,2015, Kato et al., 2008, Li et al., 2015, Liu et al., 2009,
Zhang and Schneider, 2010, Zhang and Yeung, 2010a,
Zhao et al., 2017, 2018, 2019a,b). One of the underlying
assumptions behind many multitask learning algorithms
is that the tasks are related to each other. Hence, a key
question is how to define the notion of task relatedness,and how to capture it in the learning formulation. A common assumption is that tasks can be described by weightvectors, and that they are sampled from a shared priordistribution over their space (Liu et al., 2009, Zhang and
Yeung, 2010a,b). Another strand of work assumes common feature representations to be shared among multiple
tasks, and the goal is to learn the shared representation as
well as task-specific parameters simultaneously (Argyriou
et al., 2008, Caruana, 1997, Evgeniou and Pontil, 2007,
Thrun, 1996). Moreover, when structure about multiple
tasks is available, e.g., task-specific descriptors (Bonilla
et al., 2007) or a task similarity graph (Evgeniou and Pontil, 2004), regularizers can often be incorporated into thelearning formulation to penalize hypotheses that are notconsistent with the given structure. Very recently, Senerand Koltun (2018) tackle the problem of multitask learning from the perspective of multi-objective optimization.
Specifically, this work aims to find a Pareto-optimal solution for the multi-objective function defined by multiple
tasks, and proposes to use Frank-Wolfe algorithm to findgradient update for shared parameters.
There have been several attempts to improve predictions
by either learning the relationships between different
tasks (Zhang and Yeung, 2010a), or by exploiting therelationships between different features (Argyriou et al.,2008). In this paper we consider a multiconvex framework for multitask learning that improves predictions overtabula rasa learning by assuming that all the task vectors

are sampled from a common matrix-variate normal prior.
The framework, known as MTFRL (Zhang and Schneider, 2010), learns the relationships both between tasks
and between features simultaneously via two covariancematrices, i.e., the feature covariance matrix and the task
covariance matrix. In this context, learning multiple tasks
corresponds to estimating a matrix of model parameters,
and learning feature/task relationships corresponds to estimating the row/column covariance matrices of model
parameters, respectively. This property is favorable forapplications where we not only aim for better generalization, but also seek to have a clear understanding about therelationships among different tasks.
The goal of MTFRL is to optimize over both the task vectors, as well as the two covariance matrices in the prior.
When the loss function is convex, the regularized problem
of MTFRL is multiconvex. Previous approaches (Longet al., 2017, Zhang, 2011, Zhang and Schneider, 2010)
for solving this problem hinge on the classic flip-flop algorithm (Dutilleul, 1999) to estimate the two covariancematrices. However, as we point out in Section 3, the flipflop algorithm cannot be directly applied as the maximum
likelihood estimation (MLE) formulation of the multitask
learning problem under this setting is ill-posed. As a result, in practice, heuristics have to be invented and appliedin the algorithm to ensure the positive-definiteness of bothcovariance matrices. However, it is not clear whether sucha fixed algorithm still converges or not.
In this paper we propose a well-defined variant of the MTFRL framework, and design a block coordinate-wise minimization algorithm to solve this problem. We term our
new formulation FEaTure and Relation learning (FETR).
By design, FETR is free of the nonpositive-definite problem in MTFRL. To solve FETR, we propose efficient andanalytic solutions for each of the subproblems, which allows us to get rid of the expensive iterative procedure to
optimize the covariance matrices. Specifically, we achieve
this by reducing an underlying matrix optimization problem with positive definite constraints into a minimum
weight perfect matching problem on a complete bipartite graph, where we are able to solve analytically using
combinatorial techniques. To solve the weight learning
subproblem, we propose three different strategies, including a closed form solution, a gradient descent methodwith linear convergence guarantees when the instances
are not shared by multiple tasks, and a numerical solution
based on Sylvester equation when instances are shared.
We demonstrate the efficiency of the proposed optimization algorithm by comparing it with an off-the-shelf projected gradient descent algorithm and the classic flip-flop
algorithm, on both synthetic and real-world data. Experiments show that the proposed optimization method isorders of magnitude faster than its competitors, and it

often converges to better solutions. Lastly, we extend
FETR to nonlinear setting by combining its regularization
scheme with rich nonlinear transformations using neuralnetworks. This combined approach is able to achievesignificantly better generalizations than existing methods
on real-world datasets.
To summarize, our contributions are three-fold:• We point out an ill-posed MLE problem of the existing multitask learning formulations and propose awell-defined variant, termed as FETR.
• To optimize FETR, we design an efficient block
coordinate-wise minimization algorithm and derive
analytic solutions for each of the subproblems.
• We extend our FETR formulation to nonlinear settings and empirically demonstrate its better generalizations on real-world datasets.

2

PRELIMINARY

We start by introducing notations used throughout the paper and briefly discussing the MTFRL framework (Zhang
and Schneider, 2010).
2.1

NOTATION AND SETUP

We use lowercase letters, such as y, to represent scalars,
and lowercase bold letters, such as x, to denote vectors.
Capital letters are reserved for matrices. We use Sm+ andSm++ to denote the m-dimensional symmetric positive
semidefinite cone and the m-dimensional symmetric positive definite cone, respectively. We write tr(A) for the
trace of a matrix A, and N (m, Σ) for the multivariate
normal distribution with mean m and covariance matrix
Σ. Finally, G = (A, B, E; w) is a weighted bipartite
graph with vertex sets A, B, edge set E and weight function w : E → R+ . For a matrix W ∈ Rd×m , we use
vec(W ) ∈ Rdm to denote its vectorization. We consider the following setup. Suppose we are given m learning tasks {Ti }m
i=1 , where for each learning task Ti wehave access to a training set Di with ni data instances(xji , yij ), j ∈ [ni ]. For the simplicity of discussion, here
we focus on the regression setting where xji ∈ Xi ⊆ Rd
and yij ∈ R. Extension to classification setting is straightforward. Let fi (wi , ·) : Xi → R be our model withparameter wi . In what follows, we will assume ourmodel for each task Ti to be a linear regression, i.e.,fi (wi , x) = wiT x.
2.2

MATRIX-VARIATE NORMAL
DISTRIBUTION

A matrix-variate normal distribution (Gupta and Nagar,
1999) W ∼ MN d×m (M, A, B) with mean M ∈ Rd×m ,row covariance matrix A ∈ Sd++ and column covariancematrix B ∈ Sm++ can be understood as a multivariate normal distribution with vec(W ) ∼ N (vec(M ), A ⊗ B).1
One advantage of the matrix-variate normal distribution over its equivalent multivariate counterpart is thatby imposing structure on the row and column covariance matrices, the former admits a much more compact representation than the latter (O(m2 + d2 ) versusO(m2 d2 )). The MLE of the matrix-variate normal distribution has been well studied in the literature (Dutilleul, 1999). Specifically, given an i.i.d. sample{Wi }ni=1
Pnfrom MN d×m (M, A, B), the MLE of M is
W = i=1 Wi /n. The MLE of A and B are solutions
to the following system:
(
Pn
1
−1
A = nm(Wi − W )T
i=1 (Wi − W )B
P
(1)
n
1
T −1
B = nd i=1 (Wi − W ) A (Wi − W )The above system of equations does not have a closed
form solution as the two covariance estimates depend oneach other. Hence, their estimates must be computed in an
iterative fashion until convergence, which is known as the“flip-flop” algorithm (Dutilleul, 1999, Glanz and Carvalho,2013). Furthermore, Dutilleul (1999) showed that theflip-flop algorithm is guaranteed to converge to positive
definite covariance matrices iff n ≥ max(d/m, m/d)+1.
More properties of the MLE of the matrix-variate normal
distribution can be found in (Roś et al., 2016).
2.3

MULTITASK FEATURE AND
RELATIONSHIP LEARNING

In linear regression, the likelihood function for task
i is given by: yij | xji , wi , i ∼ N (wiT x, 2i ). Let
W = (w1 , . . . , wm ) ∈ Rd×m be the model parameterfor m different tasks drawn from the matrix-variate nor−1
mal distribution MN d×m (W | 0d×m , Σ−1
1 , Σ2 ). By
maximizing the joint distribution and optimize over boththe model parameters, as well as the two covariance matrices in the prior, we reach the following optimization
problem:minimize
W,Σ1 ,Σ2nim X
X(yij − wiT xji )2 + η tr(Σ1 W Σ2 W T )i=1 j=1− η (m log |Σ1 | + d log |Σ2 |)subject to Σ1  0, Σ2  0(2)where Σ1 ∈ Sd++ , Σ2 ∈ Sm++ are the row and column
precision matrices of the matrix normal prior distribution,
respectively, and η is a constant that does not depend on
the optimization variables. It is not hard to see that the optimization problem in (2) is not convex due to the coupling
Probability density: p(X) = exp(− 12 tr(A−1 (X −M )B −1 (X − M )T ))/(2π)md/2 |A|m/2 |B|d/2 .
1

between W, Σ1 and Σ2 in the trace term. On the other
hand, since the log | · | function is concave in the positive
definite cone (Boyd and Vandenberghe, 2004), and thetrace is linear in terms of its components, it follows that(2) is multiconvex. Zhang and Schneider (2010) propose
to use the flip-flop algorithm to solve the matrix subproblem in (2), and this approach has also been widely appliedin following publications on multitask learning (Li et al.,2014, Long et al., 2017, Zhang, 2011).

3

ILL-POSED OPTIMIZATIONIn this section we first point out an important issue in theliterature on the application of the flip-flop algorithm tosolve the matrix subproblem in (2). We then proceed topropose a well-defined variant of (2) to fix the problem.
Interestingly, the variant we propose admits a closed form
solution for each block variable that can be computed
efficiently without any iterative procedure, which we will
describe and derive in more detail in Section 4.
As proved by Dutilleul (1999), one sufficient and necessary condition for the flip-flop algorithm to converge
to positive definite matrices is that the number of samples from the matrix-variate normal distribution should
satisfy n > max(d/m, m/d). However, in the context
of multitask learning, we are essentially dealing withan inference problem, where the goal is to estimate thevalue of W , which is assumed to be an unknown but
unique model parameter from the prior. This means thatin this case we have n = 1, hence the condition for theconvergence of the algorithm is violated. Technically,
for any W ∈ Rd×m where d 6= m, following the iterative update formula of the flip-flop algorithm in (1), for(0)(0)
any feasible initialization of Σ1 ∈ Sd++ and Σ2 ∈
(1)
(0)
(1)
Sm
= W (Σ2 )−1 W T /m, Σ2 =++ , we will have Σ1(0)
W T (Σ1 )−1 W/d. Now since d 6= m, we know thatrank(W ) ≤ min{d, m} < max{d, m}. As a result, after
one iteration, we will have(1)rank(Σ1 ) ≤ rank(W ) < max{d, m}(1)

rank(Σ2 ) ≤ rank(W T ) < max{d, m}
(1)

(1)

i.e., at least one of Σ1 and Σ2 is going to be rank
deficient, and in the next iteration the inverse operation is
not well-defined on at least one of them. As a fix, Zhang
and Schneider (2010) proposed to use an artificial fudge
factor to ensure that both covariance matrices stay positive
definite after each update:(t+1)= W (Σ2 )−1 W T /m + Id(t+1)= W T (Σ1 )−1 W/d + Im

Σ1

Σ2

(t)

(t)

where  > 0 is a fixed, small constant. However, since
the fudge factor  is a fixed constant which does not de-crease to 0 in the limit, it introduces extra biases into theestimation, and thus it is not clear whether or not the fixedalgorithm converges.
Perhaps what is more surprising is that (2) is not even welldefined as an optimization problem. As a counterexample,
we can fix W = 0d×m and let Σ1,σ = σId , Σ2,σ = σImwith σ > 0 so that both Σ1,σ and Σ2,σ are feasible. Now
let σ → ∞, and it is easy to verify that in this casethe objective function goes to −∞. Although we only
provide one counterexample, there is no reason to believethat the one we find is the only case where (2) fails. In
fact, Roś et al. (2016) have recently shown that the MLEof Σ1 ⊗ Σ2 does not exist if n ≤ max{d/m, m/d}, andthe only nontrivial sufficient condition known so far toguarantee the existence of the MLE is n > md. However,in the context of multitask learning, the unknown model
parameter W is unique and hence we have n = 1  md.
Given the wide applications of the above multitask learning framework in the literature, as well as the flip-flop
algorithm in this setting, we feel it important and urgent
to solve the above ill-posed and nonpositive definite problem. To this end, for some positive constants 0 < l < u,we propose a variant of (2) as follows:
minimizeΣ1 ,Σ2 ,Wnim X
X(yij − wiT xji )2 + η tr(Σ1 W Σ2 W T )i=1 j=1− η (m log |Σ1 | + d log |Σ2 |)subject to lId  Σ1  uId , lIm  Σ2  uIm(3)

The bounded constraints make the feasible set compact.
Since the objective function is continuous, by the extreme
value theorem, we know that the matrix subproblem of(3) becomes well-defined and can achieve finite lower
and upper bounds within the feasible set. Alternatively,
one can also understand this constraint as specifying a
truncated matrix normal prior over the compact set. As
we will see shortly, technically the bounded constraint
also allows us to develop an optimization procedure for
W with linear convergence rate, which is an exponential
acceleration over the unbounded case.

4

MULTICONVEX OPTIMIZATIONIn this section we propose a block coordinate-wise minimization algorithm to optimize the objective given in (3).
In each iteration, we alternatively minimize over W with
Σ1 and Σ2 fixed, then minimize over Σ1 with W and Σ2
fixed, and lastly minimize Σ2 with W and Σ1 fixed. The
whole procedure is repeated until a stationary point isfound. Due to space limit, we defer all the proofs andderivations to appendix. To simplify the notation, we assume n = ni , ∀i ∈ [m]. Let Y = (y1 , . . . , ym ) ∈ Rn×m
be the target matrix and X ∈ Rn×d be the feature matrix

shared by all the tasks. Using this notation, the objective
can be equivalently expressed in matrix form as:
1/2

1/2

minimize||Y − XW ||2F + η ||Σ1 W Σ2 ||2F

subject to− η (m log |Σ1 | + d log |Σ2 |)lId  Σ1  uId , lIm  Σ2  uIm

Σ1 ,Σ2 ,W(4)

OPTIMIZATION OF W

4.1

In order to minimize over W when both Σ1 and Σ2 are
fixed, we solve the following subproblem:
1/2

1/2

minimize h(W ) , ||Y − XW ||2F + η ||Σ1 W Σ2 ||2FW

(5)
As shown in the last section, this is an unconstrained
convex optimization problem. We present three different
algorithms to find the optimal solution of this subproblem.
The first one guarantees to find an exact solution in closed
form in O(m3 d3 ) time. The second one does gradientdescent with fixed step size to iteratively refine the solution, and we show that in our case a linear convergence
rate can be guaranteed. The third one finds the optimal
solution by solving the Sylvester equation (Bartels andStewart, 1972) characterized by the first-order optimality
condition, after a proper transformation.
A closed form solution. It is worth noting that it is
not obvious how to obtain a closed form solution directly from the formulation in (5). An application of
the first order optimality condition to (5) will lead to:
X T XW + η Σ1 W Σ2 = X T Y . Hence except for the
special case where Σ2 = cIm with c > 0 a constant,the above equation does not admit an easy closed form
solution in its matrix representation. The workaround isbased on the fact that the d × m dimensional matrix spaceis isomorphic to the dm dimensional vector space, withthe vec(·) operator implementing the isomorphism from
Rd×m to Rdm . Using this property, we have:
Proposition
4.1. (5)canbesolvedinclosed form in O(m3 d3 + mnd2 ) time;
the optimal solution W ∗ is:vec(W ∗ )=−1TT
Im ⊗ (X X) + η Σ2 ⊗ Σ1
vec(X Y ).
The computational bottleneck in the above procedure isin solving an md × md system of equations, which scalesas O(m3 d3 ) if no further sparsity structure is available.
Gradient descent. The closed form solution shown
above scales cubically in both m and d, and requires
us to explicitly form a matrix of size md × md. Thiscan be intractable even for moderate m and d. In suchcases, instead of computing an exact solution to (5), we
can use gradient descent with fixed step size to obtain an
approximate solution. The objective function h(W ) in

(5) is differentiable and its gradient can be obtained inO(m2 d + md2 ) time as ∇W h(W ) = X T (Y − XW ) +η Σ1 W Σ2 . Note that we can compute in advance both
X T Y and X T X in O(nd2 ) time, and cache them so that
we do not need to recompute them in each gradient update step. Let λi (A) be the ith largest eigenvalue of a real
symmetric matrix A. Adapted from Nesterov (2013), we
provide a linear convergence guarantee for the gradientmethod in the following proposition:
Proposition 4.2. Let λl = λd (X T X) + ηl2 and λu =2
λ1 (X T X) + ηu2 . Choose 0 < t ≤ λu +λ. For all ε > 0,lgradient descent with step size t converges to the optima
within O(log(1/ε)) steps.
The computational complexity to achieve an ε approximate solution using gradient descent is O(nd2 +log(1/ε)(m2 d + md2 )). Compared with the O(m3 d3 +mnd2 ) complexity for the exact solution, the gradient descent algorithm scales much better provided the condition
number κ , λu /λl is not too large. As a side note, when
the condition
number is large, we can effectively reduce√
it to κ by using conjugate gradient method (Shewchuk
et al., 1994).
Sylvester equation. In the field of control theory, aSylvester equation (Bhatia and Rosenthal, 1997) is a matrix equation of the form AX + XB = C, where the goal
is to find a solution matrix X given A, B and C. For thisproblem, there are efficient numerical algorithms withhighly optimized implementations that can obtain a solution within cubic time. For example, the Bartels-Stewart
algorithm (Bartels and Stewart, 1972) solves the Sylvester
equation by first transforming A and B into Schur forms
by QR factorization, and then solves the resulting triangular system via back-substitution. Our third approach isbased on the observation that we can equivalently transform the first-order optimality equation into a Sylvester
equation by multiplying both sides of the equation by
−1 T
−1 T
Σ−1
1 : Σ1 X XW + η W Σ2 = Σ1 X Y . As a result,finding the optimal solution of the subproblem amounts
to solving the above Sylvester equation. Specifically, thesolution to the above equation can be obtained using theBartels-Stewart algorithm in O(m3 + d3 + nd2 ).
Remark. Both the gradient descent and the BartelsStewart algorithm find the optimal solution in cubic time.
However, gradient descent is more widely applicable thanthe Bartels-Stewart algorithm: the Bartels-Stewart algorithm only applies to the case where all the tasks sharethe same instances, so that we can write down the matrix
equation explicitly, while gradient descent can be appliedin the case where each task has different number of inputs
and those inputs are not shared among tasks. On the other
hand, as we will see in the experiments, in practice theBartels-Stewart algorithm is faster than gradient descent,and provides a more numerically stable solution.

4.2

OPTIMIZATION OF Σ1 AND Σ2

Algorithm 1 Minimize Σ1
Input: W , Σ2 and l, u.
1: [V, ν] ← SVD(W Σ2 W T ).
2: λ ← T[l,u] (m/ν).
3: Σ1 ← V diag(λ)V T .
Algorithm 2 Minimize Σ2
Input: W , Σ1 and l, u.
1: [V, ν] ← SVD(W T Σ1 W ).
2: λ ← T[l,u] (d/ν).
3: Σ2 ← V diag(λ)V T .
Before we delve into the detailed analysis below, we firstlist the final algorithms used to optimize Σ1 and Σ2 in
Algorithm 1 and Algorithm 2, respectively. The hardthresholding function used in Line 2 of Algorithm 1 andAlgorithm 2 is defined as follows:
T[l,u] (x) = max{l, min{u, x}}(6)

The hard-thresholding function essentially keeps the valueof its argument x if l ≤ x ≤ u, otherwise it truncatesthe value of x to l(u) if x < l(x > u) respectively. Bothalgorithms are remarkably simple: each algorithm only
involves one SVD, one truncation and two matrix multiplications. The computational complexities of Algorithm 1
and Algorithm 2 are bounded by O(m2 d + md2 + d3 )and O(m2 d + md2 + m3 ), respectively.
In what follows we focus on analyzing the optimization
w.r.t. Σ1 . A symmetric analysis can be applied to solve
Σ2 as well. In order to minimize over Σ1 when W andΣ2 are fixed, we solve the following subproblem:
minimizelId Σ1 uIdtr(Σ1 W Σ2 W T ) − m log |Σ1 |(7)

Although (7) is a convex optimization problem, it is computationally expensive to solve using off-the-shelf algorithms, e.g., the interior point method, because of theconstraints, as well as the non-linearity of the objectivefunction. However, as we will show shortly, we can find aclosed form optimal solution to this problem, using tools
from the theory of doubly stochastic matrices (Dufossé
and Uçar, 2016) and perfect bipartite graph matching.
Due to space limit, we defer the detailed derivation andproof to appendix, and only show a sketch below.
Without loss of generality, for any feasible Σ1 , using
spectral decomposition, we can reparametrize Σ1 as
Σ1 = U ΛU T ,Λ = diag(λ1 , . . . , λd )(8)where u ≥ λ1 ≥ λ2 · · · ≥ λd ≥ l. Similarly, we canrepresent
W Σ2 W T = V N V T ,N = diag(ν1 , . . . , νd )(9)where 0 ≤ ν1 ≤ · · · ≤ νd . Let λ = (λ1 , · · · , λd )T andν = (ν1 , · · · , νd )T . Set K = U T V and define P to be
the Hadamard product of K, i.e., P = K ◦ K. Sinceboth U and V T are orthonormal matrices, it immediately
follows that K is also an orthonormal matrix. As a result,
we have the following two equations hold:
d
X

Pij =

d
X

j=1

j=1

d
X

d
X

Pij =

i=1

2
Kij
= 1,

∀i ∈ [d]

2
Kij
= 1,

∀j ∈ [d]i=1which implies that P is a doubly stochastic matrix. Given
U being an orthonormal matrix, we have log |Σ1 | =log |U ΛU T | = log |Λ|. On the other hand, it can bereadily verified that the following equality holds:
tr(ΛKN K T ) =

d X
d
X

2
λi Kijνj = λT P ν(10)• For each λi , construct a vertex vλi ∈ Vλ , ∀i.
• For each νj , construct a vertex vνj ∈ Vν , ∀j.
• For each pair (vλi , vνj ), construct an edgee(vλi , vνj ) with weight w(e(vλi , vνj )) = λi νj .
The following theorem relates the solution of the minimum weight matching to the partial solution of (11) w.r.t.
P:
Theorem 4.1. Let λ = (λ1 , . . . , λd ) and ν =(ν1 , . . . , νd ) with λ1 ≥ · · · ≥ λd and ν1 ≤ · · · ≤ νd .
The minimum-weight perfect matching on G is the setof edges π ∗ = {(vλi , vνi ) : 1 ≤ i ≤ d} with the miniPd
mum weight w(π ∗ ) = i=1 λi νi . Furthermore, it equals
minP λT P ν.
Proof sketch. The full proof of Theorem 4.1 is deferredto the appendix, and here we only show a sketch of thehigh-level idea. Basically, given any matching in thegraph, if there is an inverse pair (a cross) in the matching,then we can improve the matching by re-matching theinverse pair, as shown in Figure 1. Now since there are

λi

λj

λi

λj

νk

νl

νk

νl

i=1 j=1

By combining all the transformations in (8), (9) and (10)and plug them in (7), we have the following equivalent
optimization problem:
minimizeP,λsubject toT

λ Pν − mdX

Figure 1: Re-matching an inverse pair (λi , λj , νk , νl ) ={(vλi , vνl ), (vλj , vνk )} on the left side to a match withsmaller weight {(vλi , vνk ), (vλj , vνl )}. Red color is used
to highlight edges in the perfect matching.

log λii=1

l1d ≤ λ ≤ u1d(11)

where 1d denotes a vector of all ones with dimension d.
To solve (11), we make the following key observations:
1. The minimization is decomposable in terms of P
and λ. Furthermore, the optimization over P is alinear program (LP).
2. For any bounded LP, there exists at least one extreme
point that achieves the optimal solution.
3. The set of d × d doubly stochastic matrices, denoted as Bd , forms a convex polytope, known as the
Birkhoff polytope.
4. By the Birkhoff-von Neumann theorem, Bd is theconvex hull of the set of permutation matrices, i.e.,every extreme point of Bd is a permutation matrix.
Combining all the analysis above, it is clear to see that the
optimal solution P must be a permutation matrix. This
motivates us to reduce (11) to a minimum-weight perfect matching problem on a weighted complete bipartite
graph as follows: for any λ, ν ∈ Rd+ , we can constructa weighted d × d bipartite graph G = (Vλ , Vν , E; w) as
follows:

only at most finitely many number of inverse pairs, an
inductive argument shows that the optimal matching isachieved when there is no inverse pair, i.e., vλi is matchedto vνi , ∀i ∈ [d] (Figure 2).
The optimal matching in Theorem 4.1 suggests that the
optimal doubly stochastic matrix is given by P ∗ = Id ,which also implies K ∗ = P ∗ = Id and U ∗ = V . Now
plug in the P ∗ = Id into (11). The optimization w.r.t.
λ decomposes into d independent scalar optimization
problems, which can be easily solved. Using the hardthresholding function defined in (6), we can express the
optimal solution λ∗i as λ∗i = T[l,u] (m/νi ). Combine all
the analysis given above, we get the algorithms listed
at the beginning of this section to optimize Σ1 and Σ2 .
Interestingly, they have close connection to the proximal
method proposed in the literature to solve matrix completion (Cai et al., 2010), or Euclidean projection under trace
norm constraint (Chen et al., 2011, 2012). To the best of
our knowledge, this is the first algorithm that solves linear function over matrices with negative log-determinant
regularization (e.g. (7)) efficiently.

λ1

λ2

···

λd

λ1

λ2

···

λd

λ1

λ2

···

λd

ν1

ν2

···

νd

ν1

ν2

···

νd

ν1

ν2

···

νd

Figure 2: The inductive proof works by recursively removing inverse pairs from (λd , νd ) to (λ1 , ν1 ). The process stops
until there is no inverse pair in the matching. Red color is used to highlight edges in the perfect matching.

5

NONLINEAR EXTENSION

So far we discuss our FETR framework under the linear
regression model, but it can be readily extended to any
nonlinear regression/classification settings. One straightforward way to do so is to apply the (orthogonal) random
Fourier transformation (Felix et al., 2016, Rahimi and
Recht, 2008) to generate high-dimensional random features so that linear FETR in the transformed space corresponds to nonlinear models in the original feature space.
However, depending on the dimension of the random
features, this approach might lead to a huge covariancematrix Σ1 that is expensive to optimize.
Another more natural and expressive approach is to combine our regularization scheme and optimization methodwith parametrized nonlinear feature transformations, suchas neural networks. More specifically, let g(x; θ) : Rd →Rp be a neural network with learnable parameter θ thatdefines a nonlinear transformation of the input featuresfrom Rd to Rp . Essentially we can replace the feature matrix X in (4) with g(x; θ) to create a regularized multitask
neural network (Caruana, 1997) where we add one morelayer defined by the matrix W on top of the nonlinear
mapping given by g(x; θ). To train the model, we canuse backpropagation to optimize W , θ and our proposed
approach to optimize the two covariance matrices. We
will further explore this nonlinear extension in Section 6
to demonstrate its power in statistical modeling.

6

EXPERIMENTS

6.1

CONVERGENCE ANALYSIS ANDCOMPUTATIONAL EFFICIENCY

We first investigate the efficiency and scalability of thethree different algorithms for minimizing w.r.t. W on
synthetic data sets. For each experiment, we generate asynthetic data set which consists of n = 104 instancesthat are shared among all the tasks. All the instances are
randomly sampled uniformly from [0, 1]d . We gradually
increase the dimension of features, d, and the number oftasks, m to test scalability.
The first algorithm implements the closed form solution
by explicitly computing the md × md matrix product andthen solving the linear system. The second one is the proposed gradient descent, and the last one uses the BartelsStewart algorithm to solve the equivalent Sylvester equation to compute W . We use open source toolkit scipy
whose backend implementation uses highly optimized
Fortran code. For all the synthetic experiments we set
l = 0.01 and u = 100, which corresponds to a condition
number of 104 . We fix the coefficients η = 1.0. We repeat each experiment for 10 times to show both the mean
and the variance.
The experimental results are shown in Figure 3a. As
expected, the closed form solution does not scale to problems of even moderate size due to its large memory requirement. In practice the Bartels-Stewart algorithm isabout one order of magnitude faster than the gradient descent method when either m or d is large. It is also worth
pointing out here that the Bartels-Stewart algorithm is themost numerically stable algorithm among the three basedon our observations.
We compare our proposed coordinate minimization algorithm with an off-the-shelf projected gradient method andthe flip-flop algorithm to solve the optimization problem(4). Specifically, the projected gradient method updates
W, Σ1 and Σ2 in each iteration and then projects Σ1 andΣ2 onto the corresponding feasible regions. The flipflop algorithm is implemented as suggested in Zhang andSchneider (2010) and we use a fudge factor of 10−3 toavoid the nonpositive definite problem. In each iteration,
both covariance matrices are projected onto the feasible
region as well. In the SARCOS dataset all the instances
are shared among all the tasks, so that the Sylvester solver
is used to optimize W in coordinate minimization. We
repeat the experiments 10 times and report the mean andstandard deviation of the log function values versus thetime used by all three algorithms (Figure 3b). It is clear
from Figure 3b that our proposed algorithm not only converges much faster than the other two competitors, butalso achieves better results. In fact, as we observe in our
experiments, the proposed algorithm usually converges inless than 10 iterations.
6.2

REAL-WORLD DATASETSIn this section we apply FETR to two real-world datasets
to demonstrate its statistical efficiency.

(a) The mean run time (seconds) under each experimental configu- (b) The convergence speed of coordinate minimization versus
ration. The closed form solution does not scale when md ≥ 104 . projected gradient descent and the flip-flop algorithm on the
SARCOS dataset. All the experiments are repeated 10 times.

Figure 3: Experimental results of the convergence analysis on synthetic data.
Table 1: Mean squared error on the SARCOS data and the mean of normalized mean squared error (NMSE) on theschool dataset across 10-fold cross-validation.
DATASETSM ETHOD
1 ST

2 ND

S ARCOS
3 RD
4 TH

STL

31.40

22.90

9.13

10.30

0.14

0.84

0.46

0.9882 ± 0.0196

MTFL
MTRL
MTFRL

31.41
31.09
31.13

22.91
22.69
22.60

9.13
9.08
9.10

10.33
9.74
9.74

0.14
0.14
0.13

0.83
0.83
0.83

0.45
0.44
0.45

0.8891 ± 0.0380
0.9007 ± 0.0407
0.8451 ± 0.0197

FETR

31.08

22.68

9.08

9.73

0.13

0.83

0.43

0.8134 ± 0.0253

STL-NN
MT-NN
MTFRL-NN
FETR-NN

24.81
12.01
11.02
10.77

17.20
10.54
9.51
9.34

8.97
5.02
4.99
4.95

8.36
7.15
7.11
7.01

0.13
0.09
0.08
0.08

0.72
0.70
0.62
0.59

0.34
0.27
0.27
0.24

−
−
−
−

Robot Inverse Dynamics This data relates to an inverse
dynamics problem for a seven degree-of-freedom (DOF)
SARCOS anthropomorphic robot arm (Vijayakumar andSchaal, 2000). The goal is to map from a 21-dimensional
input space (7 joint positions, 7 joint velocities, 7 joint
accelerations) to the corresponding 7 joint torques. Hence
there are 7 tasks and the inputs are shared among all
the tasks. The training set and test set contain 44,484
and 4,449 examples, respectively. We further partition
the training set into a training set and a validation set,containing 31,138 and 13,346 instances, respectively.
School Data This dataset consists of the examination scores of 15,362 students from 139 secondary
schools (Goldstein, 1991). It has 27 input features, andcontains 139 tasks. Since the train/test splits are not provided, we use a 10-fold cross-validation procedure to
generate the training and test datasets.

5 TH

6 TH

7 TH

6.3

S CHOOL

SETUP AND RESULTS

We compare FETR with multitask feature learning (Evgeniou and Pontil, 2007) (MTFL), multitask relationshiplearning (Zhang and Yeung, 2010a) (MTRL), and the
MTFRL framework. We also use ridge regression as ourbaseline model, denoted as single task learning (STL).
The results reported for FETR on the SARCOS dataset
are obtained using the Sylvester equation solver, while
for the School dataset the inputs are not shared among different tasks and hence we use our gradient descent solver
for W instead. To evaluate, we compute the mean ofnormalized mean squared error (NMSE) over the output
tasks (e.g., 139 tasks for the School data). The NMSE isdefined as the ratio of the MSE and the variance on a task.
For the School dataset, we show the mean NMSE and its
standard deviation across 10 cross-validation folds, sinceno train/test splits are provided.

the instances are not shared between different tasks for
the School dataset, MT-NN, MTFRL-NN and FETR-NN
cannot be directly applied. For each method, the bestmodel on the validation set is selected.

(a) Feature covariance matrix.

The results are summarized in Table 1 (the smaller the better). Among all the methods, FETR consistently achieves
lower test set NMSEs. Moreover, we observe a significant
improvement of both MT-NN, MTFRL-NN and FETR-NN
over all the linear baselines and STL-NN. FETR-NN further improves over MT-NN and MTFRL-NN on all thetasks. The experimental results confirm that multitask
learning usually improves over single task learning when
the dataset is small and tasks are related. Furthermore,among all the competitors, we observe that nonlinear
models combined with our FETR framework give theoverall best results, demonstrating the effectiveness of the
proposed approach in both linear and nonlinear settings.
One by-product of FETR is that we also have access tothe estimated row and column covariance matrices. In
Figure 4 we plot the feature and task covariance matrices
respectively, where we can clearly observe a block diagonal structure: the first 4 tasks are negatively correlatedwith the rest 3, and the 5th and 6th task are positively correlated. Intuitively, these correlations are consistent with
the SARCOS dataset where several joints move jointly.

(b) Task covariance matrix.

Figure 4: Estimated feature and task covariance matrices
on the SARCOS dataset.
To show the power of the nonlinear extension, we also run
experiments using a single task neural network (STL-NN)
and the multitask neural network (MT-NN), based on
which we propose our FETR-NN, which incorporates theregularization scheme into the last layer of the MT-NN.
STL-NN is a model where we use a separate networkfor each task, while in MT-NN all layers except the last
output layer are shared among different tasks. As anotherbaseline, we also compare our method with the multilinear relationship network (Long et al., 2017), which canbe understood as an extension of the MTFRL method
using neural networks (MTFRL-NN). In all experiments,
STL-NN, MT-NN, MTFRL-NN and FETR-NN share exactly the same network structure: an input layer with 21
dimensions, followed by two hidden layers with 256 and100 hidden units. The output of the network in MT-NN,
MTFRL-NN and FETR-NN is a multitask layer that contains 7 output units, while in STL-NN, the output only
contains a single unit. All the methods share the same
experimental setting, including model selection. In all theexperiments we fix l = 10−3 and u = 103 . The hyperparameters range from η ∈ {10−5 , . . . , 103 }, and we use
the validation set for model selection. Note that because

7

CONCLUSIONSIn this paper we point out a common flaw in the existing
multitask feature and relationship learning frameworks,
and propose a constrained variant to fix it. Our framework admits a multiconvex formulation, which allows us
to design an efficient block coordinate-wise algorithm tooptimize. To solve the weight learning subproblem, we
propose three different strategies that can be used no matter whether the instances are shared by multiple tasks ornot. To learn the covariance matrices, we reduce the underlying matrix optimization subproblem to a minimum
weight perfect matching problem, and solve it exactly
in closed form. To the best of our knowledge, all theprevious methods have to resort to expensive iterative procedures to solve this problem. At the end, we also discuss
several possible extensions of the proposed framework
to nonlinear settings. Experimental results show that our
method is orders of magnitude faster than its competitors, and it demonstrates significantly improved statistical
performance on two real-world datasets.
Acknowledgements
HZ and GG would like to acknowledge the support fromDARPA XAI project, contract FA87501720152. OS is
supported by a CMLH Fellowship in Digital Health and
by NIH under grant U01NS098969.

References
T. Adel, H. Zhao, and A. Wong. Unsupervised domain
adaptation with a relaxed covariate shift assumption.
In AAAI, pages 1691–1697, 2017.
A. Argyriou, M. Pontil, Y. Ying, and C. A. Micchelli. A
spectral regularization framework for multi-task structure learning. In Advances in Neural Information Processing Systems, 2007.
A. Argyriou, T. Evgeniou, and M. Pontil. Convex multitask feature learning. Machine Learning, 73(3):243–
272, 2008.
R. H. Bartels and G. Stewart. Solution of the matrix
equation AX+ XB= C [F4]. Communications of theACM, 15(9):820–826, 1972.
R. Bhatia and P. Rosenthal. How and why to solve theoperator equation ax- xb= y. Bulletin of the London
Mathematical Society, 29(01):1–21, 1997.
E. V. Bonilla, F. V. Agakov, and C. Williams. Kernel
multi-task learning using task-specific features. In International Conference on Artificial Intelligence andStatistics, pages 43–50, 2007.
S. Boyd and L. Vandenberghe. Convex optimization.
2004.
J.-F. Cai, E. J. Candès, and Z. Shen. A singular value
thresholding algorithm for matrix completion. SIAM
Journal on Optimization, 20(4):1956–1982, 2010.
R. Caruana. Multitask learning. Machine learning, 28(1):41–75, 1997.
J. Chen, J. Zhou, and J. Ye. Integrating low-rank andgroup-sparse structures for robust multi-task learning.
In Proceedings of the 17th ACM SIGKDD international
conference on Knowledge discovery and data mining,
pages 42–50. ACM, 2011.
J. Chen, J. Liu, and J. Ye. Learning incoherent sparse andlow-rank patterns from multiple tasks. ACM Transactions on Knowledge Discovery from Data (TKDD), 5(4):22, 2012.
F. Dufossé and B. Uçar. Notes on birkhoff–von neumann
decomposition of doubly stochastic matrices. Linear
Algebra and its Applications, 497:108–115, 2016.
P. Dutilleul. The MLE algorithm for the matrix normal
distribution. Journal of Statistical Computation and
Simulation, 64(2):105–123, 1999.
A. Evgeniou and M. Pontil. Multi-task feature learning.
Advances in neural information processing systems, 19:
41, 2007.
T. Evgeniou and M. Pontil. Regularized multi–task learning. In Proceedings of the tenth ACM SIGKDD international conference on Knowledge discovery and data
mining, pages 109–117. ACM, 2004.

X. Y. Felix, A. T. Suresh, K. M. Choromanski, D. N.
Holtmann-Rice, and S. Kumar. Orthogonal random
features. In Advances in Neural Information Processing
Systems, pages 1975–1983, 2016.
H. Glanz and L. Carvalho. An expectation-maximization
algorithm for the matrix normal distribution. arXiv
preprint arXiv:1309.6609, 2013.
H. Goldstein. Multilevel modelling of survey data. Journal of the Royal Statistical Society. Series D (The Statistician), 40(2):235–244, 1991.
A. K. Gupta and D. K. Nagar. Matrix variate distributions,volume 104. CRC Press, 1999.
L. Jacob, J.-p. Vert, and F. R. Bach. Clustered multi-task
learning: A convex formulation. In Advances in neural
information processing systems, pages 745–752, 2009.
P. Jawanpuria, M. Lapin, M. Hein, and B. Schiele. Efficient output kernel learning for multiple tasks. In
Advances in Neural Information Processing Systems,
pages 1189–1197, 2015.
T. Kato, H. Kashima, M. Sugiyama, and K. Asai. Multitask learning via conic programming. In Advancesin Neural Information Processing Systems, pages 737–
744, 2008.
C. Li, J. Zhu, and J. Chen. Bayesian max-margin multitask learning with data augmentation. In International
Conference on Machine Learning, pages 415–423,
2014.
Y. Li, X. Tian, T. Liu, and D. Tao. Multi-task model andfeature joint learning. In Proceedings of the 24th International Conference on Artificial Intelligence, pages
3643–3649. AAAI Press, 2015.
J. Liu, S. Ji, and J. Ye. Multi-task feature learning via
efficient l 2, 1-norm minimization. In Proceedings ofthe twenty-fifth conference on uncertainty in artificial
intelligence, pages 339–348. AUAI Press, 2009.
M. Long, Z. Cao, J. Wang, and S. Y. Philip. Learning
multiple tasks with multilinear relationship networks.
In Advances in Neural Information Processing Systems,
pages 1593–1602, 2017.
Y. Nesterov. Introductory lectures on convex optimization: A basic course, volume 87. Springer Science &Business Media, 2013.
A. Rahimi and B. Recht. Random features for large-scale
kernel machines. In Advances in neural information
processing systems, pages 1177–1184, 2008.
B. Roś, F. Bijma, J. C. de Munck, and M. C. de Gunst.
Existence and uniqueness of the maximum likelihood
estimator for models with a kronecker product covariance structure. Journal of Multivariate Analysis, 143:
345–361, 2016.

O. Sener and V. Koltun. Multi-task learning as multiobjective optimization. In Advances in Neural Information Processing Systems, pages 527–538, 2018.
J. R. Shewchuk et al. An introduction to the conjugate
gradient method without the agonizing pain, 1994.
S. Thrun. Explanation-based neural network learning. In
Explanation-Based Neural Network Learning, pages
19–48. Springer, 1996.
S. Vijayakumar and S. Schaal. Locally weighted projection regression: Incremental real time learning in highdimensional space. In Proceedings of the SeventeenthInternational Conference on Machine Learning, pages
1079–1086. Morgan Kaufmann Publishers Inc., 2000.
Y. Zhang. Supervision reduction by encoding extra information about models, features and labels. 2011.
Y. Zhang and J. G. Schneider. Learning multiple tasks
with a sparse matrix-normal penalty. In Advances in
Neural Information Processing Systems, pages 2550–
2558, 2010.
Y. Zhang and D.-Y. Yeung. A convex formulation
for learning task relationships in multi-task learning.
2010a.
Y. Zhang and D.-Y. Yeung. Multi-task learning using
generalized t process. In AISTATS, pages 964–971,
2010b.
H. Zhao, S. Zhang, G. Wu, J. P. Costeira, J. M. Moura,and G. J. Gordon. Multiple source domain adaptation
with adversarial training of neural networks. arXiv
preprint arXiv:1705.09684, 2017.
H. Zhao, S. Zhang, G. Wu, J. M. Moura, J. P. Costeira,and G. J. Gordon. Adversarial multiple source domain
adaptation. In Advances in Neural Information Processing Systems, pages 8559–8570, 2018.
H. Zhao, R. T. d. Combes, K. Zhang, and G. J. Gordon. On
learning invariant representation for domain adaptation.
In Proceedings of the 36th International Conference
on Machine Learning. ACM, 2019a.
H. Zhao, J. Hu, Z. Zhu, A. Coates, and G. Gordon. Deep
generative and discriminative domain adaptation. In
Proceedings of the 18th International Conference on
Autonomous Agents and MultiAgent Systems, pages
2315–2317. International Foundation for Autonomous
Agents and Multiagent Systems, 2019b.