11institutetext: Diparimento di Elettronica, Informazione e Bioingegneria,
Politecnico di Milano, Milan, Italy 11email:
<EMAIL_ADDRESS>
# Interpetable Target-Feature Aggregation for Multi-Task Learning based on
Bias-Variance Analysis

Paolo Bonetti 11 Alberto Maria Metelli 11 Marcello Restelli 11

###### Abstract

Multi-task learning (MTL) is a powerful machine learning paradigm designed to
leverage shared knowledge across tasks to improve generalization and
performance. Previous works have proposed approaches to MTL that can be
divided into feature learning, focused on the identification of a common
feature representation, and task clustering, where similar tasks are grouped
together. In this paper, we propose an MTL approach at the intersection
between task clustering and feature transformation based on a two-phase
iterative aggregation of targets and features. First, we propose a bias-
variance analysis for regression models with additive Gaussian noise, where we
provide a general expression of the asymptotic bias and variance of a task,
considering a linear regression trained on aggregated input features and an
aggregated target. Then, we exploit this analysis to provide a two-phase MTL
algorithm (NonLinCTFA). Firstly, this method partitions the tasks into
clusters and aggregates each obtained group of targets with their mean. Then,
for each aggregated task, it aggregates subsets of features with their mean in
a dimensionality reduction fashion. In both phases, a key aspect is to
preserve the interpretability of the reduced targets and features through the
aggregation with the mean, which is further motivated by applications to Earth
science. Finally, we validate the algorithms on synthetic data, showing the
effect of different parameters and real-world datasets, exploring the validity
of the proposed methodology on classical datasets, recent baselines, and Earth
science applications.

###### Keywords:

Multi-Task Learning Variable Aggregation Bias-Variance.

##  1 Introduction

_Machine Learning_ [5, ML] approaches usually consider an individual learning
problem, eventually decomposing complex problems into independent tasks.
Inspired by the possibility of exploiting their interconnections, _Multi-task
learning_ [9, MTL] methods are designed to simultaneously learn multiple
related tasks, leveraging shared knowledge to improve performance and
generalization. In recent years, MTL has gained significant attention across
various domains, including natural language processing [11], computer vision
[3], and healthcare [32].

Focusing on supervised MTL, _task-clustering_ methods are designed to
aggregate different tasks into clusters, exploiting their relationships to
learn groups of tasks with the same model. On the other hand, _feature-based_
MTL approaches are focused on the identification of a subset of relevant
common input features (feature selection), or the extraction of a combination
of relevant original inputs (feature transformation). In the feature
transformation context, a _dimensionality reduction_ approach may be
considered to extract a set of relevant features common to all the tasks,
reducing the dimension of the feature space.

In this paper, we propose an MTL approach (NonLinCTFA) at the intersection
between task clustering and feature transformation with dimensionality
reduction. Firstly, we introduce a specific task clustering, based on
partitioning the targets, aggregating each obtained group of targets with
their mean. Then, for each aggregated task (the mean of the targets of a
cluster), we aggregate subsets of features with their mean. This way, we first
learn a single model for each group of tasks. Then, we aggregate subsets of
features with their mean in each aggregated task, producing a set of reduced
features in a dimensionality reduction fashion. In both cases, we provide
theoretical guarantees on the improvement (or not worsening) of the mean
squared error on each of the original tasks. A schematic view of this
methodology can be found in Figure 1(a). Additionally, Figure 1(b) shows a
variant that will be discussed in the next sections.

(a) Block scheme of the algorithm: partition of tasks, aggregation of each set
of targets with the mean, replication of the procedure to the $D$ features,
for each reduced task.

(b) Block scheme of a variant of the algorithm with $D$ features individually
measured for each task. In the aggregation phase, both targets and features
are averaged.

Figure 1: Block scheme of main algorithm and a homogeneous case variant.

The choice to consider the mean as an aggregation function is to preserve
interpretability, following the definition of interpretability as the property
of an ML pipeline to be understood by domain experts, without explanations by
ML experts (see [18] for a broader discussion). Indeed, we aggregate targets
and features, reducing the dimension and the variance of the models and
preserving the meaning of each aggregation, which is a mean of variables. On
the other hand, the aggregations induce an increase in bias that will be
controlled in the analysis.

(a) Motivational example.

(b) Basins reduced by NonLinCTFA.

Figure 2: Figure 2(a) represents some European sub-regions, identified with
different colors, where a target variable and meteorological features are
available. Figure 2(b) shows the aggregated regions after the application of
NonLinCTFA.

A motivational example is related to Earth science, where we may be interested
in the prediction of a target variable at different (correlated) locations,
given a set of meteorological features measured in each of them (see Figure
2(a)). In this context, the proposed algorithm aggregates highly correlated
targets and learns an individual model for their mean, resulting in a reduced
number of models, simplifying data representation (through further aggregating
highly correlated features) and enhancing the performances without loss of
interpretability, mitigating overfitting, and limiting the computational
complexity (see Figure 2(b), further discussed in the experimental section).
In this example, preserving the interpretability is essential for
climatologists, filling the gap between advanced ML methods (considered as
black-box algorithms) and their physical meaning.

Contributions After introducing notation, problem formulation, and related
works (Section 2), the contributions of this paper can be summarized as:

$1)$ General theoretical analysis of the asymptotic bias and variance,
estimating a target variable with a linear model, trained considering as
target an average of targets and as features a reduced set of basis functions
(Section 3).

$2)$ Introduction of a novel MTL approach, with theoretical guarantees,
aggregating targets and features with the mean for interpretability (Section
4).

$3)$ Extension of recent dimensionality reduction methods [7, 8], both to
multiple targets and generalizing the theoretical results in the single-task
case.

$4)$ Validation of the proposed algorithm on synthetic data, benchmark MTL
datasets and methods, and a multi-task dataset from Earth science (Section 5).

Technical proofs and additional results can be found in the Appendix.

##  2 Preliminaries

Notation Given $L$ learning tasks $\\{\mathcal{T}_{i}\\}_{i=1}^{L}$, each task
$\mathcal{T}_{i}$ is a supervised learning problem with training set
$S_{i}=\\{(\mathbf{x}^{i},y_{i})_{j}\\}_{j=1}^{n_{i}}$, where $n_{i}$ samples
for the $i$-th task are available,
$(\mathbf{x}^{i})_{j}=(x_{1}^{i},\dots,x_{D_{i}}^{i})_{j}\in\mathbb{R}^{D_{i}}$
is the $j$-th training sample associated to task $i$ composed of $D_{i}$
features, and $(y_{i})_{j}$ is the corresponding target. The random variable
$y_{i}$ represents the $i$-th target, $x_{k}^{i}$ denotes the $k$-th feature
of task $\mathcal{T}_{i}$, and the random vector $\mathbf{x^{i}}$ is the full
set of features of the $i$-th task. The symbols $\sigma^{2}_{a}$, $cov(a,b)$,
$\rho_{ab}$ and $\hat{\sigma}^{2}_{a}$, $\hat{cov}(a,b)$, $\hat{\rho}_{ab}$
denote the variance of a random variable $a$, its covariance, and correlation
with the random variable $b$ and their estimators, respectively. Finally,
$\mathbb{E}_{a}[f(a)]$ is the expected value of a function $f(\cdot)$ of the
random variable $a$ w.r.t. its distribution.

###### Remark 1

To simplify the notation, in this paper, we mainly consider a shared set of
$D$ features
$\bm{x}:=\\{x_{1},\dots,x_{D}\\}:=\\{x_{1}^{1},\dots,x_{D_{1}}^{1},\dots,x_{1}^{L},\dots,x_{D_{L}}^{L}\\}$
and the same number of training samples $n$ for each task. This simplification
motivates the workflow of Figure 1(a). We will also discuss a variant with $D$
task-specific features (workflow in Figure 1(b)). Recalling the motivational
example, this means that we will consider the full set of measurements of
meteorological features as a shared set of $D$ features. Then, we will discuss
a variant with a set of meteorological features measured in each individual
sub-region associated with its target. Additionally, this homogeneous-feature
variant could be applied in the heterogeneous settings, aggregating
corresponding features (e.g., temperature), without changing specific features
associated to a sub-region (e.g., some snow-related features may be available
only for a mountainous region).

Data generation process We consider the general relationship between the full
set of features and each target with additive Gaussian noise. The noise
variables $\\{\epsilon_{k}\\}_{k=1}^{L}$ can be correlated with each other,
but independent from features.

In the theoretical analysis, we consider a partition $\mathcal{P}$ of the
tasks. Therefore, each task $\mathcal{T}_{k}$ belongs to a set of the
partition $\mathcal{P}_{\iota(k)}$. Without loss of generality, in the
analysis we will focus on the $i$-th task $y_{i}=f_{i}(x_{1}\dots
x_{D})+\epsilon_{i},\ \epsilon_{i}\sim\mathcal{N}(0,\sigma_{i}^{2})$, and we
will denote with $\mathcal{P_{\iota}}:=\mathcal{P}_{\iota(i)}$ the set of
$\mathcal{P}$ that contains $\mathcal{T}_{i}$, and with $K_{\iota}$ its
cardinality. The aggregated (mean) target that contains $y_{i}$ is therefore
$\psi_{\iota}:=\frac{1}{K_{\iota}}\sum_{k:y_{k}\in\mathcal{P}_{\iota}}y_{k}$.
Finally, to simplify the computations, we assume the expected values of
features and targets to be zero:
$\mathbb{E}[x_{k}]=\mathbb{E}[y_{i}]=\mathbb{E}[f_{i}(x_{1}\dots x_{D})]=0,\
\forall k\in\\{1,\dots,D\\},\ \forall i\in\\{1,\dots,L\\}$.

Loss Measure As a natural performance measure for regression and in line with
related works [7, 8], in the theoretical analysis we evaluate the Mean Squared
Error (MSE), focusing on its bias and variance components (bias-variance
decomposition [13], see Equation 8 of Appendix 0.A).

Multi-Task Learning via aggregations In the theoretical analysis of Section 3,
we consider a task $y_{i}=f_{i}(x_{1}\dots x_{D})+\epsilon_{i},\
\epsilon_{i}\sim\mathcal{N}(0,\sigma_{i}^{2}),$ and we evaluate the effect (in
terms of MSE) of learning a linear model trained on a target $\psi_{\iota}$,
the mean of a given subset of original targets which contains $y_{i}$.
Additionally, $d$ given transformations $\\{\phi_{1},\dots,\phi_{d}\\}$ of the
$D$ original features are the inputs. They could be any zero-mean
transformations (basis functions), but we focus on the means of $d$ subsets of
original features. Indeed, through Algorithm 1, we exploit the theoretical
results to identify iteratively, and aggregate with the mean, convenient
partitions of features and targets in polynomial time.

######  Remark 2 (Limitations)

In the analysis, we focus on linear models to derive closed-form results.
However, the NonLinCTFA algorithm can be applied with any supervised learning
algorithm. In this case, its convenience becomes heuristic, and its
effectiveness should be tested empirically. The second main limitation is the
asymptoticity of the analysis, to identify expressions of variance and bias
without considering confidence intervals of each statistical estimator.

###  2.1 Related Works: Dimensionality Reduction, Multi-Task Learning

Dimensionality Reduction Considering each individual task, the aggregation of
its input features can be considered a dimensionality reduction method. These
approaches map $D$ features into a reduced dataset of dimension $d<D$, with
transformation functions aimed to maximize specific performance measures.
Linear, non-linear, supervised, and unsupervised methods exist. For broader
discussions and algorithms on dimensionality reduction, we refer to [27, 2,
15]. The specific idea to aggregate groups of features with their mean
resembles group regularization methods, such as Cluster LASSO [25] and OSCAR
[6], although the method proposed in this paper aggregates features
independently from the training of the supervised model. This is in line with
some recent approaches [7, 8], that we seek to generalize to a multi-task
context and by considering the full set of features to identify the
aggregations, rather than their bivariate analyses.

Multi-Task Learning The main algorithm of this paper is an MTL approach that
identifies subsets of targets, learning a single model for their mean, that we
claim to be convenient for the original individual tasks. A classical
description of MTL is [9], and more recent algorithms are in [31]. Following
their taxonomy, we briefly revise parameter-based and feature-based MTL
approaches. A few instance-based approaches also exist, designed to share
samples across tasks.

Parameter-based methods can be based on the assumption that similar tasks have
similar model parameters, forcing the coefficient matrix to be low-rank (low-
rank approaches, e.g., [30]); they can estimate and exploit task relationships
such as covariances (task-relation approaches, e.g., [34]); they can share
sets of parameters (decomposition approaches, e.g., [19]); or they can assume
that different tasks form several clusters (task-clustering approaches, e.g.,
[22]). Some recent parameter-based MTL algorithms aim to identify groups of
related tasks exploiting different metrics, such as the effect of the gradient
of the loss associated with one task on the loss of another task [12], or
summary statistics [17]. The approach of this paper can be considered a task-
clustering method since we identify disjoint groups of tasks, averaging the
corresponding targets.

Feature-based methods are based on the assumption that tasks share a common
feature representation. Some approaches learn this feature representation by
maximizing the information to each task (feature transformation). Among them,
deep learning-based approaches are the most commonly used (e.g., [29]). Other
approaches select subsets of original features, maximizing their relevance to
the targets (feature selection). They usually optimize a loss function that
both penalizes the selection of different features across tasks and minimizes
the number of important features (e.g., [33]). The proposed algorithm combines
the task-clustering phase with a feature transformation phase, where for each
cluster of tasks, we identify and aggregate subsets of features.

##  3 Bias-Variance Analysis: Theoretical Results

In this section, we present the main theoretical results, deriving the
asymptotic bias and variance of aggregated tasks. We consider a generic task
$\mathcal{T}_{i}$:

$y_{i}=f_{i}(x_{1}\dots x_{D})+\epsilon_{i},\
\epsilon_{i}\sim\mathcal{N}(0,\sigma_{i}^{2}),$ (1)

and the set of targets $\mathcal{P_{\iota}}$, that contains $y_{i}$. In this
setting, we estimate the target $y_{i}$ with a linear model, trained with the
mean of the elements of $\mathcal{P_{\iota}}$ as target, and $d$ inputs
$\\{\phi_{1},\dots,\phi_{d}\\}$. Recalling that
$\psi_{\iota}:=\frac{1}{K_{\iota}}\sum_{k:y_{k}\in\mathcal{P}_{\iota}}y_{k}$,
we estimate:

$\hat{y}_{i}=\hat{\psi}_{\iota}=\hat{w}_{1}^{\iota}\phi_{1}+\dots+\hat{w}_{d}^{\iota}\phi_{d},$
(2)

with $\hat{w}_{j}^{\iota}$ least squares estimates. In particular
$\\{\phi_{1},\dots,\phi_{d}\\}=\\{x_{1},\dots,x_{D}\\}$ is the case with
original features, and $\mathcal{P}_{\iota}=\\{y_{i}\\}$ identifies the
single-task case.

Variance We firstly derive the variance of the linear model
$\hat{\psi}_{\iota}$ (Equation 2), for the $i$-th target $y_{i}$. Proofs and
additional discussions are in Appendix 0.A.

###### Theorem 3.1

Let the relationship between the features and the target of a task
$\mathcal{T}_{i}$ be defined as Equation 1. Let also each estimator converge
in probability to the quantity that it estimates. In the asymptotic case, let
$var_{d}^{\iota}$ be the variance of a linear regression trained with the
basis functions $\\{\phi_{1},\dots,\phi_{d}\\}$ as inputs and the mean of a
cluster of targets $\psi_{\iota}$, defined in Equation 2, as output. It holds:

$var_{d}^{\iota}=\frac{\bar{\sigma}_{\iota}^{2}}{(n-1)}\cdot d,$ (3)

where
$\bar{\sigma}_{\iota}^{2}:=var(\frac{1}{K_{\iota}}\sum_{k:y_{k}\in\mathcal{P}_{\iota}}\epsilon_{k})$
is the variance of the mean of noises of the targets in $\mathcal{P}_{\iota}$,
$n$ is the number of training samples, and $d$ is the number of inputs.

###### Remark 3

The theorem follows the intuition that the asymptotic variance is proportional
to the number of inputs $d$ (i.e., the number of estimated coefficients) and
the variance gets smaller as the number of samples increases, since the
estimate is more accurate. Finally, the asymptotic variance increases with
$\bar{\sigma}_{\iota}^{2}$, which is the distortion associated with the
averaged target $\psi_{\iota}$.

###### Remark 4

Some specific cases are particularly relevant.

$1)$ When the inputs are the $D$ original features, we get
$var_{D}^{\iota}=\frac{\bar{\sigma}_{\iota}^{2}}{(n-1)}\cdot D$. For a single-
task approach, the variance is
$var_{d}^{i}=\frac{{\sigma}_{i}^{2}}{(n-1)}\cdot d$. Combining them, ($D$
features and a single-task), the variance is
$var_{D}^{i}=\frac{{\sigma}^{2}_{i}}{(n-1)}\cdot D$.

$2)$ When the variance of the noises is constant ($\forall
i:\sigma^{2}_{i}=\sigma^{2}$), we get:
$\bar{\sigma}_{\iota}^{2}=\frac{\sigma^{2}}{K_{\iota}}(1+(K_{\iota}-1)\bar{\rho}_{\iota})$,
with $\bar{\rho}_{\iota}$ average correlation among noises of the targets in
$\mathcal{P}_{\iota}$. Therefore, if $\bar{\rho}_{\iota}=1$, the task
aggregation does not reduce the variance since the variance of the averaged
noise is equal to the individual ($\sigma^{2}$). On the contrary, since
$\bar{\rho}_{\iota}\geq\frac{1}{1-K_{\iota}}$ by non-negativity of variance,
maximum variance gain corresponds to the minimum average correlation among
noises $\bar{\rho}_{\iota}=\frac{1}{1-K_{\iota}}$. Finally, when the noises
are independent ($\bar{\rho}_{\iota}=0$), the aggregated noise variance is
reduced by a factor $K_{\iota}$ ($\frac{1}{K_{\iota}}\sigma^{2}$).
Intuitively, when the noises are less correlated, the different tasks are
better exploited to refine the knowledge about each task.

$3)$ In [7, 8] a similar asymptotic result is provided in the single-task
setting, with an asymptotic variance equal to $\frac{2\sigma^{2}_{i}}{(n-1)}$
in the bivariate case and $\frac{\sigma^{2}_{i}}{(n-1)}$ in the univariate
one. Theorem 3.1 generalizes that result to $d$ input features, revealing a
cost of $\frac{\sigma^{2}}{n-1}$ for each feature considered. Therefore, by
reducing the features from $D$ to $d<D$, the asymptotic variance decrease is
$\frac{\sigma^{2}}{n-1}\cdot(D-d)$.

Bias We now focus on the bias of the linear model $\hat{\psi}_{\iota}$
(Equation 2), w.r.t. the $i$-th target $y_{i}$. Proofs and more discussions
can be found in Appendix 0.B.

###### Theorem 3.2

Let the relationship between the features and the target of a task
$\mathcal{T}_{i}$ be defined as Equation 1. Let also each estimator converge
in probability to the quantity that it estimates. In the asymptotic case, let
$bias_{d}^{\iota}$ be the bias of a linear regression trained with the basis
functions $\\{\phi_{1},\dots,\phi_{d}\\}$ as inputs and the mean of a cluster
of targets $\psi_{\iota}$, defined in Equation 2, as output. It is equal to:

$\begin{split}bias_{d}^{\iota}&=\sigma^{2}_{f_{i}}-\sigma^{2}_{\psi_{\iota}}R^{2}_{d,\iota}+2(cov(\psi_{\iota},f_{i}-\psi_{\iota}|\mathbf{\Phi})-cov(\psi_{\iota},f_{i}-\psi_{\iota})),\end{split}$
(4)

with $\sigma^{2}_{f_{i}}$ variance of the function of inputs defining the
$i$-th task, $\sigma^{2}_{\psi_{\iota}}$ variance of the aggregated target
$\psi_{\iota}$, $R^{2}_{d,\iota}$ squared coefficient of multiple correlation
between the $d$ inputs and the aggregated target $\psi_{\iota}$, and
$cov(\psi_{\iota},f_{i}-\psi_{\iota}|\mathbf{\Phi})$ partial covariance
between the two random variables $\psi_{\iota},f_{i}-\psi_{\iota}$, given the
inputs $\mathbf{\Phi}$.

###### Corollary 1

Considering the single-task setting, i.e., assuming $\psi_{\iota}=y_{i}$, the
asymptotic bias of Theorem 3.2 reduces to:

$bias_{d}^{i}=\sigma^{2}_{f_{i}}(1-R^{2}_{d,i}).$ (5)

###### Remark 5

Intuitively, Corollary 1 shows that, in single-task problems, the bias is
proportional to the information, measured through the coefficient of multiple
correlation, that the inputs share with the target. Additionally, the bias is
proportional to the variance of the function $f_{i}$ that regulates the $i$-th
data generation process. More generally, Theorem 3.2 shows that estimating the
$i$-th target with a linear model trained with the aggregated target
$\psi_{\iota}$, the bias is still proportional to the variance of the original
$i$-th target, which is an irreducible cost. Then, the bias decreases
proportionally to the coefficient of multiple correlation between aggregated
target and features, representing the skill of the linear model to predict the
aggregated target, weighted by the variance of the aggregated target itself,
which accounts for the simplification introduced by aggregating. Finally, the
third term of Equation 4 is the difference of the covariance, with and without
conditioning on the inputs, between the aggregated target and the gap between
the $i$-th task and the aggregated target itself. This way, if the features
reduce the information shared between $\psi_{\iota}$ and its gap with $f_{i}$,
the bias reduces since we are exploiting the features to improve the learning
of the $i$-th task.

###### Remark 6

Some specific cases can also be derived from the results on the bias.

$1)$ Considering the $D$ original features, the results of Theorem 3.2 and
Corollary 1 hold, evaluating the quantities of the expressions with them
(e.g., $R^{2}_{D,\iota}$).

$2)$ In [7, 8] a similar (bivariate) asymptotic analysis is provided in the
single-task setting. Corollary 1 extends those findings to a general case with
$d$ inputs. Therefore, if we perform a single-task dimensionality reduction,
aggregating $d$ sets of features with their mean, the asymptotic bias
variation is $\sigma^{2}_{f_{i}}(R^{2}_{D,i}-R^{2}_{d,i})$.

Theoretical bounds for aggregations We conclude this Section by showing a
condition for a convenient aggregation of two tasks and another for features
in single-task problems. Both results can be deduced from Theorem 3.1 and 3.2,
and they will be exploited in the two phases of the main algorithm,
respectively.

###### Corollary 2

Considering two tasks $\mathcal{T}_{i},\mathcal{T}_{j}$ regulated by Equation
1, training an individual linear model with the mean of the two targets as
output is profitable w.r.t. the individual single-task models, in terms of
MSE, if and only if:

$\begin{cases}\frac{{\sigma}^{2}_{i}}{(n-1)}\cdot
D+\sigma^{2}_{\psi_{\iota}}R^{2}_{D,\iota}\geq\frac{{\sigma}^{2}_{\iota}}{(n-1)}\cdot
D+\frac{1}{2}[\sigma^{2}_{f_{i}}R^{2}_{D,i}+\sigma^{2}_{f_{j}}R^{2}_{D,j}]\\\
\frac{{\sigma}^{2}_{j}}{(n-1)}\cdot
D+\sigma^{2}_{\psi_{\iota}}R^{2}_{D,\iota}\geq\frac{{\sigma}^{2}_{\iota}}{(n-1)}\cdot
D+\frac{1}{2}[\sigma^{2}_{f_{i}}R^{2}_{D,i}+\sigma^{2}_{f_{j}}R^{2}_{D,j}].\end{cases}$
(6)

###### Proof

We compute variance and bias from Theorem 3.1 and 3.2, substituting
$\psi_{\iota}=y_{i}$, $\psi_{\iota}=y_{j}$, and
$\psi_{\iota}=\frac{y_{i}+y_{j}}{2}$. Then, we impose that the sum of variance
and bias of the aggregated case is not worse than the individual one for both
models.

Following the intuition, Equation 6 shows the convenience of aggregating two
targets if the variance of the noise is reduced or when the predictive
capability of the aggregated model is better than the average of the single-
task ones.

###### Corollary 3

Considering a task $\mathcal{T}_{i}$, the aggregation of $D$ features into
$d<D$ aggregated ones is profitable, in terms of MSE of a linear model, if and
only if:

$\frac{\sigma^{2}_{i}}{(n-1)}\cdot(D-d)\geq\sigma^{2}_{f_{i}}(R^{2}_{D,i}-R^{2}_{d,i}).$
(7)

###### Proof

The result follows comparing the single-task bias and variance that are
particular cases of the general results of Theorem 3.1 and 3.2.

Following the intuition, Equation 7 shows the convenience in terms of variance
to reduce the number of features or in terms of bias when the predictive
capability of the linear regression is better when the features are
aggregated.

##  4 Multi-Task Learning via Aggregations: Algorithms

In this section, we present the NonLinCTFA algorithm, assuming $L$ tasks and
$D$ shared features. The algorithm is depicted in Figure 1(a), while a multi-
input variant, discussed in Remark 1, is depicted in Figure 1(b). Algorithm 1
reports the pseudo-code of NonLinCTFA, exploiting the loop of Algorithm 2 (in
Appendix 0.D) to iteratively aggregate targets firstly, and features
subsequently.

Phase I: task-aggregation. Firstly, the algorithm iteratively adds targets to
a set, forming a partition until no aggregation is convenient. At any
iteration, Equation 6 is exploited to test the convenience of adding a target
variable $y_{j}$ into a set of the partition, initialized as a singleton of a
random target. Additionally, a hyperparameter $\epsilon_{1}$ regulates the
propensity of the algorithm to aggregate.

Equation 6 regulates the convenience of an aggregation of two targets.
Therefore, the algorithm starts with $L$ individual tasks, and it identifies
two tasks that benefit from the aggregation (e.g.,
$\mathcal{T}_{1},\mathcal{T}_{2}$), moving to $L-1$ linear models. Then, the
algorithm further aggregates a third task (e.g., $\mathcal{T}_{3}$) with them,
moving to $L-2$ linear models, if the new aggregation
($\frac{y_{1}+y_{2}+y_{3}}{3}$ in the example) is convenient w.r.t. the
aggregate two-task model ($\frac{y_{1}+y_{2}}{2}$) and the univariate one
($y_{3}$). The new aggregation is therefore convenient w.r.t. the original
individual tasks, and it is further convenient w.r.t. the previous two-task
aggregation.

Similar to forward feature selection, we add a target in the current set, if
convenient, without inspecting all the possible combinations, which would be
combinatorial. This way, we do not identify the _optimal_ partition of
targets, but a convenient one w.r.t. the single-tasks, with a quadratic number
of comparisons in the number of tasks ($\mathcal{O}(L^{2})$) in the worst
case, i.e., with no aggregations.

Additionally, the aggregation depends on the ordering of targets. For this
reason, we randomize it to avoid systemic biases. A possible variation could
be to introduce a heuristic (e.g., the correlation between couples of targets)
to rank them and test for the aggregations based on this ranking.

Phase II: cluster-level feature-aggregation. In the second phase, NonLinCTFA
identifies groups of features to aggregate with their mean for each of the $l$
reduced tasks. This way, for each reduced task $\iota$, we identify
$d_{\iota}$ reduced features, reducing the dimension and improving the
performance.

Specifically, we exploit Equation 7 to iteratively identify couples of
features that are convenient to average, following the same iterative
procedure of the target aggregation phase. This is quadratic w.r.t. the number
of comparisons for each aggregated task ($\mathcal{O}(l\cdot D^{2})$).
Additionally, since the terms $\sigma^{2}_{i},\sigma^{2}_{f_{i}},n$ are
constant across different comparisons, we include them in the hyperparameter
$\epsilon_{2}$, which regulates the propensity of the algorithm to aggregate
features.

Algorithm 1 NonLinCTFA:Non-Linear Correlated Target-Feature Aggregation

features $\bm{x}=\\{x_{1}\dots x_{D}\\}$; targets $\bm{y}=\\{y_{1}\dots
y_{L}\\}$; $n$ samples, tolerances $\epsilon_{1},\epsilon_{2}$

reduced tasks $\\{\psi_{1},\dots,\psi_{l}\\}$, reduced features
$\\{\phi_{1},\dots,\phi_{d_{\iota}}\\}_{\iota=1}^{l}$

function
Compute_threshold_features$(\bm{x}_{curr},y,z_{\mathcal{P}},z_{j},\epsilon)$$\triangleright$
From Equation 7

$R_{sep}\leftarrow\text{R2score}(\bm{x}_{curr},y)$

$R_{aggr}\leftarrow\text{R2score}((\bm{x}_{curr}\setminus\\{z_{\mathcal{P}},z_{j}\\})\cup\\{mean(z_{\mathcal{P}},z_{j})\\},y)$

return $R_{sep}-R_{aggr}\leq\epsilon$

end function

function
Compute_threshold_targets$(\bm{x},y_{\mathcal{P}},y_{j},\epsilon)$$\triangleright$
From Equation 6

$y_{ag}\leftarrow mean(y_{\mathcal{P}},y_{j})$

$R_{\mathcal{P}},\sigma^{2}_{\mathcal{P}},\sigma^{2}_{f_{\mathcal{P}}}\leftarrow\text{R2score}(\bm{x},y_{\mathcal{P}}),\text{var\\_res}(\bm{x},y_{\mathcal{P}}),\text{var}(y_{\mathcal{P}})-\text{var\\_res}(\bm{x},y_{\mathcal{P}})$

$R_{j},\sigma^{2}_{j},\sigma^{2}_{f_{j}}\leftarrow\text{R2score}(\bm{x},y_{j}),\text{var\\_res}(\bm{x},y_{j}),\text{var}(y_{j})-\text{var\\_res}(\bm{x},y_{j})$

$R_{ag},\sigma^{2}_{ag},\sigma^{2}_{f_{ag}}\leftarrow\text{R2score}(\bm{x},y_{ag}),\text{var\\_res}(\bm{x},y_{ag}),\text{var}(y_{ag})-\text{var\\_res}(\bm{x},y_{ag})$

$threshold_{1}=\frac{D}{(n-1)}(\sigma^{2}_{ag}-\sigma^{2}_{\mathcal{P}})+\frac{1}{2}(R_{\mathcal{P}}\sigma^{2}_{f_{\mathcal{P}}}+R_{j}\sigma^{2}_{f_{j}})-R_{ag}\sigma^{2}_{f_{ag}}$

$threshold_{2}=\frac{D}{(n-1)}(\sigma^{2}_{ag}-\sigma^{2}_{j})+\frac{1}{2}(R_{\mathcal{P}}\sigma^{2}_{f_{\mathcal{P}}}+R_{j}\sigma^{2}_{f_{j}})-R_{ag}\sigma^{2}_{f_{ag}}$

return $(threshold_{1}\leq\epsilon)\ \mathbf{AND}\
(threshold_{2}\leq\epsilon)$

end function

function NonLinCTFA$(Input)$$\triangleright$ Main function

PHASE I: task aggregations

$\\{\psi_{1},\dots,\psi_{l}\\}\leftarrow$
Aggregation$(\bm{z}=\bm{y},phase=1,\epsilon=\epsilon_{1},\bm{x}=\bm{x},y=None)$

PHASE II: feature aggregation for each task

for each $\iota\in\\{1,\dots,l\\}$ do

$\\{\phi_{1},\dots,\phi_{d_{\iota}}\\}\leftarrow$
Aggregation$(\bm{z}=\bm{x},phase=2,\epsilon=\epsilon_{2},\bm{x}=None,y=\psi_{\iota})$

end for

return
$\\{\psi_{1},\dots,\psi_{l}\\},\\{\phi_{1},\dots,\phi_{d_{\iota}}\\}_{\iota=1}^{l}$

end function

###### Remark 7

Algorithm 1 outputs a set of reduced tasks and the associated sets of reduced
features. We do not include a final regression to decouple the algorithm from
the regression model, which can be run independently, with theoretical
guarantees in linear regression. In this sense, we propose a _filter_ method.

###### Remark 8

Recalling Remark 1, we may want to consider a multi-input setting, with $D$
features for each task individually
$\\{x_{1}^{i},\dots,x_{D}^{i}\\}_{i=1}^{L}$. In our example, they can be $D$
meteorological features measured at each location and associated with its
specific target. A variant of Algorithm 1 (Figure 1(b)) compares, at each
iteration, the model trained with the features associated with each individual
target and the model for the aggregated task, trained on averaged features.
This way, in a single phase, we aggregate couples of targets and pairwise
couples of features. In our example, we would compare individual measurements
of temperature and precipitation for single-task problems and the means of
temperatures and precipitations as the two inputs of the averaged model. This
variant can be extended to a heterogeneous case where, for example, the first
task has a snow-related feature not available for the other. In this case, we
consider it when the first task appears, keeping unchanged the aggregation of
targets and other features.

##  5 Experimental Validation

In this section, we show synthetic experiments, validating the proposed
algorithm w.r.t. single-task regressions, showing its behavior varying
parameters, and with an ablation study of its two phases. Then, applications
to real-world data show the competitiveness of the method w.r.t. single task
and benchmark MTL approaches. Code can be found at:
https://github.com/PaoloBonettiPolimi/NonLinCTFA.

###  5.1 Synthetic Experiments and Ablation Study

We start with synthetic experiments, validating the NonLinCTFA against single
tasks. Specifically, we consider $L=10$ tasks, $D=100$ features shared across
all tasks, $n=250$ training samples (same number for testing), the standard
deviation of (independent) noises $\sigma=10$, and hyperparameters
$\epsilon_{1}=0,\epsilon_{2}=0.0001$. Each target is obtained as a linear
combination of all the features, with additive Gaussian noise, randomly
sampling each coefficient from a uniform distribution in the interval
$[-1,-0.5]$ or $[0.5,1]$, obtaining two groups of tasks similar among
themselves. We perform linear regression on individual tasks and on aggregated
tasks after the first phase of the algorithm or fully applying NonLinCTFA. We
repeat the experiment $10$ times to produce confidence intervals, and we
consider as metrics the MSE and the coefficient of determination ($R^{2}$),
both in terms of absolute values and of percentage increase w.r.t. single-
task. We obtain a single-task average $R^{2}$ score of $0.48\pm 0.02$, which
increases to $0.64\pm 0.01(+33.45\pm 3.31\%)$ considering the aggregated tasks
and $0.67\pm 0.01(+39.82\pm 3.42\%)$ adding the feature aggregation.
Similarly, the MSE is equal to $9.34\pm 0.04$, reducing to $6.54\pm
0.03(-29.44\pm 1.45\%)$ and $5.98\pm 0.06(-35.36\pm 1.27\%)$. The $L=10$ tasks
become $l=2.5\pm 0.6$, and the $D=100$ features reduce to $d=3.43\pm 1.76$.
These results empirically validate the improvement provided by the proposed
algorithm w.r.t. the single-task counterparts in linear regression.
Additionally, they show a significant benefit with task-aggregation, with a
subsequent feature-aggregation phase refining the performances, with the added
value of simplifying models.

Figure 3: Test MSE (first row) and corresponding percentage decrease w.r.t.
single-task (second row), varying one parameter at a time, only aggregating
targets (Phase I) or adding the feature-aggregation phase (Phase I+II).

We further tested the proposed approach, varying one parameter at a time and
fixing the others. Figure 3 reports the results, in terms of MSE and of
percentage increase of MSE w.r.t. the single-task, showing the improvement
with task-aggregation and then the combination with feature aggregation. The
first column of the figure shows that, with a small number of samples, both
phases provide significant improvement, which is mitigated by a large number
of samples. In the second column, when the number of features increases, the
feature aggregation phase is more relevant. Similarly, the third column shows
that the task-aggregation phase is more relevant when the number of tasks
increases, assuming a constant number of features. Finally, the fourth column
shows that the problem becomes more difficult when the noise of the targets
increases, with the task aggregations that become more relevant. Figure 4 in
Appendix 0.E.1 also shows the behavior of the MSE when the hyperparameters
$\epsilon_{1},\epsilon_{2}$ are varied. In particular, $\epsilon_{1}$
represents the propensity to aggregate tasks, spanning from the single-task
case to a single aggregated task, with a value equal to $0$ that balances the
aggregation. Similarly, $\epsilon_{2}$ regulates the propensity to aggregate
features, spanning from preserving all original inputs to a single
aggregation.

###  5.2 Real World Datasets

As a first investigation on real data, we consider two classical datasets,
comparing single-task regression performance, applying NonLinCTFA, and with
some standard MTL approaches. In particular, the SARCOS dataset [28] is
composed of $7$ regression tasks, $21$ shared features, and $n=48933$ samples.
We consider $n=1000$ and the full set of $n=48933$ samples to investigate the
effect of the number of samples. Then, the _School_ dataset [4] consists of
$15362$ samples, exam records, distributed across $139$ tasks (schools), with
$27$ features. This is an example where different measures of the same
features are associated to their tasks, allowing to apply the variant of
NonLinCTFA discussed in Remark 8. In this case, we consider $27$ tasks and
$27$ features, as well as the full set of $139$ tasks.

We apply NonLinCTFA, combined with linear regression (LR), support vector
regression (SVR), and multi-layer perceptron (MLP), in comparison with the
corresponding single-task models. Additionally, we consider a random
prediction as a trivial baseline and Alternating Structural Optimization [1,
ASO], its convex relaxation [10, cASO], and Convex Clustered MTL [14, CMTL],
as benchmarks representing feature-based and clustered-based MTL methods
(adapting the implementation of https://github.com/chcorbi/MultiTaskLearning).
In line with the implementation of these baselines, we evaluate test
performance with normalized root mean squared error (NRMSE), randomly
selecting the $30\%$ of samples, with five different seeds, to produce test
confidence intervals.

Table 1: Experiments on SARCOS and _School_ datasets, considering 70% of data
for training and 30% for testing, randomizing over $5$ seeds for confidence
intervals. NRMSE is the performance measure (lower is better, best result in
bold).

| SARCOS_1000samples | SARCOS_full | School_27tasks | School_full
---|---|---|---|---
# samples $n$ | $1000$ | $48933$ | $15362$ | $15362$
# tasks $L$ | $7$ | $7$ | $27$ | $193$
Reduced # tasks (ours) | $6.0\pm 0.0$ | $6.0\pm 0.0$ | $6.0\pm 0.0$ | $23.2\pm 2.7$
# features $D$ | $21$ | $21$ |  $27$ (for each task) |  $27$ (for each task)
Reduced # features (ours) | $7.1\pm 0.3$ | $6.9\pm 0.4$ |  $27$ (for each aggregation) |  $27$ (for each aggregation)
NRMSE random | $0.363\pm 0.015$ | $0.235\pm 0.011$ | $0.634\pm 0.123$ | $0.641\pm 0.294$
NRMSE single-task LR | $0.085\pm 0.001$ | $0.069\pm 0.002$ | $0.183\pm 0.002$ | $0.165\pm 0.003$
NRMSE single-task SVR | $0.142\pm 0.010$ | $0.092\pm 0.008$ | $0.181\pm 0.009$ | $0.162\pm 0.005$
NRMSE single-task MLP | $0.069\pm 0.002$ | $0.045\pm 0.003$ | $0.182\pm 0.003$ | $0.174\pm 0.008$
NRMSE ASO | $0.075\pm 0.002$ | $0.049\pm 0.001$ | $0.172\pm 0.003$ | $\mathbf{0.152\pm 0.006}$
NRMSE cASO | $0.068\pm 0.001$ | $0.048\pm 0.001$ | $0.173\pm 0.002$ | $0.154\pm 0.006$
NRMSE CMTL | $0.111\pm 0.001$ | $0.067\pm 0.002$ | $0.843\pm 0.316$ | $1.187\pm 0.562$
NRMSE NonLinCTFA + LR (ours) | $0.054\pm 0.003$ | $0.035\pm 0.002$ | $0.162\pm 0.005$ | $0.159\pm 0.002$
NRMSE NonLinCTFA + SVR (ours) | $0.154\pm 0.021$ | $0.115\pm 0.024$ | $\mathbf{0.159\pm 0.003}$ | $0.155\pm 0.004$
NRMSE NonLinCTFA + MLP (ours) | $\mathbf{0.049\pm 0.009}$ | $\mathbf{0.031\pm 0.001}$ | $0.160\pm 0.004$ | $0.167\pm 0.003$

The results of Table 1 show that NonLinCTFA outperforms single-task models and
is competitive w.r.t. MTL baselines, with the advantage of reducing the number
of models and parameters, still preserving the interpretability.

In a second real-world application, we consider the QM9 dataset [24], a
challenging quantum chemistry dataset with $L=19$ tasks (properties of
molecules), with $139000$ graph inputs (molecules structures). We averaged
node features, position, and edge attributes, obtaining $D=19$ features.
Following the experimental setup of [23], we retrieved the dataset from
PyTorch Geometric, randomly selecting $10000$ samples for testing and the
others for training, normalizing each task and repeating the experiments three
times to produce confidence intervals.

Table 2: Experiments on QM9 dataset. $10000$ random samples are used for
testing ($\sim 138000$ for training), $3$ different seeds for confidence
intervals. MSE is the performance measure (lower is better). Training time is
also reported.

QM9 Test Results | # reduced tasks | # reduced features | MSE | Time (minutes)
---|---|---|---|---
Single-task LR | $19$ | $20$ | $0.969\pm 0.049$ | $\sim 1$
Single-task MLP | $19$ | $20$ | $0.518\pm 0.032$ | $\sim 11$
Single-task baseline of [26] | $19$ | $13$ | $0.533\pm 0.041$ | NA
HPS GNN + RLW | $19$ |  $11+$graph | $0.619\pm 0.254$ |  $\sim 1\times 300$epochs
Best GNN of [26] | $19$ |  $11+$graph | $0.216\pm 0.009$ | NA
NonLinCTFA + LR (ours) | $12\pm 1.2$ | $10.05\pm 2.96$ | $0.955\pm 0.038$ | $\sim 2$
NonLinCTFA + MLP (ours) | $12\pm 1.2$ | $10.05\pm 2.96$ | $0.469\pm 0.024$ | $\sim 14$

Table 2 shows the performance, in terms of test MSE, of the single-task linear
regression (LR) and multi-layer perceptron (MLP), together with the MSE
associated to the same models, applying NonLinCTFA. Additionally, we exploited
the implementation of LibMTL [21] of some state-of-the-art MTL methods and its
integration with the QM9 dataset for further benchmarking. In particular, the
architecture of the library compatible with graph neural networks is the hard
parameter sharing (HPS GNN) that we trained with all the $14$ weighting
strategies implemented in the library (we refer to it for details), up to
$300$ epochs, identifying the Random Loss Weighting strategy [20, RLW] as best
performing. For further comparison with the literature, in the table, we also
show the MSE of the baseline and the best-performing graph neural network
(GNN) proposed in [26]. Together with the MSE, Table 2 shows the computational
time. We can conclude that the NonLinCTFA provides significant aggregations,
improving single task performances, especially combined with the MLP.
Additionally, it is competitive w.r.t. GNN-based approaches, without
outperforming all of them given its much simpler tabular methodology, as also
highlighted by the computational time.

A final experimental setup shows an application of the NonLinCTFA Algorithm to
meteorological data, as depicted in Figure 2, where $L=29934$ European
hydrological basins are considered as tasks, each with a satellite signal as
target and $D=16$ meteorological measurements and climate indices as inputs.
In this context, we apply the variant of the algorithm described in Remark 8,
with linear regression (given the small amount of monthly measurements $n=102$
and the necessity to preserve the interpertability of the entire workflow).

Table 3: Experiments on climate dataset with NonLinCTFA, varying the
hyperparameter $\epsilon$. The number of aggregations and MSE are obtained
cross-validating.

NonLinCTFA |  $\epsilon=1$ (single-task) | $\epsilon=0.5$ | $\epsilon=0.1$ | $\epsilon=0.05$ | $\epsilon=0.01$
---|---|---|---|---|---
# reduced tasks | $29934$ | $28024$ | $2354$ | $1345$ | $969$
MSE | $1.058\pm 0.224$ | $1.045\pm 0.221$ | $0.758\pm 0.148$ | $0.755\pm 0.143$ | $0.753\pm 0.143$
| $\epsilon=0$ | $\epsilon=-0.01$ | $\epsilon=-0.05$ | $\epsilon=-0.1$ | $\epsilon=-1$
# reduced tasks | $944$ | $844$ | $680$ | $252$ | $1$
MSE | $0.750\pm 0.140$ | $0.746\pm 0.142$ | $0.756\pm 0.146$ | $0.909\pm 0.189$ | $1.132\pm 0.193$

Table 3 shows confidence intervals, in terms of MSE, associated with different
values of the hyperparameter $\epsilon$. As expected, the MSE on the original
tasks reduces when reducing the value of $\epsilon$ since the aggregations of
tasks are convenient. However, when the hyperparameter is too small, the
algorithm aggregates too many tasks, becoming detrimental to the understanding
of the behavior of the original tasks until the limit case of a single
aggregated task. The average number of reduced tasks is also reported in the
table to confirm this behavior.

##  6 Conclusions and Future Developments

In this paper, we introduced a two-phase MTL approach that aggregates sets of
targets and features with their mean, motivated by meteorological
applications, and aimed to improve the final regression performance,
preserving interpretability. We provided a bias-variance analysis, considering
linear regression, and we empirically validated the approach with synthetic
and real-world datasets, showing promising results also outside the context of
linear regression. A future development can be an extension of the analysis to
general ML models. Additionally, an applicative work with meteorological data
is under development.

## References

  * [1] Ando, R.K., Zhang, T., Bartlett, P.: A framework for learning predictive structures from multiple tasks and unlabeled data. Journal of machine learning research 6(11) (2005)
  * [2] Ayesha, S., Hanif, M.K., Talib, R.: Overview and comparative study of dimensionality reduction techniques for high dimensional data. Information Fusion 59, 44–58 (2020)
  * [3] Bachmann, R., Mizrahi, D., Atanov, A., Zamir, A.: Multimae: Multi-modal multi-task masked autoencoders. In: European Conference on Computer Vision. pp. 348–367. Springer (2022)
  * [4] Bakker, B., Heskes, T.: Task clustering and gating for bayesian multitask learning (2003)
  * [5] Bishop, C.M.: Pattern Recognition and Machine Learning (Information Science and Statistics). Springer-Verlag, Berlin, Heidelberg (2006)
  * [6] Bondell, H.D., Reich, B.J.: Simultaneous regression shrinkage, variable selection, and supervised clustering of predictors with oscar. Biometrics 64(1), 115–123 (2008)
  * [7] Bonetti, P., Metelli, A.M., Restelli, M.: Interpretable linear dimensionality reduction based on bias-variance analysis. arXiv preprint arXiv:2303.14734 (2023)
  * [8] Bonetti, P., Metelli, A.M., Restelli, M.: Nonlinear feature aggregation: Two algorithms driven by theory. arXiv preprint arXiv:2306.11143 (2023)
  * [9] Caruana, R.: Multitask learning. Machine learning 28, 41–75 (1997)
  * [10] Chen, J., Tang, L., Liu, J., Ye, J.: A convex formulation for learning shared structures from multiple tasks. In: Proceedings of the 26th annual international conference on machine learning. pp. 137–144 (2009)
  * [11] Chen, Q., Zhuo, Z., Wang, W.: Bert for joint intent classification and slot filling. arXiv preprint arXiv:1902.10909 (2019)
  * [12] Fifty, C., Amid, E., Zhao, Z., Yu, T., Anil, R., Finn, C.: Efficiently identifying task groupings for multi-task learning. Advances in Neural Information Processing Systems 34, 27503–27516 (2021)
  * [13] Hastie, T., Tibshirani, R., Friedman, J.: The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer series in statistics (2009)
  * [14] Jacob, L., Vert, J.p., Bach, F.: Clustered multi-task learning: A convex formulation. Advances in neural information processing systems 21 (2008)
  * [15] Jia, W., Sun, M., Lian, J., Hou, S.: Feature dimensionality reduction: a review. Complex & Intelligent Systems 8, 2663 – 2693 (2022)
  * [16] Johnson, R., Wichern, D.: Applied Multivariate Statistical Analysis. Pearson Prentice Hall (2007)
  * [17] Knight, P., Duan, R.: Multi-task learning with summary statistics. Advances in Neural Information Processing Systems 36 (2024)
  * [18] Kovalerchuk, B., Ahmad, M.A., Teredesai, A.: Survey of explainable machine learning with visual and granular methods beyond quasi-explanations. Interpretable artificial intelligence: A perspective of granular computing pp. 217–267 (2021)
  * [19] Li, Y., Yan, H., Jin, R.: Multi-task learning with latent variation decomposition for multivariate responses in a manufacturing network. IEEE Transactions on Automation Science and Engineering 20, 285–295 (2023)
  * [20] Lin, B., Ye, F., Zhang, Y., Tsang, I.W.: Reasonable effectiveness of random weighting: A litmus test for multi-task learning. arXiv preprint arXiv:2111.10603 (2021)
  * [21] Lin, B., Zhang, Y.: Libmtl: A python library for multi-task learning. ArXiv abs/2203.14338 (2022)
  * [22] Liu, A., Su, Y., Nie, W., Kankanhalli, M.: Hierarchical clustering multi-task learning for joint human action grouping and recognition. IEEE Transactions on Pattern Analysis and Machine Intelligence 39, 102–114 (2017)
  * [23] Navon, A., Shamsian, A., Achituve, I., Maron, H., Kawaguchi, K., Chechik, G., Fetaya, E.: Multi-task learning as a bargaining game. arXiv preprint arXiv:2202.01017 (2022)
  * [24] Ramakrishnan, R., Dral, P.O., Rupp, M., Von Lilienfeld, O.A.: Quantum chemistry structures and properties of 134 kilo molecules. Scientific data 1(1), 1–7 (2014)
  * [25] She, Y.: Sparse regression with exact clustering. Stanford University (2008)
  * [26] Tong, A., Wenkel, F., Bhaskar, D., MacDonald, K., Grady, J.D., Perlmutter, M., Krishnaswamy, S., Wolf, G.: Learnable filters for geometric scattering modules. ArXiv abs/2208.07458 (2022)
  * [27] Van Der Maaten, L., Postma, E., Van den Herik, J.: Dimensionality reduction: a comparative review. J Mach Learn Res 10, 66–71 (2009)
  * [28] Vijayakumar, S., D’souza, A., Schaal, S.: Incremental online learning in high dimensions. Neural computation 17(12), 2602–2634 (2005)
  * [29] Ye, H., Xu, D.: Taskprompter: Spatial-channel multi-task prompting for dense scene understanding. In: ICLR (2022)
  * [30] Zhang, Y., Zhang, Y., Wang, W.: Learning linear and nonlinear low-rank structure in multi-task learning. IEEE Trans. Knowl. Data Eng. 35, 8157–8170 (2023)
  * [31] Zhang, Y., Yang, Q.: A survey on multi-task learning. IEEE Transactions on Knowledge and Data Engineering 34(12), 5586–5609 (2021)
  * [32] Zhao, Y., Wang, X., Che, T., Bao, G., Li, S.: Multi-task deep learning for medical image computing and analysis: A review. Computers in Biology and Medicine 153, 106496 (2023)
  * [33] Zhong, Y., Xu, W., Gao, X.: Heterogeneous multi-task feature learning with mixed $l_{2,1}$ regularization. Mach. Learn. 113, 891–932 (2023)
  * [34] Zhou, M., Yang, P.: Automatic temporal relation in multi-task learning. Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining (2023)

##  Appendix 0.A Analysis of Variance and proof of Theorem 3.1

In this Section we firstly report the bias-variance decomposition equation,
that will be exploited in this section to prove the variance expression
reported in the main paper, and the bias in the next section.

Considering a task $\mathcal{T}_{i}$, the Mean Squared Error (MSE), can be
decomposed into three terms (bias-variance decomposition [13]):

$\begin{gathered}\underbrace{\mathbb{E}_{\mathbf{x},y_{i},\mathcal{S}_{i}}[(\mathcal{M}_{\mathcal{S}_{i}}(\mathbf{x})-y_{i})^{2}]}_{\text{MSE}}=\underbrace{\mathbb{E}_{\mathbf{x},\mathcal{S}_{i}}[(\mathcal{M}_{\mathcal{S}_{i}}(\mathbf{x})-\bar{\mathcal{M}}_{i}(\mathbf{x}))^{2}]}_{\text{variance}}\\\
\qquad+\underbrace{\mathbb{E}_{\mathbf{x}}[(\bar{\mathcal{M}_{i}}(\mathbf{x})-\bar{y}_{i}(\mathbf{x}))^{2}]}_{\text{bias}}+\underbrace{\mathbb{E}_{\mathbf{x},y_{i}}[(\bar{y}_{i}(\mathbf{x})-y_{i})^{2}]}_{\text{noise}},\end{gathered}$
(8)

where $\mathbf{x},y_{i}$ are respectively features and the target of a test
sample associated with task $\mathcal{T}_{i}$ and $\mathcal{S}_{i}$ is the
related training set. $\mathcal{M}_{\mathcal{S}_{i}}(\cdot)$ is the model
trained with $\mathcal{S}_{i}$, and $\bar{\mathcal{M}}_{i}(\cdot)$,
$\bar{y}_{i}$, are its expected value w.r.t. $\mathcal{S}_{i}$ and the
expected value of the test target $y_{i}$ w.r.t. the input features
$\mathbf{x}$.

In this section, recalling the data generation process underlying the
relationship between features and target of a generic task $i$ as:
$y_{i}=f_{i}(x_{1}\dots x_{D})+\epsilon_{i},\
\epsilon_{i}\sim\mathcal{N}(0,\sigma_{i}^{2})$, and the general linear
regression model under analysis:
$\hat{y}_{i}=\hat{\psi}_{\iota}=\widehat{\frac{1}{K_{\iota}}\sum_{k:y_{k}\in\mathcal{P}_{\iota}}y_{k}}=\hat{w}_{1}^{\iota}\phi_{1}+\dots+\hat{w}_{d}^{\iota}\phi_{d}$,
we will firstly focus on a simplified version and subsequently extend the
analysis to this more general case.

###  0.A.1 Linear model of the original features, single-task

Firstly we estimate the output $y_{i}$ with a multivariate linear regression
on the $D$ original features. Each model is predicted as:

$\hat{y}_{i}=\hat{w}_{1}^{i}x_{1}+\dots+\hat{w}_{D}^{i}x_{D}.$ (9)

To compute the variance of the linear model, we need the variance and the
expected value of the coefficients, which are reported in the following Lemma.

###### Lemma 1

In the asymptotic case, the variance and the expected value (conditioned to
the training set feature matrix $\mathbf{X}$) of the linear regression
coefficients of the model in Equation 9, trained to estimate the $i$-th task,
are respectively:

$\displaystyle\begin{split}var_{\mathcal{T}}&(\mathbf{\hat{w}^{i}}\lvert\mathbf{X})=\frac{\sigma^{2}_{i}}{(n-1)}\mathcal{P_{D}}\\\
&=\frac{\sigma^{2}_{i}}{(n-1)}\begin{bmatrix}\frac{1}{\sigma^{2}_{x_{1}|x^{-1}}}&\frac{-\rho_{x_{1},x_{2}|x^{-1,2}}}{\sigma_{x_{1}|x^{-1}}\cdot\sigma_{x_{2}|x^{-2}}}&\dots&\frac{-\rho_{x_{1},x_{D}|x^{-1,D}}}{\sigma_{x_{1}|x^{-1}}\cdot\sigma_{x_{D}|x^{-D}}}\\\
\dots&\dots&\dots&\dots\\\
\frac{-\rho_{x_{1},x_{D}|x^{-1,D}}}{\sigma_{x_{1}|x^{-1}}\cdot\sigma_{x_{D}|x^{-D}}}&\frac{-\rho_{x_{2},x_{D}|x^{-2,D}}}{\sigma_{x_{2}|x^{-2}}\cdot\sigma_{x_{D}|x^{-D}}}&\dots&\frac{1}{\sigma^{2}_{x_{D}|x^{-D}}}\end{bmatrix},\end{split}$
(10)
$\displaystyle\begin{split}\operatorname{\mathbb{E}}_{\mathcal{T}}[\mathbf{\hat{w}^{i}}\lvert\mathbf{X}]&=\frac{1}{(n-1)}\mathcal{P}_{D}\mathbf{X}^{\intercal}\mathbf{y}^{i}\\\
&=\begin{bmatrix}\frac{cov(x_{1},f_{i}(x))}{\sigma^{2}_{x_{1}\lvert
x^{-1}}}-\sum_{\begin{subarray}{c}k=1\\\ k\neq
1\end{subarray}}^{D}\Big{(}\frac{\rho_{x_{1},x_{k}\lvert
x^{-1,k}}}{\sigma_{x_{1}\lvert x^{-1}}\cdot\sigma_{x_{k}\lvert
x^{-k}}}\Big{)}cov(x_{k},f_{i}(x))\\\ \dots\\\
\frac{cov(x_{D},f_{i}(x))}{\sigma^{2}_{x_{D}\lvert
x^{-D}}}-\sum_{\begin{subarray}{c}k=1\\\ k\neq
D\end{subarray}}^{D}\Big{(}\frac{\rho_{x_{D},x_{k}\lvert
x^{-D,k}}}{\sigma_{x_{D}\lvert x^{-D}}\cdot\sigma_{x_{k}\lvert
x^{-k}}}\Big{)}cov(x_{k},f_{i}(x))\end{bmatrix},\end{split}$ (11)

where $\mathcal{P_{D}}$ is the inverse of the covariance matrix, i.e., the
precision matrix, $\sigma_{x_{j}|x^{-j}}$ is the partial variance of a feature
$x_{j}$ given the others, and $\rho_{x_{j},x_{k}|x^{-j,k}}$ is the partial
correlation between the features $x_{j},x_{k}$, given the others.

###### Proof

Recalling that in general, for a linear regression model [16]:

$var_{\mathcal{T}}(\hat{w}\lvert\mathbf{X})=(\mathbf{X}^{T}\mathbf{X})^{-1}\sigma^{2},$

we compute the variance of the estimated coefficients conditioned on the
training set:

$\displaystyle var_{\mathcal{T}}(\mathbf{\hat{w}^{i}}\lvert\mathbf{X})$
$\displaystyle=\Bigg{(}\begin{bmatrix}x^{1}_{1}&\dots&x^{n}_{1}\\\ &\dots&\\\
x_{D}^{1}&\dots&x_{D}^{n}\end{bmatrix}\begin{bmatrix}x^{1}_{1}&\dots&x^{1}_{D}\\\
&\dots&\\\ x_{1}^{n}&\dots&x_{D}^{n}\end{bmatrix}\Bigg{)}^{-1}\sigma^{2}_{i}$
$\displaystyle=\Bigg{(}\begin{bmatrix}\sum_{i=1}^{n}(x_{1}^{i})^{2}&\sum_{i=1}^{n}(x_{1}^{i}x_{2}^{i})&\dots&\sum_{i=1}^{n}(x_{1}^{i}x_{D}^{i})\\\
\sum_{i=1}^{n}(x_{1}^{i}x_{2}^{i})&\sum_{i=1}^{n}(x_{2}^{i})^{2}&\dots&\sum_{i=1}^{n}(x_{2}^{i}x_{D}^{i})\\\
\dots&\dots&\dots&\dots\\\
\sum_{i=1}^{n}(x_{1}^{i}x_{D}^{i})&\sum_{i=1}^{n}(x_{2}^{i}x_{D}^{i})&\dots&\sum_{i=1}^{n}(x_{D}^{i})^{2}\end{bmatrix}\Bigg{)}^{-1}\sigma^{2}_{i}$
$\displaystyle=\frac{\sigma^{2}_{i}}{(n-1)}\Bigg{(}\begin{bmatrix}\hat{\sigma}^{2}_{x_{1}}&\hat{cov}(x_{1},x_{2})&\dots&\hat{cov}(x_{1},x_{D})\\\
\dots&\dots&\dots&\dots\\\
\hat{cov}(x_{1},x_{D})&\hat{cov}(x_{2},x_{D})&\dots&\hat{\sigma}^{2}_{x_{D}}\end{bmatrix}\Bigg{)}^{-1}$
$\displaystyle=\frac{\sigma^{2}}{(n-1)}\mathcal{P_{D}}.$

The last equality holds since the inverse of the covariance matrix of the $D$
original features is the precision matrix $\mathcal{P_{D}}$ (or concentration
matrix), which can be rewritten as follows (directly considering the
asymptotic case to substitute the estimators with the quantities they
estimate):

$\mathcal{P_{D}}=\begin{bmatrix}\frac{1}{\sigma^{2}_{x_{1}|x^{-1}}}&\frac{-\rho_{x_{1},x_{2}|x^{-1,2}}}{\sigma_{x_{1}|x^{-1}}\cdot\sigma_{x_{2}|x^{-2}}}&\dots&\frac{-\rho_{x_{1},x_{D}|x^{-1,D}}}{\sigma_{x_{1}|x^{-1}}\cdot\sigma_{x_{D}|x^{-D}}}\\\
\dots&\dots&\dots&\dots\\\
\frac{-\rho_{x_{1},x_{D}|x^{-1,D}}}{\sigma_{x_{1}|x^{-1}}\cdot\sigma_{x_{D}|x^{-D}}}&\frac{-\rho_{x_{2},x_{D}|x^{-2,D}}}{\sigma_{x_{2}|x^{-2}}\cdot\sigma_{x_{D}|x^{-D}}}&\dots&\frac{1}{\sigma^{2}_{x_{D}|x^{-D}}}\end{bmatrix}.$

In the equation, the partial variance is the variance of the residual of the
regression of the other features on the current one, while the partial
correlation is the correlation of the residual of the linear regression of the
other features on them.

This proves the expression of variance of the lemma. Similarly, with the
expected values, we have:

$\displaystyle\operatorname{\mathbb{E}}_{\mathcal{T}}[\mathbf{\hat{w}^{i}}\lvert\mathbf{X}]$
$\displaystyle=(\mathbf{X}^{\intercal}\mathbf{X})^{-1}\mathbf{X}^{\intercal}\mathbf{y_{i}}=\frac{1}{(n-1)}\mathcal{P}_{D}\mathbf{X}^{\intercal}\mathbf{y_{i}}$
$\displaystyle=\frac{1}{(n-1)}\begin{bmatrix}\sum_{j=1}^{n}\Big{[}\frac{x_{1}^{j}}{\sigma^{2}_{x_{1}\lvert
x^{-1}}}-\sum_{\begin{subarray}{c}k=1\\\ k\neq
1\end{subarray}}^{D}\frac{\rho_{x_{1},x_{k}\lvert x^{-1,k}}\cdot
x_{k}^{j}}{\sigma_{x_{1}\lvert x^{-1}}\cdot\sigma_{x_{k}\lvert
x^{-k}}}\Big{]}\cdot f_{i}(x)^{j}\\\ \dots\\\
\sum_{j=1}^{n}\Big{[}\frac{x_{D}^{j}}{\sigma^{2}_{x_{D}\lvert
x^{-D}}}-\sum_{\begin{subarray}{c}k=1\\\ k\neq
D\end{subarray}}^{D}\frac{\rho_{x_{k},x_{D}\lvert x^{-k,D}}\cdot
x_{k}^{j}}{\sigma_{x_{k}\lvert x^{-k}}\cdot\sigma_{x_{D}\lvert
x^{-D}}}\Big{]}\cdot f_{i}(x)^{j}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\frac{cov(x_{1},y)}{\sigma^{2}_{x_{1}\lvert
x^{-1}}}-\sum_{\begin{subarray}{c}k=1\\\ k\neq
1\end{subarray}}^{D}\Big{(}\frac{\rho_{x_{1},x_{k}\lvert
x^{-1,k}}}{\sigma_{x_{1}\lvert x^{-1}}\cdot\sigma_{x_{k}\lvert
x^{-k}}}\Big{)}cov(x_{k},f_{i}(x))\\\ \dots\\\
\frac{cov(x_{D},y)}{\sigma^{2}_{x_{D}\lvert
x^{-D}}}-\sum_{\begin{subarray}{c}k=1\\\ k\neq
D\end{subarray}}^{D}\Big{(}\frac{\rho_{x_{D},x_{k}\lvert
x^{-D,k}}}{\sigma_{x_{D}\lvert x^{-D}}\cdot\sigma_{x_{k}\lvert
x^{-k}}}\Big{)}cov(x_{k},f_{i}(x))\end{bmatrix}$

This concludes the proofs, showing the second equality reported in the lemma.

###### Remark 9

In the case of $D$-dimensional linear data generation process, the expected
value of the coefficients is the vector of the coefficients of the data
generation process itself, and the model is unbiased, as expected.

We are now ready to show the variance of the model in this single-task
setting.

###### Theorem 0.A.1

In the asymptotic case, the variance of the linear regression model of
Equation 9, trained to estimate the $i$-th task, is related to the
characteristics of the features and with the coefficients of the linear
regression model as:

$\displaystyle var_{D}^{i}$
$\displaystyle=\frac{\sigma^{2}_{i}}{(n-1)}\cdot\sum_{k=1}^{D}\Big{\\{}\frac{\sigma^{2}_{x_{k}}}{\sigma^{2}_{x_{k}\lvert
x^{-k}}}-\sum_{j=1,j\neq
k}^{D}\frac{cov(x_{k},x_{j})\cdot\rho_{x_{k},x_{j}\lvert
x^{-k,j}}}{\sigma_{x_{k}\lvert x^{-k}}\cdot\sigma_{x_{j}\lvert
x^{-j}}}\Big{\\}}$
$\displaystyle=\frac{\sigma^{2}_{i}}{(n-1)}\cdot\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}^{k}_{k}\lvert\mathbf{X}],$

where
$\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}^{k}_{k}\lvert\mathbf{X}]$ is
the expected value of the coefficient associated with the variable $x_{k}$ of
the regression of the full set of features (including the variable $x_{k}$
itself) on $x_{k}$, which is equal to $1$ (and the other coefficients equal to
$0$) when the features are independent.

Furthermore, it is equal to:

$\displaystyle var_{D}^{i}=$ $\displaystyle\frac{\sigma^{2}_{i}}{(n-1)}\cdot
D.$

###### Proof

By definition of variance of the model, in this setting it is defined as
$var_{D}^{i}:=\operatorname{\mathbb{E}}_{x}\operatorname{\mathbb{E}}_{\mathcal{T}}[(\hat{w}_{1}^{i}x_{1}+\dots+\hat{w}_{D}^{i}x_{D}-\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}^{i}_{1}x_{1}+\dots+\hat{w}^{i}_{D}x_{D}])^{2}]$.
Recalling the expression of variance of the coefficients of the previous
lemma, in the asymptotic case we get:

$\displaystyle var_{D}^{i}=$
$\displaystyle\operatorname{\mathbb{E}}_{x}\operatorname{\mathbb{E}}_{\mathcal{T}}[(\hat{w}_{1}^{i}x_{1}+\dots+\hat{w}_{D}^{i}x_{D}-\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}^{i}_{1}x_{1}+\dots+\hat{w}^{i}_{D}x_{D}])^{2}]$
$\displaystyle=$ $\displaystyle\sum_{k=1}^{D}\sigma^{2}_{x_{k}}\cdot
var_{\mathcal{T}}(\hat{w}^{i}_{k})+2\sum_{k=1}^{D-1}\sum_{j=k+1}^{D}cov_{x}(x_{k},x_{j})\cdot
cov_{\mathcal{T}}(\hat{w}_{k}^{i},\hat{w}_{j}^{i})$ $\displaystyle=$
$\displaystyle\frac{\sigma^{2}_{i}}{(n-1)}\cdot\Big{\\{}\sum_{k=1}^{D}\frac{\sigma^{2}_{x_{k}}}{\sigma^{2}_{x_{k}\lvert
x^{-k}}}-2\sum_{k=1}^{D-1}\sum_{j=k+1}^{D}\frac{cov(x_{k},x_{j})\cdot\rho_{x_{k},x_{j}\lvert
x^{-k,j}}}{\sigma_{x_{k}\lvert x^{-k}}\cdot\sigma_{x_{j}\lvert
x^{-j}}}\Big{\\}}.$

Firstly, exploiting the expression of expected values of coefficients adapted
from the previous lemma, we prove the first expression of the theorem:

$\displaystyle var_{D}^{i}=$
$\displaystyle\frac{\sigma^{2}_{i}}{(n-1)}\cdot\Big{\\{}\sum_{k=1}^{D}\frac{\sigma^{2}_{x_{k}}}{\sigma^{2}_{x_{k}\lvert
x^{-k}}}-2\sum_{k=1}^{D-1}\sum_{j=k+1}^{D}\frac{cov(x_{k},x_{j})\cdot\rho_{x_{k},x_{j}\lvert
x^{-k,j}}}{\sigma_{x_{k}\lvert x^{-k}}\cdot\sigma_{x_{j}\lvert
x^{-j}}}\Big{\\}}$ $\displaystyle=$
$\displaystyle\frac{\sigma^{2}_{i}}{(n-1)}\cdot\sum_{k=1}^{D}\Big{\\{}\frac{\sigma^{2}_{x_{k}}}{\sigma^{2}_{x_{k}\lvert
x^{-k}}}-\sum_{j=1,j\neq
k}^{D}\frac{cov(x_{k},x_{j})\cdot\rho_{x_{k},x_{j}\lvert
x^{-k,j}}}{\sigma_{x_{k}\lvert x^{-k}}\cdot\sigma_{x_{j}\lvert
x^{-j}}}\Big{\\}}$
$\displaystyle=\frac{\sigma^{2}_{i}}{(n-1)}\cdot\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}^{k}_{k}\lvert\mathbf{X}].$

For the second result, recalling the expression of the variance, The following
considerations hold:

$\displaystyle var_{D}^{i}=$
$\displaystyle\sum_{k=1}^{D}\sum_{j=1}^{D}cov_{x}(x_{k},x_{j})\cdot
cov_{\mathcal{T}}(\hat{w}_{k}^{i},\hat{w}_{j}^{i}).$

This is the sum of $D\times D$ elements, basically the sum of the elements of
the Hadamard product between the elements of the covariance matrix of the
coefficient and the covariance matrix of the features. Recalling that
$var_{\mathcal{T}}(\hat{w}^{i}_{D}\lvert\mathbf{X})=\frac{\sigma^{2}_{i}}{(n-1)}(\mathbf{X}^{\intercal}\mathbf{X})^{-1}$
and $cov(\mathbf{x})=\mathbf{X}^{\intercal}\mathbf{X}$, the following
equivalences hold, proving the result:

$\displaystyle
var_{D}^{i}=\sum_{k=1}^{D}\sum_{j=1}^{D}cov_{x}(x_{k},x_{j})\cdot
cov_{\mathcal{T}}(\hat{w}_{k}^{i},\hat{w}_{j}^{i})=$
$\displaystyle=\frac{\sigma^{2}_{i}}{(n-1)}\begin{bmatrix}1&\dots&1\end{bmatrix}(\mathbf{X}^{\intercal}\mathbf{X})^{-1}(\mathbf{X}^{\intercal}\mathbf{X})\begin{bmatrix}1\\\
\dots\\\ 1\end{bmatrix}=\frac{\sigma^{2}_{i}}{(n-1)}D.$

###### Remark 10

The result can be trivially confirmed in the settings that consider two and
three input features, where the inverse
$(\mathbf{X}^{\intercal}\mathbf{X})^{-1}$ can be computed explicitly,
exploiting the closed form of the inverses of $2\times 2,3\times 3$ matrices.

###  0.A.2 Extension to transformed features and aggregated targets

In this section, we elaborate on the results of the previous subsection.

Firstly, given the original targets $y_{1}\dots y_{L}$, we apply a
multivariate linear regression after having aggregated clusters of targets,
resulting in $l$ aggregated targets $\psi_{1}\dots\psi_{l}$. We still consider
the $D$ original features. Focusing again on the $i$-th task with target
$y_{i}$ and assuming that it is associated with the $\iota$ cluster
$\psi_{\iota}=\frac{1}{K}\sum_{k:y_{k}\in\mathcal{P}_{\iota}}y_{k}$, we
estimate the target $y_{i}$ with:

$\hat{y}_{i}=\hat{\psi}_{\iota}=\widehat{\frac{1}{K_{\iota}}\sum_{k:y_{k}\in\mathcal{P}_{\iota}}y_{k}}=\hat{w}_{1}^{\iota}x_{1}+\dots+\hat{w}_{D}^{\iota}x_{D}.$

Defining
$\bar{\sigma}_{\iota}^{2}:=var(\frac{1}{K_{\iota}}\sum_{k:y_{k}\in\mathcal{P}_{\iota}}\epsilon_{k})$,
recalling that we have the same features of the previous subsection, we can
write the variance of the coefficients as:

$\displaystyle var_{\mathcal{T}}(\hat{w}^{\iota}_{D}\lvert\mathbf{X})$
$\displaystyle=\frac{\bar{\sigma}_{\iota}^{2}}{(n-1)}\mathcal{P_{D}}$
$\displaystyle=\frac{\bar{\sigma}_{\iota}^{2}}{(n-1)}\begin{bmatrix}\frac{1}{\sigma^{2}_{x_{1}|x^{-1}}}&\frac{-\rho_{x_{1},x_{2}|x^{-1,2}}}{\sigma_{x_{1}|x^{-1}}\cdot\sigma_{x_{2}|x^{-2}}}&\dots&\frac{-\rho_{x_{1},x_{D}|x^{-1,D}}}{\sigma_{x_{1}|x^{-1}}\cdot\sigma_{x_{D}|x^{-D}}}\\\
\dots&\dots&\dots&\dots\\\
\frac{-\rho_{x_{1},x_{D}|x^{-1,D}}}{\sigma_{x_{1}|x^{-1}}\cdot\sigma_{x_{D}|x^{-D}}}&\frac{-\rho_{x_{2},x_{D}|x^{-2,D}}}{\sigma_{x_{2}|x^{-2}}\cdot\sigma_{x_{D}|x^{-D}}}&\dots&\frac{1}{\sigma^{2}_{x_{D}|x^{-D}}}\end{bmatrix}.$

Then, we estimate each task with a single-task multivariate linear regression
on the $d$ reduced features $\phi_{1}\dots\phi_{d}$. The $i$-th model is
therefore predicted as:

$\hat{y}_{i}=\hat{w}_{1}^{i}\phi_{1}+\dots+\hat{w}_{d}^{i}\phi_{d}.$

The variance of the coefficients, therefore, becomes:

$\displaystyle var_{\mathcal{T}}(\hat{w}^{i}_{d}\lvert\mathbf{X})$
$\displaystyle=\frac{\sigma^{2}_{i}}{(n-1)}\mathcal{P}_{d}(\phi)$
$\displaystyle=\frac{\sigma^{2}_{i}}{(n-1)}\begin{bmatrix}\frac{1}{\sigma^{2}_{\phi_{1}|\phi^{-1}}}&\frac{-\rho_{\phi_{1},\phi_{2}|\phi^{-1,2}}}{\sigma_{\phi_{1}|\phi^{-1}}\cdot\sigma_{\phi_{2}|\phi^{-2}}}&\dots&\frac{-\rho_{\phi_{1},\phi_{D}|\phi^{-1,d}}}{\sigma_{\phi_{1}|\phi^{-1}}\cdot\sigma_{\phi_{d}|\phi^{-d}}}\\\
\dots&\dots&\dots&\dots\\\
\frac{-\rho_{\phi_{1},\phi_{d}|\phi^{-1,d}}}{\sigma_{\phi_{1}|\phi^{-1}}\cdot\sigma_{\phi_{d}|\phi^{-d}}}&\frac{-\rho_{\phi_{2},\phi_{d}|\phi^{-2,d}}}{\sigma_{\phi_{2}|\phi^{-2}}\cdot\sigma_{\phi_{d}|\phi^{-d}}}&\dots&\frac{1}{\sigma^{2}_{\phi_{d}|\phi^{-d}}}\end{bmatrix}.$

Finally, combining both aggregations of targets and features, we get the
following variance of the coefficients:

$\displaystyle var_{\mathcal{T}}(\hat{w}^{\iota}_{d}\lvert\mathbf{X})$
$\displaystyle=\frac{\bar{\sigma}_{\iota}^{2}}{(n-1)}\mathcal{P}_{d}(\phi)$
$\displaystyle=\frac{\bar{\sigma}_{\iota}^{2}}{(n-1)}\begin{bmatrix}\frac{1}{\sigma^{2}_{\phi_{1}|\phi^{-1}}}&\frac{-\rho_{\phi_{1},\phi_{2}|\phi^{-1,2}}}{\sigma_{\phi_{1}|\phi^{-1}}\cdot\sigma_{\phi_{2}|\phi^{-2}}}&\dots&\frac{-\rho_{\phi_{1},\phi_{D}|\phi^{-1,d}}}{\sigma_{\phi_{1}|\phi^{-1}}\cdot\sigma_{\phi_{d}|\phi^{-d}}}\\\
\dots&\dots&\dots&\dots\\\
\frac{-\rho_{\phi_{1},\phi_{d}|\phi^{-1,d}}}{\sigma_{\phi_{1}|\phi^{-1}}\cdot\sigma_{\phi_{d}|\phi^{-d}}}&\frac{-\rho_{\phi_{2},\phi_{d}|\phi^{-2,d}}}{\sigma_{\phi_{2}|\phi^{-2}}\cdot\sigma_{\phi_{d}|\phi^{-d}}}&\dots&\frac{1}{\sigma^{2}_{\phi_{d}|\phi^{-d}}}\end{bmatrix}.$

Therefore, we can exploit these results to obtain the expression of the
variance in the three configurations. Note that the last configuration is the
most general. It contains all the previous ones as particular cases, and it
proves Theorem 3.1 of the main paper.

###### Theorem 0.A.2

In the asymptotic case, estimating the $i$-th task with the average
$\psi_{\iota}$ of a set of tasks that contains $y_{i}$, the variance of the
model is:

$\displaystyle var_{D}^{\iota}=$
$\displaystyle\frac{\bar{\sigma}_{\iota}^{2}}{(n-1)}\cdot D.$

On the other hand, estimating the $i$-th task in a single-task setting, with a
set of transformed features $\phi_{1}\dots\phi_{d}$, the variance of the model
is:

$var_{d}^{i}=\frac{\sigma^{2}_{i}}{(n-1)}\cdot d.$

Finally, combining the aggregation of targets and the transformation of
features, we get the expression reported in Theorem 3.1 in the main paper:

$var_{d}^{\iota}=\frac{\bar{\sigma}_{\iota}^{2}}{(n-1)}\cdot d.$

###### Proof

Exploiting the expression of variance of the coefficients, with similar
arguments of the previous subsection, in the first case, we get:

$\displaystyle
var_{D}^{\iota}=\sum_{k=1}^{D}\sum_{j=1}^{D}cov_{x}(x_{k},x_{j})\cdot
cov_{\mathcal{T}}(\hat{w}_{k}^{\iota},\hat{w}_{j}^{\iota})=$
$\displaystyle=\frac{\sigma^{2}_{\iota}}{(n-1)}\begin{bmatrix}1&\dots&1\end{bmatrix}(\mathbf{X}^{\intercal}\mathbf{X})^{-1}(\mathbf{X}^{\intercal}\mathbf{X})\begin{bmatrix}1\\\
\dots\\\ 1\end{bmatrix}=\frac{\sigma^{2}_{i}}{(n-1)}D.$

In the second case we get:

$\displaystyle
var_{d}^{i}=\sum_{k=1}^{d}\sum_{j=1}^{d}cov_{x}(\phi_{k},\phi_{j})\cdot
cov_{\mathcal{T}}(\hat{w}_{k}^{i},\hat{w}_{j}^{i})=$
$\displaystyle=\frac{\sigma^{2}_{i}}{(n-1)}\begin{bmatrix}1&\dots&1\end{bmatrix}(\mathbf{\Phi}^{\intercal}\mathbf{\Phi})^{-1}(\mathbf{\Phi}^{\intercal}\mathbf{\Phi})\begin{bmatrix}1\\\
\dots\\\ 1\end{bmatrix}=\frac{\sigma^{2}_{i}}{(n-1)}d.$

Finally, by combining the two cases, we get:

$\displaystyle
var_{d}^{\iota}=\sum_{k=1}^{d}\sum_{j=1}^{d}cov_{x}(\phi_{k},\phi_{j})\cdot
cov_{\mathcal{T}}(\hat{w}_{k}^{\iota},\hat{w}_{j}^{\iota})=$
$\displaystyle=\frac{\sigma^{2}_{\iota}}{(n-1)}\begin{bmatrix}1&\dots&1\end{bmatrix}(\mathbf{\Phi}^{\intercal}\mathbf{\Phi})^{-1}(\mathbf{\Phi}^{\intercal}\mathbf{\Phi})\begin{bmatrix}1\\\
\dots\\\ 1\end{bmatrix}=\frac{\sigma^{2}_{\iota}}{(n-1)}d.$

##  Appendix 0.B Analysis of Bias and proof of Theorem 3.2

As in the previous section, to show the bias of the model, we will first focus
on a single-task linear regression of the original features, subsequently
extending to the general case of an aggregated target $\phi_{\iota}$ and a
reduced set of inputs $\phi_{1},\dots,\phi_{d}$.

###  0.B.1 Linear model of the original features, single-task

Firstly we estimate the output $y_{i}$ with a multivariate linear regression
on the $D$ original features. Each model is predicted as shown in Equation 9.
To compute the bias of this linear model, we will exploit the expression of
the expected value of the coefficients that can be found in Lemma 1.
Additionally, we will exploit the following Lemma, that shows the expression
of the coefficient o multiple correlation in a linear regression setting.

###### Lemma 2

Considering the linear regression on the target $y_{i}$ of the features
$x_{1},\dots,x_{D}$, the coefficient of multiple correlation $R^{2}_{D,i}$ is:

$\displaystyle R^{2}_{D,i}$
$\displaystyle=\sum_{k=1}^{D}\Big{\\{}\frac{\sigma^{2}_{x_{k}}}{\sigma^{2}_{x_{k}\lvert
x^{-k}}}\rho^{2}_{x_{k},y_{i}}-\rho_{x_{k},y_{i}}\sum_{\begin{subarray}{c}h=1\\\
h\neq k\end{subarray}}^{D}\Big{(}\frac{\rho_{x_{h},x_{k}\lvert
x^{-h,k}}}{\sigma_{x_{h}\lvert x^{-h}}\cdot\sigma_{x_{k}\lvert
x^{-k}}}\Big{)}\rho_{x_{h},y_{i}}\sigma_{x_{h}}\sigma_{x_{k}}\Big{\\}}$
$\displaystyle=\frac{1}{\sigma^{2}_{y_{i}}}\sum_{k=1}^{D}\Big{\\{}\frac{cov(x_{k},y_{i})^{2}}{\sigma^{2}_{x_{k}\lvert
x^{-k}}}-cov(x_{k},y_{i})\sum_{\begin{subarray}{c}h=1\\\ h\neq
k\end{subarray}}^{D}\Big{(}\frac{\rho_{x_{h},x_{k}\lvert
x^{-h,k}}}{\sigma_{x_{h}\lvert x^{-h}}\cdot\sigma_{x_{k}\lvert
x^{-k}}}\Big{)}cov(x_{h},y_{i})\Big{\\}}.$

###### Proof

Considering the definition of coefficient of determination (i.e., squared
coefficient of multiple correlation):

$\displaystyle R^{2}_{D,i}$
$\displaystyle=\begin{bmatrix}\rho_{x_{1},y_{i}}&\dots&\rho_{x_{D},y_{i}}\end{bmatrix}corr(\mathbf{X})^{-1}\begin{bmatrix}\rho_{x_{1},y_{i}}\\\
\dots\\\ \rho_{x_{D},y_{i}}\end{bmatrix},$

we need to compute the inverse of the correlation matrix of the features.

Given the inverse of the covariance matrix:

$cov(\mathbf{X})^{-1}=\begin{bmatrix}\frac{1}{\sigma^{2}_{x_{1}}|x^{-1}}&\frac{\rho_{x_{1},x_{2}\lvert
x^{-1,2}}}{\sigma_{x_{1}\lvert x^{-1}}\cdot\sigma_{x_{2}\lvert
x^{-2}}}&\dots&\frac{\rho_{x_{1},x_{D}\lvert x^{-1,D}}}{\sigma_{x_{1}\lvert
x^{-1}}\cdot\sigma_{x_{D}\lvert x^{-D}}}\\\ \dots&\dots&\dots&\dots\\\
\frac{\rho_{x_{1},x_{D}\lvert x^{-1,D}}}{\sigma_{x_{1}\lvert
x^{-1}}\cdot\sigma_{x_{D}\lvert x^{-D}}}&\frac{\rho_{x_{D},x_{2}\lvert
x^{-D,2}}}{\sigma_{x_{D}\lvert x^{-D}}\cdot\sigma_{x_{2}\lvert
x^{-2}}}&\dots&\frac{1}{\sigma^{2}_{x_{D}}|x^{-D}}\\\ \end{bmatrix}$

decomposing it as follows:

$cov(\mathbf{x})=\begin{bmatrix}\sigma_{1}&0&\dots&0\\\
\dots&\dots&\dots&\dots\\\
0&\dots&0&\sigma_{D}\end{bmatrix}\begin{bmatrix}1&\rho_{x_{1},x_{2}}&\dots&\rho_{x_{1},x_{D}}\\\
\dots&\dots&\dots&\dots\\\
\rho_{x_{1},x_{D}}&\rho_{x_{2},x_{D}}&\dots&1\end{bmatrix}\begin{bmatrix}\sigma_{1}&0&\dots&0\\\
\dots&\dots&\dots&\dots\\\ 0&\dots&0&\sigma_{D}\end{bmatrix},$

and recalling that the inverse of a product is the product of the inverses
$(ABA)^{-1}=A^{-1}B^{-1}A^{-1}$, we get:

$\displaystyle corr(\mathbf{X})^{-1}=\begin{bmatrix}\sigma_{1}&0&\dots&0\\\
\dots&\dots&\dots&\dots\\\
0&\dots&0&\sigma_{D}\end{bmatrix}cov(\mathbf{X})^{-1}\begin{bmatrix}\sigma_{1}&0&\dots&0\\\
\dots&\dots&\dots&\dots\\\ 0&\dots&0&\sigma_{D}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\sigma_{1}&0&\dots&0\\\
\dots&\dots&\dots&\dots\\\
0&\dots&0&\sigma_{D}\end{bmatrix}\begin{bmatrix}\frac{1}{\sigma^{2}_{x_{1}}|x^{-1}}&\dots&\frac{\rho_{x_{1},x_{D}\lvert
x^{-1,D}}}{\sigma_{x_{1}\lvert x^{-1}}\cdot\sigma_{x_{D}\lvert x^{-D}}}\\\
\dots&\dots&\dots\\\ \frac{\rho_{x_{1},x_{D}\lvert
x^{-1,D}}}{\sigma_{x_{1}\lvert x^{-1}}\cdot\sigma_{x_{D}\lvert
x^{-D}}}&\dots&\frac{1}{\sigma^{2}_{x_{D}}|x^{-D}}\\\
\end{bmatrix}\begin{bmatrix}\sigma_{1}&0&\dots&0\\\ \dots&\dots&\dots&\dots\\\
0&\dots&0&\sigma_{D}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\frac{\sigma^{2}_{x_{1}}}{\sigma^{2}_{x_{1}}|x^{-1}}&\frac{\sigma_{x_{1}}\sigma_{x_{2}}\rho_{x_{1},x_{2}\lvert
x^{-1,2}}}{\sigma_{x_{1}\lvert x^{-1}}\cdot\sigma_{x_{2}\lvert
x^{-2}}}&\dots&\frac{\sigma_{x_{1}}\sigma_{x_{D}}\rho_{x_{1},x_{D}\lvert
x^{-1,D}}}{\sigma_{x_{1}\lvert x^{-1}}\cdot\sigma_{x_{D}\lvert x^{-D}}}\\\
\dots&\dots&\dots&\dots\\\
\frac{\sigma_{x_{1}}\sigma_{x_{D}}\rho_{x_{1},x_{D}\lvert
x^{-1,D}}}{\sigma_{x_{1}\lvert x^{-1}}\cdot\sigma_{x_{D}\lvert
x^{-D}}}&\frac{\sigma_{x_{D}}\sigma_{x_{2}}\rho_{x_{D},x_{2}\lvert
x^{-D,2}}}{\sigma_{x_{D}\lvert x^{-D}}\cdot\sigma_{x_{2}\lvert
x^{-2}}}&\dots&\frac{\sigma_{x_{D}}^{2}}{\sigma^{2}_{x_{D}}|x^{-D}}\\\
\end{bmatrix}.$

Therefore:

$\displaystyle
R^{2}_{D,i}=\begin{bmatrix}\rho_{x_{1},y_{i}}&\dots&\rho_{x_{D},y_{i}}\end{bmatrix}\begin{bmatrix}\frac{\sigma^{2}_{x_{1}}}{\sigma^{2}_{x_{1}}|x^{-1}}&\dots&\frac{\sigma_{x_{1}}\sigma_{x_{D}}\rho_{x_{1},x_{D}\lvert
x^{-1,D}}}{\sigma_{x_{1}\lvert x^{-1}}\cdot\sigma_{x_{D}\lvert x^{-D}}}\\\
\dots&\dots&\dots\\\ \frac{\sigma_{x_{1}}\sigma_{x_{D}}\rho_{x_{1},x_{D}\lvert
x^{-1,D}}}{\sigma_{x_{1}\lvert x^{-1}}\cdot\sigma_{x_{D}\lvert
x^{-D}}}&\dots&\frac{\sigma_{x_{D}}^{2}}{\sigma^{2}_{x_{D}}|x^{-D}}\\\
\end{bmatrix}\begin{bmatrix}\rho_{x_{1},y_{i}}\\\ \dots\\\
\rho_{x_{D},y_{i}}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}{\scriptstyle\frac{\sigma^{2}_{x_{1}}}{\sigma^{2}_{x_{1}\lvert
x^{-1}}}\rho_{x_{1},y_{i}}-\sum_{\begin{subarray}{c}h=1\\\ h\neq
1\end{subarray}}^{D}\Big{(}\frac{\rho_{x_{h},x_{1}\lvert
x^{-h,1}}}{\sigma_{x_{h}\lvert x^{-h}}\cdot\sigma_{x_{1}\lvert
x^{-1}}}\Big{)}\rho_{x_{h},y_{i}}\sigma_{x_{h}}\sigma_{x_{1}}}&\dots\end{bmatrix}\begin{bmatrix}\rho_{x_{1},y_{i}}\\\
\dots\\\ \rho_{x_{D},y_{i}}\end{bmatrix}$
$\displaystyle=\sum_{k=1}^{D}\Big{\\{}\frac{\sigma^{2}_{x_{k}}}{\sigma^{2}_{x_{k}\lvert
x^{-k}}}\rho^{2}_{x_{k},y_{i}}-\rho_{x_{k},y_{i}}\sum_{\begin{subarray}{c}h=1\\\
h\neq k\end{subarray}}^{D}\Big{(}\frac{\rho_{x_{h},x_{k}\lvert
x^{-h,k}}}{\sigma_{x_{h}\lvert x^{-h}}\cdot\sigma_{x_{k}\lvert
x^{-k}}}\Big{)}\rho_{x_{h},y_{i}}\sigma_{x_{h}}\sigma_{x_{k}}\Big{\\}}.$

###### Remark 11

The same result of the Lemma can be computed by considering the definition of
the coefficient of determination in terms of the residual and total sum of
squares:

$R^{2}_{D,i}=\frac{SS_{reg}}{SS_{tot}}=\frac{\mathbf{\hat{w}}^{\intercal}X^{\intercal}\mathbf{y}_{i}-\frac{1}{N}(\mathbf{1}^{\intercal}\mathbf{y}_{i})^{2}}{\mathbf{y}_{i}^{\intercal}\mathbf{y}_{i}-\frac{1}{N}(\mathbf{1}^{\intercal}\mathbf{y}_{i})^{2}}.$

The following theorem finally shows the expression of bias of the model in
this linear single-task setting.

###### Theorem 0.B.1

In the asymptotic case, the bias of the linear regression model of Equation 9,
trained to estimate the $i$-th task, is related to the characteristics of the
features and with the coefficients of the linear regression model as:

$\displaystyle bias_{D}^{i}$
$\displaystyle=var_{x}(\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]x_{k})+\sigma^{2}_{f_{i}}-2cov(\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]x_{k},f_{i}(x)),$

where the expected values of the regression coefficients
$\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]$ can be computed
with the results of Lemma 1.

Furthermore, it is equal to:

$\displaystyle bias_{D}^{i}=\sigma^{2}_{f_{i}}(1-R^{2}_{D,i}).$

###### Proof

The first expression can be obtained by definition of bias of the model and
exploiting the hypothesis of null expected values of $f_{i}$ and
$x_{1},\dots,x_{D}$, as:

$\displaystyle bias_{D}^{i}$
$\displaystyle=\operatorname{\mathbb{E}}_{x}[(\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]x_{k}-f_{i}(x))^{2}]=var_{x}(\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]x_{k}-f_{i}(x))$
$\displaystyle=var_{x}(\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]x_{k})+var_{x}(f_{i}(x))-2cov(\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]x_{k},f_{i}(x)).$

This already shows that the bias depends on the variance of the predicted and
the real model, decreased by how much the real and the predicted model are
related (through covariance). However, by exploiting the coefficient of
multiple correlation and the partial correlation, we can refine this
expression to estimate it directly from the data.

Exploiting the expression of the expected value of the coefficients, the last
term of the equation above can be written as:

$\displaystyle-2cov(\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]x_{k},f_{i}(x))=-2\sum_{k=1}^{D}cov(x_{k},f_{i})\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]$
$\displaystyle=-2\sum_{k=1}^{D}\Big{\\{}\frac{cov(x_{k},f_{i})^{2}}{\sigma^{2}_{x_{k}\lvert
x^{-k}}}-\sum_{\begin{subarray}{c}h=1\\\ h\neq
k\end{subarray}}^{D}\Big{(}\frac{\rho_{x_{h},x_{k}\lvert
x^{-h,k}}}{\sigma_{x_{h}\lvert x^{-h}}\cdot\sigma_{x_{k}\lvert
x^{-k}}}\Big{)}cov(x_{h},f_{i})cov(x_{k},f_{i})\Big{\\}}$
$\displaystyle=-2\sigma^{2}_{f_{i}}R^{2}_{D,i}.$

The bias is, therefore, equal to:

$\displaystyle
bias_{D}^{i}=var_{x}(\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]x_{k})+\sigma^{2}_{f_{i}}-2\sigma^{2}_{f_{i}}R^{2}_{D,i}.$

We finally characterize the first term, i.e. the variance of the linear
estimator, concluding the proof of the theorem.

###### Lemma 3

The following equality holds:

$\sigma^{2}_{f_{i}}R^{2}_{D,i}=var_{x}(\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]x_{k}).$

###### Proof

Exploiting the properties of variance:

$\displaystyle
var_{x}(\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]x_{k})$
$\displaystyle=\sum_{k=1}^{D}\sum_{j=1}^{D}cov(x_{k},x_{j})\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{j}^{i}]$
$\displaystyle=\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]\Big{(}\sum_{j=1}^{D}cov(x_{k},x_{j})\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{j}^{i}]\Big{)},$

and recalling that, from above:

$\displaystyle\sigma^{2}_{f_{i}}R^{2}_{D,i}=\sum_{k=1}^{D}cov(x_{k},f_{i})\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}],$

we need to show that, for each $k\in\\{1,\dots,D\\}$, the following equality
holds:

$\sum_{j=1}^{D}cov(x_{k},x_{j})\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{j}^{i}]=cov(x_{k},f_{i}).$

This equality can be shown by exploiting the definition of partial covariance:

$0=cov(x_{k},f_{i}|\mathbf{X})=cov(x_{k},f_{i})-cov(x_{k},\mathbf{X})cov(\mathbf{X})^{-1}cov(\mathbf{X},f_{i})$

where the first equality holds because the residual of the regression of
$\mathbf{X}$ on $x_{k}$ is zero since $x_{k}$ is among the regressors, and the
second equality holds by definition of partial covariance. Rewriting the
equality that we want to show as:

$\displaystyle
cov(x_{k},f_{i})=\sum_{j=1}^{D}cov(x_{k},x_{j})\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{j}^{i}]=cov(x_{k},\mathbf{X})\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{\mathbf{w}}^{i}]$
$\displaystyle=cov(x_{k},\mathbf{X})(\mathbf{X}^{\intercal}\mathbf{X})^{-1}\mathbf{X}^{\intercal}\mathbf{f_{i}}=cov(x_{k},\mathbf{X})cov(\mathbf{X})^{-1}cov(\mathbf{X},f_{i}),$

proves the statement of the lemma and concludes the proof of the theorem.

###### Remark 12

The result is confirmed in the two- and three-dimensional case and in the case
of $D$ mutually independent features, where all the quantities can be
explicitly computed.

###  0.B.2 Extension to transformed features and aggregated targets

As done for the variance, in this section, we elaborate on the results of the
previous subsection and extend them.

Firstly, given the original targets $y_{1}\dots y_{L}$, we apply a
multivariate linear regression after having aggregated clusters of targets,
resulting in $l$ aggregated targets $\psi_{1}\dots\psi_{l}$. We still consider
the $D$ original features. Focusing again on the $i$-th task, defining
$\bar{f}_{\iota}(x):=\frac{1}{K}\sum_{k:y_{k}\in\mathcal{P}_{\iota}}f_{k}(x)$,
the expected values of coefficients are:

$\displaystyle\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}^{\iota}_{D}\lvert\mathbf{X}]=\begin{bmatrix}\frac{cov(x_{1},\bar{f}_{\iota}(x))}{\sigma^{2}_{x_{1}\lvert
x^{-1}}}-\sum_{\begin{subarray}{c}k=1\\\ k\neq
1\end{subarray}}^{D}\Big{(}\frac{\rho_{x_{1},x_{k}\lvert
x^{-1,k}}}{\sigma_{x_{1}\lvert x^{-1}}\cdot\sigma_{x_{k}\lvert
x^{-k}}}\Big{)}cov(x_{k},\bar{f}_{\iota}(x))\\\ \dots\\\
\frac{cov(x_{D},\bar{f}_{\iota}(x))}{\sigma^{2}_{x_{D}\lvert
x^{-D}}}-\sum_{\begin{subarray}{c}k=1\\\ k\neq
D\end{subarray}}^{D}\Big{(}\frac{\rho_{x_{D},x_{k}\lvert
x^{-D,k}}}{\sigma_{x_{D}\lvert x^{-D}}\cdot\sigma_{x_{k}\lvert
x^{-k}}}\Big{)}cov(x_{k},\bar{f}_{\iota}(x))\end{bmatrix}.$

###### Theorem 0.B.2

In the asymptotic case, estimating the $i$-th task with the average
$\psi_{\iota}$ of a set of tasks that contains $y_{i}$, the bias of the model
is:

$\begin{split}bias_{d}^{\iota}&=\sigma^{2}_{f_{i}}+\sigma^{2}_{\psi_{\iota}}R^{2}_{d,\iota}-2(cov(f_{i},\psi_{\iota})-cov(f_{i},\psi_{\iota}|\mathbf{\Phi}))\\\
&=\sigma^{2}_{f_{i}}-\sigma^{2}_{\psi_{\iota}}R^{2}_{d,\iota}+2(cov(\psi_{\iota},f_{i}-\psi_{\iota}|\mathbf{\Phi})-cov(\psi_{\iota},f_{i}-\psi_{\iota})).\end{split}$
(12)

###### Proof

In this case, the bias w.r.t. the $i$-th task is:

$\displaystyle\operatorname{\mathbb{E}}_{x}$
$\displaystyle[(\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{\iota}]x_{k}-f_{i}(x))^{2}]=var_{x}(\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{\iota}]x_{k}-f_{i}(x))$
$\displaystyle=var_{x}(\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{\iota}]x_{k})+var_{x}(f_{i}(x))-2cov(\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{\iota}]x_{k},f_{i}(x)).$

Considering the linear regression on the target $\psi_{\iota}$ of the features
$x_{1},\dots,x_{D}$ the coefficient of multiple correlation $R^{2}_{D,\iota}$
is:

$\displaystyle R^{2}_{D,\iota}$
$\displaystyle=\sum_{k=1}^{D}\Big{\\{}\frac{\sigma^{2}_{x_{k}}}{\sigma^{2}_{x_{k}\lvert
x^{-k}}}\rho^{2}_{x_{k},\psi_{\iota}}-\rho_{x_{k},\psi_{\iota}}\sum_{\begin{subarray}{c}h=1\\\
h\neq k\end{subarray}}^{D}\Big{(}\frac{\rho_{x_{h},x_{k}\lvert
x^{-h,k}}}{\sigma_{x_{h}\lvert x^{-h}}\cdot\sigma_{x_{k}\lvert
x^{-k}}}\Big{)}\rho_{x_{h},\psi_{\iota}}\sigma_{x_{h}}\sigma_{x_{k}}\Big{\\}}$
$\displaystyle=\frac{1}{\sigma^{2}_{\psi_{\iota}}}\sum_{k=1}^{D}\Big{\\{}\frac{cov(x_{k},\psi_{\iota})^{2}}{\sigma^{2}_{x_{k}\lvert
x^{-k}}}-cov(x_{k},\psi_{\iota})\sum_{\begin{subarray}{c}h=1\\\ h\neq
k\end{subarray}}^{D}\Big{(}\frac{\rho_{x_{h},x_{k}\lvert
x^{-h,k}}}{\sigma_{x_{h}\lvert x^{-h}}\cdot\sigma_{x_{k}\lvert
x^{-k}}}\Big{)}cov(x_{h},\psi_{\iota})\Big{\\}}.$

Recalling the expression of the expected value of each coefficient:

$\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}^{\iota}_{k}\lvert\mathbf{X}]=\frac{cov(x_{k},\psi_{\iota})}{\sigma^{2}_{x_{k}\lvert
x^{-k}}}-\sum_{\begin{subarray}{c}h=1\\\ h\neq
k\end{subarray}}^{D}\Big{(}\frac{\rho_{x_{h},x_{k}\lvert
x^{-h,k}}}{\sigma_{x_{h}\lvert x^{-h}}\cdot\sigma_{x_{k}\lvert
x^{-k}}}\Big{)}cov(x_{h},\psi_{\iota}),$

the bias can be rewritten as:

$bias_{D}^{\iota}=\sigma^{2}_{\psi_{\iota}}R^{2}_{D,\iota}+\sigma^{2}_{f_{i}}-2cov(\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{\iota}]x_{k},f_{i}(x)).$

It remains to explicit the last term:

$\displaystyle-2cov(\sum_{k=1}^{D}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{\iota}]x_{k},f_{i}(x))=-2\sum_{k=1}^{D}cov(x_{k},f_{i})\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{\iota}]$
$\displaystyle=-2\sum_{k=1}^{D}\Big{\\{}\frac{cov(x_{k},f_{i})cov(x_{k},\psi_{\iota})}{\sigma^{2}_{x_{k}\lvert
x^{-k}}}-cov(x_{k},f_{i})\sum_{\begin{subarray}{c}h=1\\\ h\neq
k\end{subarray}}^{D}\Big{(}\frac{\rho_{x_{h},x_{k}\lvert
x^{-h,k}}}{\sigma_{x_{h}\lvert x^{-h}}\cdot\sigma_{x_{k}\lvert
x^{-k}}}\Big{)}cov(x_{h},\psi_{\iota})\Big{\\}}$
$\displaystyle=-2cov(f_{i},\mathbf{X})cov(\mathbf{X},\mathbf{X})^{-1}cov(\mathbf{X},\psi_{\iota})$
$\displaystyle=-2(cov(f_{i},\psi_{\iota})-cov(f_{i},\psi_{\iota}|\mathbf{X})),$

which concludes the proof.

The single-task and multi-task settings with transformed inputs
$\phi_{1},\dots,\phi_{d}$ can be proved with the same passages used in this
section to prove the results related to the regression with original features.
Therefore, the following theorem concludes the overview of results related to
biases, proving Theorem 3.2 of the main paper.

###### Theorem 0.B.3

In the asymptotic case, estimating the $i$-th task in a single-task setting,
with a set of transformed features $\phi_{1}\dots\phi_{d}$, the bias of the
model is:

$bias_{d}^{i}=\sigma^{2}_{f_{i}}(1-R^{2}_{d,i}).$

Finally, combining the aggregation of targets and the transformation of
features, we get the expression reported in Theorem 3.2 in the main paper:

$\begin{split}bias_{d}^{\iota}&=\sigma^{2}_{f_{i}}+\sigma^{2}_{\psi_{\iota}}R^{2}_{d,\iota}-2(cov(f_{i},\psi_{\iota})-cov(f_{i},\psi_{\iota}|\mathbf{\Phi}))\\\
&=\sigma^{2}_{f_{i}}-\sigma^{2}_{\psi_{\iota}}R^{2}_{d,\iota}+2(cov(\psi_{\iota},f_{i}-\psi_{\iota}|\mathbf{\Phi})-cov(\psi_{\iota},f_{i}-\psi_{\iota})).\end{split}$
(13)

###### Proof

In the first case, with the same considerations of the previous subsection, we
get:

$\displaystyle bias_{d}^{i}$
$\displaystyle=var_{x}(\sum_{k=1}^{d}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]\phi_{k})+var_{x}(f_{i}(x))-2cov(\sum_{k=1}^{d}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]\phi_{k},f_{i}(x))$
$\displaystyle=var_{x}(\sum_{k=1}^{d}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]\phi_{k})+var_{x}(f_{i}(x))-2cov(\sum_{k=1}^{d}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]\phi_{k},f_{i}(x))$
$\displaystyle=var_{x}(\sum_{k=1}^{d}\operatorname{\mathbb{E}}_{\mathcal{T}}[\hat{w}_{k}^{i}]\phi_{k})+\sigma^{2}_{f_{i}}-2\sigma^{2}_{f_{i}}R^{2}_{d,i}$
$\displaystyle=\sigma^{2}_{f_{i}}(1-R^{2}_{d,i}).$

Similarly, in the second case the theorem follows with the same argumentations
of Theorem , proving also Theorem 3.2 in the main paper.

##  Appendix 0.C Additional results: bias-variance comparisons between models

Single vs. Multi task model Firstly, focusing on the comparison between the
two settings with the original set of $D$ features and comparing the case with
a single-task with a multi-task model that learns on an aggregated target that
is the average of a cluster of features $\psi_{\iota}$, the asymptotic
decrease of variance when the target is aggregated is equal to:

$\displaystyle\Delta var_{D}^{i-\iota}$
$\displaystyle:=var_{D}^{i}-var_{D}^{\iota}$
$\displaystyle=\frac{\sigma^{2}_{i}-\bar{\sigma}_{\iota}^{2}}{(n-1)}\cdot\Big{\\{}\sum_{k=1}^{D}\frac{\sigma^{2}_{x_{k}}}{\sigma^{2}_{x_{k}\lvert
x^{-k}}}-2\sum_{k=1}^{D-1}\sum_{j=k+1}^{D}\frac{cov(x_{k},x_{j})\cdot\rho_{x_{k},x_{j}\lvert
x^{-k,j}}}{\sigma_{x_{k}\lvert x^{-k}}\cdot\sigma_{x_{j}\lvert
x^{-j}}}\Big{\\}}$
$\displaystyle=\frac{\sigma^{2}_{i}-\bar{\sigma}_{\iota}^{2}}{(n-1)}\cdot\mathrm{D}\overset{\sigma^{2}_{i}=\sigma^{2}_{k}\forall
k\in\mathcal{P}_{\iota}}{=}(\frac{K-1}{K})\cdot\frac{\sigma^{2}_{i}(1-\bar{\rho}_{i})}{(n-1)}\mathrm{D},$

Where $\bar{\rho}_{i}:=\frac{1}{\frac{K(K-1)}{2}}\cdot\sum_{h\neq
k}\rho_{h,k}$ is the average correlation among noises. As discussed in the
main paper, remembering that by non-negativity of variance
$\bar{\rho}_{i}\geq\frac{1}{1-K}$, no variance gain corresponds to
$\bar{\rho}_{i}=1$, while the maximum one corresponds to the minimum average
correlation $\bar{\rho}_{i}=\frac{1}{1-K}$. The gain is satisfactory when the
noises are independent, $\bar{\rho}_{i}=0$, leading to
$\frac{K-1}{K}\cdot\frac{\sigma^{2}_{i}}{(n-1)}\mathrm{D}$. This is in line
with the intuition that the learning problem is simpler when the target is an
average of less correlated tasks. However, this is balanced by the asymptotic
difference of biases, which increases if we aggregate tasks that are less:

$\displaystyle\Delta bias_{D}^{\iota-i}$
$\displaystyle=\sigma^{2}_{f_{i}}-\sigma^{2}_{\psi_{\iota}}R^{2}_{D,\iota}-2(cov(\psi_{\iota},f_{i}-\psi_{\iota})-cov(\psi_{\iota},f_{i}-\psi_{\iota}|\mathbf{X}))-$
$\displaystyle\sigma^{2}_{f_{i}}+\sigma^{2}_{f_{i}}R^{2}_{D,i}$
$\displaystyle=\sigma^{2}_{f_{i}}R^{2}_{D,i}-\sigma^{2}_{\psi_{\iota}}R^{2}_{D,\iota}-2(cov(\psi_{\iota},f_{i}-\psi_{\iota})-cov(\psi_{\iota},f_{i}-\psi_{\iota}|\mathbf{X}))$
$\displaystyle=\sigma^{2}_{\psi_{\iota}}R^{2}_{D,\iota}+\sigma^{2}_{f_{i}}R^{2}_{D,i}-2(cov(f_{i},\psi_{\iota})-cov(f_{i},\psi_{\iota}|\mathbf{X}))$

Full vs. Aggregated features Focusing on the single-task case for simplicity,
we can inspect the effect in terms of bias and variance due to the
transformation of features from the $D$ original ones to a set of $d<D$ basis
functions. The results are valid for general zero-mean transformation of
features, although the main interest of our approach is in the aggregation
with the mean of $d$ groups of features that form a partition of the original
ones.

The asymptotic decrease of variance in this single-task setting is therefore:

$\displaystyle\Delta var_{D-d}^{i}$
$\displaystyle=\frac{\sigma^{2}_{i}}{(n-1)}\cdot\Big{\\{}\Big{[}\sum_{k=1}^{D}\frac{\sigma^{2}_{x_{k}}}{\sigma^{2}_{x_{k}\lvert
x^{-k}}}-\sum_{k=1}^{d}\frac{\sigma^{2}_{\phi_{k}}}{\sigma^{2}_{\phi_{k}\lvert\phi^{-k}}}\Big{]}$
$\displaystyle-2\Big{[}\sum_{k=1}^{D-1}\sum_{j=k+1}^{D}\frac{cov(x_{k},x_{j})\cdot\rho_{x_{k},x_{j}\lvert
x^{-k,j}}}{\sigma_{x_{k}\lvert x^{-k}}\cdot\sigma_{x_{j}\lvert x^{-j}}}$
$\displaystyle-\sum_{k=1}^{d-1}\sum_{j=k+1}^{d}\frac{cov(\phi_{k},\phi_{j})\cdot\rho_{\phi_{k},\phi_{j}\lvert\phi^{-k,j}}}{\sigma_{\phi_{k}\lvert\phi^{-k}}\cdot\sigma_{\phi_{j}\lvert\phi^{-j}}}\Big{]}\Big{\\}}$
$\displaystyle=\frac{\sigma^{2}_{i}}{(n-1)}\cdot(D-d).$

On the other hand, the increase of bias is equal to:

$\Delta
bias_{d-D}^{i}=\sigma^{2}_{f_{i}}(1-R^{2}_{d,i})-\sigma^{2}_{f_{i}}(1-R^{2}_{D,i})=\sigma^{2}_{f_{i}}(R^{2}_{D,i}-R^{2}_{d,i}).$

From the variance comparisons, we clearly see that the asymptotic advantage of
dimensionality reduction resides in a reduced number of coefficients that need
to be estimated. Additionally, the bias shows that it is profitable to reduce
features as long as the amount of information shared with the target, measured
in terms of the (squared) coefficient of multiple correlation, remains
satisfactory.

##  Appendix 0.D Algorithm

Algorithm 2 reports the pseudo-code of the auxiliary aggregation loop
considered both for the target and the task aggregation phase. In particular,
in the first phase, it receives as input the full set of targets, the set of
features, the hyperparameter, and an indication of the phase. The loop
iterates over targets, forming the partition through the exploitation of the
threshold function defined in Algorithm 1. Similarly, in the second phase, the
algorithm receives the current aggregated target and the full set of features,
iteratively identifying the feature aggregations exploiting the threshold
function for features defined in Algorithm 1.

Algorithm 2 Auxiliary loop for the iterative aggregation of a set of variables
with the mean, used both for targets and features, holds for generic random
vector.

Input variables $\bm{z}=\\{z_{1},\dots,z_{\lvert\bm{z}\rvert}\\}$; $D$ input
features $\bm{x}=\\{x_{1},\dots,x_{D}\\}$ if we are aggregating targets, one
aggregated target $y=\psi_{\iota}$ if we are aggregating features for an
aggregated task in the second phase, flag _phase_ identifies it; tolerance
$\epsilon$.

reduced variables (tasks or features):
$\\{z_{\mathcal{P}_{1}},\dots,z_{\mathcal{P}_{K}}\\}$

function Aggregation$(\bm{z},phase,\epsilon,\bm{x},y)$

$\triangleright$ Iterative aggregation of sets of variables with the mean,
used both for targets and features, holds for generic random vector $\bm{z}$

$\bm{\mathcal{P}}\leftarrow\\{\\}$

$\mathcal{V}\leftarrow\\{\\}$ $\triangleright$ Set of already considered
elements

for each $i\in\\{1,\dots,\lvert\bm{z}\rvert\\}$ do

if $i\not\in\mathcal{V}$ then

$\mathcal{P}\leftarrow\\{i\\}$

$\mathcal{V}\leftarrow\mathcal{V}\cup\\{i\\}$

for each $j\in\\{i+1,\dots,\lvert\bm{z}\rvert\\}$ do

$z_{\mathcal{P}}\leftarrow mean(z_{k}:k\in\mathcal{P})$ $\triangleright$ Mean
of elements in current set

if $phase==1$ then$\triangleright$ First phase, we are aggregating targets

$threshold\leftarrow$
Compute_threshold_targets$(\bm{x},z_{\mathcal{P}},z_{j},\epsilon)$

else

$\bm{x}_{curr}\leftarrow\\{x_{k}:k\notin\mathcal{V}\\}\cup\\{mean(x_{h}:h\in\mathcal{P}_{k}),k\in\\{1,\dots,\lvert\bm{\mathcal{P}}\rvert\\}\\}$

$threshold\leftarrow$
Compute_threshold_features$(\bm{x}_{curr},y,z_{\mathcal{P}},z_{j},\epsilon)$

end if

if $threshold==True$ then $\triangleright$ Add $j$-th element to current set

$\mathcal{P}\leftarrow\mathcal{P}\cup\\{j\\}$

$\mathcal{V}\leftarrow\mathcal{V}\cup\\{j\\}$

end if

end for

$\bm{\mathcal{P}}\leftarrow\bm{\mathcal{P}}\cup\\{\mathcal{P}\\}$

end if

end for

$K\leftarrow\lvert\bm{\mathcal{P}}\rvert$

for each $k\in\\{1,\dots,K\\}$ do

$z_{\mathcal{P}_{k}}=mean(z_{h}:h\in\mathcal{P}_{k})$

end for

return $\\{z_{\mathcal{P}_{1}},\dots,z_{\mathcal{P}_{K}}\\}$

end function

##  Appendix 0.E Experimental Validation

###  0.E.1 Synthetic Experiments

This section includes additional results related to subsection 5.1 of the main
paper, where synthetic experiments are introduced. Considering the
configuration described in the main paper, where the effect of varying
different parameters has been described, Figure 4 shows the test performances
varying the hyperparameters of the NonLinCTFA approach. From the first column
of the figure, we can see that $\epsilon_{1}$ regulates the propensity to
aggregate targets, with a value close to zero that balances the aggregations,
large negative value that tends to produce a singleton with all tasks averaged
(which becomes detrimental for the final performance on the original tasks
when too many aggregations are performed) and large positive values that tend
to the single-task case, without any aggregation. In this case, the subsequent
feature aggregation phase is beneficial when sufficient information is
preserved by the task aggregation. The second column of the figure shows the
behavior of the algorithm varying the hyperparameter $\epsilon_{2}$, which is
stable for values close to zero, while it tends to aggregate all features when
the value is large, becoming detrimental.

Figure 4: Test Mean Squared Error (first row) and corresponding percentage
decrease w.r.t. the single task model (second row), varying one hyperparameter
at a time, only aggregating targets (Phase I) or adding the feature
aggregation (Phase I+II). The first parameter $\epsilon_{1}$ refers to the
propensity of the user to aggregate tasks. On the other hand, the second
parameter $\epsilon_{2}$ refers to the propensity of the user to aggregate
features.

###  0.E.2 Real World Experiments

In this subsection we provide additional details and results on the real-world
experiments performed. All experiments have been conducted with an HPC-system,
specifically a BullSequana XH2000 supercomputer using the 3rd generation of
AMD EPYC CPUs (Milan), NVIDIA A100 GPUs, and a 130 Petabyte DDN filesystem. On
this machine, we exploited $4$ nodes of CPU and $2$ GPUs for training neural
networks.

In the first set of experiments, we consider School and SARCOS data, as
described in the main paper, comparing the regression scores obtained with LR,
SVR, and MLP after the application of the NonLinCTFA. We selected these
datasets to test the algorithm on classical data, exploiting the
implementation of some MTL baselines available for benchmarking, as described
in Section 5. The SARCOS dataset considers $D$ shared features, therefore the
main algorithm has been applied on it. On the contrary, the School dataset
presents the same features, measured differently for each school (task).
Therefore, we applied the variant of the algorithm discussed in 8.
Additionally, given the different amounts of samples available for each
school, in the aggregated case, we randomly extracted a number of samples from
the school with more samples equal to the number of samples of the other
school, ordering both samples depending on the target and aggregating
consequently. The idea is that in this way, we are aggregating a sample of the
best student of one school with a sample of the best student of the other
school, then the second best student of the first is aggregated with the
second best of the other school, and so on. We considered the NRMSE as a
performance score, being a valid elaboration of the MSE, the metric that we
are optimizing, and exploiting the fact that it was already implemented for
the benchmark MTL approaches considered.

In the second real-world framework, we consider the QM9 dataset, as discussed
in Section 5. This is a recent, challenging graph-based dataset, where
thousands of molecule structures are provided as inputs, with their
characteristics being the targets. In this case, we could therefore apply
Algorithm 1. We had to adapt the input to be a tabular set of features,
therefore we averaged the features of nodes and edges, as done in other
literature approaches. To provide a fair comparison, we exploited a state-of-
the-art library, LibMTL, as described in the main paper, producing $14$
results with a hard parameter sharing graph neural network, trained
considering $14$ different weighting strategies (all described in the paper of
the library [21] and implemented from recent papers. Additionally, we added
some results, in terms of MSE, from a recent paper that considers the QM9
dataset similarly, to show the best performance, to our knowledge, of a GNN
architecture on these data. From the results, we conclude that our approach is
able to outperform the single-task counterparts and to be competitive with
graph neural network architectures with hard parameter sharing. However, it
does not achieve the performance of the best GNN architecture to the best of
our knowledge. This is reasonable since we are transforming the graph inputs
to tabular data, losing the graph information. However, we consider this
result an instructive example of a tradeoff between reaching the best possible
performance with a completely black box approach, such as a complex graph
neural network, and an approach like NonLinCTFA, where the entire pipeline is
interpretable, aimed to produce a simple yet effective ML workflow.

Finally, we conclude the experimental section, providing an example of the
application of our method to Earth science data. In particular, we consider
about $30000$ hydrological sub-basins of Europe, for which we have a satellite
target variable and a set of meteorological features and climate indices at
each location. In this case, we show the performance of our approach, w.r.t.
the single-task counterparts (which corresponds to $\epsilon=1$, showing a
clear improvement of the MSE when a balanced aggregation is performed and a
deterioration of the performances when the algorithm aggregates too much or
when it does not aggregate enough. We did not consider additional MTL
baselines in this context since its main purpose is to conclude the
experimental validation section providing a reliable example of application in
line with the motivational example that has been mentioned throughout the
paper. Additionally, the main applicative interest is to provide an
improvement w.r.t. single task, preserving the interpretability of the
pipeline and improving existing climate indices, as will be further described
in a future applicative work.