# Learning Multitask Gaussian Bayesian Networks

Shuai Liulabel=e1 [<EMAIL_ADDRESS>Yixuan
<EMAIL_ADDRESS>[ Baojuan<EMAIL_ADDRESS>[
Huaning<EMAIL_ADDRESS>[ Xiangyu Changlabel=e5 [
<EMAIL_ADDRESS>School of Management, Xi’an Jiaotong
University,  School of Statistics and Management, Shanghai University of
Finance and Economics,  School of Biomedical Engineering, Fourth Military
Medical University,  Department of Psychiatry, Xijing Hospital, Air Force
Medical University,

###### Abstract

Major depressive disorder (MDD) is closely related to the brain functional
connectivity alterations in patients, which can be uncovered by resting-state
functional magnetic resonance imaging (rs-fMRI) data. We consider the problem
of identifying alterations of the brain functional connectivity for a single
MDD patient, which is particularly difficult as the amount of data collected
during an fMRI scan is too limited to provide sufficient information for
individual analysis. Additionally, rs-fMRI data usually have the
characteristics of incompleteness, sparsity, variability, high dimensions, and
large noise.

To address the challenges above, especially learning from insufficient rs-fMRI
data, we propose a multitask Gaussian Bayesian network framework to identify
disorder-induced alterations for individual MDD patients. We achieve this goal
by treating each patient in a class of observations as a task, and leveraging
data from others. Since the patients have the same type of disorder, they
share similar connectivity alterations in brain functional networks.
Therefore, we introduce a shared prior covariance matrix among all tasks and
utilize the hierarchical model structure to estimate multiple Gaussian
Bayesian networks. For computing, we derive closed-form expressions for the
complete likelihood function and its gradient and use the Monte Carlo
expectation-maximization algorithm to efficiently search for the approximately
best network structures. We assess the performance of the proposed method with
extensive and systematic numerical experiments and highlight its usefulness in
analyzing rs-fMRI data and facilitating the study of MDD.

Bayesian network,

multitask learning,

major depressive disorder,

Hamiltonian Monte Carlo,

expectation-maximization algorithm,

###### keywords:

, , , , and

##  1 Introduction

Major depressive disorder (MDD) is a common mental disorder characterized by
low mood, worthlessness, anxiety, and cognitive impairments (Fu et al., 2008).
MDD causes more than 800,000 deaths worldwide every year and is the second
leading cause of disability worldwide, with point prevalence exceeding 4% (Yan
et al., 2019). Therefore, understanding the mechanism of MDD is crucial for
effective prevention, diagnosis, and treatment. The existing knowledge about
MDD is very limited, but recent studies have found that MDD is closely related
to the alterations of brain functional connectivity (Wu et al., 2011; Liu et
al., 2013; Yamashita et al., 2020). Resting-state functional magnetic
resonance imaging (rs-fMRI) is an effective modality for understanding the
biological and cognitive information of mental disorders (Xia and He, 2017).
It can be used for the noninvasive investigation of MDD by quantifying
functional connections of the whole brain that are related to spontaneous
blood-oxygen-level-dependent signal fluctuations (Smith et al., 2013).

Bayesian network (BN), as an approach to studying the cause of interdependence
between variables in a system, has drawn much attention (Koller and Friedman,
2009). BNs visualize the way of information propagation in an intuitive
graphical representation (Koski and Noble, 2012). Such an interpretable
structure makes them widely applied in scientific domains such as
bioinformatics and neuroscience (Bielza and Larranaga, 2014). Many BN models
have been developed for disease screening, diagnosis, treatment, and
prognosis. For example, Moreira et al. (2016) builds a BN model for the
diagnosis of preeclampsia, and Zarikas, Papageorgiou and Regner (2015)
utilizes the BN model to predict the future severity of acute lung disease.

However, BNs have not been widely used in clinical practice because of the gap
between developing an accurate BN and proving its clinical availability. This
is more obvious when using BNs to study brain functional connectivity of MDD
using rs-fMRI data. The main barriers are: (1) In addition to high
dimensionality, rs-fMRI data usually have the characteristics of
incompleteness, sparsity, variability, high noise, and significant error. When
the number of variables exceeds 50 to 100, the number of possible network
structures increases super-exponentially. Although some progress has been made
in learning large BNs (Gámez, Mateo and Puerta, 2011; Aragam, Gu and Zhou,
2017), learning and interpreting BNs from rs-fMRI data is still challenging.
(2) Due to the difficult and expensive data acquisition, the amount of rs-fMRI
data is relatively small. In the Strategic Research Program for Brain Sciences
database, only 255 MDD patients’ rs-fMRI data were collected (Yamashita et
al., 2021). Moreover, because of the considerable heterogeneity of depression,
the brain functional connectivity in patients with MDD varies from person to
person (Young et al., 2016), and there are different subtypes of depression
itself (Beijers et al., 2019). Therefore, the available rs-fMRI datasets are
too small to simultaneously learn robust BNs for individuals or subtype groups
of MDD. Learning more helpful information from limited data is still one of
the research hotspots of BN learning.

To learn robust BNs efficiently via insufficient rs-fMRI data, we propose a
multitask learning (MTL) method to learn multiple Gaussian Bayesian networks
(GBNs) simultaneously, which we call the multitask Gaussian Bayesian network
(MTGBN). Our original interest in developing this method arises from the
noninvasive biomarker discovery in rs-fMRI data. For example, in the study of
brain functional connectivity alterations for a single MDD patient, data
derived from an individual’s neuroimage do not deliver enough knowledge due to
the complex data structure and high data dimensionality. However, since the
patients are related in the sense that they have the same disorder, they share
disorder-induced connectivity alteration in brain functional networks.
Therefore, in the BN learning, we also hope to leverage data from other
patients with the assumption that such alterations show some degrees of
similarity. Such models would help practitioners understand how systems are
structured jointly from related tasks.

Few works exist on using fMRI data to learn network structures from related
tasks to our best knowledge. Chen et al. (2021) considers estimating multiple
differential networks with fMRI datasets from multiple research centers and
exploits sparse group minimax concave penalty to induce sparse structures of
each differential network. He et al. (2021) proposes an algorithm called
Trans-Copula-CLIME for estimating undirected graphical models from fMRI
datasets, which characterizes the similarity between tasks by the sparsity of
a divergence matrix. Different from the work by Chen et al. (2021) and He et
al. (2021), we focus on GBN learning. The benefit that GBN brings to studying
network structures is that GBN can characterize the direct dependencies
between variables, which can provide more information than differential
networks and undirected graphical models. This is the first work to use the
Bayesian approach to learn multiple GBNs from related tasks on MTL, as far as
we know.

The proposed MTGBN method improves the learning ability of GBNs via the MTL
technique. Specifically, we treat each subject in a class of observations as a
task and assume that the tasks share some common prior knowledge. A Bayesian
hierarchical model thus models the observed data. This setting can help us
learn useful information from limited data. Remarkably, in the structure
learning of BNs, we derive a closed-form formula for the complete likelihood
function and its gradient and use the Monte Carlo expectation-maximization
algorithm (MCEM, Wei and Tanner, 1990) to efficiently search for the
approximately best BN structures. The contributions of this article are
threefold.

  * •

The proposed MTGBN method is able to learn BNs from insufficient high
dimensional data by jointly learning multiple BN structures from related
observations. It addresses one common problem in pattern recognition: how to
learn network structures from related observations such that the
representations can share related information while also maintaining subject-
specific characteristics.

  * •

MTGBN overcomes the computational intractability issue of the existing MTL
framework for GBNs. We derive a closed-form expression for the complete
likelihood function and provide exact gradient evaluations to the Hamiltonian
Monte Carlo (HMC, Neal et al., 2011) sampler that is used by the MCEM
algorithm.

  * •

We apply MTGBN to rs-fMRI data of MDD patients and show that MTGBN can
effectively identify the increase of functional connectivity in MDD patients,
especially the connectivity within the dorsal medial prefrontal cortex.
Doctors can use it to understand the potential abnormalities of brain
functional connections of individual MDD patients and carry out personalized
treatments.

This article is organized as follows. In Section 2 we begin by introducing the
rs-fMRI dataset analyzed in this work. In Section 3, we briefly review BN
learning and MTL methods. We present the formulation of MTGBN in Section 4,
and introduce its structure learning in Section 5. An efficient computing
algorithm for MTGBN is given in Section 6. We conduct simulations in Section
7, and analyze a real rs-fMRI dataset and assess the advantages of MTGBN in
Section 8. Finally, we conclude this article in Section 9.

##  2 Real-world rs-fMRI Data

###  2.1 Data acquisition

The data analyzed in this study were obtained from the DecNef Project Brain
Data Repository111https://bicr-resource.atr.jp/srpbsopen/, and are a subset of
the SRPBS multi-disorder MRI dataset (Tanaka et al., 2021). The rs-fMRI data
for each participant at resting state were obtained using 3-Tesla Siemens
whole body MRI systems at Center of Innovation in Hiroshima University. During
the fMRI scan, participants were instructed to be relax, stay awake and
refrain from thinking anything particular. The rs-fMRI data have the following
parameters: TR = 2.5s, TE = 30ms, flip angle = $80^{\circ}$, matrix size =
$64\times 64$, field of view = $212$mm, in-plane resolution = $3.3\times 3.3$,
slice gap = 0.8mm, slice thickness = 3.2mm, number of slices = 40, number of
volumes = 240 + 4 (dummy), total scan time = 10min + 10s (dummy). The
T1-weighted images were then acquired with the following parameters: voxel
size = $1\times 1\times 1$, TR = 2300ms, TE = 2.98ms, TI = 900ms; flip angle =
$90^{\circ}$, FOV = $256$, matrix = $256\times 256$.

###  2.2 Samples

We use the rs-fMRI data of 15 MDD patients and 15 healthy control (HC)
participants for the experiment. The selection criteria are as follows. The
MDD patients are aged between 40 and 50, right-handed, free of other diseases,
and their BDI scores are between 20 and 35. BDI is a clinical instrument for
the severity of depression. The higher the score, the more severe the
depression. The HC participants are selected to purposely match the MDD
patients in terms of other selection criteria except for BDI score, _i.e._ ,
they are in the same age cohort as the MDD patients, right-handed, free of
other diseases, but with a BDI score between 0 and 7 (_i.e._ , they do not
have depression).

We applied the following pre-processing steps to the data. For each
participant, 230 fMRI volumes were remained after removing the first ten
volumes. Slice acquisition time and head motion were corrected on the
remaining images. Participants with excessive head motion were not enrolled.
The T1-weighted images use unified segment and diffeomorphic anatomical
registration to guide fMRI registration through exponentiated Lie algebra as
implemented in the SPM12
software222https://www.fil.ion.ucl.ac.uk/spm/software/spm12/. The rs-fMRI data
were smoothed with a 6mm full-width at half-maximum Gaussian kernel and
further denoised by regressing out several nuisance signals, including the
Friston-24 head motion parameters and signals from cerebrospinal fluid and
white matter, before linear detrending and temporal band-pass filtering
($0.01-0.08$ Hz).

After pre-processing, we select 30 brain regions (each viewed as a region of
interest, ROI) based on three different network templates related to the
default mode network (DMN, Raichle and Snyder, 2007), including the Andrews-
Hanna DMN (aDMN, Buckner, Andrews-Hanna and Schacter, 2008), the dorsal DMN
(dDMN, Chen et al., 2017), and the ventral DMN (vDMN, Lee, Parthasarathi and
Kable, 2021). Details of brain regions are presented in Table 5. The voxel-
wise fMRI time courses are averaged into one regional average time course
within each selected region. Then each participant’s dataset has the form of a
$230\times 30$ matrix (230 sampling points in each regional average fMRI time
course and 30 selected brain regions). The dataset obtained through the steps
above can be reasonably considered to follow Gaussian distributions because
each brain region includes at least several hundred voxels and the measurement
data of each region is the regional average over the belonging voxels.
Therefore, the central limit theorem is applicable. Also, such hypotheses are
common for fMRI studies (Valdes-Sosa et al., 2011).

##  3 Background and related work

###  3.1 Learning Bayesian networks

BN is a formal method to analyze decision-making strategies under uncertain
conditions (Pearl, 1988). Given a random vector
$\boldsymbol{X}=(X_{1},\ldots,X_{p})^{\top}$, we use $\mathbb{P}(\mathbf{x})$
to denote its joint probability mass function (for the discrete case) or the
joint density function (for the continuous case), where $\mathbf{x}$ is an
observation of $\boldsymbol{X}$. From the chain rule of probability, we can
write $\mathbb{P}(\mathbf{x})$ as

$\mathbb{P}(\mathbf{x})=\mathbb{P}(x_{1})\prod_{i=2}^{p}\mathbb{P}\left(x_{i}\mid
x_{1},\ldots,x_{i-1}\right).$ (3.1)

For each $X_{i}$, let $\Pi_{i}\subseteq\left\\{X_{1},\ldots,X_{i-1}\right\\}$
be a set of variables such that

$\mathbb{P}\left(x_{i}\mid
x_{1},\ldots,x_{i-1}\right)=\mathbb{P}\left(x_{i}\mid\Pi_{i}\right).$ (3.2)

A BN $B=(\mathcal{G},\Theta)$ can represent the joint distribution of
$\boldsymbol{X}$ through a directed network model, where $\mathcal{G}$ is a
network structure that encodes the assertions of conditional independence in
equation (3.2), and $\Theta$ includes the parameters of the BN model
corresponding to $\mathcal{G}$ (Koller and Friedman, 2009). In particular,
$\mathcal{G}$ is a directed acyclic graph (DAG) such that (1) each variable in
$\boldsymbol{X}$ corresponds to a node in $\mathcal{G}$, and (2) the parents
of the node corresponding to $X_{i}$ are the nodes corresponding to the
variables in $\Pi_{i}$. Throughout this article, we use $X_{i}$ to refer to
both the variable and its corresponding node in the graph. Combining equations
(3.1) and (3.2), we can see that the joint distribution for $\boldsymbol{X}$
can be represented graphically by a BN, namely,

$\mathbb{P}(\mathbf{x})=\mathbb{P}(x_{1})\prod_{i=2}^{p}\mathbb{P}\left(x_{i}\mid\Pi_{i}\right).$

GBNs are a class of important and widely-used BNs (Geiger and Heckerman,
1994). In GBNs, $\boldsymbol{X}$ follows a multivariate normal distribution
$\mathcal{N}\left(\mu,\Sigma\right)$, with the joint density function

$\mathbb{P}(\mathbf{x})=(2\pi)^{-p/2}\left|\Sigma\right|^{-1/2}\exp\left(-\frac{1}{2}(\mathbf{x}-\mu)^{\top}\Sigma^{-1}(\mathbf{x}-\mu)\right),$

where $\mu$ is the expectation of $\boldsymbol{X}$, $\Sigma$ is the covariance
matrix, and $|\Sigma|$ is the determinant of $\Sigma$. Clearly, the parameters
of a GBN are $\Theta=\\{\mu,\Sigma\\}$.

The structure $\mathcal{G}$ and the parameters $\Theta$ can be learned from
the observed data, which we call the BN learning procedure. The primary
purpose of BN learning is to estimate a BN $B=(\mathcal{G},\Theta)$ that can
best explain the $n$ independent observations
$D=\\{\mathbf{x}_{i}=(x_{i1},\dots,x_{ip})^{\top}\\}_{i=1}^{n}$, leading to
two major tasks: the structure learning, which attempts to learn a DAG
$\mathcal{G}$ that encodes the structural dependence of the variables in
$\boldsymbol{X}$, and the parameter learning, which estimates the parameters
$\Theta$ associated with $\mathcal{G}$.

#### Structure learning

Structure learning is an essential part of BNs, which can mainly be divided
into two categories: (1) constraint-based methods; (2) score-based methods.
Constraint-based methods try to identify dependent and independent areas
within the data (Koller and Friedman, 2009). One widely-used method for
structure learning is the SGS algorithm (Spirtes et al., 2000) that identifies
the dependencies of variable pairs by evaluating their independence. This
method is computationally expensive because the number of tests overgrows as
the network size increases. Compared with SGS, the PC algorithm (Spirtes and
Glymour, 1991) is more efficient because it has a size threshold. The PC
algorithm assumes that all nodes in the undirected graph are connected and
uses hypothesis tests to delete the connections.

Score-based methods attempt to compute a score measuring all possible network
structures, using a scoring function such as BIC (Schwarz et al., 1978), K2
(Cooper and Herskovits, 1992), and BDeu (Heckerman, Geiger and Chickering,
1995). Because learning BN structures from data is an NP-hard problem,
existing score-based methods usually use heuristic search techniques to reduce
the computational cost. For large networks, the search space is extremely
large, and ordinary score-based methods also become intractable. Therefore,
Friedman, Nachman and Pe’er (2013) proposes to improve the learning efficiency
of score-based methods by restricting the search space.

#### Parameter learning

There are two main methods for parameter learning: (1) maximum likelihood
estimation (MLE, see Ji, Xia and Meng, 2015); (2) Bayesian estimation (Geiger
et al., 1997). MLE method estimates the parameters by maximizing the
likelihood of the data. However, MLE becomes problematic when the data are
sparse, as the dataset does not contain all possible combinations of the
involved variables. In such cases, the Bayesian estimation is preferred, which
places a prior distribution on the variables (Altendorf, Restificar and
Dietterich, 2005). It should be pointed out that parameter learning itself is
not sufficient, as it is based on the assumption that the network structure is
known in advance, or at least one structure is held fixed during parameter
learning.

In this work, we focus on the score-based structure learning. The score
functions for BN are typically derived via a Bayesian approach. For example,
one Bayesian measure of the goodness of a BN structure $\mathcal{G}$ is its
posterior probability given the data $D$:

$\mathbb{P}(\mathcal{G}\mid D)=\frac{1}{c}\
\mathbb{P}(\mathcal{G})\times\mathbb{P}(D\mid\mathcal{G}),$ (3.3)

where $\mathbb{P}(\mathcal{G})$ is the prior of $\mathcal{G}$, and
$c=\mathbb{P}(D)=\sum_{\mathcal{G}}\mathbb{P}(\mathcal{G})\times\mathbb{P}(D\mid\mathcal{G})$
is a normalizing constant. Computing $c$ is in general intractable, as there
may be a huge number of network structures to sum over, even in a small
domain. But we never need to explicitly compute $c$, since maximizing the
posterior probability is equivalent to maximizing the joint distribution
$\mathbb{P}(\mathcal{G})\times\mathbb{P}(D\mid\mathcal{G})=\mathbb{P}(D,\mathcal{G})$,
which we use as the score function in what follows.

Once the score function is determined, it can then be used to find the best
$\mathcal{G}$ among the set of all network structures. Unfortunately, this
optimization problem is NP-hard. Therefore, heuristic algorithms, including
greedy search, greedy search with restarts, best-first search, and Monte Carlo
methods, are usually required (Heckerman, Geiger and Chickering, 1995).

###  3.2 Multitask learning

General machine learning methods usually require many samples to learn a
robust model. In the case of insufficient data, MTL (Caruana, 1997) is a
general approach that incorporates task relatedness during learning and can
learn robust models (Zhang and Yang, 2021). Given $m$ learning tasks
$\mathcal{T}=\left\\{T_{1},\ldots,T_{m}\right\\}$, where at least one of them
is related, the aim of the MTL is to improve the performance and
generalization of a model $T_{i}$ by using the knowledge contained in all or
some of related $m$ tasks. The tasks in MTL can be general learning tasks such
as supervised tasks (Schwab et al., 2018), unsupervised tasks (Cao et al.,
2019), semi-supervised tasks (Rei, 2017), reinforcement learning tasks (Teh et
al., 2017), or graphical tasks (Honorio and Samaras, 2010).

According to Zhang and Yang (2021), the current research of MTL focuses on
characterizing the relatedness between tasks, which can be divided into three
classes, including instance-based MTL (Bickel et al., 2008), feature-based MTL
(Evgeniou and Pontil, 2007), and parameter-based MTL (Parameswaran and
Weinberger, 2010). In this paper, we focus on the parameter-based MTL. The
most common approach in parameter-based MTL is the task relation learning
approach, which quantifies the relatedness among tasks in three ways: task
similarity, task correlation, and task covariance. Early work such as Evgeniou
and Pontil (2004) considers the task relations as prior information. Several
studies have taken the relatedness into account by placing a common prior on
model parameters of the tasks (Evgeniou et al., 2005), e.g., with Gaussian
copula models (Gonçalves, Von Zuben and Banerjee, 2016) and matrix generalized
inverse Gaussian distribution (Li et al., 2015). However, when task relations
are not available in applications, it is necessary to automatically learn task
relationships from data. For instance, Williams, Bonilla and Chai (2007)
proposes a multitask Gaussian process to learn task covariance from data, and
Zhang and Yeung (2010) proposes a multitask generalized $t$ process by placing
an inverse-Wishart prior.

##  4 The proposed MTGBN framework

The MTGBN framework is primarily motivated by various real-world problems.
Consider a medical task that studies a specific disease with $m$ subtypes, and
suppose that the collected variables for each subtype are modeled by a BN.
However, the data sizes for some subtypes are small, making it challenging to
learn robust network structures for all subtypes. Also, it is inappropriate to
pool the data and only learn a single BN shared by all the subtypes, as it
ignores the heterogeneity between subtypes.

In the above scenario, the MTL method that learns multiple BNs simultaneously
is much preferred, where we treat the modeling of each subtype as one task. In
this article, we focus on GBNs and propose the MTGBN model that can address
the issue of insufficient sample sizes in individual GBN learning tasks. MTGBN
achieves this target by assuming that the $m$ subtypes share some common
information, compensating for the lack of data for minor subtypes.

Formally, let $\mathcal{D}=\left\\{D_{1},\ldots,D_{m}\right\\}$ be $m$ related
observational datasets, where
$D_{i}=\\{\mathbf{x}_{i1},\ldots,\mathbf{x}_{in_{i}}\\}$ represents the
collected data for task $i$. Assume that each observation
$\mathbf{x}_{ij}\in\mathbb{R}^{p}$ is the realization of a multivariate normal
distribution with zero mean and covariance matrix $\Sigma_{i}$.

In many applications, we have prior beliefs about the relatedness of tasks
based on metadata or domain knowledge. In MTGBN, we define a shared prior
parameter, $\Sigma_{h}$, to represent the relatedness of all tasks. Figure 1
shows the generative framework of the MTGBN method, where the covariance
matrices $\left\\{\Sigma_{i}\right\\}_{i=1}^{m}$ are sampled from a common
prior distribution parameterized by $\Sigma_{h}$, and each $\Sigma_{i}$
generates a dataset $D_{i}$ for task $i$. $\Sigma_{h}$ is inferred from all
the tasks, which helps transfer the learned knowledge from one task to
another.

Figure 1: Illustration of the generative framework of MTGBN

We formalize the proposed MTGBN framework as follows. First, the parameter
$\Sigma_{h}$ is assigned a uniform prior distribution over all $p\times p$
positive definite matrices. Then given $\Sigma_{h}$, the $m$ covariance
matrices $\boldsymbol{\Sigma}=\left\\{\Sigma_{i}\right\\}_{i=1}^{m}$ are
independently sampled from an inverse Wishart distribution
$IW\left(\Sigma_{h},\nu_{0}\right)$ (Robert, 2014), where $\nu_{0}>p-1$ and
with the following joint density function:

$\displaystyle\mathbb{P}\left(\boldsymbol{\Sigma}\mid\Sigma_{h}\right)=\prod_{i=1}^{m}\mathbb{P}\left(\Sigma_{i}\mid\Sigma_{h}\right),$
(4.1)
$\displaystyle\mathbb{P}\left(\Sigma_{i}\mid\Sigma_{h}\right)=\frac{\left|\Sigma_{h}\right|^{\nu_{0}/2}}{2^{\nu_{0}p/2}\Gamma_{p}\left(\nu_{0}/2\right)}\left|\Sigma_{i}\right|^{-\left(\nu_{0}+p+1\right)/2}\cdot\exp\left\\{-\frac{1}{2}\operatorname{tr}\left(\Sigma_{h}\Sigma_{i}^{-1}\right)\right\\},$

where for $\alpha>(m-1)/2$,

$\Gamma_{m}(\alpha)=\pi^{m(m-1)/4}\prod_{i=1}^{m}\Gamma\left(\alpha-(i-1)/2\right)$

is the multivariate gamma function (Muirhead, 2009), $\Gamma(\cdot)$ is the
gamma function, $p$ denotes the dimension of $\Sigma_{i}$, and $\nu_{0}$ is
the degree of freedom in the inverse Wishart distribution, viewed as a fixed
hyperparameter.

Finally, we assume that the observed data $D_{i}$ in class $i$ follow a zero
mean multivariate normal distribution, and the $m$ datasets
$\mathcal{D}=\left\\{D_{i}\right\\}_{i=1}^{m}$ are mutually independent given
the covariance matrices
$\boldsymbol{\Sigma}=\left\\{\Sigma_{i}\right\\}_{i=1}^{m}$. Clearly, the
joint density function of $\mathcal{D}$ given $\boldsymbol{\Sigma}$ is

$\displaystyle\mathbb{P}\left(\mathcal{D}\mid\boldsymbol{\Sigma}\right)=\prod_{i=1}^{m}\mathbb{P}\left(D_{i}\mid\Sigma_{i}\right),$
(4.2)
$\displaystyle\mathbb{P}\left(D_{i}\mid\Sigma_{i}\right)=(2\pi)^{-n_{i}p/2}\left|\Sigma_{i}\right|^{-n_{i}/2}\exp\left\\{-\frac{1}{2}\operatorname{tr}\left(n_{i}\Sigma_{i}^{-1}S_{i}\right)\right\\},$

where $n_{i}$ denotes the sample size of $D_{i}$, and $S_{i}$ is the sample
covariance matrix of $D_{i}$.

The parameter learning of MTGBN is concerned with inferring all the parameters
$\Sigma_{h}$ and $\boldsymbol{\Sigma}$ from
$\mathcal{D}=\left\\{D_{1},\ldots,D_{m}\right\\}$. However, one major
advantage of GBN is that it can reveal the essential causal relationship
between target variables, and hence in this work our main interest is on the
structure learning of MTGBN. For convenience, in what follows we omit the
parameter component $\Theta$ in a BN $B=(\mathcal{G},\Theta)$, and use $B_{i}$
to represent the BN structure of the collected dataset for class $i$.
Therefore, the learning goal of MTGBN is to simultaneously infer the $m$ GBNs
$\mathcal{B}=\left\\{B_{1},\ldots,B_{m}\right\\}$ from $m$ related datasets
$\mathcal{D}=\left\\{D_{1},\ldots,D_{m}\right\\}$, with $B_{i}$ corresponding
to $D_{i}$.

##  5 Structure learning of MTGBN

As is described in Section 3.1, to learn the BN structures from observations,
we need to first construct an appropriate score function, and then use
heuristic methods to search for the best network structure based on the score
function. Similar to equation (3.3), we follow the likelihood principle,
_i.e._ , taking $\mathbb{P}\left(\mathcal{B},\mathcal{D}\right)$ as the score
function. However, one major computational challenge of MTGBN stems from the
hypothesis prior covariance matrix $\Sigma_{h}$. Although the introduction of
$\Sigma_{h}$ allows the sharing of information among tasks, it also brings
difficulties to maximizing the score function. The hierarchical model
illustrated by Figure 1 defines the joint distribution of $\Sigma_{h}$,
$\mathcal{B}$, and $\mathcal{D}$, and we treat the parameter $\Sigma_{h}$ as a
latent variable. Therefore, $\mathbb{P}\left(\mathcal{B},\mathcal{D}\right)$
is computed by integrating out $\Sigma_{h}$ in the joint density function
$\mathbb{P}\left(\Sigma_{h},\mathcal{B},\mathcal{D}\right)$,

$\mathbb{P}\left(\mathcal{B},\mathcal{D}\right)=\int\mathbb{P}\left(\Sigma_{h},\mathcal{B},\mathcal{D}\right)\mathrm{d}\Sigma_{h},$
(5.1)

and we typically work on the log-likelihood
$\mathcal{L}\left(\mathcal{B},\mathcal{D}\right)=\log\mathbb{P}\left(\mathcal{B},\mathcal{D}\right)$
for convenience.

Due to the existence of the latent variable $\Sigma_{h}$, direct evaluation
and maximization of $\mathcal{L}\left(\mathcal{B},\mathcal{D}\right)$ is
intractable. Fortunately, the expectation-maximization (EM, Dempster, Laird
and Rubin, 1977) algorithm provides an effective way to tackle this problem.
Specifically, we treat the latent variable $\Sigma_{h}$ as missing data, and
the EM algorithm proceeds by utilizing the following two ingredients: (1) the
complete log-likelihood function
$\log\mathbb{P}\left(\Sigma_{h},\mathcal{B},\mathcal{D}\right)$; (2)
expectations with respect to the conditional distribution
$\mathbb{P}\left(\Sigma_{h}\mid\mathcal{B},\mathcal{D}\right)$. We first
derive the closed-form expressions for these two quantities in Sections 5.1
and 5.2, respectively, and then outline the EM algorithm in Section 5.3.

###  5.1 The complete likelihood

The complete likelihood admits the decomposition

$\mathbb{P}\left(\Sigma_{h},\mathcal{B},\mathcal{D}\right)=\mathbb{P}\left(\Sigma_{h}\right)\prod_{i=1}^{m}\mathbb{P}\left(B_{i}\mid\Sigma_{h}\right)\mathbb{P}\left(D_{i}\mid
B_{i},\Sigma_{h}\right).$ (5.2)

In Section 4 we have assumed that $\mathbb{P}\left(\Sigma_{h}\right)$ follows
a uniform distribution, and to show the expressions for the remaining two
parts, $\mathbb{P}\left(B_{i}\mid\Sigma_{h}\right)$ and
$\mathbb{P}\left(D_{i}\mid B_{i},\Sigma_{h}\right)$, it is necessary to
introduce a number of definitions and notations.

######  Definition 1.

Let $G=(V,E)$ be an undirected graph, where $V=\\{1,\dots,d\\}$ is the set of
nodes and $E$ is the set of undirected edges. Any subset of edges
$E^{\prime}\subset E$ can be used to define a _subgraph_
$G^{\prime}=(V,E^{\prime})$. A subgraph is called _complete_ if every pair of
nodes are adjacent. Define a _clique_ of $G$ as a maximal complete subgrah of
$G$. By maximal, we mean that a clique is not contained in a lager complete
subgraph. Define a _decompsoition_ of $G=(V,E)$ as a pair $(A,B)$ of subsets
of $V$ such that $V=A\cup B$ and $A\cap B$ is complete. An undirected graph is
_decomposable_ if it is complete, or if there exists a proper decomposition
$(A,B)$ into complete subgraphs $G^{A}$ and $G^{B}$ (Dawid and Lauritzen,
1993). In other words, a decomposable graph can be decomposed into several
cliques, and each clique is a complete subgraph of the graph.

Let $\widetilde{G}$ be a decomposable undirected graph consisting of $k$
cliques in the sequence $\left\\{C_{1},\ldots,C_{k}\right\\}$. Define the
_histories_ of the sequence as the sets $H_{j}=\bigcup_{i=1}^{j}C_{i}$ for
$j=1,\ldots,k$. Similarly, define the _separators_ of the sequence as
$S_{j}=H_{j-1}\cap C_{j}$ for $j=2,\ldots,k$. The merit of this decomposition
is that we can find the corresponding low-dimensional factorization for a
high-dimensional graph. For any two subsets $A$ and $B$, $S$ is said to
separate $A$ and $B$, denoted by $A\perp B\mid S$, if for every node pair
$a\in A$ and $b\in B$, each path connecting $a$ and $b$ includes at least one
element of $S$. A probability density $\mathbb{P}$ is said to be globally
Markov, if $X=(X_{A},X_{B},X_{D})\sim\mathbb{P}$ and $X_{A}\perp X_{B}\mid
X_{D}$ for any triple subsets $A$, $B$, $D$ such that $A\perp B\mid D$.
$X_{A}\perp X_{B}\mid X_{D}$ means that $X_{A}$ and $X_{B}$ are independent
conditional on $X_{D}$. Let $\mathbb{P}_{V}$ be a joint density that is
globally Markov with respect to a decomposable graph $\widetilde{G}=(V,E)$.
Let $\mathcal{C}$ and $\mathcal{S}$ be the cliques and separators of
$\widetilde{G}$, with corresponding marginal densities
$\left\\{\mathbb{P}_{C}:C\in\mathcal{C}\right\\}$ and
$\left\\{\mathbb{P}_{S}:S\in\mathcal{S}\right\\}$. Then $\mathbb{P}_{V}$ can
be factorized as (Lauritzen, 1996)

$\mathbb{P}_{V}\left(x_{V}\right)=\frac{\prod_{C\in\mathcal{C}}\mathbb{P}_{C}\left(x_{C}\right)}{\prod_{S\in\mathcal{S}}\mathbb{P}_{S}\left(x_{S}\right)}.$

A BN $B$ can be transformed into a unique undirected graph $G$, also known as
a moral graph. The undirected graph $G$ is generated from $B$, where if a node
has more than one parent node, edges between these parent nodes are inserted,
and all directed edges become undirected. As an illustration, Figure 2(a) and
(b) show the transformation from a BN to an undirected graph, which is
achieved by (1) adding an edge between the nodes $V_{B}$ and $V_{E}$, as both
$V_{B}$ and $V_{E}$ are the parent nodes of $V_{F}$, and (2) deleting all the
directions in $B$. Additionally, Figure 2(c) demonstrates that $G$ can be
decomposed into four cliques and three separators.

Figure 2: An example of graph tranformation and decomposition.

For Gaussian graphical models, the graph structure is usually related to the
inverse of the covariance matrix. Let $M$ be a $p\times p$ positive definite
matrix and $A$ is an index set with $A\subset\\{1,\dots,p\\}$. Then we use
$M_{AA}$ to denote the corresponding $|A|\times|A|$ submatrix of $M$, which is
also a positive definite matrix. We are now ready to present the first main
theorem based on the definitions above.

######  Theorem 5.1.

Suppose that $\tilde{G}_{i}$ is the decomposable undirected graph generated by
$B_{i}$, and $\mathcal{C}$ and $\mathcal{S}$ are the cliques and separators of
$\widetilde{G}_{i}$, respectively. Then
$\mathbb{P}\left(B_{i}\mid\Sigma_{h}\right)$ has the closed-form expression

$\mathbb{P}\left(B_{i}\mid\Sigma_{h}\right)=\frac{\prod_{C\in\mathcal{C}}\left|2^{-1}\Sigma_{hCC}\right|^{\left(\nu_{0}+|C|-1\right)/2}\Gamma_{|C|}^{-1}\left(\frac{\nu_{0}+|C|-1}{2}\right)}{\prod_{S\in\mathcal{S}}\left|2^{-1}\Sigma_{hSS}\right|^{\left(\nu_{0}+|S|-1\right)/2}\Gamma_{|S|}^{-1}\left(\frac{\nu_{0}+|S|-1}{2}\right)}\coloneqq
h(\Sigma_{h},\widetilde{G}_{i}),$ (5.3)

where $|C|$ and $|S|$ are the cardinality of the clique $C$ and the separator
$S$, respectively, $\Sigma_{hCC}$ and $\Sigma_{hSS}$ denote the submatrices of
$\Sigma_{h}$ according to the graph decomposation.

We make the following remark to further explain the implications of Theorem
5.1.

######  Remark 1.

$\mathbb{P}\left(B_{i}\mid\Sigma_{h}\right)$ can be expressed as
$\mathbb{P}\left(B_{i}\mid\Sigma_{h}\right)=\int\mathbb{P}\left(B_{i}\mid\Sigma_{i}\right)\mathbb{P}\left(\Sigma_{i}\mid\Sigma_{h}\right)\mathrm{d}\Sigma_{i}$.
As we defined in equation (4.1),
$\mathbb{P}\left(\Sigma_{i}\mid\Sigma_{h}\right)$ is the density function of
an inverse Wishart distribution $IW\left(\Sigma_{h},\nu_{0}\right)$.
Therefore, the most challenging part in deriving Theorem 5.1 is to measure the
density of a sparse BN $B_{i}$ given a covariance matrix $\Sigma_{i}$, _i.e._
, $\mathbb{P}\left(B_{i}\mid\Sigma_{i}\right)$. Fortunately, Dawid and
Lauritzen (1993) introduces the hyper inverse Wishart (HIW) distribution to
address this issue in the context of decomposable undirected graphs. Given a
decomposable undirected graph $\widetilde{G}_{i}$, the role of HIW
distribution is to limit the parameter space of an inverse Wishart
distribution according to this graph. In our case,
$\mathbb{P}\left(\widetilde{G}_{i}\mid\Sigma_{i}\right)\mathbb{P}\left(\Sigma_{i}\mid\Sigma_{h}\right)$
is an HIW distribution, denoted as
$HIW(\Sigma_{i},\Sigma_{h},\widetilde{G}_{i})$. Marginalizing out
$\Sigma_{i}$, we obtain the normalizing constant

$h(\Sigma_{h},\widetilde{G}_{i})=\int\mathbb{P}\left(\widetilde{G}_{i}\mid\Sigma_{i}\right)\mathbb{P}\left(\Sigma_{i}\mid\Sigma_{h}\right)\mathrm{d}\Sigma_{i}$

and have the identity
$\mathbb{P}\left(B_{i}\mid\Sigma_{h}\right)=\mathbb{P}\left(\widetilde{G}_{i}\mid\Sigma_{h}\right)=h(\Sigma_{h},\widetilde{G}_{i})$.

######  Theorem 5.2.

For the target structure $B_{i}$ corresponding to dataset $D_{i}$,

$\mathbb{P}\left(D_{i}\mid
B_{i},\Sigma_{h}\right)=\prod_{k=1}^{p}\frac{\mathbb{P}\left(D_{i}^{\left\\{X_{k}\right\\}\cup\Pi_{k}}\mid\Sigma_{h}\right)}{\mathbb{P}\left(D_{i}^{\Pi_{k}}\mid\Sigma_{h}\right)},$
(5.4)

where $\Pi_{k}$ are the parents of $X_{k}$,
$D_{i}^{\left\\{X_{k}\right\\}\cup\Pi_{k}}$ is the dataset $D_{i}$ restricted
to $X_{k}$ and $\Pi_{k}$, and $D_{i}^{\Pi_{k}}$ is the dataset $D_{i}$
restricted to $\Pi_{k}$. Each term in equation (5.4) is of the form

$\mathbb{P}\left(D_{i}^{C}\mid\Sigma_{h}\right)=\frac{\left|\Sigma_{hCC}\right|^{\nu_{0}/2}\Gamma_{p}\left(\frac{\nu_{0}+n_{i}}{2}\right)}{\pi^{n_{i}p/2}\left|\Sigma_{hCC}+n_{i}S_{iCC}\right|^{\left(\nu_{0}+n_{i}\right)/2}\Gamma_{p}\left(\frac{\nu_{0}}{2}\right)},$
(5.5)

where $\Sigma_{hCC}$ is the sub-matrix of $\Sigma_{h}$, and $S_{iCC}$ is the
sub-matrix of the sample covariance matrix $S_{i}$.

Finally, multiplying $\mathbb{P}\left(D_{i}\mid B_{i},\Sigma^{h}\right)$,
$\mathbb{P}\left(B_{i}\mid\Sigma_{h}\right)$, and
$\mathbb{P}\left(\Sigma_{h}\right)$, we obtain the complete likelihood
function $\mathbb{P}\left(\Sigma_{h},\mathcal{D},\mathcal{B}\right)$ in (5.2).

###  5.2 The conditional distribution

The EM algorithm requires computing the expectation of the complete log-
likelihood function with respect to the conditional distribution
$\mathbb{P}\left(\Sigma_{h}\mid\mathcal{B},\mathcal{D}\right)$. As a first
step, Theorem 5.3 gives its density function in a closed form.

######  Theorem 5.3.

Suppose that $\tilde{G}_{i}$ is the decomposable undirected graph generated by
$B_{i}$, and then we have

$\mathbb{P}\left(\Sigma_{h}\mid\mathcal{B},\mathcal{D}\right)=\prod_{i=1}^{m}\frac{\left|\Sigma_{h}\right|^{\nu_{0}/2}\cdot
h\left(\Sigma_{h}+n_{i}S_{i},\widetilde{G}_{i}\right)}{\left|\Sigma_{h}+n_{i}S_{i}\right|\left(\nu_{0}+n_{i}\right)/2},$

where the term $h(\cdot,\widetilde{G}_{i})$ is given in equation (5.3).

###  5.3 The EM algorithm

To maximize $\mathcal{L}\left(\mathcal{B},\mathcal{D}\right)$ with respect to
$\mathcal{B}$, the EM algorithm starts with an initial value
$\mathcal{B}^{(0)}$, and then repeats the following two steps until
convergence:

  * •

E-step: Define the function

$\mathcal{Q}\left(\mathcal{B}\mid\mathcal{B}^{(t-1)},\mathcal{D}\right)=\int\log\mathbb{P}\left(\Sigma_{h},\mathcal{B},\mathcal{D}\right)\mathbb{P}\left(\Sigma_{h}\mid\mathcal{B}^{(t-1)},\mathcal{D}\right)\mathrm{d}\Sigma_{h}.$
(5.6)

  * •

M-step: Solve the optimization problem below to generate the next iterate
$\mathcal{B}^{(t)}$:

$\mathcal{B}^{(t)}=\underset{\mathcal{B}}{\arg\max}\
\mathcal{Q}\left(\mathcal{B}\mid\mathcal{B}^{(t-1)},\mathcal{D}\right).$

The basic convergence properties of the EM algorithm were established by
Boyles (1983) and Wu (1983), and Friedman (2013) has proved that the algorithm
converges when there is no further improvement in the objective score.
However, in our setting, the $\mathcal{Q}$-function in (5.6) is still
intractable due to the high-dimensional integration. Instead, we view the
integral in (5.6) as an expectation, and then use Monte Carlo samples to
approximate it. This leads to the MCEM algorithm, which is a practical
computing method for the EM algorithm.

It should be pointed out that the convergence properties of MCEM is more
complicated than ordinary EM. Some efforts have been made to establish the
convergence properties of MCEM; see for example Chan and Ledolter (1995) and
Fort and Moulines (2003). We describe our implementation of the MCEM algorithm
for MTGBN in Section 6.

##  6 Computing method

###  6.1 MC E-step

The basic idea of MCEM is to replace the exact
$\mathcal{Q}\left(\mathcal{B}\mid\mathcal{B}^{(t-1)},\mathcal{D}\right)$
function in (5.6) by the following unbiased estimator:

$\widetilde{\mathcal{Q}}\left(\mathcal{B}\mid\mathcal{B}^{(t-1)},\mathcal{D}\right)=\frac{1}{N}\sum_{l=1}^{N}\log\mathbb{P}\left(\Sigma_{h}^{l},\mathcal{B},\mathcal{D}\right),$
(6.1)

where $\\{\Sigma_{h}^{l}\\}_{l=1}^{N}$ is a Monte Carlo sample of
$\mathbb{P}\left(\Sigma_{h}\mid\mathcal{B}^{(t-1)},\mathcal{D}\right)$.
Therefore, the key to the MC E-step is a proper sampling algorithm for
$\mathbb{P}\left(\Sigma_{h}\mid\mathcal{B}^{(t-1)},\mathcal{D}\right)$.

There are two main challenges that need to be addressed: (1) the positive-
definiteness of $\Sigma_{h}$ must be maintained; (2) the sampling scheme needs
to be fast and efficient. For the first challenge, we use the Cholesky
decomposition to transform $\Sigma_{h}$ into an unrestricted vector $V$, and
implement the sampling on the unrestricted space. For the second problem, we
employ the HMC sampler (Neal et al., 2011) for efficient sampling, which is a
powerful tool for Bayesian computation, especially in high-dimensional and
other complicated scenarios. It has been empirically justified that HMC can
offer greater computational efficiency than the traditional Gibbs sampling and
Metropolis–Hastings algorithm.

Specifically, we factorize $\Sigma_{h}$ using the Cholesky decomposition
$\Sigma_{h}=LL^{T}$, where $L$ is a lower-triangular matrix with positive
diagonal elements (_i.e._ , $l_{ii}>0$ for $i=1,2,\ldots,p$ and $l_{ij}=0$ for
$i<j$). To further remove the constraints on $l_{ii}$, we work on the
reparameterized variables333The similar method is also used by the Stan
software: https://mc-stan.org/docs/2_28/reference-manual/covariance-
matrices-1.html., $\log(l_{ii})$, and let $V=(v_{ij})$ be a lower-triangular
matrix such that $v_{ii}=\log(l_{ii})$ for $i=1,2,\ldots,p$ and
$v_{ij}=l_{ij}$ for $i>j$. In this way, $V$ collects all unrestricted
variables for $\Sigma_{h}$. HMC requires evaluating both
$\log\mathbb{P}(V|\mathcal{B}^{(t-1)},\mathcal{D})$ and
$\nabla_{V}\log\mathbb{P}(V|\mathcal{B}^{(t-1)},\mathcal{D})$, and we
summarize the results in Theorem 6.1.

######  Theorem 6.1.

Let $\mathbf{diag}\\{A\\}$ be the diagonal of a matrix $A$, and
$\mathbf{lower}\\{A\\}$ be the lower-triangular part of $A$ excluding the
diagonal. Denote by $I$ the identity matrix, and $I_{C}$ a matrix constructed
by selecting rows of $I$ indexed by $C$. Let $u\circ v$ be the elementwise
multiplication of two vectors $u$ and $v$. Then the density function of $V$
given $\mathcal{B}^{(t-1)}$ and $\mathcal{D}$ is

$\mathbb{P}(V|\mathcal{B}^{(t-1)},\mathcal{D})=\prod_{i=1}^{m}\frac{\left|\Sigma_{h}\right|^{\nu_{0}/2}\cdot
h\left(\Sigma_{h}+n_{i}S_{i},\mathcal{B}_{i}^{(t-1)}\right)}{\left|\Sigma_{h}+n_{i}S_{i}\right|\left(\nu_{0}+n_{i}\right)/2}\cdot
2^{p}\prod_{i=1}^{p}l_{ii}^{p-i+2},$ (6.2)

and the gradient of the log-density function is

$\displaystyle\mathbf{diag}\\{\nabla_{V}\log\mathbb{P}(V|\mathcal{B}^{(t-1)},\mathcal{D})\\}$
$\displaystyle=\mathbf{diag}\\{\widetilde{L}+q(m,\nu_{0},p)\cdot
I\\}\circ\mathbf{diag}\\{L\\},$ (6.3)
$\displaystyle\mathbf{lower}\\{\nabla_{V}\log\mathbb{P}(V|\mathcal{B}^{(t-1)},\mathcal{D})\\}$
$\displaystyle=\mathbf{lower}\\{\widetilde{L}\\},$

where

$\displaystyle\widetilde{L}$
$\displaystyle=g(\nu_{0}+n_{i},\Sigma_{h}+n_{i}S_{i},\mathcal{B}^{(t-1)})-f(\nu_{0}+n_{i},\Sigma_{h}+n_{i}S_{i}),$
$\displaystyle g(\nu,\Sigma,\mathcal{B})$
$\displaystyle=2\nu\left\\{\sum_{C\in\mathcal{C}}I_{C}^{\mathrm{T}}(\Sigma_{CC})^{-1}I_{C}L-\sum_{S\in\mathcal{S}}I_{S}^{\mathrm{T}}(\Sigma_{SS})^{-1}I_{S}L\right\\},$
$\displaystyle f(\nu,\Sigma)$ $\displaystyle=\nu\cdot\Sigma^{-1}L,$
$\displaystyle q(m,\nu,p)$ $\displaystyle=\frac{p(p+2m\nu+3)}{2}.$

Based on equations (6.2) and (6.3), we can use the HMC algorithm to obtain a
set of sample $V$ through iterations, denoted as $\\{V^{(l)}\\}_{l=1}^{N}$.
Then by transfroming $V$ back to $\Sigma_{h}$, we can obtain
$\\{\Sigma_{h}^{(l)}\\}_{l=1}^{N}$. Details of the derivation and the
algorithm are given in Appendix D.

###  6.2 MC M-step

Finally in the MC M-step, we maximize the approximate $\mathcal{Q}$ function
$\widetilde{\mathcal{Q}}\left(\mathcal{B}\mid\mathcal{B}^{(t-1)},\mathcal{D}\right)$
as in (6.1), and take $\mathcal{B}^{(t)}$ to be the maximizer. We use the max-
min hill-climbing algorithm (MMHC, Tsamardinos, Brown and Aliferis, 2006) to
search the BN structures, which combines ideas from local learning,
constraint-based methods, and search-and-score techniques. It has been
empirically justified that on average MMHC outperforms several prototypical
algorithms (Wang, 2018), such as PC (Kalisch and Bühlman, 2007), sparse
candidate (Friedman, Nachman and Pe’er, 2013), greedy equivalence search
(Ramsey et al., 2017), and greedy search (Chickering, 2002).

It is observed that $\\{B^{(t-1)}_{i}\\}_{i=1}^{m}$ are independent with each
other in the objective function
$\widetilde{\mathcal{Q}}\left(\mathcal{B}\mid\mathcal{B}^{(t-1)},\mathcal{D}\right)$.
Therefore, we can apply MMHC $m$ times on the $m$ independent sub-problems:

$B_{i}^{(t)}=\underset{B_{i}}{\arg\max}\
\sum_{l=1}^{N}\left\\{\log\mathbb{P}(D_{i}\mid
B_{i},\Sigma_{h}^{l})+\log\mathbb{P}(B_{i}\mid\Sigma_{h}^{l})\right\\}.$

The overall MCEM algorithm of MTGBN is thereby given in Algorithm 1.

Algorithm 1 MCEM algorithm of MTGBN

0: $\mathcal{D}$: Observational datasets; $\mathcal{B}^{(0)}$: Initial BN
structures; $\varepsilon$: Tolerance.

0: $\hat{\mathcal{B}}$: Estimated BN structures.

1: for $t=1,\ldots,T$ do

2: MC E-step:

3: $\\{\Sigma_{h}^{(l)}\\}_{l=1}^{N}$ $\leftarrow$ HMC ($\mathcal{D}$,
$\mathcal{B}^{(t-1)}$), see Appendix D for details.

4: Define function
$\widetilde{\mathcal{Q}}\left(\mathcal{B}\mid\mathcal{B}^{(t-1)},\mathcal{D}\right)$
by equation (6.1).

5: MC M-step:

6: for $i=1,\ldots,m$ do

7: $B_{i}^{(t)}$ $\leftarrow$ MMHC ($D_{i},B_{i}^{(t-1)}$,
$\\{\Sigma_{h}^{(l)}\\}_{l=1}^{N}$).

8: MMHC is a heuristic algorithm to detect BN structures (Tsamardinos, Brown
and Aliferis, 2006).

9: end for

10: Set $\mathcal{B}^{(t)}=\\{B_{i}^{(t)}\\}_{i=1}^{m}$.

11: break if
$\left|\widetilde{\mathcal{Q}}\left(\mathcal{B}^{(t)}\mid\mathcal{B}^{(t-1)},\mathcal{D}\right)-\widetilde{\mathcal{Q}}\left(\mathcal{B}^{(t-1)}\mid\mathcal{B}^{(t-1)},\mathcal{D}\right)\right|\leq\varepsilon.$

12: end for

13: return $\hat{\mathcal{B}}=\mathcal{B}^{(t)}$.

##  7 Simulations

In this section, we empirically evaluate the finite sample performance of the
proposed MTGBN method on synthetic data and benchmark network data and compare
it with several competing methods.

###  7.1 Simulation design of synthetic data

We simulate synthetic datasets with our data generation mechanism (see Figure
1). First, we random generate $\Sigma_{h}$ from the uniform distribution over
all $p\times p$ positive definite matrices with unit trace, using an approach
similar to Mittelbach, Matthiesen and Jorswieck (2012). Then given
$\Sigma_{h}$, we sample $m$ intermediate matrices
$\\{\tilde{\Sigma}_{l}\\}_{l=1}^{m}$ from $IW\left(\Sigma_{h},\nu_{0}\right)$,
so each $\tilde{\Sigma}_{l}$ is a $p\times p$ positive definite matrix.

The $m$ undirected graph structures $\left\\{G_{l}\right\\}_{l=1}^{m}$ are
defined as follows. For each $l=1,\dots,m$, let
$\tilde{\Omega}_{l}=(\tilde{\omega}_{ij,l})=\tilde{\Sigma}_{l}^{-1}\in\mathbb{R}^{p\times
p}$, $V=\\{1,\ldots,p\\}$, and define $E_{l}=\\{(i,j):i\neq
j,|\tilde{\omega}_{ij,l}|>t_{l}\\}$, where $t_{l}>0$ is a threshold to control
the density of $E_{l}$. Let $\Omega_{l}=(\omega_{ij,l})$ be a matrix such that
$\omega_{ii,l}=\tilde{\omega}_{ii,l}$, $\omega_{ij,l}=\tilde{\omega}_{ij,l}$
if $(i,j)\in E_{l}$, and $\omega_{ij,l}=0$ otherwise. We assume that
$\Omega_{l}$ is still positive definite, which can be achieved by repeatedly
sampling $\tilde{\Sigma}_{l}$ until the condition is met. Then the graphs are
given by $G_{l}=(V,E_{l})$.

Denote by $|E_{l}|$ the cardinality of $E_{l}$, so $d_{l}=|E_{l}|/p^{2}$
represents the density of $G_{l}$, and is dependent on the threshold $t_{l}$.
By adjusting $t_{l}$, we can obtain $G_{l}$ with different density levels. In
fact, the zero elements in $\Omega_{l}$ imply the conditional independence
structure of the variables in $G_{l}$ (Dawid and Lauritzen, 1993).

Finally, we generate $m$ independent datasets
$\left\\{D_{l}\right\\}_{l=1}^{m}$ by sampling each
$D_{l}=\left\\{\mathbf{x}_{l1},\ldots,\mathbf{x}_{ln_{l}}\right\\}$ from a
multivariate normal distribution $\mathcal{N}\left(0,\Sigma_{l}\right)$, where
$\Sigma_{l}=\Omega_{l}^{-1}$, and $n_{l}$ is the sample size for the $l$-th
dataset.

To evaluate and compare the performance of the MTGBN method, we also learn
each network independently with single-task learning (SIG for short), and
learn a common network structure for all contexts (AVG for short), _i.e._ , an
“average” network that treats all subjects as samples in one task. AVG was
previously used as a baseline in Oyen and Lane (2012). To compare the
performance of MTGBN, SIG, and AVG, especially on the differences between the
structures, we use edge error, edge precision, edge F-score, and edge distance
as the primary evaluation measures.

Let $\left\\{B_{output,l}\right\\}_{l=1}^{m}$ represent the $m$ BN structures
obtained from MTGBN, SIG, or AVG method, and
$\left\\{G_{output,l}\right\\}_{l=1}^{m}$ be the undirected graph structures
corresponding to $\left\\{B_{output,l}\right\\}_{l=1}^{m}$. In the experiment
of synthetic data, since the generated networks
$\left\\{G_{l}\right\\}_{l=1}^{m}$ are undirected graphs, we mainly evaluate
the network structures without considering the edge orientation, which we call
the adjacency measure. To compute the evaluation metrics, we first calculate
four basic statistics: true positives (TP), the number of edges that are both
in $G_{output,1}$ and $G_{1}$; false positives (FP), the number of edges that
are present in $G_{output,1}$ but not in $G_{1}$; false negatives (FN), the
number of edges that are present in $G_{1}$ but not in $G_{output,1}$; and
true negatives (TN), the number of vertex pairs that are neither edges in
$G_{1}$ nor in $G_{output,1}$. Then the edge adjacency error, edge adjacency
precision, edge adjacency F-score, and edge adjacency distance are defined as:

  * •

Error:

$\text{\emph{Error}}=\frac{FP+FN}{TP+TN+FP+FN}.$ (7.1)

  * •

Precision:

$\text{ \emph{Precision} }=\frac{TP}{TP+FP}.$

  * •

F-score:

$\text{ \emph{F-score}}=2\cdot\frac{\text{ \emph{Precision}
}\cdot\text{\emph{Recall} }}{\text{ \emph{Precision} }+\text{\emph{Recall}}},$

where $\text{ \emph{Recall}}=TP/(TP+FN)$.

  * •

Edge distance: Given two structures, this measure is the minimum number of
edges needed to convert the learned graph into the true one.

We compare the performance of MTGBN, SIG, and AVG on the following three
aspects.

Performance v.s. sample size: We fix the task number $m=10$, number of nodes
$p=15$, density level $d_{l}=0.3$, and compare the model performance of
network structure learning for different sample size
$n_{l}=50,100,150,200,250,300$.

Performance v.s. task number: We fix the number of nodes $p=15$, density level
$d_{l}=0.3$, sample size $n_{l}=250$, and consider different task number
$m=2,4,6,8,10$.

Performance v.s. density level: We fix the task number $m=10$, number of nodes
$p=15$, and consider a combination of different sample size
$n_{l}=100,200,300$ and density levels ranging from $0.2$ to $0.6$.

The evaluation metrics for one simulation run are averaged over the $m$
datasets, and we repeat each experiment $10$ times. The final evaluation
metrics are thereby averaged over the $10$ simulation runs.

###  7.2 Simulation design of benchmark network data

In this experiment, we select benchmark networks from the Bayesian Network
Repository (Constantinou et al., 2020). We choose three medium-sized benchmark
networks and generate low sample size data to validate our assumption that the
proposed MTGBN method can improve the understanding of the underlying
structure and generalize better even at small sample sizes. The information of
the selected benchmark networks is provided in Table 1.

Table 1: Facts of benchmark networks from the Bayesian Network Repository Benchmark Networks | Nodes | Arcs
---|---|---
ECOLI70 | 46 | 70
MAGIC-NIAB | 44 | 66
MAGIC-IRRI | 64 | 102

To obtain a set of related networks, we apply perturbations to the given
benchmark network to create similar but different networks. Denote by
$B_{truth}$ the benchmark BN with an adjacency matrix
$A_{t}=(a^{t}_{ij})\in\mathbb{R}^{p\times p}$, where $a^{t}_{ij}=1$ if there
is an edge from node $i$ to node $j$, and $a^{t}_{ij}=0$ otherwise,
$i,j=1,\ldots,p$. Then the generated BN, denoted by $B_{gen}$, is constructed
in the following way. Let $A_{g}$ be the adjacency matrix of $B_{gen}$, and
define the perturbation level as the ratio of the number of changed edges in
$A_{g}$ to the number of all possible edges. For example, if the perturbation
level is $30\%$, we first set $A_{g}=A_{t}$, and then randomly select $30\%$
of the off-diagonal elements $a^{t}_{ij}(i\neq j)$ in $A_{t}$. For each
selected element $a^{t}_{ij}$, if $a^{t}_{ij}=0$, then we set $a^{g}_{ij}=1$
(adding edge); if $a^{t}_{ij}=1$, then with probability $1/2$ we set
$a^{g}_{ij}=0$ (deleting edge), otherwise let $a^{g}_{ji}=1$ (reversing edge).
For each benchmark BN, we generate $m$ networks by such a perturbation scheme,
and in our experiment we set $m=10$. Each set of generated networks can be
viewed as a class of tasks sharing the same prior, since they are generated
from the same benchmark BN. Then we simulate data for each generated BN with a
sample size $n_{l}=250,l=1,\dots,m$.

We mainly compare the edge arrowhead error under different network
perturbation levels in the benchmark network experiment. The edge arrowhead
error is the same as equation (7.1), except that TP, FP, FN, and TN are
calculated considering the edge orientation. Same as Section 7.1, we repeat
the experiment 10 times and compute the averaged errors for MTGBN, SIG, and
AVG.

###  7.3 Simulation results

In the synthetic data experiment, we first evaluate the performance of the
proposed MTGBN method under different sample sizes, with the results shown in
Figure 3. It is clear that MTGBN significantly outperforms the other two
methods. The average edge adjacency error of MTGBN under different sample
sizes is in the range of 0.03 - 0.06, whereas for SIG, even using five times
the number of training data (sample size $=250$), its average edge adjacency
error is still larger than MTGBN (sample size $=50$). The same trend is
observed in the edge adjacency distance. For the other two metrics, MTGBN
consistently results in higher average edge adjacency precision and F-score
than AVG and SIG. Compared with SIG, the average edge adjacency precision of
MTGBN is also more stable under different sample sizes, and the best result is
achieved at the sample size of 200. Such findings show that the proposed MTGBN
method can produce satisfactory results using a small amount of data, thereby
demonstrating its usefulness in cases where data are expensive or hard to
collect.

(a)  (b)

(c)  (d)

Figure 3: Simulation results for synthetic data under different sample sizes.
(a) Average edge adjacency error (lower is better); (b) Average edge adjacency
distance (lower is better); (c) Average edge adjacency precision (higher is
better); (d) Average edge adjacency F-score (higher is better).

(a)  (b)

(c)  (d)

Figure 4: Simulation results for synthetic data under different task numbers.
(a) Average edge adjacency error (lower is better); (b) Average edge adjacency
distance (lower is better); (c) Average edge adjacency precision (higher is
better); (d) Average edge adjacency F-score (higher is better).

We also compare the performance of MTGBN, AVG, and SIG under different task
numbers $m$. Figure 4 indicates that MTGBN performs better than AVG and SIG
for all the test cases. As expected, SIG does not have better performance at a
larger value of $m$, since it only utilizes the information of a single task.
In contrast, MTGBN can use the information of related tasks, thus leading to
increased performance. We find that MTGBN can already get a good result when
$m=4$, and with the increase of $m$, the performance tends to be stable. These
results are also consistent with the existing literature (Oyen and Lane,
2012).

In addition to the sample size and task number, we find that the network
density also affects the model’s performance. We compare the average edge
adjacency error of MTGBN, AVG, and SIG under different sample sizes and
density levels, with results shown in Figure 5. It can be seen that with the
increase in sample size and decrease in density, the edge adjacency error of
the model gets smaller. Such a phenomenon is more obvious in MTGBN. Compared
with AVG and SIG, the proposed MTGBN method is especially advantageous with a
small sample size ($n_{l}=100$) and a medium density ($d_{l}=0.4$).

Figure 5: Average edge adjacency errors for synthetic data for different
sample sizes and density levels.

In the experiment of the benchmark network data, we compare the average
arrowhead errors for MTGBN, AVG, and SIG, and show the results in Table 2. We
consider four perturbation levels, $5\%$, $10\%$, $20\%$, and $30\%$,
representing the proportion of possible edges that are randomly modified
during simulation. Again, the comparison demonstrates that MTGBN is better
than SIG and AVG all the time. Even at a high perturbation level of $30\%$,
MTGBN still has decent error rates, whereas SIG and AVG do not provide good
results. We also find that the more related the tasks are, the better results
MTGBN can achieve.

Table 2: Average edge arrowhead errors of the learned BNs under different network perturbation levels. |  | 5% Perturbation | 10% Perturbation | 20% Perturbation | 30% Perturbation
---|---|---|---|---|---
| SIG | 0.048 | 0.079 | 0.085 | 0.098
ECOLI70 | AVG | 0.047 | 0.073 | 0.082 | 0.095
| MTGBN | 0.039 | 0.060 | 0.074 | 0.083
| SIG | 0.053 | 0.079 | 0.093 | 0.110
MAGIC-NIAB | AVG | 0.049 | 0.075 | 0.091 | 0.105
| MTGBN | 0.045 | 0.062 | 0.079 | 0.089
| SIG | 0.091 | 0.110 | 0.124 | 0.132
MAGIC-IRRI | AVG | 0.089 | 0.091 | 0.110 | 0.120
| MTGBN | 0.065 | 0.080 | 0.090 | 0.110

##  8 Applications to rs-fMRI data

Recent studies have found that MDD is closely related to the changes in the
brain connection network (Yamashita et al., 2020), especially the changes in
DMN. Several previous studies have reported abnormalities of the DMN in MDD
patients and found some altered functional connections (Greicius et al., 2007;
Bluhm et al., 2009; Sheline et al., 2010).

Such disorder-related DMN brain functional connectivity alteration can
elucidate neural mechanisms of MDD and further facilitate the diagnosis of
MDD. Therefore, we can use rs-fMRI data to identify DMN brain functional
connectivity patterns in MDD patients. This section will conduct such
experiments and apply the proposed MTGBN method to explore connectivity
alteration patterns in DMN that can facilitate MDD diagnosis.

Following the data acquisition and selection procedure introduced in Section
2, we analyze the rs-fMRI data of 15 MDD patients, which are regarded as 15
related tasks. Then we can obtain a brain connection network for each patient.
To study the differences in brain functional connectivity, the brain
connection networks of 15 HC participants, which are viewed as another set of
related tasks, are also obtained using the proposed MTGBN method. For
reference, we compare MTGBN with SIG and existing standard analysis methods of
rs-fMRI data: Granger causal analysis (GC, Seth, 2010) and dynamic causal
model (DCM, Marreiros, Kiebel and Friston, 2008).

We use MTGBN to estimate the connections between ROI nodes in three different
DMN network templates: aDMN, dDMN, and vDMN. Table 3 shows the difference in
connection degrees of $30$ ROI nodes (including in, out, and total connection
degree). Here, the connection degree represents the amount of functional
connections. The greater the connection degree, the higher the effective
connectivity. In aDMN, the total connection degrees are 292 for MDD and 232
for HC, indicating that MDD has 25.9% more effective connectivity than HC.
Similarly, In dDMN, MDD has 29.6% more effective connectivity than HC with
total connection degrees of 254 and 196, respectively. This increase in
effective connectivity in MDD patients is also consistent with the findings of
previous studies (Greicius et al., 2007; Lemogne et al., 2009; Belleau,
Taubitz and Larson, 2015). The rise of DMN connectivity is an apparent robust
biomarker of MDD (Zamoscik et al., 2014), which can also be found in our
results.

Table 3: Connection degrees of MDD patients and HC participants in three DMN network templates Node | Total degree | In | Out | Node | Total degree | In | Out
---|---|---|---|---|---|---|---
aDMN (MDD) |  |  |  | aDMN (HC) |  |  |
aDMN1 | 9 | 4 | 5 | aDMN1 | 2 | 2 | 0
aDMN2 | 34 | 13 | 21 | aDMN2 | 35 | 14 | 21
aDMN3 | 29 | 10 | 19 | aDMN3 | 7 | 3 | 4
aDMN4 | 42 | 18 | 24 | aDMN4 | 19 | 9 | 10
aDMN5 | 30 | 14 | 16 | aDMN5 | 25 | 9 | 16
aDMN6 | 12 | 10 | 2 | aDMN6 | 9 | 4 | 5
aDMN7 | 21 | 7 | 14 | aDMN7 | 24 | 13 | 11
aDMN8 | 33 | 19 | 14 | aDMN8 | 22 | 12 | 10
aDMN9 | 30 | 17 | 13 | aDMN9 | 29 | 14 | 15
aDMN10 | 31 | 21 | 10 | aDMN10 | 26 | 18 | 8
aDMN11 | 21 | 14 | 7 | aDMN11 | 34 | 18 | 16
Total | 292 |  |  | Total | 232 |  |
dDMN (MDD) |  |  |  | dDMN (HC) |  |  |
dDMN1 | 30 | 10 | 20 | dDMN1 | 22 | 2 | 20
dDMN2 | 20 | 11 | 9 | dDMN2 | 12 | 10 | 2
dDMN3 | 25 | 11 | 14 | dDMN3 | 16 | 5 | 11
dDMN4 | 28 | 15 | 13 | dDMN4 | 27 | 17 | 10
dDMN5 | 31 | 13 | 18 | dDMN5 | 18 | 12 | 6
dDMN6 | 24 | 13 | 11 | dDMN6 | 23 | 13 | 10
dDMN7 | 28 | 13 | 15 | dDMN7 | 15 | 9 | 6
dDMN8 | 44 | 24 | 20 | dDMN8 | 36 | 17 | 19
dDMN9 | 24 | 17 | 7 | dDMN9 | 27 | 13 | 14
Total | 254 |  |  | Total | 196 |  |
vDMN (MDD) |  |  |  | vDMN (HC) |  |  |
vDMN1 | 7 | 2 | 5 | vDMN1 | 2 | 0 | 2
vDMN2 | 29 | 11 | 18 | vDMN2 | 21 | 8 | 13
vDMN3 | 23 | 14 | 9 | vDMN3 | 28 | 11 | 17
vDMN4 | 26 | 13 | 13 | vDMN4 | 24 | 15 | 9
vDMN5 | 18 | 8 | 10 | vDMN5 | 20 | 11 | 9
vDMN6 | 26 | 10 | 16 | vDMN6 | 31 | 12 | 19
vDMN7 | 30 | 16 | 14 | vDMN7 | 28 | 12 | 16
vDMN8 | 19 | 9 | 10 | vDMN8 | 21 | 16 | 5
vDMN9 | 38 | 21 | 17 | vDMN9 | 39 | 19 | 20
vDMN10 | 22 | 15 | 7 | vDMN10 | 26 | 16 | 10
Total | 238 |  |  | Total | 240 |  |

From Table 3 we can find some ROI nodes with high connection degrees, such as
aDMN4 (temporal-parietal junction), aDMN2 (posterior cingulate cortex), aDMN8
(posterior inferior parietal lobule), and aDMN10 (parahippocampal cortex) in
aDMN. Specifically, aDMN4 and aDMN2 have more outgoing connections than
incoming connections, whereas aDMN8 and aDMN10 have more incoming connections
than outgoing connections. It indicates that aDMN4 and aDMN2 have a more
essential impact on other ROI nodes, and aDMN8 is mainly affected by other ROI
nodes. These findings are also consistent with the physiological discoveries
of MDD (Li et al., 2018).

Figure 6 depicts the network connection counts (the number of subjects with
the same connection) of MDD patients in three DMN networks. ROI connections
with high connection counts are generally considered universal connections
among MDD patients. Many such ROI connections have been noted as
neurophysiological diagnostic biomarkers of MDD. From Figure 6, we can find
some important functional connections with high connection counts: aDMN3
(dorsal medial prefrontal cortex) $\rightarrow$ aDMN4 (temporal-parietal
junction), aDMN4 (temporal-parietal junction) $\rightarrow$ aDMN5 (lateral
temporal cortex), aDMN3 (dorsal medial prefrontal cortex) $\rightarrow$ aDMN5
(lateral temporal cortex ), aDMN4 (temporal-parietal junction) $\rightarrow$
aDMN10 (parahippocampal cortex), and aDMN9 (retrosplenial cortex)
$\rightarrow$ aDMN2 (posterior cingulate cortex). These functional connections
are also consistent with the MDD imaging biomarkers identified in the existing
literature for distinguishing MDD from HC (Khundakar and Thomas, 2009;
Andrews-Hanna et al., 2010; Zhao et al., 2014).

Figure 6: The network connection counts in three DMN networks of MDD patients. Table 4: Two-sample proportion tests for published differential connections. Template | Connection | MTGBN | SIG | GC | DCM
---|---|---|---|---|---
aDMN | aDMN3 $\rightarrow$ aDMN4 | 0.0031 | 0.0339 | 0.0127 | 0.2280
aDMN4 $\rightarrow$ aDMN5 | 0.0155 | 0.2193 | 0.1367 | 0.3452
aDMN3 $\rightarrow$ aDMN5 | 0.0339 | 0.0716 | 0.0607 | 0.1345
aDMN4 $\rightarrow$ aDMN10 | 0.0708 | 0.2714 | 0.1279 | 0.2280
aDMN9 $\rightarrow$ aDMN2 | 0.0708 | 0.4721 | 0.3563 | 0.3575
dDMN | dDMN1 $\rightarrow$ dDMN7 | 0.0014 | 0.0680 | 0.2321 | 0.1345
dDMN2 $\rightarrow$ dDMN4 | 0.0127 | 0.1279 | 0.2321 | 0.1345
dDMN7 $\rightarrow$ dDMN8 | 0.0493 | 0.3452 | 0.2045 | 0.3563
dDMN6 $\rightarrow$ dDMN2 | 0.1345 | 0.6437 | 0.3524 | 0.6437
vDMN | vDMN2 $\rightarrow$ vDMN4 | 0.0127 | 0.1279 | 0.0327 | 0.3452
vDMN9 $\rightarrow$ vDMN4 | 0.0607 | 0.1367 | 0.1367 | 0.6548
vDMN7 $\rightarrow$ vDMN9 | 0.1345 | 0.6437 | 0.2280 | 0.7807
vDMN8 $\rightarrow$ vDMN3 | 0.1345 | 0.7807 | 0.3563 | 0.6437

To further test whether MTGBN is able to learn the difference between MDD
patients and HC participants, we do a two-sample proportion test for each
published differential connection, and report the $p$-values in Table 4. For
example, in aDMN there are five functional connections that are known to
effectively distinguish MDD from HC (the “Connection” column in Table 4), and
a small $p$-value indicates that MTGBN is more likely to identify such a
connection in the MDD group than in the HC group. For comparison, the testing
results for SIG, GC, and DCM are also included. To ensure the consistency of
measurement standards, we control these methods to learn the same density of
the network structure.

From Table 4 we can find that overall MTGBN produces the most significant
results, whereas in some cases, for example aDMN4 $\rightarrow$ aDMN5, all
other methods have insignificant outcomes even at the $\alpha=0.1$ level. This
finding validates that MTGBN is good at detecting correlation structures that
reflect the actual biological functions.

As an interpretation of the results in Table 4, we find that within the aDMN
template, MDD shows significantly increased connectivity for aDMN3
$\rightarrow$ aDMN4, aDMN4 $\rightarrow$ aDMN5, and aDMN3 $\rightarrow$ aDMN5
compared to HC. According to the definition of DMN subsystem (Andrews-Hanna et
al., 2010), these are within-system connections in the dorsal medial
prefrontal cortex (dMPFC) subsystem, which is considered to play more social
reflective roles (Andrews-Hanna et al., 2010; Andrews-Hanna, Smallwood and
Spreng, 2014), and is related to various cognitive functions, evaluation and
expression of negative emotions (Etkin, Egner and Kalisch, 2011). These
observations may reflect the resting-state hyperconnectivity in the dMPFC
subsystem in MDD patients. The neuropathology studies in MDD also report
abnormalities in the dMPFC subsystem (Khundakar and Thomas, 2009; Zhao et al.,
2014).

##  9 Conclusion

Learning Bayesian network structures has received lots of attention, but also
encounters significant challenges with insufficient data. This study proposes
a multitask learning method to estimate multiple networks simultaneously from
related observations, which effectively compensates for the small sample
sizes. We achieve this goal by introducing a shared prior covariance matrix
among all tasks and utilizing the hierarchical model structure for estimation.
Notably, we develop a tractable computing method for the structure learning of
networks, powered by analytic likelihood and gradient evaluations that lead to
efficient HMC sampling. Extensive and systematic numerical experiments
validate the excellent performance of the proposed method, and we highlight
its usefulness in analyzing rs-fMRI data and facilitating the study of the
major depressive disorder.

##  Appendix A Proof of Theorem 5.1

Here we provide the detailed proof of Theorem 5.1. The conditional
distribution $\mathbb{P}\left(B_{i}\mid\Sigma_{h}\right)$ can be expressed as
following formula

$\mathbb{P}\left(B_{i}\mid\Sigma_{h}\right)=\int\mathbb{P}\left(B_{i}\mid\Sigma_{i}\right)\mathbb{P}\left(\Sigma_{i}\mid\Sigma_{h}\right)d\Sigma_{i},$

where the conditional distribution
$\mathbb{P}\left(\Sigma_{i}\mid\Sigma_{h}\right)$ is the density function of
an inverse Wishart distribution $IW\left(\Sigma_{h},\nu_{0}\right)$.
Therefore, the most challenging part in deriving Theorem 5.1 is to measure the
density of a sparse BN $B_{i}$ given a covariance matrix $\Sigma_{i}$, _i.e._
, $\mathbb{P}\left(B_{i}\mid\Sigma_{i}\right)$.

Fortunately, Dawid and Lauritzen (1993) constructs the HIW distribution to
address this issue in the context of decomposable undirected graph. Given a
decomposable undirected graph $\widetilde{G}_{i}$ (consider
$\widetilde{G}_{i}$ that is the corresponding undirected graph of $B_{i}$),
the role of HIW distribution is to limit the parameter space of inverse
Wishart distribution according to this graph. In the context of undirected
graph, we can regard
$\mathbb{P}\left(\widetilde{G}_{i}\mid\Sigma_{i}\right)\mathbb{P}\left(\Sigma_{i}\mid\Sigma_{h}\right)$
as a HIW distribution (denote as
$HIW(\Sigma_{i},\Sigma_{h},\widetilde{G}_{i})$). The role of
$\mathbb{P}\left(\widetilde{G}_{i}\mid\Sigma_{i}\right)$ here is to limit the
parameter space of $\mathbb{P}\left(\Sigma_{i}\mid\Sigma_{h}\right)$ in a
subspace $\Psi_{i}$ in which $\Sigma_{i}$ is consistent with
$\widetilde{G}_{i}$. Therefore, the HIW distribution
$HIW(\Sigma_{i},\Sigma_{h},\widetilde{G}_{i})$ can be written as follows:

$HIW(\Sigma_{i},\Sigma_{h},\widetilde{G}_{i})=\mathbb{P}\left(\widetilde{G}_{i}\mid\Sigma_{i}\right)\mathbb{P}\left(\Sigma_{i}\mid\Sigma_{h}\right)=\mathbb{P}\left(\Sigma_{i}\mid\Sigma_{h}\right),\Sigma_{i}\in\Psi_{i},$

where $\Psi_{i}$ is the subspace of $\Sigma_{i}$. Let
$\Omega_{i}=\Sigma_{i}^{-1}$, $\Psi_{i}$ be the set of $p\times p$ positive
definite symmetric matrices with $\Omega_{ij}=0$ if and only if
$(i,j)\notin\widetilde{G}_{i}$. Therefore, we further derive the following
expression:

$\small\mathbb{P}\left(\widetilde{G}_{i}\mid\Sigma_{h}\right)=\int_{\Psi_{i}}\mathbb{P}\left(\Sigma_{i}\mid\Sigma_{h}\right)d\Sigma_{i},\quad\Sigma_{i}\in\Psi_{i}.$
(A.1)

In fact, equation (A.1) tells the truth that after marginalizing out
$\Sigma_{i}$ from $HIW(\Sigma_{i},\Sigma_{h},\widetilde{G}_{i})$, we can
obtain $\mathbb{P}\left(\widetilde{G}_{i}\mid\Sigma_{h}\right)$. Actually,
$\mathbb{P}\left(\widetilde{G}_{i}\mid\Sigma_{h}\right)$ is the normalizing
constant of $HIW(\Sigma_{i},\Sigma_{h},\widetilde{G}_{i})$, denote as
$h(\Sigma_{h},\widetilde{G}_{i})$, which can be derived to following equation:

$\mathbb{P}\left(\widetilde{G}_{i}\mid\Sigma^{h}\right)=h(\Sigma_{h},\widetilde{G}_{i})=\frac{\prod_{C\in\mathcal{C}}\left|2^{-1}\Sigma_{hCC}\right|^{\left(\nu_{0}+|C|-1\right)/2}\Gamma_{|C|}^{-1}\left(\frac{\nu_{0}+|C|-1}{2}\right)}{\prod_{S\in\mathcal{S}}\left|2^{-1}\Sigma_{hSS}\right|^{\left(\nu_{0}+|S|-1\right)/2}\Gamma_{|S|}^{-1}\left(\frac{\nu_{0}+|S|-1}{2}\right)},$
(A.2)

where $\Gamma(\cdot)$ denotes Gamma function, $\mathcal{C}$ and $\mathcal{S}$
are the cliques and separators of the decomposable undirected graph,
$\Sigma_{hCC}$ and $\Sigma_{hSS}$ denote the submatrix of $\Sigma_{h}$
according to graph decomposation, $|C|$ and $|S|$ denote the cardinality of
the clique $C$ and separator $S$ respectively.

Fortunately, a BN $B_{i}$ has an unique undirected graph $\widetilde{G}_{i}$
corresponding to it, thus we have the following relationship

$\mathbb{P}\left(B_{i}\mid\Sigma_{h}\right)=\mathbb{P}\left(\widetilde{G}_{i}\mid\Sigma_{h}\right)=h(\Sigma_{h},\widetilde{G}_{i}).$

##  Appendix B Proof of Theorem 5.2

Here we provide the detailed proof of Theorem 5.2. Let $\boldsymbol{X}$ be a
random vector of $p$ continuous variables $(X_{1},\ldots,X_{p})^{\top}$, and
$D=\left\\{\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\right\\}$ be $n$ observations
of $\boldsymbol{X}$. Assume $\boldsymbol{X}$ follows a multivariate normal
distribution, and $B$ is the GBN associated with $\boldsymbol{X}$. Then,
Geiger and Heckerman (1994) shows that the conditional distribution of $D$
given the network structure $B$ can be factorized as

$\mathbb{P}\left(D\mid
B\right)=\prod_{k=1}^{p}\frac{\mathbb{P}\left(D^{\left\\{X_{k}\right\\}\cup\Pi_{k}}\right)}{\mathbb{P}\left(D^{\Pi_{k}}\right)},$
(B.1)

where $D^{\left\\{X_{k}\right\\}\cup\Pi_{k}}$ is the dataset $D$ restricted to
the variable $X_{k}$ and its parents $\Pi_{k}$. Each term in equation (B.1)
follows a multivariate normal distribution.

Similar to equation (B.1), we can derive $\mathbb{P}\left(D_{i}\mid
B_{i},\Sigma_{h}\right)$ by encoding the BN structure $B_{i}$ into
$\mathbb{P}(D_{i}|\Sigma_{h})$, note that:

$\mathbb{P}\left(D_{i}\mid
B_{i},\Sigma_{h}\right)=\prod_{k=1}^{p}\frac{\mathbb{P}\left(D_{i}^{\left\\{X_{k}\right\\}\cup\Pi_{k}}\mid\Sigma_{h}\right)}{\mathbb{P}\left(D_{i}^{\Pi_{k}}\mid\Sigma_{h}\right)}.$

Therefore, the main fact we need to do is to derive
$\mathbb{P}(D_{i}|\Sigma_{h})$. From equation (4.1) and (4.2), we have

$\displaystyle\mathbb{P}(D_{i},\Sigma_{i}|\Sigma_{h})$
$\displaystyle=\mathbb{P}(D_{i}|\Sigma_{i})\mathbb{P}(\Sigma_{i}|\Sigma_{h})$
$\displaystyle=\frac{(2\pi)^{-n_{i}p/2}|\Sigma_{h}|^{\nu_{0}/2}}{2^{\nu_{0}p/2}\Gamma_{p}(\nu_{0}/2)}|\Sigma_{i}|^{-(\nu_{0}+p+n_{i}+1)/2}\exp\left\\{-\frac{1}{2}\mathrm{tr}\left[\left(\Sigma_{h}+n_{i}S_{i}\right)\Sigma_{i}^{-1}\right]\right\\}.$
(B.2)

Moreover, since $\mathbb{P}(D_{i}|\Sigma_{i}){\sim}\mathcal{N}(0,\Sigma_{i})$
and $\mathbb{P}(\Sigma_{i}|\Sigma_{h})\sim IW(\Sigma_{h},\nu_{0})$. Due to the
conjugacy, we have $\mathbb{P}(\Sigma_{i}|D_{i},\Sigma_{h})\sim
IW(\Sigma_{h}+n_{i}S_{i},\nu_{0}+n_{i})$, which can be written as

$\mathbb{P}(\Sigma_{i}|D_{i},\Sigma_{h})=\frac{|\Sigma_{h}+n_{i}S_{i}|^{(\nu_{0}+n_{i})/2}}{2^{(\nu_{0}+n_{i})p/2}\Gamma_{p}\left(\frac{\nu_{0}+n_{i}}{2}\right)}|\Sigma_{i}|^{-(\nu_{0}+p+n_{i}+1)/2}\exp\left\\{-\frac{1}{2}\mathrm{tr}\left[\left(\Sigma_{h}+n_{i}S_{i}\right)\Sigma_{i}^{-1}\right]\right\\}.$

and as a result, we obtain

$\mathbb{P}(D_{i}|\Sigma_{h})=\frac{\mathbb{P}(D_{i},\Sigma_{i}|\Sigma_{h})}{\mathbb{P}(\Sigma_{i}|D_{i},\Sigma_{h})}=\frac{|\Sigma_{h}|^{\nu_{0}/2}\Gamma_{p}\left(\frac{\nu_{0}+n_{i}}{2}\right)}{\pi^{n_{i}p/2}|\Sigma_{h}+n_{i}S_{i}|^{(\nu_{0}+n_{i})/2}\Gamma_{p}\left(\frac{\nu_{0}}{2}\right)}.$

##  Appendix C Proof of Theorem 5.3

In this section we provide the detailed proof of Theorem 5.3. To derive
$\mathbb{P}(\Sigma_{h}|\mathcal{B},\mathcal{D})$, where
$\mathcal{B}=\left\\{B_{i}\right\\}_{i=1}^{m}$, note that

$\mathbb{P}(\Sigma_{h}|\mathcal{B},\mathcal{D})\propto\mathbb{P}(\Sigma_{h},\mathcal{B},\mathcal{D})\propto\mathbb{P}(\mathcal{B},\mathcal{D}|\Sigma_{h})\mathbb{P}(\Sigma_{h})\propto\mathbb{P}(\mathcal{B},\mathcal{D}|\Sigma_{h}).$

By conditioning on $\\{\Sigma_{i}\\}_{i=1}^{m}$, we obtain

$\displaystyle\mathbb{P}(\mathcal{B},\mathcal{D}|\Sigma_{h})$
$\displaystyle=\int\mathbb{P}(\mathcal{B}|\\{\Sigma_{i}\\}_{i=1}^{m})\mathbb{P}(\mathcal{D}|\\{\Sigma_{i}\\}_{i=1}^{m})\mathbb{P}(\\{\Sigma_{i}\\}_{i=1}^{m}|\Sigma_{h})d\Sigma_{1}\cdots
d\Sigma_{m}$
$\displaystyle=\prod_{i=1}^{m}\int\mathbb{P}(B_{i}|\Sigma_{i})\mathbb{P}(D_{i}|\Sigma_{i})\mathbb{P}(\Sigma_{i}|\Sigma_{h})d\Sigma_{i}.$
(C.1)

The challenging of derivation in equation (C.1) is how to integrate
$\Sigma_{i}$ and generate $B_{i}$ and $D_{i}$ from $\Sigma_{h}$. Fortunately,
Theorem 5.1 provides us a bridge so that we can derive
$\mathbb{P}(\mathcal{B},\mathcal{D}|\Sigma_{h})$ without actually integrating
$\Sigma_{i}$. Thus, we can re-write equation (C.1) as :

$\mathbb{P}(\mathcal{B},\mathcal{D}|\Sigma_{h})=\prod_{i=1}^{m}\int_{\Psi_{i}}\mathbb{P}(D_{i}|\Sigma_{i})\mathbb{P}(\Sigma_{i}|\Sigma_{h})d\Sigma_{i},\quad\Sigma_{i}\in\Psi_{i}.$

In equation (B.2), we have already shown that

$\mathbb{P}(D_{i}|\Sigma_{i})\mathbb{P}(\Sigma_{i}|\Sigma_{h})=\frac{(2\pi)^{-n_{i}p/2}|\Sigma_{h}|^{\nu_{0}/2}}{2^{\nu_{0}p/2}\Gamma_{p}(\nu_{0}/2)}|\Sigma_{i}|^{-(\nu_{0}+p+n_{i}+1)/2}\exp\left\\{-\frac{1}{2}\mathrm{tr}\left[\left(\Sigma_{h}+n_{i}S_{i}\right)\Sigma_{i}^{-1}\right]\right\\}.$

So we need to compute the integral

$\int_{\Psi_{i}}|\Sigma_{i}|^{-(\nu_{0}+p+n_{i}+1)/2}\exp\left\\{-\frac{1}{2}\mathrm{tr}\left[\left(\Sigma_{h}+n_{i}S_{i}\right)\Sigma_{i}^{-1}\right]\right\\}d\Sigma_{i}.$
(C.2)

Recall that $\mathbb{P}(\Sigma_{i}|D_{i},\Sigma_{h})\sim
IW(\Sigma_{h}+n_{i}S_{i},\nu_{0}+n_{i})$, so

$\mathbb{P}(\Sigma_{i}|D_{i},\Sigma_{h})=\frac{|\Sigma_{h}+n_{i}S_{i}|^{(\nu_{0}+n_{i})/2}}{2^{(\nu_{0}+n_{i})p/2}\Gamma_{p}\left(\frac{\nu_{0}+n_{i}}{2}\right)}|\Sigma_{i}|^{-(\nu_{0}+p+n_{i}+1)/2}\exp\left\\{-\frac{1}{2}\mathrm{tr}\left[\left(\Sigma_{h}+n_{i}S_{i}\right)\Sigma_{i}^{-1}\right]\right\\},$

which is proportional to the integral in equation (C.2). Therefore,

$\displaystyle\int_{\Psi_{i}}|\Sigma_{i}|^{-(\nu_{0}+p+n_{i}+1)/2}\exp\left\\{-\frac{1}{2}\mathrm{tr}\left[\left(\Sigma_{h}+n_{i}S_{i}\right)\Sigma_{i}^{-1}\right]\right\\}d\Sigma_{i}$
$\displaystyle=$
$\displaystyle\frac{2^{(\nu_{0}+n_{i})p/2}\Gamma_{p}\left(\frac{\nu_{0}+n_{i}}{2}\right)}{|\Sigma_{h}+n_{i}S_{i}|^{(\nu_{0}+n_{i})/2}}\int_{\Psi_{i}}\mathbb{P}(\Sigma_{i}|\mathcal{D}_{i},\Sigma_{h})d\Sigma_{i}.$
(C.3)

Comparing $\mathbb{P}(\Sigma_{i}|\Sigma_{h})$ in the equation (A.1) and
$\mathbb{P}(\Sigma_{i}|\mathcal{D}_{i},\Sigma_{h})d\Sigma_{i}$ in equation
(C.3), we find that they only differ in the parameters of the inverse Wishart
distribution. Then we conclude that
$\int_{\Psi_{i}}\mathbb{P}(\Sigma_{i}|D_{i},\Sigma_{h})d\Sigma_{i}$ has a
similar closed-form expression to (A.2), given by

$\int_{\Psi_{i}}\mathbb{P}(\Sigma_{i}|D_{i},\Sigma_{h})\mathrm{d}\Sigma_{i}=h(\Sigma_{h}+n_{i}S_{i},\widetilde{G}_{i}).$

As a result,

$\displaystyle\int_{\Psi_{i}}|\Sigma_{i}|^{-(\nu_{0}+p+n_{i}+1)/2}\exp\left\\{-\frac{1}{2}\mathrm{tr}\left[\left(\Sigma_{h}+n_{i}S_{i}\right)\Sigma_{i}^{-1}\right]\right\\}d\Sigma_{i}$
$\displaystyle=$
$\displaystyle\frac{2^{(\nu_{0}+n_{i})p/2}\Gamma_{p}\left(\frac{\nu_{0}+n_{i}}{2}\right)}{|\Sigma_{h}+n_{i}S_{i}|^{(\nu_{0}+n_{i})/2}}\int_{\Psi_{i}}\mathbb{P}(\Sigma_{i}|\mathcal{D}_{i},\Sigma_{h},\nu_{0})d\Sigma_{i},$

and

$\mathbb{P}(\Sigma_{h}|\mathcal{B},\mathcal{D})=\prod_{i=1}^{m}\frac{|\Sigma_{h}|^{v_{0}/2}\cdot
h(\Sigma_{h}+n_{i}S_{i},\widetilde{G}_{i})}{|\Sigma_{h}+n_{i}S_{i}|^{(\nu_{0}+n_{i})/2}}.$

##  Appendix D Details of HMC Algorithm

In this section we continue the discussion on the HMC algorithm in Step 5 of
Algorithm 1. The trickiest thing about sampling $\Sigma_{h}$ from
$\mathbb{P}\left(\Sigma_{h}\mid\mathcal{B},\mathcal{D}\right)$ is that
$\Sigma_{h}$ is constrained to the positive semidefinite cone. Next, we
discuss how to transform $\Sigma_{h}$ to a unconstrained lower-triangular
matrix $V$ (Transform) and inverse-transform from $V$ to $\Sigma_{h}$
(Inverse-Transform).

Transform: For the positive semidefinite matrix $\Sigma_{h}$, there is a
unique lower-triangular matrix $L=chol(\Sigma_{h})$ with positive diagonal
entries, call a Cholesky factor, such that:

$\Sigma_{h}=LL^{\mathrm{T}}.$

The off-diagonal entries of the Cholesky factor $L$ are unconstrained, but the
diagonal entries $l_{ii}$ must be positive. To complete the transform, the
diagonal is log-transformed to produce a fully unconstrained lower-triangular
matrix $V$ defined by

$v_{ij}=\left\\{\begin{array}[]{cl}0&\text{ if }i<j,\\\ \log l_{ii}&\text{ if
}i=j,\\\ l_{ij}&\text{ if }i>j.\end{array}\right.$

Inverse-Transform: The inverse-transform reverses the two steps of the
transform. Given an unconstrained lower-triangular matrix $V$, the first step
is to recover the intermediate matrix $L$ by reversing the log trandfrom,

$l_{ij}=\left\\{\begin{array}[]{cl}0&\text{ if }i<j,\\\ \exp v_{ii}&\text{ if
}i=j,\\\ v_{ij}&\text{ if }i>j.\end{array}\right.$

In Section 6, we derived the expression of
$\mathbb{P}\left(V\mid\mathcal{B},\mathcal{D}\right)$. The core of HMC
algorithm is the calculation of logarithmic density function and its gradient.
Thus, next we will compute the logarithm of
$\mathbb{P}\left(V\mid\mathcal{B},\mathcal{D}\right)$ and its gradient. The
logarithm of $\mathbb{P}\left(V\mid\mathcal{B},\mathcal{D}\right)$, which
denotes as $\log\\{\mathbb{P}(V|\mathcal{B},\mathcal{D})\\}$ that can be
expressed as follows:

$\displaystyle\log\\{\mathbb{P}(V|\mathcal{B},\mathcal{D})\\}$
$\displaystyle=\sum_{i=1}^{m}\log
h(\Sigma_{h}+n_{i}S_{i},B_{i})-\sum_{i=1}^{m}\frac{\nu_{0}+n_{i}}{2}\log|\Sigma_{h}+n_{i}S_{i}|$
$\displaystyle\quad+\sum_{i=1}^{p}(m\nu_{0}+p-i+2)v_{ii}+C.$ (D.1)

After deriving $\log\\{\mathbb{P}(V|\mathcal{B},\mathcal{D})\\}$ we also need
to derive the gradient of $\log\\{\mathbb{P}(V|\mathcal{B},\mathcal{D})\\}$
with respect to $V$, which denote as
$\nabla_{V}\log\\{\mathbb{P}(V|\mathcal{B},\mathcal{D})\\}$. The expression of
$\log\\{\mathbb{P}(V|\mathcal{B},\mathcal{D})\\}$ mainly includes three parts,
we will compute the gradient for each term separately.

The first term consists of functions of the form
$g(\Sigma_{h})=\log|(\Sigma_{h}+A)_{CC}|$, where $C$ is an index set. Let
$I_{C}$ be a matrix constructed by selecting rows of the identity matrix
indexed by $C$, and then
$\Sigma_{hCC}=I_{C}\Sigma_{h}I_{C}^{\mathrm{T}}=I_{C}LL^{\mathrm{T}}I_{C}^{\mathrm{T}}$.
It is well known that $\partial
g/\partial\Sigma_{hCC}=\\{(\Sigma_{h}+A)_{CC}^{-1}\\}^{\mathrm{T}}=(\Sigma_{h}+A)_{CC}^{-1}$,
so we can write
$\mathrm{d}g=\mathrm{tr}\\{(\Sigma_{h}+A)_{CC}^{-1}(\mathrm{d}\Sigma_{hCC})\\}$.
Since

$\mathrm{d}\Sigma_{hCC}=I_{C}(\mathrm{d}\Sigma_{h})I_{C}^{\mathrm{T}}=I_{C}L(\mathrm{d}L)^{\mathrm{T}}I_{C}^{\mathrm{T}}+I_{C}(\mathrm{d}L)L^{\mathrm{T}}I_{C}^{\mathrm{T}}.$

Then we have

$\displaystyle\mathrm{d}g$
$\displaystyle=\mathrm{tr}\\{(\mathrm{d}L)^{\mathrm{T}}I_{C}^{\mathrm{T}}(\Sigma_{h}+A)_{CC}^{-1}I_{C}L\\}+\mathrm{tr}\\{L^{\mathrm{T}}I_{C}^{\mathrm{T}}(\Sigma^{h}+A)_{CC}^{-1}I_{C}(\mathrm{d}L)\\}$
$\displaystyle=2\mathrm{tr}\\{L^{\mathrm{T}}I_{C}^{\mathrm{T}}(\Sigma^{h}+A)_{CC}^{-1}I_{C}(\mathrm{d}L)\\}.$

Recall that $L$ is lower-triangular, so we get $\partial g/\partial
L=2\mathbf{lower}\\{I_{C}^{\mathrm{T}}(\Sigma_{h}+A)_{CC}^{-1}I_{C}L\\}$,
where $\mathbf{lower}\\{A\\}$ stands for the lower-triangular part of a matrix
$A$. Finally, it is easy to see that

$\frac{\partial g}{\partial v_{ij}}=\begin{cases}\frac{\partial g}{\partial
l_{ii}}\cdot l_{ii},&i=j,\\\ \frac{\partial g}{\partial
l_{ij}},&i>j.\end{cases}$

The second term involves $f(\Sigma_{h})=\log|\Sigma_{h}+A|$, where $A$ is a
positive definite matrix. Similar to $\partial g/\partial\Sigma_{hCC}$, we
have $\partial f/\partial\Sigma_{h}=(\Sigma_{h}+A)^{-1}$, and
$\mathrm{d}f=\mathrm{tr}\\{(\Sigma_{h}+A)^{-1}(\mathrm{d}\Sigma_{h})\\}$.
Again using the relation
$\mathrm{d}\Sigma_{h}=L(\mathrm{d}L)^{\mathrm{T}}+(\mathrm{d}L)L^{\mathrm{T}}$,
we obtain

$\displaystyle\mathrm{d}f$
$\displaystyle=\mathrm{tr}\\{(\Sigma_{h}+A)^{-1}(L(\mathrm{d}L)^{\mathrm{T}}+(\mathrm{d}L)L^{\mathrm{T}})\\}$
$\displaystyle=\mathrm{tr}\\{(\mathrm{d}L)^{\mathrm{T}}(\Sigma_{h}+A)^{-1}L\\}+\mathrm{tr}\\{L^{\mathrm{T}}(\Sigma_{h}+A)^{-1}(\mathrm{d}L)\\}$
$\displaystyle=2\mathrm{tr}\\{L^{\mathrm{T}}(\Sigma_{h}+A)^{-1}(\mathrm{d}L)\\}.$

So $\partial f/\partial L=2\mathbf{lower}\\{(\Sigma_{h}+A)^{-1}L\\}$. Similar
to $\partial g/\partial v_{ij}$ we have of a matrix $A$. Finally, it is easy
to see that

$\frac{\partial f}{\partial v_{ij}}=\begin{cases}\frac{\partial f}{\partial
l_{ii}}\cdot l_{ii},&i=j,\\\ \frac{\partial f}{\partial
l_{ij}},&i>j.\end{cases}$

Thus $\nabla_{V}\log\mathbb{P}(V|\mathcal{B},\mathcal{D})$ can be expressed as
follows:

$\displaystyle\mathbf{diag}\\{\nabla_{V}\log\mathbb{P}(V|\mathcal{B},\mathcal{D})\\}$
$\displaystyle=\mathbf{diag}\\{\widetilde{L}+q(m,\nu_{0},p)\cdot
I\\}\circ\mathbf{diag}\\{L\\},$ (D.2)
$\displaystyle\mathbf{lower}\\{\nabla_{V}\log\mathbb{P}(V|\mathcal{B},\mathcal{D})\\}$
$\displaystyle=\mathbf{lower}\\{\widetilde{L}\\},$

where

$\displaystyle\widetilde{L}$
$\displaystyle=g(\nu_{0}+n_{i},\Sigma_{h}+n_{i}S_{i},\mathcal{B})-f(\nu_{0}+n_{i},\Sigma_{h}+n_{i}S_{i}),$
$\displaystyle g(\nu,\Sigma,\mathcal{B})$
$\displaystyle=2\nu\left\\{\sum_{C\in\mathcal{C}}I_{C}^{\mathrm{T}}(\Sigma_{CC})^{-1}I_{C}L-\sum_{S\in\mathcal{S}}I_{S}^{\mathrm{T}}(\Sigma_{SS})^{-1}I_{S}L\right\\},$
$\displaystyle f(\nu,\Sigma)$ $\displaystyle=\nu\cdot\Sigma^{-1}L,$
$\displaystyle q(m,\nu,p)$ $\displaystyle=\frac{p(p+2m\nu+3)}{2}.$

According to the above definition, the specific solution process and iterative
steps of HMC algorithm are shown in Algorithm 2.

Algorithm 2 HMC for Sampling $\Sigma_{h}$

0: Observational datasets $\mathcal{D}$, BN structures $\mathcal{B}$ and
initial $\Sigma_{h}^{(0)}$.

0: $\\{\Sigma_{h}^{(l)}\\}_{l=1}^{N}$: $l$ samples of $\Sigma_{h}$.

1: Initialization: The number of iterations $N$, the number of steps $K$, the
step length $\epsilon$, and a specified covariance matrix to simulate momentum
$M$.

2: $V^{(0)}$ $\leftarrow$ Transform ($\Sigma_{h}^{(0)}$).

3: $\log\\{\mathbb{P}(V|\mathcal{B},\mathcal{D})\\}$ = equation (D.1).

4: Calculate $\log\\{\mathbb{P}(V^{(0)}|\mathcal{B},\mathcal{D})\\}$.

5: for $t=1,\ldots,N$ do

6: $\mathbf{p}\leftarrow\mathcal{N}(0,\mathbf{M})$,

7: $V^{(t)}\leftarrow V^{(t-1)},\tilde{V}\leftarrow
V^{(t-1)},\tilde{\mathbf{p}}\leftarrow\mathbf{p}$,

8: for $i=1,\ldots,K$ do

9: $\tilde{V},\tilde{\mathbf{p}}\leftarrow$ Leapfrog
$(\tilde{V},\tilde{\mathbf{p}},\epsilon,\mathbf{M})$.

10: end for

11: $\alpha\leftarrow\min\left(1,\frac{\exp\left(\log
f(\tilde{V})-\frac{1}{2}\tilde{\mathbf{p}}^{T}\mathbf{M}^{-1}\tilde{\mathbf{p}}\right)}{\exp\left(\log
f\left(\tilde{V}^{(t-1)}\right)-\frac{1}{2}\mathbf{p}^{T}\mathbf{M}^{-1}\mathbf{p}\right)}\right)$,

12: With probability $\alpha,V^{(t)}\leftarrow\tilde{V}\text{ and
}\mathbf{p}^{(t)}\leftarrow-\tilde{\mathbf{p}}$,

13: $\Sigma_{h}^{(t)}$ $\leftarrow$ Inverse-Transform ($V^{(t)}$).

14: end for

15: return $\Sigma_{h}^{(1)},\ldots,\Sigma_{h}^{(N)}$.

16: function Leapfrog $\left(V^{*},\mathbf{p}^{*},\epsilon,\mathbf{M}\right)$

17: $\tilde{\mathbf{p}}\leftarrow\mathbf{p}^{*}+(\epsilon/2)\nabla_{V}\log
f\left(V^{*}\right)$,

18: $\tilde{V}\leftarrow V^{*}+\epsilon\mathbf{M}^{-1}\tilde{\mathbf{p}}$,

19: $\tilde{\mathbf{p}}\leftarrow\tilde{\mathbf{p}}+(\epsilon/2)\nabla_{V}\log
f(\tilde{V})$,

20: return $\tilde{V},\tilde{\mathbf{p}}$.

21: end function

##  Appendix E ROI information of three DMN networks

See Table 5.

Table 5: Description of ROI nodes in the three DMN networks and their Montreal Neurological Institute (MNI) coordinates. Region name | Abbreviation | Broadman Areas | x | y | z |
---|---|---|---|---|---|---
aDMN (Andrews-Hanna Default Mode Network) |  |  |  |  |  |
P (PCC-aMPFC Core) |  |  |  |  |  |
Anterior medial prefrontal cortex | aDMN1 | 10, 32 | -6 | 52 | -2 |
Posterior cingulate cortex | aDMN2 | 23, 31 | -8 | -56 | 26 |
D (dMPFC Subsystem) |  |  |  |  |  |
Dorsal medial prefrontal cortex | aDMN3 | 9, 32 | 0 | 52 | 26 |
Temporal parietal junction | aDMN4 | 40, 39 | -54 | -54 | 28 |
Lateral temporal cortex | aDMN5 | 21, 22 | -60 | -24 | -18 |
Temporal pole | aDMN6 | 21 | -50 | 14 | -40 |
M (MTL Subsystem) |  |  |  |  |  |
Ventral medial prefrontal cortex | aDMN7 | 11, 24, 25, 32 | 0 | 26 | -18 |
Posterior inferior parietal lobule | aDMN8 | 39 | -44 | -74 | 32 |
Retrosplenial cortex | aDMN9 | 29, 30, 19 | -14 | -52 | 8 |
Parahippocampal cortex | aDMN10 | 20, 36, 19 | -28 | -40 | -12 |
Hippocampal formation | aDMN11 | 20, 36 | -22 | -20 | -26 |
dDMN (Dorsal Default Mode Network) |  |  |  |  |  |
Medical anterior cingulate cortex | dDMN1 | 9, 10, 24, 32, 11 | -3 | 49 | 14 |
Left angular gyrus | dDMN2 | 39 | -48 | -68 | 35 |
Right superior frontal gyrus | dDMN3 | 9 | 19 | 38 | 47 |
Posterior cingulate cortex, precuneus | dDMN4 | 23, 30 | 1 | -53 | 28 |
Middle cingulate cortex | dDMN5 | 23 | 2 | -15 | 36 |
Right angular gyrus | dDMN6 | 29 | 56 | -66 | 36 |
Left and right thalamus | dDMN7 | N/A | 6 | -6 | 2 |
Left hippocampus | dDMN8 | 20, 36, 30 | -24 | -29 | -13 |
Right hippocampus | dDMN9 | 20, 36, 30 | 27 | -23 | -17 |
vDMN (Ventral Default Mode Network) |  |  |  |  |  |
Left posterior cingulate cortex | vDMN1 | 29, 30, 23 | -12 | -58 | 15 |
Left middle frontal gyrus | vDMN2 | 8, 6 | -24 | 12 | 55 |
Left parahippocampal gyrus | vDMN3 | 37, 20 | -28 | -37 | -15 |
Left middle occipital gyrus | vDMN4 | 19, 39 | -44 | -74 | 30 |
Right posterior cingulate cortex | vDMN5 | 30, 23 | 13 | -53 | 14 |
Precuneus | vDMN6 | 7, 5 | 1 | -57 | 54 |
Right superior frontal gyrus, | vDMN7 | 9, 8 | 26 | 26 | 45 |
Right parahippocampal gyrus | vDMN8 | 37, 30 | 28 | -33 | -19 |
Right angular gyrus | vDMN9 | 39, 19 | 43 | -74 | 32 |
Right lobule | vDMN10 | N/A | 8 | -26 | 70 |

## References

  * Altendorf, Restificar and Dietterich (2005) [author] Altendorf, Eric EE. E., Restificar, Angelo CA. C. and Dietterich, Thomas GT. G. (2005). Learning from sparse data by exploiting monotonicity constraints. In Proceedings of the Twenty-First Conference on Uncertainty in Artificial Intelligence 18–26.
  * Andrews-Hanna, Smallwood and Spreng (2014) [author] Andrews-Hanna, Jessica RJ. R., Smallwood, JonathanJ. and Spreng, R NathanR. N. (2014). The default network and self-generated thought: component processes, dynamic control, and clinical relevance. Annals of the New York Academy of Sciences 1316 29–52.
  * Andrews-Hanna et al. (2010) [author] Andrews-Hanna, Jessica RJ. R., Reidler, Jay SJ. S., Sepulcre, JorgeJ., Poulin, ReneeR. and Buckner, Randy LR. L. (2010). Functional-anatomic fractionation of the brain’s default network. Neuron 65 550–562.
  * Aragam, Gu and Zhou (2017) [author] Aragam, BryonB., Gu, JiayingJ. and Zhou, QingQ. (2017). Learning large-scale Bayesian networks with the sparsebn package. arXiv preprint arXiv:1703.04025.
  * Beijers et al. (2019) [author] Beijers, LianL., Wardenaar, Klaas JK. J., van Loo, Hanna MH. M. and Schoevers, Robert AR. A. (2019). Data-driven biological subtypes of depression: systematic review of biological approaches to depression subtyping. Molecular Psychiatry 24 888–900.
  * Belleau, Taubitz and Larson (2015) [author] Belleau, Emily LE. L., Taubitz, Lauren EL. E. and Larson, Christine LC. L. (2015). Imbalance of default mode and regulatory networks during externally focused processing in depression. Social Cognitive and Affective Neuroscience 10 744–751.
  * Bickel et al. (2008) [author] Bickel, SteffenS., Bogojeska, JasminaJ., Lengauer, ThomasT. and Scheffer, TobiasT. (2008). Multi-task learning for HIV therapy screening. In Proceedings of the 25th International Conference on Machine Learning 56–63.
  * Bielza and Larranaga (2014) [author] Bielza, ConchaC. and Larranaga, PedroP. (2014). Bayesian networks in neuroscience: a survey. Frontiers in Computational Neuroscience 8 131–153.
  * Bluhm et al. (2009) [author] Bluhm, RobynR., Williamson, PeterP., Lanius, RuthR., Théberge, JeanJ., Densmore, MariaM., Bartha, RobertR., Neufeld, RichardR. and Osuch, ElizabethE. (2009). Resting state default-mode network connectivity in early depression using a seed region-of-interest analysis: Decreased connectivity with caudate nucleus. Psychiatry and Clinical Neurosciences 63 754–761.
  * Boyles (1983) [author] Boyles, Russell AR. A. (1983). On the convergence of the EM algorithm. Journal of the Royal Statistical Society: Series B (Methodological) 45 47–50.
  * Buckner, Andrews-Hanna and Schacter (2008) [author] Buckner, Randy LR. L., Andrews-Hanna, Jessica RJ. R. and Schacter, Daniel LD. L. (2008). The brain’s default network: anatomy, function, and relevance to disease. Annals of the New York Academy of Sciences 1124 1–38.
  * Cao et al. (2019) [author] Cao, WenmingW., Qian, ShengS., Wu, SiS. and Wong, Hau-SanH.-S. (2019). Unsupervised multi-task learning with hierarchical data structure. Pattern Recognition 86 248–264.
  * Caruana (1997) [author] Caruana, RichR. (1997). Multitask learning. Machine learning 28 41–75.
  * Chan and Ledolter (1995) [author] Chan, KSK. and Ledolter, JohannesJ. (1995). Monte Carlo EM estimation for time series models involving counts. Journal of the American Statistical Association 90 242–252.
  * Chen et al. (2017) [author] Chen, Jingyuan EJ. E., Glover, Gary HG. H., Greicius, Michael DM. D. and Chang, CatieC. (2017). Dissociated patterns of anti-correlations with dorsal and ventral default-mode networks at rest. Human Brain Mapping 38 2454–2465.
  * Chen et al. (2021) [author] Chen, HaoH., Guo, YingY., He, YongY., Liu, DongD., Liu, LeiL. and Zhou, Xiao-HuaX.-H. (2021). Joint Learning of Multiple Differential Networks with fMRI data for Brain Connectivity Alteration Detection. arXiv preprint arXiv:2106.03334.
  * Chickering (2002) [author] Chickering, David MaxwellD. M. (2002). Optimal structure identification with greedy search. Journal of Machine Learning Research 3 507–554.
  * Constantinou et al. (2020) [author] Constantinou, ACA., Liu, YY., Chobtham, KK., Guo, ZZ. and Kitson, NKN. (2020). The Bayesys data and Bayesian network repository. Queen Mary University of London: London, UK.
  * Cooper and Herskovits (1992) [author] Cooper, Gregory FG. F. and Herskovits, EdwardE. (1992). A Bayesian method for the induction of probabilistic networks from data. Machine Learning 9 309–347.
  * Dawid and Lauritzen (1993) [author] Dawid, A PhilipA. P. and Lauritzen, Steffen LS. L. (1993). Hyper Markov laws in the statistical analysis of decomposable graphical models. The Annals of Statistics 21 1272–1317.
  * Dempster, Laird and Rubin (1977) [author] Dempster, Arthur PA. P., Laird, Nan MN. M. and Rubin, Donald BD. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B (Methodological) 39 1–22.
  * Etkin, Egner and Kalisch (2011) [author] Etkin, AmitA., Egner, TobiasT. and Kalisch, RaffaelR. (2011). Emotional processing in anterior cingulate and medial prefrontal cortex. Trends in Cognitive Sciences 15 85–93.
  * Evgeniou and Pontil (2004) [author] Evgeniou, TheodorosT. and Pontil, MassimilianoM. (2004). Regularized multi–task learning. In Proceedings of the Tenth ACM SIGKDD International Conference on Knowledge Discovery and Data mining 109–117.
  * Evgeniou and Pontil (2007) [author] Evgeniou, AnA. and Pontil, MassimilianoM. (2007). Multi-task feature learning. Advances in Neural Information Processing Systems 19\.
  * Evgeniou et al. (2005) [author] Evgeniou, TheodorosT., Micchelli, Charles AC. A., Pontil, MassimilianoM. and Shawe-Taylor, JohnJ. (2005). Learning multiple tasks with kernel methods. Journal of Machine Learning Research 6 615–637.
  * Fort and Moulines (2003) [author] Fort, GersendeG. and Moulines, EricE. (2003). Convergence of the Monte Carlo expectation maximization for curved exponential families. The Annals of Statistics 31 1220–1259.
  * Friedman (2013) [author] Friedman, NirN. (2013). The Bayesian structural EM algorithm. arXiv preprint arXiv:1301.7373.
  * Friedman, Nachman and Pe’er (2013) [author] Friedman, NirN., Nachman, IftachI. and Pe’er, DanaD. (2013). Learning Bayesian network structure from massive datasets: The” sparse candidate” algorithm. arXiv preprint arXiv:1301.6696.
  * Fu et al. (2008) [author] Fu, Cynthia HYC. H., Mourao-Miranda, JanainaJ., Costafreda, Sergi GS. G., Khanna, AkashA., Marquand, Andre FA. F., Williams, Steve CRS. C. and Brammer, Michael JM. J. (2008). Pattern classification of sad facial processing: toward the development of neurobiological markers in depression. Biological Psychiatry 63 656–662.
  * Gámez, Mateo and Puerta (2011) [author] Gámez, José AJ. A., Mateo, Juan LJ. L. and Puerta, José MJ. M. (2011). Learning Bayesian networks by hill climbing: efficient methods based on progressive restriction of the neighborhood. Data Mining and Knowledge Discovery 22 106–148.
  * Geiger and Heckerman (1994) [author] Geiger, DanD. and Heckerman, DavidD. (1994). Learning gaussian networks. In Uncertainty Proceedings 1994 235–243. Elsevier.
  * Geiger et al. (1997) [author] Geiger, DanD., Heckerman, DavidD. et al. (1997). A characterization of the Dirichlet distribution through global and local parameter independence. Annals of Statistics 25 1344–1369.
  * Gonçalves, Von Zuben and Banerjee (2016) [author] Gonçalves, André RA. R., Von Zuben, Fernando JF. J. and Banerjee, ArindamA. (2016). Multi-task sparse structure learning with gaussian copula models. The Journal of Machine Learning Research 17 1205–1234.
  * Greicius et al. (2007) [author] Greicius, Michael DM. D., Flores, Benjamin HB. H., Menon, VinodV., Glover, Gary HG. H., Solvason, Hugh BH. B., Kenna, HeatherH., Reiss, Allan LA. L. and Schatzberg, Alan FA. F. (2007). Resting-state functional connectivity in major depression: abnormally increased contributions from subgenual cingulate cortex and thalamus. Biological Psychiatry 62 429–437.
  * He et al. (2021) [author] He, YongY., Li, QiushiQ., Hu, QinqinQ. and Liu, LeiL. (2021). Transfer Learning in High-dimensional Semi-parametric Graphical Models with Application to Brain Connectivity Analysis. arXiv preprint arXiv:2112.13356.
  * Heckerman, Geiger and Chickering (1995) [author] Heckerman, DavidD., Geiger, DanD. and Chickering, David MD. M. (1995). Learning Bayesian networks: The combination of knowledge and statistical data. Machine Learning 20 197–243.
  * Honorio and Samaras (2010) [author] Honorio, JeanJ. and Samaras, DimitrisD. (2010). Multi-task learning of Gaussian graphical models. In ICML.
  * Ji, Xia and Meng (2015) [author] Ji, ZhiweiZ., Xia, QibiaoQ. and Meng, GuanminG. (2015). A review of parameter learning methods in Bayesian network. In International Conference on Intelligent Computing 3–12. Springer.
  * Kalisch and Bühlman (2007) [author] Kalisch, MarkusM. and Bühlman, PeterP. (2007). Estimating high-dimensional directed acyclic graphs with the PC-algorithm. Journal of Machine Learning Research 8 613–636.
  * Khundakar and Thomas (2009) [author] Khundakar, Ahmad AA. A. and Thomas, Alan JA. J. (2009). Morphometric changes in early-and late-life major depressive disorder: evidence from postmortem studies. International Psychogeriatrics 21 844–854.
  * Koller and Friedman (2009) [author] Koller, DaphneD. and Friedman, NirN. (2009). Probabilistic graphical models: principles and techniques. MIT press.
  * Koski and Noble (2012) [author] Koski, Timo JTT. J. and Noble, JohnJ. (2012). A review of Bayesian networks and structure learning. Mathematica Applicanda 40 51–103.
  * Lauritzen (1996) [author] Lauritzen, Steffen LS. L. (1996). Graphical models 17\. Clarendon Press.
  * Lee, Parthasarathi and Kable (2021) [author] Lee, SangilS., Parthasarathi, TrishalaT. and Kable, Joseph WJ. W. (2021). The ventral and dorsal default mode networks are dissociably modulated by the vividness and valence of imagined events. Journal of Neuroscience 41 5243–5250.
  * Lemogne et al. (2009) [author] Lemogne, CédricC., le Bastard, GuillaumeG., Mayberg, HelenH., Volle, EmmanuelleE., Bergouignan, LoretxuL., Lehéricy, StéphaneS., Allilaire, Jean-FrançoisJ.-F. and Fossati, PhilippeP. (2009). In search of the depressive self: extended medial prefrontal network during self-referential processing in major depression. Social Cognitive and Affective Neuroscience 4 305–312.
  * Li et al. (2015) [author] Li, YingmingY., Yang, MingM., Qi, ZhongangZ. and Zhang, Zhongfei MarkZ. M. (2015). Bayesian multi-task relationship learning with link structure. IEEE Transactions on Knowledge and Data Engineering 28 873–887.
  * Li et al. (2018) [author] Li, Bao-JuanB.-J., Friston, KarlK., Mody, MariaM., Wang, Hua-NingH.-N., Lu, Hong-BingH.-B. and Hu, De-WenD.-W. (2018). A brain network model for depression: From symptom understanding to disease intervention. CNS Neuroscience & Therapeutics 24 1004–1019.
  * Liu et al. (2013) [author] Liu, FengF., Guo, WenbinW., Liu, LingL., Long, ZhiliangZ., Ma, ChaoqiongC., Xue, ZhiminZ., Wang, YifengY., Li, JunJ., Hu, MaorongM., Zhang, JianweiJ. et al. (2013). Abnormal amplitude low-frequency oscillations in medication-naive, first-episode patients with major depressive disorder: a resting-state fMRI study. Journal of Affective Disorders 146 401–406.
  * Marreiros, Kiebel and Friston (2008) [author] Marreiros, André CA. C., Kiebel, Stefan JS. J. and Friston, Karl JK. J. (2008). Dynamic causal modelling for fMRI: a two-state model. Neuroimage 39 269–278.
  * Mittelbach, Matthiesen and Jorswieck (2012) [author] Mittelbach, MartinM., Matthiesen, BhoB. and Jorswieck, Eduard AE. A. (2012). Sampling uniformly from the set of positive definite matrices with trace constraint. IEEE Transactions on Signal Processing 60 2167–2179.
  * Moreira et al. (2016) [author] Moreira, Mário WLM. W., Rodrigues, Joel JPCJ. J., Oliveira, Antonio MBA. M., Ramos, Ronaldo FR. F. and Saleem, KashifK. (2016). A preeclampsia diagnosis approach using Bayesian networks. In 2016 IEEE International Conference on Communications (ICC) 1–5. IEEE.
  * Muirhead (2009) [author] Muirhead, Robb JR. J. (2009). Aspects of multivariate statistical theory. John Wiley & Sons.
  * Neal et al. (2011) [author] Neal, Radford MR. M. et al. (2011). MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo 11\.
  * Oyen and Lane (2012) [author] Oyen, DianeD. and Lane, TerranT. (2012). Leveraging domain knowledge in multitask Bayesian network structure learning. In Proceedings of the AAAI Conference on Artificial Intelligence 26\.
  * Parameswaran and Weinberger (2010) [author] Parameswaran, ShibinS. and Weinberger, Kilian QK. Q. (2010). Large margin multi-task metric learning. Advances in Neural Information Processing Systems 23 1867–1875.
  * Pearl (1988) [author] Pearl, JudeaJ. (1988). Probabilistic reasoning in intelligent systems: networks of plausible inference. Morgan kaufmann.
  * Raichle and Snyder (2007) [author] Raichle, Marcus EM. E. and Snyder, Abraham ZA. Z. (2007). A default mode of brain function: a brief history of an evolving idea. Neuroimage 37 1083–1090.
  * Ramsey et al. (2017) [author] Ramsey, JosephJ., Glymour, MadelynM., Sanchez-Romero, RubenR. and Glymour, ClarkC. (2017). A million variables and more: the fast greedy equivalence search algorithm for learning high-dimensional graphical causal models, with an application to functional magnetic resonance images. International Journal of Data Science and Analytics 3 121–129.
  * Rei (2017) [author] Rei, MarekM. (2017). Semi-supervised Multitask Learning for Sequence Labeling. In Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers) 2121–2130.
  * Robert (2014) [author] Robert, ChristianC. (2014). Machine learning, a probabilistic perspective.
  * Schwab et al. (2018) [author] Schwab, PatrickP., Keller, EmanuelaE., Muroi, CarlC., Mack, David JD. J., Strässle, ChristianC. and Karlen, WalterW. (2018). Not to cry wolf: Distantly supervised multitask learning in critical care. In International Conference on Machine Learning 4518–4527. PMLR.
  * Schwarz et al. (1978) [author] Schwarz, GideonG. et al. (1978). Estimating the dimension of a model. Annals of Statistics 6 461–464.
  * Seth (2010) [author] Seth, Anil KA. K. (2010). A MATLAB toolbox for Granger causal connectivity analysis. Journal of Neuroscience Methods 186 262–273.
  * Sheline et al. (2010) [author] Sheline, Yvette IY. I., Price, Joseph LJ. L., Yan, ZhiziZ. and Mintun, Mark AM. A. (2010). Resting-state functional MRI in depression unmasks increased connectivity between networks via the dorsal nexus. Proceedings of the National Academy of Sciences 107 11020–11025.
  * Smith et al. (2013) [author] Smith, Stephen MS. M., Vidaurre, DiegoD., Beckmann, Christian FC. F., Glasser, Matthew FM. F., Jenkinson, MarkM., Miller, Karla LK. L., Nichols, Thomas ET. E., Robinson, Emma CE. C., Salimi-Khorshidi, GholamrezaG., Woolrich, Mark WM. W. et al. (2013). Functional connectomics from resting-state fMRI. Trends in Cognitive Sciences 17 666–682.
  * Spirtes and Glymour (1991) [author] Spirtes, PeterP. and Glymour, ClarkC. (1991). An algorithm for fast recovery of sparse causal graphs. Social Science Computer Review 9 62–72.
  * Spirtes et al. (2000) [author] Spirtes, PeterP., Glymour, Clark NC. N., Scheines, RichardR. and Heckerman, DavidD. (2000). Causation, Prediction, and Search. MIT press.
  * Tanaka et al. (2021) [author] Tanaka, Saori CS. C., Yamashita, AyumuA., Yahata, NoriakiN., Itahashi, TakashiT., Lisi, GiuseppeG., Yamada, TakashiT., Ichikawa, NahoN., Takamura, MasahiroM., Yoshihara, YujiroY., Kunimatsu, AkiraA. et al. (2021). A multi-site, multi-disorder resting-state magnetic resonance image database. Scientific data 8 1–15.
  * Teh et al. (2017) [author] Teh, Yee WhyeY. W., Bapst, VictorV., Czarnecki, Wojciech MarianW. M., Quan, JohnJ., Kirkpatrick, JamesJ., Hadsell, RaiaR., Heess, NicolasN. and Pascanu, RazvanR. (2017). Distral: Robust multitask reinforcement learning. arXiv preprint arXiv:1707.04175.
  * Tsamardinos, Brown and Aliferis (2006) [author] Tsamardinos, IoannisI., Brown, Laura EL. E. and Aliferis, Constantin FC. F. (2006). The max-min hill-climbing Bayesian network structure learning algorithm. Machine Learning 65 31–78.
  * Valdes-Sosa et al. (2011) [author] Valdes-Sosa, Pedro AP. A., Roebroeck, AlardA., Daunizeau, JeanJ. and Friston, KarlK. (2011). Effective connectivity: influence, causality and biophysical modeling. Neuroimage 58 339–361.
  * Wang (2018) [author] Wang, YueY. (2018). Analysis of the Max-Min Hill-Climbing Algorithm. In 2018 International Conference on Transportation & Logistics, Information & Communication, Smart City (TLICSC 2018) 509–511. Atlantis Press.
  * Wei and Tanner (1990) [author] Wei, Greg CGG. C. and Tanner, Martin AM. A. (1990). A Monte Carlo implementation of the EM algorithm and the poor man’s data augmentation algorithms. Journal of the American statistical Association 85 699–704.
  * Williams, Bonilla and Chai (2007) [author] Williams, ChrisC., Bonilla, Edwin VE. V. and Chai, Kian MK. M. (2007). Multi-task Gaussian process prediction. Advances in Neural Information Processing Systems 153–160.
  * Wu (1983) [author] Wu, CF JeffC. J. (1983). On the convergence properties of the EM algorithm. The Annals of Statistics 95–103.
  * Wu et al. (2011) [author] Wu, Qi-ZhuQ.-Z., Li, Dong-MingD.-M., Kuang, Wei-HongW.-H., Zhang, Ti-JiangT.-J., Lui, SuS., Huang, Xiao-QiX.-Q., Chan, Raymond CKR. C., Kemp, Graham JG. J. and Gong, Qi-YongQ.-Y. (2011). Abnormal regional spontaneous neural activity in treatment-refractory depression revealed by resting-state fMRI. Human Brain Mapping 32 1290–1299.
  * Xia and He (2017) [author] Xia, MingruiM. and He, YongY. (2017). Functional connectomics from a “big data” perspective. Neuroimage 160 152–167.
  * Yamashita et al. (2020) [author] Yamashita, AyumuA., Sakai, YukiY., Yamada, TakashiT., Yahata, NoriakiN., Kunimatsu, AkiraA., Okada, NaohiroN., Itahashi, TakashiT., Hashimoto, RyuichiroR., Mizuta, HirotoH., Ichikawa, NahoN. et al. (2020). Generalizable brain network markers of major depressive disorder across multiple imaging sites. PLoS biology 18 e3000966.
  * Yamashita et al. (2021) [author] Yamashita, AyumuA., Sakai, YukiY., Yamada, TakashiT., Yahata, NoriakiN., Kunimatsu, AkiraA., Okada, NaohiroN., Itahashi, TakashiT., Hashimoto, RyuichiroR., Mizuta, HirotoH., Ichikawa, NahoN. et al. (2021). Common Brain Networks Between Major Depressive-Disorder Diagnosis and Symptoms of Depression That Are Validated for Independent Cohorts. Frontiers in Psychiatry 12 1–14.
  * Yan et al. (2019) [author] Yan, Chao-GanC.-G., Chen, XiaoX., Li, LeL., Castellanos, Francisco XavierF. X., Bai, Tong-JianT.-J., Bo, Qi-JingQ.-J., Cao, JunJ., Chen, Guan-MaoG.-M., Chen, Ning-XuanN.-X., Chen, WeiW. et al. (2019). Reduced default mode network functional connectivity in patients with recurrent major depressive disorder. Proceedings of the National Academy of Sciences 116 9078–9083.
  * Young et al. (2016) [author] Young, CBC., Chen, TT., Nusslock, RR., Keller, JJ., Schatzberg, AFA. and Menon, VV. (2016). Anhedonia and general distress show dissociable ventromedial prefrontal cortex connectivity in major depressive disorder. Translational Psychiatry 6 e810.
  * Zamoscik et al. (2014) [author] Zamoscik, VeraV., Huffziger, SilkeS., Ebner-Priemer, UlrichU., Kuehner, ChristineC. and Kirsch, PeterP. (2014). Increased involvement of the parahippocampal gyri in a sad mood predicts future depressive symptoms. Social Cognitive and Affective Neuroscience 9 2034–2040.
  * Zarikas, Papageorgiou and Regner (2015) [author] Zarikas, VasiliosV., Papageorgiou, ElpinikiE. and Regner, PeterP. (2015). Bayesian network construction using a fuzzy rule based approach for medical decision support. Expert Systems 32 344–369.
  * Zhang and Yang (2021) [author] Zhang, YuY. and Yang, QiangQ. (2021). A survey on multi-task learning. IEEE Transactions on Knowledge and Data Engineering.
  * Zhang and Yeung (2010) [author] Zhang, YuY. and Yeung, Dit-YanD.-Y. (2010). Multi-task learning using generalized t process. In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics 964–971. JMLR Workshop and Conference Proceedings.
  * Zhao et al. (2014) [author] Zhao, Y-JY.-J., Du, M-YM.-Y., Huang, X-QX.-Q., Lui, SS., Chen, Z-QZ.-Q., Liu, JJ., Luo, YY., Wang, X-LX.-L., Kemp, GJG. and Gong, Q-YQ.-Y. (2014). Brain grey matter abnormalities in medication-free patients with major depressive disorder: a meta-analysis. Psychological Medicine 44 2927–2937.