\section{Introduction}

Current  machine learning systems  achieve remarkable results in several challenging  tasks, but are limited by the amount of human supervision required. Leveraging similarity among different problems is widely acknowledged to be  a key approach to reduce the need for supervised data. Indeed, this idea is at the basis of multi-task learning, where the joint solution of different  problems (tasks) has the potential to exploit tasks relatedness (structure) to improve learning accuracy. This idea has motivated a variety of methods, including frequentist \cite{micchelli04,argyriou08,argyriou08b} and Bayesian methods (see e.g. \cite{alvarez12} and references therein), with connections to structured learning~\cite{bakir07,tsochantaridis04}. \\
The focus of our study is the development of a general regularization  framework  to   learn multiple tasks as well as their structure.  Following~\cite{micchelli04,evgeniou05} we consider a setting where tasks are modeled as the components of a vector-valued function and their structure corresponds to the choice of suitable functional spaces. Exploiting the theory of reproducing kernel Hilbert spaces for vector-valued functions (RKHSvv)~\cite{micchelli04},  we consider and analyze  a flexible regularization framework, within which a variety of previously proposed approaches can be recovered as special cases, see e.g.  ~\cite{jacob08,lozano10,minh11,zhang10,dinuzzo11,sindhwani13}. Our main technical contribution  is a unifying  study of the minimization problem corresponding to such a  regularization framework. More precisely, we devise an optimization approach that can efficiently compute a solution and for which we prove convergence under weak assumptions. Our approach is based on a barrier method that is combined with block coordinate descent techniques \cite{tseng01,razaviyayn13}. In this sense our analysis generalizes the results in \cite{argyriou08} for which a low-rank assumption was considered; however the extension is not straightforward, since we consider a much larger class of regularization schemes (any convex penalty). Up to our knowledge, this is the first result in multi-task learning proving the convergence of alternating minimization schemes for such a general family of problems.\\
The RKHSvv setting allows to naturally deal both with linear and non-linear models and the approach we propose provides a general computational framework for learning output kernels as formalized in \cite{dinuzzo11}.\\
The rest of the paper is organized as follows: in Sec~\ref{sec:RKHSvv} we review basic ideas of regularization in RKHSvv. In Sec.~\ref{sec:unified_perspective} we discuss the equivalence of different approaches to encode known structures among multiple tasks. In Sec.~\ref{sec:learn_joint} we discuss a general framework for learning multiple tasks and their relations where we consider a wide family of structure-inducing penalties and study an optimization strategy to solve them. This setting allows us, in Sec.~\ref{sec:examples}, to recover several previous methods as special cases. Finally in Sec.~\ref{sec:experiments} we evaluate the performance of the optimization method proposed.

\paragraph{Notation.} With $S^n_{++} \subset S^n_+ \subset S^n \subset \mathbb{R}^{n \times n}$ we denote respectively the space of positive definite, positive semidefinite (PSD) and symmetric $n \times n$ real-valued matrices. $O^n$ denotes the space of orthonormal $n \times n$ matrices. For any square matrix $M\in\mathbb{R}^{n \times n}$ and $p\geq1$, we denote by $\|M\|_p = (\sum_{i=1}^n \sigma_i(M)^p)^{1/p}$ the $p$-Schatten norm of $M$, where $\sigma_i(M)$ is the $i$-th largest singular value of $M$. For any $M\in\mathbb{R}^{n \times m}$, $M^\top$ denotes the transpose of $M$. For any PSD matrix $A\in S_+^n$, $A^\dagger$ denotes the pseudoinverse of $A$.  We denote by $I_n\in S_{++}^n$ the $n \times n$ identity matrix. The notation $\ensuremath{\text{\rm Ran}}(M)\subseteq\mathbb{R}^m$ identifies the range of columns of a matrix $M\in\mathbb{R}^{m \times n}$.

\section{Background}\label{sec:RKHSvv}

We study the problem of jointly learning multiple tasks by modeling individual task-predictors as the components of a vector-valued function. Let us assume to have $T$ supervised scalar learning problems (or tasks), each with a ``training'' set of input-output observations $\mathcal{S}_t = \{(x_{it},y_{it})\}_{i=1}^{n_t}$ with $x_{it} \in \mathcal{X}$ input space and $y_{it}\in\mathcal{Y}$ output space\footnote{To avoid clutter in the notation, we have restricted ourselves to the typical situation where all tasks share same input and output spaces, i.e. $\mathcal{X}_t = \mathcal{X}$ and $\mathcal{Y}_t \subseteq \mathbb{R}$.}. Given a loss function $\mathcal{L}:\mathbb{R}\times\mathbb{R}\to\mathbb{R}_+$ that measures the per-task prediction errors, we want to solve the following joint regularized learning problem
\begin{equation}\label{eq:learning_problem}
\underset{f\in\mathcal{H}}{\text{minimize}} \ \ \sum_{t=1}^T \frac{1}{n_t}\sum_{i=1}^{n_t} \mathcal{L}(y_{i}^{(t)},f_t(x_{i}^{(t)})) + \lambda \|f\|_\mathcal{H}^2
\end{equation}
where $\mathcal{H}$ is an Hilbert space of vector-valued functions $f:\mathcal{X} \to \mathcal{Y}^T$with scalar components $f_t:\mathcal{X}\to\mathcal{Y}$. In order to define a suitable space of hypotheses $\mathcal{H}$, in this section we briefly recall concepts from the theory of reproducing kernel Hilbert spaces for vector-valued functions (RKHSvv) and corresponding regularization theory, which plays a key role in our work. In particular, we focus on a class of reproducing kernels (known as separable kernels) that can be designed to encode specific tasks structures (see ~\cite{evgeniou05,argyriou13} and Sec.~\ref{sec:unified_perspective}). Interestingly, separable kernels are related to ideas such as defining a metric on the output space or a label encoding in multi-label problems (see Sec.~\ref{sec:unified_perspective})

\begin{remark}[Multi-task and multi-label learning]
Multi-label learning is a class of supervised learning problems in which the goal is to associate input examples with a label or a set of labels chosen from a discrete set. In general, due to discrete nature of the output space, these problems cannot be solved directly; hence, a so-called {\it surrogate} problem is often introduced, which is computationally tractable and whose solution allows to recover the solution of the original problem~\cite{steinwart08,bartlett06,mroueh12}.\\
Multi-label learning and multi-task learning are strongly related. Indeed, surrogate problems typically consist in a set of distinct supervised learning problems (or tasks) that are solved simultaneously and therefore have a natural formulation in the multi-task setting. For instance, in multi-class classification problems the ``One vs All'' strategy is often adopted, which consists in solving a set of multiple binary classification problems, one for each class.
\end{remark}


\subsection{Learning Multiple Tasks with RKHSvv}


In the scalar setting, reproducing kernel Hilbert spaces have already been proved to be a powerful tool for machine learning applications. Interestingly, the theory of RKHSvv and corresponding Tikhonov regularization scheme follow closely the derivation in the scalar case. 

\begin{definition}
Let $(\mathcal{H},\langle\cdot,\cdot\rangle_{\mathcal{H}})$ be a Hilbert space of functions from $\mathcal{X}$ to $\mathbb{R}^T$. A symmetric, positive definite, matrix-valued function $\Gamma: \mathcal{X} \times \mathcal{X} \to \mathbb{R}^{T \times T}$ is called a reproducing kernel for $\mathcal{H}$ if for all $x \in \mathcal{X}, c \in \mathbb{R}^T$ and $f \in \mathcal{H}$ we have that $\Gamma(x,\cdot)c \in \mathcal{H}$ and the following reproducing property holds
$
\langle f(x), c\rangle_{\mathbb{R}^T} = \langle f,\Gamma(x,\cdot)c\rangle_{\mathcal{H}}.
$
\end{definition}
In analogy to the scalar setting, it can be proved (see~\cite{micchelli04}) that the Representer Theorem holds also for regularization in RKHSvv. In particular we have that any solution of the learning problem introduced in Eq.~\eqref{eq:learning_problem} can be written in the form 
\begin{equation}
f(x) = \sum_{t = 1}^T \sum_{i=1}^{n_t} \Gamma(x,x_{i}^{(t)}) c_{i}^{(t)}   
\end{equation}
with $c_i^{(t)} \in \mathbb{R}^T$ coefficient vectors.\\
The choice of kernel $\Gamma$ induces a joint representation of the inputs as well as a structure among the output components~\cite{alvarez12}; In the rest of the paper we will focus on so-called separable kernels, where these two aspects are factorized. In Section~\ref{sec:learn_joint}, we will see how separable kernels provide a natural way to learn the tasks structure as well as the tasks. 





\subsection{Separable Kernels}
Separable (reproducing) kernels are functions of the form  
$\Gamma(x,x')=k(x,x')A$ \ $\forall x,x'\in\mathcal{X}$ where $k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ is a scalar reproducing kernel and $A\in S_+^T$ is a positive semi-definite (PSD) matrix. In this case, the representer theorem allows to rewrite problem~\eqref{eq:learning_problem} in a more compact matrix notation as
\begin{equation}\label{eq:learning_problem_matrix}\tag{$\mathcal{P}$}
\underset{C\in\mathbb{R}^{n \times T}}{\text{minimize}} \ \ V(Y,KCA) + \lambda \ tr(AC^\top K C).
\end{equation}
Here $Y\in\mathbb{R}^{n \times T}$ is a matrix with $n = \sum_{t=1}^T n_t$ rows containing the output  points; $K\in S_+^n$ is the empirical kernel matrix associated to $k$ and $V:\mathbb{R}^{n \times T} \times \mathbb{R}^{n \times T} \to \mathbb{R}_+$ generalizes the loss in~\eqref{eq:learning_problem} and consists in a linear combination of the entry-wise application of $\mathcal{L}$. Notice that this formulation accounts also the situation where not all training outputs $y^{(t)}$ are observed when a given input $x \in \mathcal{X}$ is provided: in this case the functional $V$ weights $0$ the loss values of those entries of $Y$ (and the associated entries of $KCA$) that are not available in training.\\
Finally, the second term in~\eqref{eq:learning_problem_matrix} follows by observing that,  for all  $f\in\mathcal{H}$ of the form $f(\cdot) =\sum_{i=1}^n k(x_i,\cdot)A c_i$, the squared norm  can be written as $\|f\|_\mathcal{H}^2=\sum_{i,j}^n k(x_i,x_j) c_i^\top A c_j = tr(AC^\top KC)$ where $C\in\mathbb{R}^{n \times T}$ is the matrix with $i$-th row corresponding to the coefficient vector $c_i\in\mathbb{R}^T$ of $f$. Notice that we have re-ordered the index $i$ to be in $\{1,\dots,n\}$ to ease the notation.

\subsection{Incorporating Known Tasks Structure}\label{sec:unified_perspective}
Separable kernels provide a natural way to incorporate the task structure when the latter is known a priori. This strategy is quite general and indeed in the following we comment on how the matrix $A$ can be chosen to recover several multi-task methods previously proposed in contexts such as regularization, coding/embeddings or output metric learning, postponing a more detailed discussion in the supplementary material. These observations motivate the extension in Sec.~\ref{sec:learn_joint} of the learning problem~\eqref{eq:learning_problem_matrix} to a setting where it is possible to infer $A$ from the data.

\paragraph{Regularizers.}
Tasks relations can be enforced by devising suitable regularizers~\cite{evgeniou05}. Interestingly, for a large class of such methods it can be shown that this is equivalent to the choice of the matrix $A$ (or rather its pseudoinverse) \cite{micchelli04}. If we consider the squared norm of a function $f = \sum_{i=1}^n k(x_i,\cdot)Ac_i \in\mathcal{H}$ we have (see~\cite{evgeniou05})
\begin{equation}
\|f\|_\mathcal{H}^2 = \sum_{t,s=1}^T A^\dagger_{ts} \langle f_t, f_s \rangle_{\mathcal{H}_k}
\end{equation}
where $A_t$ is the $t$-th column of $A$, $\mathcal{H}_k$ is the RKHS associated to the scalar kernel $k$ and $f_t= \sum_{i=1}^n k(x_i,\cdot) A_t^\top c_i \in\mathcal{H}_k$  is the $t$-th component of $f$. The above equation suggests to interpret $A^\dagger$ as the matrix that models the structural relations between tasks by directly coupling different predictors. 
For instance, by setting $A^\dagger= I_T + \gamma (\mathbf{1}\mathbf{1}^\top)/T$, with $\mathbf{1}\in\mathbb{R}^T$ the vector of all $1$s, we have that the parameter $\gamma$ controls the variance $\sum_{t=1}^T \|\bar{f} - f_t\|_{\mathcal{H}_k}^2$ of the tasks with respect to their mean $\bar{f}=\frac{1}{T} \sum_{t=1}^T f_t$. If we have access to some notion of similarity among tasks in the form of a graph with adjacency matrix $W\in S^T$, we can consider the regularizer $\sum_{t,s=1}^T W_{t,s} \|f_t - f_s\|_{\mathcal{H}_k}^2 + \gamma \sum_{t}^T \|f_t\|_{\mathcal{H}_k}^2$ which corresponds to $A^\dagger=L + \gamma I_T$ with $L$ the graph Laplacian induced by $W$.
\paragraph{Output Metric.}
A different approach to model tasks relatedness consists in choosing a suitable  metric on the output space to reflect the tasks structure~\cite{lozano10}.  Clearly a change of metric on the output space with the standard inner product $\langle y,y\prime \rangle_{\mathbb{R}^T}$ between two output points $y,y\prime\in\mathcal{Y}^T$ corresponds to the choice of a different inner product $\langle y,y\prime\rangle_\Theta = \langle y, \theta y\prime \rangle_{\mathbb{R}^T}$ for some positive definite matrix $\Theta \in S_{++}^T$. Indeed this can be direct related to the choice of a suitable separable kernel.  In particular, for the least squares loss function a direct equivalence holds between choosing a metric deformation associated to a $\Theta\in S_{++}^T$ and a separable kernel $k(\cdot,\cdot)I_T$ or use the canonical metric (i.e. with $\Theta=I_T$ the identity) and kernel $k(\cdot,\cdot)\Theta$. The details of this equivalence can be found in the supplementary material.
\paragraph{Output Representation.}
The tasks structure can also be modeled by designing an ad-hoc embedding for the output space. This approach is particularly useful for multi-label scenarios, where output embedding can be designed to encode complex structures such as (e.g. trees, strings, graphs, etc.)~\cite{fergus10,joachims09,crammer00}. Interestingly in these cases, or more generally whenever the embedding map $L:\mathcal{Y}^T\to\widetilde{\mathcal{Y}}$, from the original to the new output space, is linear, then it is possible to show that the learning problem with new code is equivalent to~\eqref{eq:learning_problem} for a suitable choice of separable kernel with $A=L^\top L$. We refer again to the supplementary material for the details of this equivalence.




 



\section{Learning the Tasks and their Structure}\label{sec:learn_joint}
Clearly, an interesting setting occurs when knowledge of the tasks structure is not available and therefore it is not possible to design a suitable separable kernel. In this case a favorable approach is to infer the tasks relations directly from the data. To this end we propose to consider the following extension of problem~\eqref{eq:learning_problem_matrix}
\begin{equation}\label{eq:nonconvex}\tag{$\mathcal{Q}$}
\begin{aligned}
\underset{C\in\mathbb{R}^{n \times T}, A \in S_+^T}{\text{minimize}} & \ \  
V(Y,KCA) + \lambda tr(AC^\top KC) + F(A),
\end{aligned}
\end{equation}
where the penalty $F:S_+^T\to\mathbb{R}_+$ is designed to learn specific tasks structures encoded in the matrix $A$. The above regularization is general enough to encompass a large number of previously proposed approaches by simply specifying a choice of the scalar kernel and the penalty $F$. A detailed discussion of these connections is postponed to Section~\ref{sec:examples}.
In this section, we focus on computational aspects. Throughout, we restrict ourselves to convex loss functions $V$ and convex (and coercive) penalties $F$.  In this case, the objective function in~\eqref{eq:nonconvex} is separately convex  in  $C$ and $A$ but not  jointly convex. Hence,   block coordinate methods, which are often used in practice, e.g. alternating minimization over $C$ and $A$,  are not guaranteed to converge to a global minimum. Our  study provides a general  framework to provably compute a solution to problem~\eqref{eq:nonconvex}. First, In Section \ref{CMBM}, we prove our main results providing a  characterization  of the solutions of Problem~\eqref{eq:nonconvex} and  studying a barrier method to cast their computation as a convex optimization problem. Second, in Section \ref{BCM}, we discuss how block coordinate methods can be naturally used to solve such a problem, analyze their convergence properties and discuss some general cases of interest.

\subsection{Characterization of Minima and A Barrier Method}\label{CMBM}


We begin, in Section \ref{QR}, providing a characterization of the solutions to Problem \eqref{eq:nonconvex} by showing that it has an equivalent formulation in terms of the minimization of a convex objective function, namely Problem \eqref{eq:convex_equivalence}.  
Depending on the behavior of the objective function on  the boundary of the optimization domain,  
Problem \eqref{eq:convex_equivalence} might not be solved using standard optimization techniques. This possible issue motivates the introduction, in Section \ref{RSd}, of a barrier method; a family of ``perturbated'' convex programs is introduced whose solutions are shown to converge to those of Problem \eqref{eq:convex_equivalence} (and hence of the original \eqref{eq:nonconvex}).


\subsubsection{An Equivalent formulation for~\eqref{eq:nonconvex}}\label{QR}

The objective functional in~\eqref{eq:nonconvex} is not convex, therefore in principle it is hard to find a global minimizer. As it turns out however, it is possible to circumvent this issue and efficiently find a global solution to~\eqref{eq:nonconvex}. The following result represents a first step in this direction. 

\begin{theorem}\label{teo:convex_equivalence}
Let $K\in S_+^n$ and consider the convex set
$$
\mathcal{C}=\left\{(C,A) \in \mathbb{R}^{n \times T} \times S^T_+ \ | \ \ensuremath{\text{\rm Ran}}(C^\top KC) \subseteq \ensuremath{\text{\rm Ran}}(A) \right\}.
$$
Then, for any $F:S_+^T\to \mathbb{R}_+$ convex and coercive, problem
\begin{equation}\label{eq:convex_equivalence}\tag{$\mathcal{R}$}
\begin{aligned}
\underset{(C,A) \ \in \ \mathcal{C}}{\mbox{\emph{minimize}}} &V(Y,KC)+\lambda tr\left(A^{\dagger}C^\top K C\right) + F(A)
\end{aligned}
\end{equation}
has convex objective function and it is equivalent to \eqref{eq:nonconvex}. In particular, the two problems achieve the same minimum value and, given a solution $(C_R,A_R)$ for \eqref{eq:convex_equivalence}, the couple $(C_RA_R^\dagger,A_R)$ is a minimizer for \eqref{eq:nonconvex}. Vice-versa, given a solution $(C_Q,A_Q)$ for \eqref{eq:nonconvex}, the couple $(C_QA_Q,A_Q)$ is a minimizer for \eqref{eq:convex_equivalence}. 
\end{theorem}
The above result highlights a remarkable connection between the problems \eqref{eq:nonconvex} (non-convex) and \eqref{eq:convex_equivalence} (convex). In particular, we have the following Corollary, which provides us with a useful characterization of the local minimizers of problem~\eqref{eq:nonconvex}. 

\begin{corollary}\label{cor:invexity} 
Let $Q:\mathbb{R}^{n \times T} \times S_+^T \to \mathbb{R}$ be the objective  function of problem~\eqref{eq:nonconvex}. Then, every local minimizer for $Q$ on the open set $\mathbb{R}^{n \times T} \times S_{++}^T$ is also a global minimizer. 
\end{corollary}
Corollary~\ref{cor:invexity} follows from Theorem~\ref{teo:convex_equivalence} and the fact that, on the restricted domain $\mathbb{R}^{n \times T} \times S_{++}^T$, the map $Q$ is the combination of the objective functional of~\eqref{eq:convex_equivalence} and the invertible function $(C,A)\longmapsto(CA,A)$.  Moreover, if $Q$ is differentiable,  i.e. $V$ and the penalty $F$ are differentiable, this is exactly the definition of a \textit{convexifiable} function, which in particular implies  {\em invexity}~\cite{craven95}. 
The latter property ensures that, in the differentiable case, all the {\em stationary} points (rather than only local minimizers) are global minimizers. This result was originally proved in~\cite{dinuzzo11} for the special case of $V$ the least-squares loss and $F(\cdot)=\|\cdot\|_F^2$ the Frobenius norm; Here we have proved its generalization to all convex losses $V$ and penalties $F$.\\
We end this section adding two comments. First, we note that, while the objective function in Problem \eqref{eq:convex_equivalence} is convex, the corresponding minimization problem might not be a convex program  (in the sense that  the feasible set~$\mathcal{C}$ is not identified by a set of linear equalities and non-linear convex inequalities~\cite{boyd04}). Second, Corollary~\eqref{cor:invexity} holds only on the interior of the minimization domain $\mathbb{R}^{n \times T} \times S_+^T$ and does not characterize the behavior of the target functional on its boundary. In fact, one can see that both issues can be tackled defining a {\em perturbed} objective functional having a suitable behavior on the 
boundary of the minimization domain. This is the key motivation for the barrier method we discuss in the next section.


\subsubsection{A Barrier Method to Optimize~\eqref{eq:convex_equivalence}}\label{RSd}
 
Here we propose a barrier approach inspired by the work in~\cite{argyriou08} by introducing a perturbation of problem~\eqref{eq:convex_equivalence} that enforces the objective functions to be equal to $+\infty$ on the boundary of $\mathbb{R}^{n \times T} \times S_+^T$. As a consequence, each perturbed problem can be solved as a convex optimization constrained on a closed cone.
The latter comment is made more precise in the following result that we prove in the supplementary material.



\begin{theorem}\label{teo:perturbation}
Consider the family of optimization problems
\begin{equation}\label{eq:perturbation}\tag{$\mathcal{S}^\delta$}
\begin{aligned}
\underset{\substack{C\in\mathbb{R}^{n \times T}, \\ A\in S_{+}^T}}{\mbox{\emph{minimize}}} & V(Y,KC) +\lambda tr(A^{-1}(C^\top KC + \delta^2 I_T) ) + F(A)
\end{aligned}
\end{equation}
with $I_T \in S_{++}^T$ the identity matrix. Then, for each $\delta>0$ the problem~\eqref{eq:perturbation} admits a minimum. Furthermore, the set of minimizers for~\eqref{eq:perturbation} converges to the set of minimizers for~\eqref{eq:convex_equivalence} as $\delta$ tends to zero. More precisely, given any sequence $\delta_m>0$ such that $\delta_m\to0$ and a sequence of minimizers $(C_m,A_m)\in\mathbb{R}^{n \times T} \times S_{+}^T$ for~\eqref{eq:perturbation}, there exists a sequence $(C^*_m,A^*_m)\in\mathbb{R}^{n \times T} \times S_{+}^T$ of minimizers for~\eqref{eq:convex_equivalence} such that $\|C_m-C^*_m\|_F + \|A_m-A^*_m\|_F \to0$ as $m\to+\infty$.
\end{theorem}
The barrier $\delta^2 tr(A^{-1})$ is fairly natural and can be seen as preconditioning of the problem leading to favorable computations. The proposed barrier method is similar in spirit to the approach developed in~\cite{argyriou08} and indeed Theorem~\ref{teo:perturbation} and next Corollary~\ref{cor:bcd} are a generalization over the two main results in~\cite{argyriou08} to any convex penalty $F$ on the cone of PSD matrices. However, notice that since we are considering a much wider family of penalties (than the trace norm as in~\cite{argyriou08}) our results cannot directly derived from those in~\cite{argyriou08}. In the next section we discuss how to compute the solution of Problem \eqref{eq:perturbation} considering a block coordinate approach.







\subsection{Block Coordinate  Descent Methods}\label{BCM}

The characteristic block variable structure of the objective function in problem \eqref{eq:perturbation}, suggests that it might be beneficial to use  block coordinate methods (BCM) (see~\cite{beck11}) to solve it. Here with  BCM we identify a large class of  methods that, in our setting, iterate steps of an optimization on $C$, with $A$ fixed, followed by an optimization of $A$, for $C$ fixed.\\
A {\em meta} block coordinate algorithm to solve~\eqref{eq:perturbation} is reported in in Algorithm~\ref{alg:bcd}. Here we interpret each optimization step over $C$ as a supervised step, and each optimization step over $A$ as a  an unsupervised step  (in the sense that it involves the inputs but not the outputs). Indeed, when the structure matrix $A$ is fixed, problem~\eqref{eq:convex_equivalence} boils down to the standard supervised multi-task learning frameworks where a priori knowledge regarding the tasks structure is available. Instead, when the coefficient matrix $C$ is fixed, the problem of learning $A$ can be interpreted as an unsupervised setting in which the goal is to actually find the underlying task structure \cite{tenenbaum10}.\\
Several optimization methods can be used as procedures for both \textsc{SupervisedStep} and \textsc{UnsupervisedStep} in Algorithm~\ref{alg:bcd}.  In particular, a first class of methods is called Block Coordinate Descent (BCD) and identifies a wide class of iterative methods that perform (typically inexact) minimization of the objective function one block of variables at the time. Different strategies to choose which direction minimize at each step have been proposed: pre-fixed cyclic order, greedy search~\cite{razaviyayn13} or randomly, according to a predetermined distribution~\cite{nesterov12}. For a review of several BCD algorithms we refer the reader to~\cite{razaviyayn13} and references therein.\\
A second class of methods is called alternating minimization and corresponds to the situation where at each step in Algorithm~\ref{alg:bcd} and exact minimization is performed.  This latter approach is favorable when a closed form solution exists for at least one block of variables (see Section \ref{CF}) and has been  studied  extensively in \cite{tseng01} in the abstract setting where an oracle provides a block-wise minimizer at each iteration.
The following Corollary describes  the convergence properties of BCD and Alternate minimization sequences provided by applying Algorithm~\ref{alg:bcd} to~\eqref{eq:perturbation}.

\begin{algorithm}[t]
   \caption{\textsc{Convex Multi-task Learning}}
   \label{alg:bcd}
\begin{algorithmic}
   \State {\bfseries Input:} $K, Y,\epsilon$ tolerance, $\delta$ perturbation parameter, $S$ objective functional of~\eqref{eq:perturbation}, $V$ loss, $F$ structure penalty.
   \State {\bfseries Initialize:} $(C,A)=(C_0,A_0), t=0$ 
   \Repeat
  
  
 
   \State $C_{t+1} \gets$ \textsc{SupervisedStep} $(V,K,Y,C_{t},A_{t})$
    \State $A_{t+1} \gets$ \textsc{UnsupervisedStep}$(F,K,\delta,C_{t+1},A_{t})$
     \State $t \gets t+1$
  
  
  
   \Until{$| S(C_{t+1},A_{t+1})-S(C_{t},A_{t}) | < \epsilon$}
\end{algorithmic}
\end{algorithm}

\begin{corollary}\label{cor:bcd}
Let the Problem~\eqref{eq:perturbation} be defined as in Theorem~\ref{teo:perturbation} then:
\begin{itemize}
\item[(a)] \textbf{Alternating Minimization:} Let the two procedures in Algorithm~\ref{alg:bcd} each provide a block-wise minimizer of the functional with the other block held fixed.  Then every limiting point of a minimization sequence provided by Algorithm~\ref{alg:bcd}, is a global minimizer for~\eqref{eq:perturbation}.
\item[(b)] \textbf{Block Coordinate Descent:} Let the two procedures in Algorithm~\ref{alg:bcd} each consist in a single step of a first order optimization method (e.g. Projected Gradient Descent, Proximal methods, etc.). Then every limiting point of a minimizing sequence provided by Algorithm~\ref{alg:bcd} is a global minimizer for~\eqref{eq:perturbation}.
\end{itemize}
\end{corollary}


Corollary~\eqref{cor:bcd} follows by applying previous results on BCD and Alternate minimization. In particular, for the proof of part $(a)$ we refer to Theorem $4.1$ in ~\cite{tseng01}, while for part $(b)$ we refer to Theorem $2$ in~\cite{razaviyayn13}.\\
In the following we discuss the actual implementation of both \textsc{Supervised} and \textsc{Unsupervised} procedures in the case where $V$ is chosen to be least-squares loss and the penalty $F$ to be a spectral $p$-Schatten norm. This should provide the reader with a practical example of how the meta-algorithm introduced in this section can be specialized to a specific multi-task learning setting.

\begin{remark}(Convergence of Block Coordinate Methods)
Several works in multi-task learning have proposed some form of BCM strategy to solve the learning problem. However, up to our knowledge, so far only the authors in~\cite{argyriou08} have considered the issue of convergence to a global optimum. Their results where proved for a specific choice of structure penalty in a framework similar to that of problem~\eqref{eq:convex_equivalence} (see Section \ref{sec:examples}) but do not extend straightforwardly to other settings. Corollary~\ref{cor:bcd} aims to fill this gap, providing convergence guarantees for block coordinate methods for a large class of multi-task learning problems.
\end{remark}





\subsubsection{Closed Form solutions for Alternating Minimization: Examples}\label{CF}

Here we focus on the alternating minimization case and discuss some settings in which it is possible to obtain a closed form solution for the procedures \textsc{SupervisedStep} and \textsc{UnsupervisedStep}.
 
\paragraph{(\textsc{SupervisedStep}) Least Square Loss.}
When the loss function $V$ is chosen to be least squares (i.e. $V(Y,Z) = \|Y-Z\|_F^2$ for any two matrices $Y,Z \in \mathbb{R}^{n \times m}$) and the structure matrix $A$ is fixed, a closed form solution for the coefficient matrix $C$ returned by the \textsc{SupervisedStep} procedure can be easily derived (see for instance~\cite{alvarez12}):
$$
vec(C)=(I_T \otimes K+\lambda A^{-1} \otimes I_n)^{-1}vec(Y).
$$
Here, the symbol $\otimes$ denotes the Kronecker product, while the notation $vec(M) \in \mathbb{R}^{nm}$ for a matrix $M\in\mathbb{R}^{n \times m}$ identifies the concatenation of its columns in a single vector. In~\cite{minh11} the authors proposed a faster approach to solve this problem in closed form based on Sylvester's method.

\paragraph{(\textsc{UnsupervisedStep}) $p$-Schatten penalties.}\label{sec:p-schatten_penalties}
We consider the case in which $F$ is chosen to be a spectral penalty of the form $F(\cdot) = \|\cdot\|_p^p$ with $p\geq1$. Also in this setting the optimization problem has a closed form solution, as shown in the following.

\begin{proposition}\label{prop:p_solution}
Let the penalty of problem~\eqref{eq:perturbation} be $F = \|\cdot\|_p^p$ with $p\geq1$. Then, for any  $C\in\mathbb{R}^{n \times T}$ fixed, the optimization problem~\eqref{eq:perturbation} in the block variable $A$ has a minimizer of the form
\begin{equation}\label{eq:p_solution}
A_C^{\delta} = \sqrt[p+1]{(C^\top K C + \delta^2 I_T)/\lambda}.
\end{equation}
\end{proposition}

Proposition~\ref{prop:p_solution} generalizes a similar result originally proved in in~\cite{argyriou08} for the special case $p=1$ and provides an explicit formula for the \textsc{UnsupervisedStep} of Algorithm~\ref{alg:bcd}. We report the proof in the supplementary material.


















\section{Previous Work: Comparison and Discussion}\label{sec:examples}


The framework introduced in problem~\eqref{eq:nonconvex} is quite general and accounts for several choices of loss function and task-structural priors. Section~\ref{sec:learn_joint} has been mainly devoted to derive efficient and generic optimization procedures; in this section we focus our analysis on the modeling aspects, investigating the impact of different structure penalties on the multi-task learning problem. In particular, we will briefly review some multi-task learning method previously proposed, discussing how they can be formulated as special cases of problem~\eqref{eq:nonconvex} (or, equivalently, \eqref{eq:convex_equivalence}).

\paragraph{Spectral Penalties.} The penalty $F = \|\cdot\|_F^2$ was considered in~\cite{dinuzzo11}, together with a least squares loss function and the non convex problem~\eqref{eq:nonconvex} is solved directly by alternating minimization.
However, as pointed out in~Sec.~\ref{sec:learn_joint}, solving the non convex problem (although invex, see the discussion on Corollary~\ref{cor:invexity}) directly could in principle become problematic when the alternating minimization sequence gets close to the boundary of $\mathbb{R}^{n \times T} \times S_{++}^T$. A related idea is that of considering $F(A) = tr(A)$ (i.e. the $1$-Schatten norm). This latter approach can shown to be equivalent to the Multi-Task Feature Learning setting of~\cite{argyriou08} (see supplementary material).



\begin{figure*}
        \begin{center}
        
            \includegraphics[width=.24\textwidth]{MTCL_d.eps}
            \includegraphics[width=.24\textwidth]{MTCL_T.eps}
            \includegraphics[width=.24\textwidth]{MTFL_d.eps}
            \includegraphics[width=.24\textwidth]{MTFL_T.eps}
        \end{center}
       \caption{Comparison of the computational performance of the alternating minimization strategy studied in this paper with respect to the optimization methods proposed for MTCL in~\cite{jacob08} and MTFL~\cite{argyriou08} in the original papers. Experiments are repeated for different number of tasks and input-space dimensions as described in Sec.~\ref{sec:speed}.}\label{fig:speed}%
\end{figure*}

\paragraph{Cluster Tasks Learning.} In \cite{jacob08}, the authors studied a multi-task setting where tasks are assumed to be organized in a fixed number $r$ of unknown disjoint clusters. While the original formulation was conceived for linear setting, it can be easily extended to non-linear kernels and cast in our framework. Let $E\in\{0,1\}^{T \times r}$ be the binary matrix whose entry $E_{st}$ has value $1$ or $0$ depending on whether task $s$ is in cluster $t$ or not. Set $M=I - E^\dagger E^\top$, and $U=\frac{1}{T}11^{\top}$. In~\cite{jacob08} the authors considered a regularization setting of the form of~\eqref{eq:convex_equivalence} where the structure matrix $A$ is parametrized by the matrix $M$ in order to reflect the cluster structure of the tasks. More precisely:
$$A^{-1}(M)=\epsilon_{M}U+\epsilon_{B}(M-U)+\epsilon_{W}(I-M)$$
where the first term characterizes a global penalty on the average of all tasks predictors, the second term penalizes the between-clusters variance, and the third term controls the tasks variance within each cluster. Clearly, it would be ideal to identify an optimal matrix $A(M)$ minimizing problem~\eqref{eq:convex_equivalence}. However, $M$ belongs to a discrete non convex set, therefore authors propose a convex relaxation by constraining $M$ to be in a convex set $\mathcal{S}_c=\{ M \in S^{T}_{+} , 0\preceq M \preceq I, tr(M)=r \}$. In our notations $F(A)$ is therefore the indicator function over the set of all matrices $A = A(M)$ such that $M \in \mathcal{S}_c$. The authors propose a pseudo gradient descent method to solve the problem jointly.


\paragraph{Convex Multi-task Relation Learning.}
Starting from a multi-task Gaussian Process setting, in~\cite{zhang10}, authors propose a model where 
the covariance among the coefficient vectors of the $T$ individual tasks is controlled by a matrix $A\in S_{++}^T$ in the form of a prior. The initial maximum likelihood estimation problem is relaxed to a convex optimization with target functional of the form
\begin{equation}
\|Y - KC \|_F^2 + \lambda_1 \ tr(C^\top K C) + \lambda_2 \ tr(A^{-1} C^\top K C) 
\end{equation}
constrained to the set $\mathcal{A}=\{A \ | \ A\in S_{++}^T, tr(A)=1)$. This setting is equivalent to problem~\eqref{eq:convex_equivalence} (by choosing $F$ to be the indicator function of $\mathcal{A}$) with the addition of the term $tr(C^\top K C)$.












\paragraph{Non-Convex Penalties.}
Often times, interesting structural assumptions cannot be cast in a convex form and indeed several works have proposed non-convex penalties to recover interpretable relations among multiple tasks. For instance~\cite{argyriou13} requires $A$ to be a graph Laplacian, or~\cite{dinuzzo13} imposes a low-rank factorization of $A$ in two smaller matrices. In~\cite{mroueh11,kumar12} different sparsity models are proposed.\\
Interestingly, most of these methods can be naturally cast in the form of problem~\eqref{eq:nonconvex} or \eqref{eq:convex_equivalence}. Unfortunately our analysis of the barrier method does not necessarily hold also in these settings and therefore Alternating Minimization is not guaranteed to lead to a stationary point.



\section{Experiments}\label{sec:experiments}

We empirically evaluated the efficacy of the block coordinate optimization strategy proposed in this paper on both artificial and real datasets. Synthetic experiments were performed to assess the computational aspects of the approach, while we evaluated the quality of solutions found by the system on realistic settings.

\subsection{Computational Times}\label{sec:speed}

As discussed in Sec.~\ref{sec:examples}, several methods previously proposed in the literature, such as Multi-task Cluster Learning (MTCL) \cite{jacob08} and Multi-task Feature Learning (MTFL~\cite{argyriou08}]), can be formulated as special cases of problem~\eqref{eq:nonconvex} or \eqref{eq:convex_equivalence}. It is natural to compare the proposed alternating minimization strategy with the optimization solution originally proposed for each method. To assess the system's performance with respect to varying dimensions of the feature space and an increasing number of tasks, we chose to perform this comparison in an artificial setting.\\
We considered a linear setting where the input data lie in $\mathbb{R}^d$ and are distributed according to a normal distribution with zero mean and identity covariance matrix. $T$ linear models $w_t \in \mathbb{R}^d$ for $t=1,\dots,T$ were then generated according to a normal distribution in order to sample $T$ distinct training sets, each comprising of $30$ examples $(x_i^{(t)},y_i^{(t)})$ such that $y_i^{(t)} = \langle w_t, x_i^{(t)}\rangle + \epsilon$ with $\epsilon$ Gaussian noise with zero mean and $0.1$ standard deviation. On these learning problems we compared the computational performance of our alternating minimization strategy and the original optimization algorithms originally proposed for MTCL and MTFL and for which the code has been made available by the authors'. In our algorithm we used $A_0 = I$ identity matrix as initialization for the alternating minimization procedure. We used a least-squares loss for all experiments.\\
Figure~\ref{fig:speed} reports the comparison of computational times of alternating minimization and the original methods to converge to the same minima (of respectively the functional of MTCL and MTFL). We considered two settings: one where the number of tasks was fixed to $T=100$ and $d$ increased from $5$ to $150$ and a second one wher $d$ was fixed to $100$ and $T$ varied bewteen $5$ and $150$. To account for statistical stability we repeated the experiments for each couple $(T,d)$ and different choices of hyperparameters while generating a new random datasets at each time. We can make two observations from these results: 1) in the setting where $T$ is kept fixed we observe a linear increase in the computational times for both original MTCL and MTFL methods, while alternating minimization is almost constant with respect to the input space dimension. 2) When $d$ is fixed and the number of tasks increases, all optimization strategies require more time to converge. This shows that in general alternating minimization is a viable option to solve these problems and in particular, when $T << min(d,n)$ -- which is often the case in non-linear settings --this method is particularly efficient.

\begin{table*}[t]
\scriptsize
\begin{center}
\hspace*{-2cm}            
    \begin{tabular}{l c >{\columncolor{gray!35}}c c >{\columncolor{gray!35}}c c >{\columncolor{gray!35}}c c >{\columncolor{gray!35}}c c}
    
        & \multicolumn{2}{c}{50 tr. samples per class} & \multicolumn{2}{c}{100 tr. samples per class} & \multicolumn{2}{c}{150 tr. samples per class} & \multicolumn{2}{c}{200 tr. samples per class} &  \tstrut \bstrut \\ 
        & {\bf nMSE ($\pm$ std)} & \cellcolor{white} {\bf nI} & {\bf nMSE ($\pm$ std)} & \cellcolor{white} {\bf nI} & {\bf nMSE ($\pm$ std)} & \cellcolor{white} {\bf nI} & {\bf nMSE ($\pm$ std)} & \cellcolor{white} {\bf nI} & \tstrut \bstrut \\ 

        \specialrule{.1em}{.05em}{.0em} 

        {\bf STL}       & $0.2436 \pm 0.0268$ & $0$ & $0.1723 \pm 0.0116$ & $0$ & $0.1483 \pm 0.0077$ & $0$ & $0.1312 \pm 0.0021$ & $0$     & \tstrut \bstrut \\ 

        {\bf MTFL}      & $0.2333 \pm 0.0213$ & $0.0416$ & $0.1658 \pm 0.0107$ & $0.0379$ & $0.1428 \pm 0.0083$ & $0.0281$ & $0.1311 \pm 0.0055$ & $0.0003$ & \tstrut \bstrut \\ 

        {\bf MTRL}      & $0.2314 \pm 0.0217$ & $0.0404$ & $0.1653 \pm 0.0112$ & $0.0401$ & $0.1421 \pm 0.0081$ & $0.0288$ & $0.1303 \pm 0.0058$ & $0.0071$ & \tstrut \bstrut \\ 

        {\bf OKL}       & $0.2284 \pm 0.0232$ & $0.0630$ & $0.1604 \pm 0.0123$ & $0.0641$ & $\mathbf{0.1410 \pm 0.0087}$ & $0.0350$ & $0.1301 \pm 0.0073$ & $0.0087$ & \tstrut \bstrut \\ 

    \end{tabular}


\end{center}
\caption{Comparison of Multi-task learning methods on the Sarcos dataset. The advantage of learning the tasks jointly decreases as more training examples became available.}
\label{tab:sarcos}
\end{table*}




\subsection{Real dataset}


We assessed the benefit of adopting multi-task learning approaches on two real dataset. In particular we considered the following algorithms: Single Task Learning (STL) as a baseline, Multi-task Feature Learning (MTFL)~\cite{argyriou08}, Multi-task Relation Learning (MTRL)~\cite{zhang10}, Output Kernel Learning (OKL) \cite{dinuzzo11}. We used least squares loss for all experiments.

\paragraph{Sarcos.}

Sarcos\footnote{url{http://www.gaussianprocess.org/gpml/data/}} is a regression dataset designed to evaluate machine learning solutions for inverse dynamics problems in robotics. It consists in a collection of $21$-dimensional inputs, i.e. the joint positions, velocities and acceleration of a robotic arm with $7$ degrees of freedom and $7$ outputs (the tasks), which report the corresponding torques measured at each joint.\\
For each task, we randomly sampled $50,100,150$ and $200$ training examples while we kept a test set of $5000$ examples in common for all tasks. We used a linear kernel and performed $5$-fold crossvalidation to find the best regularization parameter according to the normalized mean squared error (nMSE) of predicted torques. We averaged the results over $10$ repetitions of these experiments. The results, reported in Table~\ref{tab:sarcos}, show clearly that to adopt a multi-task approach in this setting is favorable; however, in order to quantify more clearly such improvement, we report in Table~\ref{tab:sarcos} also the {\it normalized improvement} (\textit{nI}) over single-task learning (STL). For each multi-task method MTL, the normalized improvement nI(MTL) is computed as the average
$$
\mbox{nI(MTL)} = \frac{1}{n_{exp}} \sum_{i=1}^{n_{exp}} \frac{\mbox{nMSE}_i(\mbox{STL})-\mbox{nMSE}_i(\mbox{MTL})}{\sqrt{\mbox{nMSE}_i(\mbox{STL})\cdot\mbox{nMSE}_i(\mbox{MTL})}}
$$
over all the $n_{exp} = 10$ experiments of the normalized differences between the nMSE achieved by respectively the STL approach and the given multi-task method MTL.


\begin{table}
\scriptsize
\begin{center}
\rowcolors{3}{}{gray!35}

\begin{tabular}{lcccccc}
    &  \multicolumn{6}{c}{\bf Accuracy (\%) per \# tr. samples per class}  \tstrut \bstrut \\
    & \multicolumn{2}{c}{$50$} & \multicolumn{2}{c}{$100$} & \multicolumn{2}{c}{$150$} \tstrut \bstrut \\
   
    \specialrule{.1em}{.05em}{.0em} 


    {\bf STL}                    & $72.23$ & $\pm 0.04$ & $76.61$ & $\pm 0.02$ & $79.23$ & $\pm 0.01$ \tstrut \bstrut \\
                                                

    {\bf MTFL}                   & $73.23$ & $\pm .08$ & $77.24$ & $\pm .05$  & $80.11$ & $\pm .03$  \tstrut \bstrut \\
                                                
    {\bf MTRL}                   & $73.13$ & $\pm 0.08$ & $77.53$ & $\pm 0.04$ & $80.21$ & $\pm 0.05$ \tstrut \bstrut \\
                                                
    {\bf OKL}                    & $72.25$ & $\pm 0.03$ & $77.06$ & $\pm 0.01$ & $80.03$ & $\pm 0.01$ \tstrut \bstrut \\
                                                
   
   
    \end{tabular}
\end{center}
\caption{Classification results on the $15$-scene dataset. Four multi-task methods and the single-task baseline are compared.}
\label{tab:15scenes}
\end{table}



\paragraph{$15$-Scenes.}
$15$-Scenes\footnote{http://www-cvr.ai.uiuc.edu/ponce\_grp/data/} is a dataset designed for scene recognition, consisting in a $15$-class classification problem. We represented images using LLC coding~\cite{wang10} and trained the system on a training set comprising $50$, $100$ and $150$ examples per class. The test set consisted in $7500$ images evenly divided with respect to the $15$ scenes. Table~\ref{tab:15scenes} reports the mean classification accuracy on $20$ repetitions of the experiments. It can be noticed that while all multi-task approach seem to achieve approximately similar performance, these are consistently outperforming the STL baseline.









\section{Conclusions}

We have studied a general multi-task learning framework where the tasks structure can be modeled compactly in a matrix. For a wide family of models, the problem of jointly learning the tasks and their relations can be cast as a convex program, generalizing previous results for special cases~\cite{argyriou08,dinuzzo11}. Such an optimization can be naturally approached by block coordinate minimization, which can be seen as alternating between supervised and unsupervised learning steps optimizing respectively the tasks or their structure. We evaluated our method real data, confirming the benefit of multi-task learning when tasks share similar properties.\\
From an optimization perspective, future work will focus on studying the theoretical properties of block coordinate methods, in particular regarding convergence rates. Indeed, the empirical evidence we report suggests that similar strategies can be remarkably efficient in the multi-task setting. From a modeling perspective, future work will focus on studying wider families of matrix-valued kernels, overcoming the limitations of separable ones. Indeed, this would allow to account also for structures in the interaction space between the input and output domains jointly, which is not the case for separable models. 






{\small
\bibliographystyle{ieee}
