text:
\section{Introduction}
\label{sec:introduction}

MDD is a common mental disease that is characterized by low mood, worthlessness, anxiety, and cognitive impairments \citep{fu2008pattern}. MDD causes more than 800,000 deaths worldwide every year and is also the leading cause of disability \citep {otte2016major}. 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 connection \citep{wu2011abnormal,liu2013abnormal,yamashita2020generalizable}. Resting-state functional magnetic resonance imaging (rs-fMRI) is an effective modality for understanding the biological and cognitive information of disease \citep{xia2017functional}. 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 (BOLD) signal fluctuations \citep{smith2013functional}.

Bayesian network (BN), as an approach to indicate the cause of interdependence between variables in a system, has drawn much attention \citep{koller2009probabilistic}. BNs visualize the way of information propagation in an intuitive graphical representation \citep{koski2012review}. Such interpretable structure makes them widely applied in scientific domains such as bioinformatics and neuroscience \citep{bielza2014bayesian}. Many BN models have been developed for disease screening, diagnosis, treatment, and prognosis. For example, \cite{moreira2016preeclampsia} build a BN model for the diagnosis of preeclampsia. \cite{zarikas2015bayesian} utilize 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 applying BNs for functional connections of MDD using rs-fMRI data. The main barriers are: (1) In addition to high dimensionality, rs-fMRI data usually has the characteristics of incompleteness, sparsity, variability, high noise, and high error. When the number of variables exceeds 50 to 100, the number of possible network structures also increases super-exponentially. Although some progress has been made in learning large BNs \citep{gamez2011learning,aragam2017learning}, learning and interpreting BNs from rs-fMRI data is still challenging. (2) Due to the difficulty and expensive data acquisition, the amount of rs-fMRI data is relatively small. In the Strategic Research Program for Brain Sciences (SRPBS) database, only 255 MDD patients' rs-fMRI data were collected \citep{yamashita2021common}. Since the heterogeneity of brain functional connections in patients with depression varies from person to person \cite{young2016anhedonia}. There are also different subtypes of depression itself \cite{beijers2019data}. Therefore, the available rs-fMRI data sets are too small to learn robust BNs for individuals or sub-type groups of MDD simultaneously. 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 GBNs simultaneously. 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 individual's neuroimage does 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 disease, their shared disease-induced connectivity alteration in brain functional networks. Therefore, in the learning, we also hope to leverage data from other patients with the assumption that such disease-induced alterations show some degrees of similarity. We assume that such disease-induced alterations show some degree of similarity with the tool to learn such network structures from observations to understand how systems are structured jointly from related tasks. Few work exists on using fMRI data to learn network structures from related tasks. \cite{chen2021joint} consider estimating multiple differential networks with fMRI data sets from multiple research centers, and exploit sparse group Minimax Concave Penalty (gMCP) to induce sparse structures of each differential network. \cite{he2021transfer} propose an algorithm called Trans-Copula-CLIME for estimating undirected graphical model from fMRI datasets, and characterizing the similarity between tasks by the sparsity of a divergence matrix. Differently from the work by \cite{chen2021joint} and \cite{he2021transfer}, we focus on GBN learning. The benefit that GBN brings to studying network structures is that GBN can represent the direct dependencies between the corresponding variable, which can provide more information than differential network and undirected graphical model. As far as we know, this is the first work to use Bayesian approach to learn multiple GBNs from related tasks on MTL.

The proposed MTGBN method enhances the learning ability of GBNs from observations and makes further improvements with the proposed multitask learning technique. Specifically, we treat each subject in a class of observation as a task and then learn the GBNs of this data class by learning from all tasks that share a default covariance matrix that encodes prior knowledge. This setting can help us to learn more information from limited data. In structure learning of BNs, we derive a closed-form formula of the complete likelihood function and use the MCEM algorithm to efficiently search for the approximately best BN structures.

Our contributions to this paper are threefold.
\begin{itemize}
    \item We propose an MTGBN method.
    To cope with the "learning BNs from insufficient high dimensional data" issue, MTGBN jointly learns multiples BN structures from related observations. The approach addresses one common problem in pattern recognition: how to learn representing network structures from related observations such that the representations can share related information while also maintaining subject-specific characteristics.
    \item We deal with the computationally intractable issue of the existing MTL framework for GBNs. A closed-form expression is derived for the complete likelihood and design Hamiltonian Monte-Carlo (HMC) with exact gradient evaluation from the positive semidefinite cone so that the MCEM algorithm can be used.
    \item We focus on MDD. The computational results from rs-fMRI data of MDD patients show that our method can produce a better explanation. Doctors can use it to understand clearly the potential brain functional connections of each MDD patient and carry out the personalized treatment.
\end{itemize}

The paper is organized as follows. In Section \ref{sec:data} we begin by presenting the rs-fMRI data sets analyzed in this work. In Section \ref{sec:background}, we review the BN learning and MTL. In Section \ref{sec:MTGBN}, we present the formulation of the proposed MTGBN method. In Section \ref{sec:structure_learning_mtgbn}, we present the structure learning of MTGBN. An efficient structure learning algorithm for the MTGBN is given in Section \ref{sec:computing}. We conduct simulations and discuss the results in Section \ref{sec:simulation}. In Section \ref{sec:applications} we analyze a real rs-fMRI data set and assess the advantages of our method. Next, we conclude this paper in Section \ref{sec:conclusion}.

\section{Real-world rs-fMRI Data}
\label{sec:data}

In this work we consider real-world  rs-fMRI data of MDD.

\subsection{Data acquisition}
\label{subsec:acquisition}

The rs-fMRI data for each participant subject at resting state was 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.5$s, TE=$30$ms, 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.8$mm, slice thickness=$3.2$mm, number of slices= $40$, number of volumes= $240 + 4$ (dummy), total scan time= $10$ min + $10$s (dummy). The T1-weighted images were then acquired with the following parameters: voxel size= $1\times 1\times 1$, TR=$2300$ms, TE=$2.98$ms, TI=$900$ms; flip angle = $90^\circ$, FOV= $256$, matrix= $256\times 256$. The data utilized in this study can be downloaded publicly from the DecNef Project Brain Data Repository. \footnote{ \url{https://bicr-resource. atr.jp/srpbsopen/}}

\subsection{Samples}
\label{subsec:samples}

We use the rs-fMRI data of $15$ patients with MDD and $15$ Healthy Controls (HCs) 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 with BDI scores between $20-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-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 (DARTEL) in SPM12. The rs-fMRI data were smoothed with a 6-mm 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 selected $30$ brain regions (regard each region as an region of interest (ROI)) based on three different default mode network (DMN) templates networks related to the DMN (including the Andrews-Hanna Default Mode Networks (aDMN, see \citealp{buckner2008brain}), the dorsal Default Mode Networks (dDMN, see \citealp{chen2017dissociated}), and the ventral Default Mode Networks (vDMN, see \citealp{lee2021ventral}). Datails of brain regions are presented in  Table \ref{tab: roi}. The voxel-wise fMRI time courses were 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 above steps 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. Thus, the central limit theorem is applicable. Also, the usual hypothesis is common for fMRI studies \citep{valdes2011effective}.

\section{Background and related work}
\label{sec:background}

\subsection{Learning Bayesian networks}
\label{subsec:BNs}

BN is a formal method to analyze decision-making strategies under uncertain conditions \citep{pearl1988probabilistic}.
Given a random vector $\boldsymbol{X}=(X_{1}, \ldots, X_{p})^\top$, we use $\mathbb{P}(\boldsymbol{x})$ to denote its joint probability mass function (for the discrete case) or the joint density function (for the continuous case), where $\boldsymbol{x}$ is an observation of $\boldsymbol{X}$. From the chain rule of probability, we can write $\mathbb{P}(\boldsymbol{x})$ as
\begin{equation}
\mathbb{P}(\boldsymbol{x})=\mathbb{P}(x_{1})\prod_{i=2}^{p} \mathbb{P}\left(x_{i} \mid x_{1}, \ldots, x_{i-1}\right).
\label{eq:Px}
\end{equation}
For each $X_{i}$, let $\Pi_{i} \subseteq\left\{X_{1}, \ldots, X_{i-1}\right\}$ be a set of variables such that
\begin{equation}
\mathbb{P}\left(x_{i} \mid x_{1}, \ldots, x_{i-1}\right)=\mathbb{P}\left(x_{i} \mid \Pi_{i}\right).
\label{eq:Pxcond}
\end{equation}

A BN $B=(\mathcal{G}, \Theta)$ can represent the joint distribution of $\boldsymbol{X}$ through a directed graphical model, where $\mathcal{G}$ is a network structure that encodes the assertions of conditional independence in equation \eqref{eq:Pxcond}, and $\Theta$ includes the parameters of the BN model corresponding to $\mathcal{G}$ \citep{koller2009probabilistic}. 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 \eqref{eq:Px} and \eqref{eq:Pxcond}, we can see that the joint distribution for $\boldsymbol{X}$ can be represented graphically by a BN, namely,
\[
\mathbb{P}(\boldsymbol{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 \citep{geiger1994learning}. In GBNs, $\boldsymbol{X}$ follows a multivariate normal distribution $\mathcal{N}\left(\mu, \Sigma \right)$, with the joint density function
\[
   \mathbb{P}(\boldsymbol{x})= (2 \pi)^{-p/2}\left|\Sigma\right|^{-1/2}\exp\left(-\frac{1}{2}(\boldsymbol{x}-\mu)^{\top} \Sigma^{-1}(\boldsymbol{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$ 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: parameter learning and structure learning.

\paragraph*{Parameter learning}

In parameter learning, there are two main methods: (1) maximum likelihood estimation (MLE, see \citealp{ji2015review}); (2) Bayesian estimation \citep{geiger1997characterization}. MLE method estimates the parameters by maximizing the likelihood of the data. However, when the data is sparse, MLE has problems because the data set does not contain all possible combinations of the involved variables. To solve this problem, a prior distribution is usually placed on the variables \citep{altendorf2005learning} for the Bayesian estimation. However, only parameter learning is not enough. Generally speaking, parameter learning is based on the assumption that we know the network structure in advance (or at least we determine one regardless of whether it is correct). Therefore, learning the network structure without knowing the network structure is more important for BNs.

\paragraph*{Structure learning}

Structure learning is an essential part of BNs, mainly divided into (1) constrained-based methods; (2) score-based methods. The strategy of learning BNs through constrained-based methods tries to identify dependent and independent areas within the data \citep{koller2009probabilistic}. One of the critical algorithms for structure learning is the SGS \citep{spirtes2000causation} 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 increases. Compared with SGS, the PC algorithm \citep{spirtes1991algorithm} 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 \citep{schwarz1978estimating}; K2 \citep{cooper1992bayesian}; BDeu \citep{heckerman1995learning}. Because searching network structures is NP-hard. For large networks, it is unfeasible to use score only. \cite{friedman2013learning} improve it by reducing the search space of score-based structure learning methods.

Given the observed dataset $D$, structure learning attempts to learn a DAG $\mathcal{G}$ that encodes the structural dependence of the variables in $\boldsymbol{X}$, whereas parameter learning estimates the parameters $\Theta$ associated with $\mathcal{G}$. In this work, we focus on 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$:
\begin{equation}
\mathbb{P}(\mathcal{G} \mid D) = \frac{1}{c}\ \mathbb{P}(\mathcal{G}) \times \mathbb{P}(D \mid \mathcal{G}),
\label{eq:Pgd}
\end{equation}
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. Note that even in a small domain, there may be a huge number of network structures to sum over, making the computing of $c$ intractable. Therefore, we use the complete likelihood $\mathbb{P}(\mathcal{G}) \times \mathbb{P}(D \mid \mathcal{G})=\mathbb{P}(D ,\mathcal{G})$ as an equivalent but tractable score function.

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 \citep{heckerman1995learning}.

\subsection{Multitask learning}
\label{sebsec:MTL}

General machine learning methods usually require many samples to learn the robust model. In the case of insufficient data, MTL \citep{caruana1997multitask} is a general approach that incorporates task relatedness during learning and can learn robust models \citep{zhang2021survey}. 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 multitask learning 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 \citep{schwab2018not}, unsupervised tasks \citep{cao2019unsupervised}, semi-supervised tasks \citep{rei2017semi}, reinforcement learning tasks \citep{teh2017distral} or graphical tasks \citep{honorio2010multi}.

According to \cite{zhang2021survey}, the current research of MTL focuses on characterizing the relatedness between tasks, which can be divided into three parts, including instance-based MTL \citep{bickel2008multi}, feature-based MTL \citep{evgeniou2007multi}, parameter-based MTL \citep{parameswaran2010large}. 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 \citep{evgeniou2004regularized} 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 \citep{evgeniou2005learning}, e.g., with Gaussian copula models \citep{gonccalves2016multi} and Matrix Generalized Inverse Gaussian distribution \citep{li2015bayesian}. However, when task relations are not available in most applications, it is necessary to automatically learn task relationships from data. \cite{williams2007multi} propose a multitask Gaussian process to learn task covariance from data. \cite{zhang2010multi} propose a multitask generalized $t$ process by placing an inverse-Wishart prior.

\section{The proposed MTGBN framework}
\label{sec:MTGBN}

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, which compensates for the lack of data for minor subtypes.

Formally, let $\mathcal{D}=\left\{D_{1}, \ldots, D_{m}\right\}$ be $m$ related observational data sets, 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 \ref{framework}  shows the generative framework of the MTGBN method, where the covariance matrices $\Sigma_i$ are sampled from a common prior distribution parameterized by $\Sigma_h$, and each $\Sigma_i$ generates a data set $D_i$ for task $i$. $\Sigma_h$ is inferred from all the tasks, so it helps transfer the learned knowledge from one task to another.

\begin{figure}
\includegraphics[width=0.7\textwidth]{MTGBN}
\caption{Illustration of the learning framework of MTGBN}
\label{framework}
\end{figure}

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 $I W\left(\Sigma_{h}, \nu_{0}\right)$ \citep{robert2014machine}, with the following joint density function:
\begin{equation}
\begin{aligned}
&\mathbb{P}\left(\boldsymbol{\Sigma} \mid \Sigma_{h}\right)=\prod_{i=1}^{m} \mathbb{P}\left(\Sigma_{i} \mid \Sigma_{h}\right),\\
&\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\},
\label{eq:Psigma}
\end{aligned}
\end{equation}
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 \citep{muirhead2009aspects}, $\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$ data sets $\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
\begin{equation}
\begin{aligned}
&\mathbb{P}\left(\mathcal{D} \mid\boldsymbol{\Sigma}\right)=\prod_{i=1}^{m} \mathbb{P}\left(D_{i} \mid \Sigma_{i}\right),\\
&\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\},
\end{aligned}
\label{eq:Pdsigmai}
\end{equation}
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. 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 data set 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 data sets $\mathcal{D}=\left\{D_{1}, \ldots, D_{m}\right\}$, with $B_i$ corresponding to $D_i$.

\section{Structure learning of MTGBN}
\label{sec:structure_learning_mtgbn}

As is described in Section \ref{subsec:BNs}, 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 \eqref{eq:Pgd}, we follow the likelihood principle, \emph{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 \ref{framework} 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)$,
\begin{equation}
\mathbb{P}\left(\mathcal{B}, \mathcal{D}\right)=\int \mathbb{P}\left(\Sigma_{h}, \mathcal{B}, \mathcal{D}\right) \mathrm{d} \Sigma_{h},
\label{eq:Pgdb}
\end{equation}
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, \citealp{dempster1977maximum}) 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 \ref{subsec:complete_likelihood} and \ref{subsec:conditional_distr}, respectively, and then outline the EM algorithm in Section \ref{subsec:em_algorithm}.

\subsection{The complete likelihood}
\label{subsec:complete_likelihood}

 The complete likelihood admits the decomposition
\begin{equation}
\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).
\label{eq:completell}
\end{equation}
In Section \ref{sec:MTGBN} we have assumed that $\mathbb{P}\left(\Sigma_{h}\right)\propto 1$, 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.

\begin{definition}
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 \emph{subgraph} $G^{\prime}=(V,E^{\prime})$. A subgraph is called \emph{complete} if every pair of nodes are adjacent. Define a \emph{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 \emph{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 \emph{decomposable} if it is complete, or if there exists a proper decomposition $(A,B)$ into complete subgraphs $G^{A}$ and $G^{B}$ \citep{dawid1993hyper}. In other words, a decomposable graph can be decomposed into several cliques, and each clique is a complete subgraph of the graph.
\end{definition}

Let $\widetilde{G}$ be a decomposable undirected graph consisting of $k$ cliques in the sequence $\left\{C_{1}, \ldots, C_{k}\right\}$. Define the \emph{histories} of the sequence as the sets $H_{j}=\bigcup_{i=1}^{j} C_{i}$ for $j=1, \ldots, k$. Similarly, define the \emph{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_{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 \citep{lauritzen1996graphical}
\[
\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 \ref{fig:graph_trans}(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~\ref{fig:graph_trans}(c) demonstrates that $G$ can be decomposed into four cliques and three separators.

\begin{figure}
\includegraphics[width=0.7\textwidth]{graph_trans}
\caption{An example of graph tranformation and decomposition.}
\label{fig:graph_trans}
\end{figure}

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.

\begin{theorem}
\label{Theorem1}
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
\begin{equation}
\mathbb{P}\left(B_{i} \mid \Sigma_{h}\right)=\frac{\prod_{C \in \mathcal{C}}\left|2^{-1} \Sigma_{h C C}\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_{h S S}\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}),
\label{eq:hfunc}
\end{equation}
where $|C|$ and $|S|$ are the cardinality of the clique $C$ and the separator $S$, respectively, and $\Sigma_{h C C}$ and $\Sigma_{h S S}$ denote the submatrices of $\Sigma_{h}$ according to the graph decomposation.
\end{theorem}

We make the following remark to further explain the implications of Theorem \ref{Theorem1}.

\begin{remark}
$\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 \eqref{eq:Psigma}, $\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 \ref{Theorem1} is to measure the density of a sparse BN $B_{i}$ given a covariance matrix $\Sigma_{i}$, \emph{i.e.}, $\mathbb{P}\left(B_{i} \mid \Sigma_{i}\right)$. Fortunately, \cite{dawid1993hyper} 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})$.
\end{remark}

\begin{theorem}
\label{Theorem2}
For the target structure $B_i$ corresponding to dataset $D_i$,
\begin{equation}
\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)},
\label{eq:Pdbsigma}
\end{equation}
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 \eqref{eq:Pdbsigma} is of the form
\begin{equation}
\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)},
\label{eq:Pdsigmah}
\end{equation}
where $\Sigma_{hCC}$ is the sub-matrix of $\Sigma_{h}$, and $S_{iCC}$ is the sub-matrix of the sample covariance matrix $S_{i}$.
\end{theorem}

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 \eqref{eq:completell}.

\subsection{The conditional distribution}
\label{subsec:conditional_distr}

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 \ref{Theorem3} gives its density function in a closed form.

\begin{theorem}
\label{Theorem3}
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 Eq.\eqref{eq:hfunc}.
\end{theorem}

\subsection{The EM algorithm}
\label{subsec: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:
\begin{itemize}
\item \textbf{E-step}: Define the function
	\begin{equation}
	\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},\label{eq:Qfunc}
	\end{equation}
\item \textbf{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).
	\]
\end{itemize}

The basic convergence properties of the EM algorithm were established by \cite{boyles1983convergence} and \cite{wu1983convergence}. \cite{wu1983convergence} discusses several convergence problems of EM algorithm in detail. He shows that when the complete data are from a curved exponential family with compact parameter space, and when the $\mathcal{Q}$-function satisfies a certain mild differentiability condition, then any EM sequence converges to a stationary point  of the likelihood function.  When applying EM algorithm to structure learning, there is no corresponding properties. Fortunately, \cite{friedman2013bayesian} have proved that the algorithm converges when there is no further improvement in the objective score.

However, the $\mathcal{Q}$-function in the E-step \eqref{eq:Qfunc} is still intractable due to the high-dimensional integration. Instead, we view the integral in \eqref{eq:Qfunc} 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 is more complicated to point out the convergence properties of MCEM than ordinary EM. Some scholars have made some efforts to establish convergence properties of MCEM. \cite{chan1995monte} treat the data as fixed data and keep the Monte Carlo sample size unchanged during MCEM iterations. Then, they set the sample size to infinity and study the asymptotic properties of MCEM sequence as the Monte Carlo approximation to the ordinary EM sequence with the same starting value. \cite{fort2003convergence} regard the data as fixed, the Monte Carlo sample size increases with the number of MCEM iterations and tends to infinity with the number of iterations, so as to establish the convergence of the sequence. We describe the details of MCEM in Section \ref{sec:computing}.

\section{Computing method}
\label{sec:computing}

\subsection{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 \eqref{eq:Qfunc} by the following unbiased estimator:
\begin{equation}
\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),
\label{eq:Qtilde}
\end{equation}
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 \citep{neal2011mcmc} 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 (\emph{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 variables\footnote{The similar method is also used by the Stan software: \url{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 \ref{thm:distr_v}.

\begin{theorem}
\label{thm:distr_v}
Let $\mathbf{diag}\{A\}$ be a vector formed by the diagonal elements of a matrix $A$, and $\mathbf{lower}\{A\}$ be the lower-triangular part 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 element-wise multiplication of two vectors $u$ and $v$. Then the conditional density of $V$ is
\begin{equation}
	\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}\cdot2^{p}\prod_{i=1}^{p}l_{ii}^{p-i+2},\label{eq:density_v}
\end{equation}
and the gradient of the log-density function is
\begin{equation}
\begin{aligned}
	\mathbf{diag}\{\nabla_{V} \log \mathbb{P}(V|\mathcal{B}^{(t-1)},\mathcal{D})\}&=\mathbf{diag}\{\widetilde{L}+q(m,\nu_{0},p)\cdot I\} \circ \mathbf{diag}\{L\},\\
	\mathbf{lower}\{\nabla_{V} \log \mathbb{P}(V|\mathcal{B}^{(t-1)},\mathcal{D})\}&=\mathbf{lower}\{\widetilde{L}\},
\end{aligned}
\label{eq:deriv_v}
\end{equation}
where
\begin{align*}
	\widetilde{L}&=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}),\\
	g(\nu,\Sigma,\mathcal{B})&=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\},\\
	f(\nu,\Sigma)&=\nu \cdot \Sigma^{-1}L,\\
	q(m,\nu,p)&=\frac{p(p+2m\nu+3)}{2}.
\end{align*}
\end{theorem}
Based on equations \eqref{eq:density_v} and \eqref{eq:deriv_v}, 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 \ref{Appendix d}.

\subsection{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 \eqref{eq:Qtilde}, and take $\mathcal{B}^{(t)}$ to be the maximizer. 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\}.
\]

We use Max-Min Hill-Climbing (MMHC, see \citealp{tsamardinos2006max}) to search the Bayesian network structures, which combines ideas from local learning, constraint-based and search-and-score techniques. It has been empirically justified that MMHC outperforms on average several prototypical algorithms \cite{wang2018analysis}, such as PC \cite{kalisch2007estimating}, Sparse Candidate \cite{friedman2013learning}, Greedy Equivalence Search \cite{ramsey2017million}, and Greedy Search \cite{chickering2002optimal}.

According to the above definition, the specific MCEM algorithm of MTGBN given in Algorithm \ref{mcem}.
\begin{figure}
		\renewcommand{\algorithmicrequire}{\textbf{Input:}}
		\renewcommand{\algorithmicensure}{\textbf{Output:}}
		\let\@latex@error\@gobble
		\begin{algorithm}[H]
			\caption{ MCEM of MTGBN}
			\begin{algorithmic}[1]\label{mcem}
				\REQUIRE $\mathcal{D}$: Observational datasets; $\mathcal{B}^{(0)}$: Initial BN structures; $\epsilon$: Tolerance.
					\ENSURE
		      $\mathcal{B}^{(T)}$.
				\STATE \textbf{loop} $t= 0,1 , \ldots, T$
					\STATE \quad \textbf{MC E-step:}
					\STATE \qquad  Computing $\widetilde{\mathcal{Q}}\left(\mathcal{B} \mid \mathcal{B}^{(t-1)}, \mathcal{D}\right)$ by equation \eqref{eq:Qtilde}.
					\STATE \quad \textbf{MC M-step:}
					\STATE \qquad $\{\Sigma_{h}^{(l)}\}_{l=1}^{N}$ $\leftarrow$ \textbf{HMC} ($\mathcal{D}$, $\mathcal{B}^{(t-1)}$), where the detail of HMC algorithm for MTGBN can be found in Appendix~\ref{Appendix d}.
					\STATE \qquad\textbf{for} $i=1, \ldots, m$ \textbf{do}
					\STATE \quad\qquad $B_{i}^{(t)}$ $\leftarrow$ \textbf{MMHC} ($D_{i}, B_{i}^{(t-1)}$, $\{\Sigma_{h}^{(l)}\}_{l=1}^{N}$),
					where the MMHC is the classical heuristic algorithm for detecting the Bayesian network structure \citep{tsamardinos2006max}.
					\STATE \qquad\textbf{end for}
					\STATE \qquad Let $\mathcal{B}^{(t)}=\left\{B_{i}^{(t)}\right\}_{i=1}^{m}$.
                    \STATE \quad \textbf{Break}: if $$|\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)|\leq \epsilon,$$ \textbf{then}
                    \STATE \quad \textbf{return} $\mathcal{B}^{(T)}$.
				\STATE \textbf{end loop}
			\end{algorithmic}
		\end{algorithm}
\end{figure}

\section{Simulations}
\label{sec:simulation}

\subsection{Simulation design}
This study empirically evaluates the finite sample performance of proposed MTGBN on synthetic data and benchmark network data compared with several alternative methods.
\paragraph*{Simulation design of synthetic data}

We generate synthetic datasets by our data generation mechanism (see Figure \ref{framework}). Thus, first we need to construct the shared prior corvariance $\Sigma_{h}$. As mentioned above, the hyper prior distribution of $\Sigma_{h}$ follows a uniform distribution over positive definite matrices. We use an approach similar to \cite{mittelbach2012sampling}, randomly generate a  $(p \times p)$ positive definite matrix with unit trace from the uniform distribution as the initial value of $\Sigma_{h}$. Then we sample $m$ covariance matrices $\left\{\Sigma_{l}\right\}_{l=1}^{m}$ from $\Sigma_{h}$, and each $\Sigma_{l}$ in $\left\{\Sigma_{l}\right\}_{l=1}^{m}$ is a $(p \times p)$ positive definite matrix.

The corresponding $m$ undirected graph structures  $\left\{G_{l}\right\}_{l=1}^{m}$ are defined as follows. Let $\Omega_l=\Sigma_l^{-1}=(\omega_{ij,l})\in\mathbb{R}^{p^2},l=1,\dots,m$. Let $V=\{1, \ldots, p\}$ and define $E_l=\{(i, j), i \neq j$ such that $\omega_{ij,l} \neq 0 \}$, then $G_l=(V,E_l)$. We further denote $|E_l|$ as the cardinality of $E_l$, so $s_l=|E_l|/p^2$ can be used to represent the sparsity of $\Omega_l$ or $G_l$. In fact, $s_l$ represents the ratio of nonzero elements to the total number of elements in $\Omega_l$. In order to control the sparsity of $\Omega_l$, we shrink the part of $\omega_{ij,l}$ with small absolute value to zero (under the premise that transformed $\Omega_l$ is a invertible matrix). In this way, we can obtain $\Omega_l$ with different sparsity $s_l$. Actually, the graph of $\Sigma_l$ gives the configuration of nonzero off-diagonal elements in $\Omega_l$ and $\omega_{ij,l} \neq 0$ implies the conditional independence structure \citep{dawid1993hyper}.

As discussed in Section \ref{subsec:BNs}, we generate $m$ observational datasets $\left\{D_{l}\right\}_{l=1}^{m}$ as follows. For each $D_{l}$ in $\left\{D_{l}\right\}_{l=1}^{m}$, we generate $D_{l}=\left\{\mathbf{x}_{l 1}, \ldots, \mathbf{x}_{l n_{i}}\right\}$ using $\Sigma_{l}$ as the true covariance matrices by sampling observation from the corresponding multivariate normal distribution $\mathcal{N}\left(0, \Sigma_{l}\right)$. By changing the value of $n_{i}$, we can obtain data with different sample size.

To evaluate and compare the performance of the MTGBN method (MTL for short), we also learn each network independently with single-task learning (STL for short) and learn a single network structure for all contexts (AVG for short), i.e., an "average" network that represents all tasks in the class. Note that AVG treats all subjects as samples in one task and learns the identical structure. This method was previously used as a baseline in \cite{oyen2012leveraging}. To compare the performance of MTL, STL, and AVG, especially the differences between the structures, we use edge adjacency error, precision, F-score, and edge distance as the primary evaluation measures.

Let $\left\{G_{l}\right\}_{l=1}^{m}$ represent the truth structures which are corresponding to $\left\{\Sigma_{l}^{-1}\right\}_{l=1}^{m}$. Let $\left\{B_{output, l}\right\}_{l=1}^{m}$ represent the $m$ BN structures that obtained from MTL, SIG and AVG method. Let $\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, because the generated network  $\left\{G_{l}\right\}_{l=1}^{m}$ are undirected graphs, we mainly evaluate the network structures without considering the edge orientation. To compute adjacency error, precision, and F-score, we first calculate four basic statitics: true positives (TP), false positives (FP), false negatives (FN) and true negatives (TN). For example, when $l=1$, TP is the number of edges that are common in both $G_{output, 1}$ and $G_{1}$ without considering the edge orientation. FP is the number of edges that are present in $G_{output, 1}$ but not in $G_{ 1}$. FN is the number of edges that are present in $G_{ 1}$ but not in $G_{output, 1}$. TN is the number of no edges that are common in both $G_{output, 1}$ and $G_{1}$. Then, the edge adjacency error, precision, F-score and edge distance can be derived as

\begin{itemize}
    \item \textbf{Error:}
     \begin{equation}
     \text {\emph{Error}}=\frac{FP+F N}{T P+T N+F P+F N}. \label{eq:error}
    \end{equation}
    \item \textbf{Precision:}
    \[
     \text { \emph{Precision} }=\frac{T P}{T P+F P}.
    \]
     \item \textbf{F-score:}
     \[
      \text { \emph{F-score}}=2 * \frac{(\text { \emph{Precision} } * \text {\emph{Recall} })}{(\text { \emph{Precision} }+\text {\emph{Recall}})},
     \]
     where $\text { \emph{Recall}}=T P/(T P+FN)$.
     \item \textbf{Edge distance:} Given two structures, this measure is the minimum number of edges needed to convert the learnt graph into the true graph.
\end{itemize}

We will compare the performance between MTL, STL, and AVG, focusing on the following different situations.

\textbf{Performance vs. Sample Size:} We use $m$ to represent the task number (the network structure learning of each dataset  $D_{l}$ in $\left\{D_{l}\right\}_{l=1}^{m}$ is regarded as a task), and use $p$ to represent the number of nodes. Then we set task number $m=10$, number of nodes $p=15$, sparsity $s=0.3$, and compare the model performance of network structure learning when the sample size $n_i=(50, 100, 150, 200, 250, 300)$. We regard the network structure learning of $m$ datasets as one learning step, and each learning step calculates the average of the evaluation measures. We repeat the above learning step $10$ times and further compute the average value of evaluation measurement over those runs.

\textbf{Performance vs. Task Number :} We set number of nodes $p=15$, sparsity $s=0.3$, sample size $n_i=250$, and compare the model performance of network structure learning when task number $m=(2, 4, 6, 8, 10)$. Similar to above, we repeat the learning step 10 times, and compute the average value of evalution measurment.

\textbf{Performance vs. Sparisy:} We also evaluate the model performance against sparsity. We set task number $m=10$, number of nodes $p=15$ and compare the model performance of network structure learning when sample size $n_i$ vary from $100$ to $300$ and sparsity vary from $0.2$ to $0.6$. Also, we repeat the learning step $10$ times and compute the average value of average errors output from each learning step. This comparison will help to discuss the applicability of MTGBN.

\paragraph*{Simulation design of benchmark networks data}

In the experiment of benchmark network data, we select benchmark networks from the Bayesian Network Repository (BNR, \citealp{constantinou2020bayesys}). We choose three medium-sized benchmark networks and generate low sample size data to validate our assumption that the proposed MTGBN can improve understanding of the underlying structure and generalization even at a low sample size. The information of the selected benchmark is provided in Table \ref{tab:benchmark}.
\begin{table}[h!]
\centering
\caption{Facts of Benchmark Networks from Bayesian Network Repository (BNR) }
\begin{tabular}{@{}llll@{}}
\toprule
Benchmark Networks & Nodes & Arcs & Parameters \\
\midrule
ECOLI70 & 46 & 70 & 162 \\
MAGIC-NIAB & 44 & 66 & 154 \\
MAGIC-IRRI & 64 & 102 & 230 \\
\bottomrule
\end{tabular}
\label{tab:benchmark}
\end{table}

To obtain a set of related networks, we modify some structures of the given benchmark network to create similar but different networks. For each benchmark BN, we simulate $m$ generated BNs by varying different perturbation levels to them. The perturbation level represents the degree of difference between the generated BN and the benchmark BN, equal to the proportion of the number of changed edges and the number of possible edges $(p \cdot(p-1))$. We denote benchmark BN as $B_{truth}$ with adjacency matrix $A_{t}=(a^{t}_{ij})\in\mathbb{R}^{p^2}$, and denote generated BN as $B_{generated}$ with the adjacency matrix $A_{g}=(a^{g}_{ij})\in\mathbb{R}^{p^2}$. For example, when the perturbation level equals $30\%$, we randomly select $30\%$ of the edges to modify. 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 randomly selected element, under the condition that $B_{generated}$ is a directed acyclic graph, we modify it as follows: if $a^{t}_{ij}=0$, then set $a^{g}_{ij}=1$ (adding edge), else if $a^{t}_{ij}=1$, set $a^{g}_{ij}=0$ (delete edge) or set $a^{g}_{ji}=1$ (reverse edge) with the probability $1/2$. For each benchmark BN, we generate ten of them ($m=10$) by such perturbation. Each set of generated networks can be considered 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_i=250$.

We will mainly compare edge arrowhead error under different network perturbation levels in benchmark network comparison. The expression of edge arrowhead error is the same as equation \eqref{eq:error}, expect that TP, FP, FN, TN are calculated considering the edge orientation. Taking the same definition as in synthetic data, we use $\left\{B_{output, l}\right\}_{l=1}^{m}$ to represent the $m$ BN structures learned from benchmark network data by using MTL, SIG, and AVG. For example, when  $l=1$, TP is the number of edges that are common in both $B_{output,1}$ and $B_{truth}$ and share the same edge orientation. For each experiment of benchmark network data, we regard the network structure learning of $10$ data sets generated from the benchmark network as one learning step, and each learning step calculates the average edge arrowhead error between MTL, STL, and AVG. We repeat the above learning step 10 times and compute the average value of average arrowhead errors output from each learning step.

\subsection{Simulation results}

In this section, we perform the simulation results for sybthetic data and benchmark network Data.
\paragraph*{Simulation results for synthetic data}

We evaluate the performance of the proposed MTGBN method based on simulation data of different sample sizes. Figure \ref{fig:preformance_datasize} show the performance of the various methods (MTL, AVG, and SIG) for the GBNs with different sample sizes. Both results show that MTL significantly outperforms STL all the time. The average adjacency error of MTL in different data sizes is in the range of $0.03-0.06$. STL, even using five times the number of training data (sample size $=250$), the average adjacency error is still far from reaching MTL (sample size $=50$). The same trend is shown in the average edge distance. As the figure indicates, MTL often results in higher average adjacency precision and average adjacency F-score than AVG and MTL. Compared with SIG, the average adjacency precision of MTL is more stable with different sample sizes, and the best results are given when the sample size $=200$. Such analysis shows that the proposed method can produce satisfactory results using a small amount of data (efficiency) and demonstrate its usefulness (applicability) in cases where data are expensive or hard to collect.

\begin{figure}[htbp]
    \centering
	  \subfloat[]{
       \includegraphics[width=0.45\linewidth]{Error_datasize.pdf}}
    \label{2a}\hfill
	  \subfloat[]{
        \includegraphics[width=0.45\linewidth]{Edgedistance_datasize.pdf}}
    \label{2b}\\
	  \subfloat[]{
        \includegraphics[width=0.45\linewidth]{Precision_datasize.pdf}}
    \label{2c}\hfill
	  \subfloat[]{
        \includegraphics[width=0.45\linewidth]{F_score_datasize.pdf}}
     \label{2d}
	  \caption{Performance against sample size.}
	  \label{fig:preformance_datasize}
	\end{figure}

In addition, we also compare the performance of MTL, AVG, and SIG when the task number is different. Figure \ref{fig:preformance_task} indicates that MTL always performs better than AVG and SIG. We provide $250$ samples per task. Although we increase task number, these samples are still not enough for good SIG learning because the increase in task number does not provide additional information. It also indirectly shows that MTL can use the information of different tasks and contribute to better learning. We find that when the task number $=4$, MTL can get a good result, and with the increase in task number, the performance tends to be stable. These results are also consistent with existing literature \citep{oyen2012leveraging}.

\begin{figure}[htbp]
    \centering
	  \subfloat[]{
       \includegraphics[width=0.45\linewidth]{Error_task.pdf}}
    \label{3a}\hfill
	  \subfloat[]{
        \includegraphics[width=0.45\linewidth]{Edgedistance_task.pdf}}
    \label{3b}\\
	  \subfloat[]{
        \includegraphics[width=0.45\linewidth]{Precision_task.pdf}}
    \label{3c}\hfill
	  \subfloat[]{
        \includegraphics[width=0.45\linewidth]{F_scoret_task.pdf}}
     \label{3d}
	  \caption{Performance against task number.}
	  \label{fig:preformance_task}
	\end{figure}

In the experiment, we find that sparsity also affects the model's performance in addition to the sample size and task number. To explore the applicability of our method, we compare the performance of MTL, AVG, and SIG under different sample sizes and sparsity. In Figure \ref{fig:preformance_sparsity}, the smaller the x-axis value, the greater the sparsity. Figure \ref{fig:preformance_sparsity} indicates that with the increase in sample size and sparsity, the average adjacency error performance of the model is better and better, which is more obvious in MTL.

\begin{figure*}[]
    \centering
	  \subfloat[]{
       \includegraphics[width=0.45\linewidth]{sig.pdf}}
    \label{4a}
	  \subfloat[]{
        \includegraphics[width=0.45\linewidth]{avg.pdf}}
    \label{4b}
	  \subfloat[]{
        \includegraphics[width=0.45\linewidth]{mtl.pdf}}
    \label{4c}
	  \caption{Performance against sparsity.}
	  \label{fig:preformance_sparsity}
	\end{figure*}

\paragraph*{Simulation results for benchmark network data}

The average arrowhead error of the estimated GBNs under different levels of edge perturbation can be found in the following Table \ref{tab:error}. Table \ref{tab:error} shows the result based on $5\%$, $10\%$, $20\%$, and $30\%$ of possible edges that are randomly modified during simulation. From the comparison, we demonstrate that MTL is better than SIG and AVG all the time. We also notice that the more related the task share, the better results the proposed method will produce based on different perturbation levels. However, each task with a $30\%$ perturbation level MTL still has decent error rates when SIG and AVG do not provide good results.

\begin{table}[]
\centering
\caption{Error under different network perturbation level }
\begin{tabular}{@{}cccccc@{}}
\toprule
 &  & 5\% Perturbation & 10\% Perturbation & 20\% Perturbation & 30\% Perturbation \\ \midrule
 & SIG & 0.048 & 0.079 & 0.085 & 0.098 \\
ECOLI70 & AVG & 0.047 & 0.073 & 0.082 & 0.095 \\
 & MTL & 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 \\
 & MTL & 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 \\
 & MTL & 0.065 & 0.080 & 0.090 & 0.110 \\ \bottomrule
\end{tabular}
\label{tab:error}
\end{table}

\section{Applications to rs-fMRI data}
\label{sec:applications}

Major depressive disorder (MDD) is a common mental disease that is characterized by low mood, worthlessness, anxiety, and cognitive impairments \citep{fu2008pattern}. The current understanding of MDD is relatively limited, but recent studies have found that MDD is closely related to the changes in the brain connection network \citep{yamashita2020generalizable}, especially the changes in the default mode network (DMN). Several previous studies have reported abnormalities of the DMN in MDD and found some altered functional connections \citep{greicius2007resting,bluhm2009resting,sheline2010resting}].

Such disease-related DMN brain 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 connectivity patterns in MDD patients. This section will conduct such experiments and apply the proposed MTGBN method to explore connectivity reduction patterns in DMN that can facilitate MDD diagnosis.

We consider each MDD participant as a related task within the MDD group, and the patterns learned will be used to contrast with the HC group, where each HC participant is considered as a related task, respectively. We compare the differences in learning performance between the MTGBN, STL, and AVG based on real-world rs-fMRI data. In addition, we also compare the MTGBN method with the existing standard analysis methods of rs-fMRI data: Granger causal analysis (GC, \citealp{seth2010matlab}) and dynamic causal model (DCM, \citealp{marreiros2008dynamic}).

\subsection{Data}
For the experiment, we use the rs-fMRI data of $15$ patients with MDD. These data are regarded as $15$ related tasks, and thus we can obtain a brain connection network for each patient. To compare the MDD connectivity network, the networks of $15$ HC, which are considered to be another set of related tasks, are also obtained using the proposed MTGBN method.  Details of data acquisition and selection criteria are shown in Section \ref{sec:data}.

\subsection{Results}

We apply the proposed MTGBN to compare the connection degree between ROI nodes in three different DMN network templates (aDMN \cite{buckner2008brain}, dDMN \cite{chen2017dissociated}, and vDMN \cite{lee2021ventral}). Table \ref{tab:degree} shows the difference in connection degree of $30$ ROI nodes (including in, out, and total connection degree). Here, the connection degree represents the amount of functional connectivity. In aDMN, the entire connection degrees are $292$ and $232$ for MDD and HC. In other words, MDD has $25.7\%$ more effective connectivity than HC. Similarly, In dDMN, MDD has $29.6\%$ more effective connectivity than HC with a total connection degree of $254$ and $196$, respectively. This increase in connectivity in MDD is also consistent with the findings of previous studies \citep{greicius2007resting,lemogne2009search,belleau2015imbalance}. The rise of DMN connectivity is an apparent robust marker of MDD \citep{zamoscik2014increased}, which can also be found in our results.

\begin{table}
\centering
\caption{Node degree (number of connections) of MDD and HC in three DMN network templates}
\begin{tabular}{@{}lccclccc@{}}
\toprule
Node & Total degree & In degree & Out degree & \multicolumn{1}{c}{Node} & Total degree & In degree & Out degree \\ \midrule
\textbf{aDMN(MDD)} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \textbf{aDMN(HC)} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} \\
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 & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & Total & 232 & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} \\
\textbf{dDMN(MDD)} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \textbf{dDMN(HC)} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} \\
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 & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & Total & 196 & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} \\
\textbf{vDMN(MDD)} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \textbf{vDMN(HC)} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} \\
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 & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & Total & 240 & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} \\ \bottomrule
\end{tabular}
\label{tab:degree}
\end{table}

We can also find some significant 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 in aDMN have more outgoing connections than incoming connections, while aDMN8 and aDMN10 have more incoming connections than outgoing connections. It indicates that aDMN4 and aDMN2 have a more significant impact on other nodes, and aDMN8 is mainly affected by other nodes. These findings are also consistent with the physiological discoveries of MDD \citep{li2018brain}.

Figure \ref{fig:subjectnumber} depicts the difference in network connection counts (the number of subjects with the same connection) in three DMN networks of MDD subjects. ROI connections with high connection counts are generally considered universal connections among MDD participants. Many such ROI connections have been noted as neurophysiological diagnostic biomarkers of MDD disease. From Figure \ref{fig:subjectnumber}, we can find some important connections: 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) in the aDMN.

\begin{figure*}
  \centering
  {\includegraphics[width=0.7\textwidth]{Subjecnnumber.pdf}}
  \caption{The network connection counts in three DMN networks of MDD subjects.}
  \label{fig:subjectnumber}
\end{figure*}

To compare the performance between MTL and SIG, we assess the pairwise statistical significance of the significant connectivity between MDD participants and HC participants by hypothesis testing. The p-values of the hypothesis testing are given in Table \ref{tab:pvalue}.

\begin{table}
\centering
\caption{The statistical significance result of the significant edges.}
\begin{tabular}{@{}lcclll@{}}
\toprule
\multicolumn{1}{c}{Conection} & MTL & SIG & Conection & \multicolumn{1}{c}{MTL} & \multicolumn{1}{c}{SIG} \\ \midrule
\textbf{aDMN} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \textbf{vDMN} &  &  \\
aDMN3 $\rightarrow$ aDMN4 & 0.0040 & 0.0186 & vDMN2 $\rightarrow$ vDMN4 & \multicolumn{1}{c}{0.0401} & \multicolumn{1}{c}{0.5830} \\
aDMN4 $\rightarrow$ aDMN5 & 0.0401 & 0.3356 & vDMN9 $\rightarrow$ vDMN4 & \multicolumn{1}{c}{0.0548} & \multicolumn{1}{c}{0.3356} \\
aDMN3 $\rightarrow$ aDMN5 & 0.0186 & 0.1894 & vDMN7 $\rightarrow$ vDMN9 & \multicolumn{1}{c}{0.1894} & \multicolumn{1}{c}{0.5000} \\
aDMN4 $\rightarrow$ aDMN10 & 0.0823 & 0.1648 & vDMN8 $\rightarrow$ vDMN3 & \multicolumn{1}{c}{0.1894} & \multicolumn{1}{c}{0.7202} \\
aDMN9 $\rightarrow$ aDMN2 & 0.1648 & 0.5830 &  &  &  \\
\textbf{dDMN} &  &  &  &  &  \\
dDMN1 $\rightarrow$ dDMN7 & 0.0032 & 0.4346 &  &  &  \\
dDMN2 $\rightarrow$ dDMN4 & 0.0186 & 0.2178 &  &  &  \\
dDMN7 $\rightarrow$ dDMN8 & 0.1894 & 0.1648 &  &  &  \\
dDMN6 $\rightarrow$ dDMN2 & 0.2722 & 0.6714 &  &  &  \\ \bottomrule
\end{tabular}
\label{tab:pvalue}
\end{table}

To be summarized, in the aDMN, the significant functional connections involve seven ROIs. According to MTL and STL methods, the connectivity between aDMN3 and aDMN4 has the most significant increase ($p=0.004$) in MDD participants compared with HC participants, while MTL produces much more significant results. The connectivity between aDMN4 and aDMN5 has the second most significant ($p=0.041$) increase in MDD participants compared with HC participants using MTL, while STL produces a non-significant outcome. In addition, compared with HC participants, MDD participants also show increased connectivity in aDMN3 $\rightarrow$ aDMN5, aDMN4 $\rightarrow$ aDMN10, and aDMN9 $\rightarrow$ aDMN2. According to the definition of DMN subsystem \citep{andrews2010functional}, we can find that aDMN3 $\rightarrow$ aDMN4, aDMN4 $\rightarrow$ aDMN5, and aDMN3 $\rightarrow$ aDMN5 are within-system connectivity in the dMPFC subsystem, whereas the others are inter-system connectivity between the dMPFC and MTL subsystems. The results show the significantly greater connectivity between aDMN3 (dorsal medial prefrontal cortex), aDMN4 (temporal-parietal junction), and aDMN5 (lateral temporal cortex). These observations may reflect the resting-state hyperconnectivity in the dMPFC subsystem in MDD participants. The dMPFC subsystem is considered to play more social reflective roles \citep{andrews2010functional,andrews2014default}, which is related to various cognitive functions, evaluation and expression of negative emotions \citep{etkin2011emotional}. The neuropathology studies in MDD also report abnormalities in the dMPFC subsystem \citep{khundakar2009morphometric,zhao2014brain}. Furthermore, the inter-system connectivity between the dMPFC and MTL subsystems increases compared with HC participants. The increased connectivity between subsystems indicates that MDD participants have stronger self-related neural activity, which may be related to the self-abnormalities in MDD \citep{zhu2017rumination}.

\begin{table}
\centering
\caption{The performance of different methods in identifying important connections.}
\begin{tabular}{@{}lccccc@{}}
\toprule
Significant edges & Participants & MTL  & SIG & GC & DCM \\ \midrule
aDMN3 $\rightarrow$ aDMN4 & MDD & 73\%  & 53\% & 60\% & 67\% \\
 & HC & 13\%  & 20\% & 20\% & 53\% \\
aDMN4 $\rightarrow$ aDMN5 & MDD & 87\%  & 67\% & 73\% & 73\% \\
 & HC & 7\%  & 27\% & 27\% & 67\% \\
dDMN1 $\rightarrow$ dDMN7 & MDD & 67\%  & 56\% & 53\% & 53\% \\
 & HC & 7\%  & 27\% & 40\% & 33\% \\
dDMN2 $\rightarrow$ dDMN4 & MDD & 63\%  & 49\% & 53\% & 53\% \\
 & HC & 17\%  & 27\% & 40\% & 33\% \\
vDMN2 $\rightarrow$ vDMN4 & MDD & 83\%  & 73\% & 57\% & 67\% \\
 & HC & 40\%  & 53\% & 27\% & 69\% \\
vDMN9 $\rightarrow$ vDMN4 & MDD & 79\% & 63\% & 57\% & 63\% \\
 & HC & 53\%  & 43\% & 37\% & 75\% \\ \bottomrule
\end{tabular}
\label{tab:performance}
\end{table}

To compare the performance differences between our method and other existing methods (GC and DCM). To ensure the consistency of measurement standards, we control various methods to learn the same sparsity of the network structure. We compare the performance of different methods in identifying significant connections. Here, the percentage represents the proportion of the number of subjects identified significant edges to the total number of subjects. For the same significant edge, a good method should determine a higher ratio in the MDD and the opposite in HC. The results in Table \ref{tab:performance} show that our MTL method is better than other methods in identifying substantial connections. When the GC and DCM methods cannot significantly distinguish the difference between MDD and HC, our method is still effective.

\section{Conclusion}
\label{sec:conclusion}
Learning Bayesian network structure has received lots of attention and proved to lack learning from insufficient observations. This study proposes a multitask learning method to learn multiple Bayesian networks simultaneously from observations. The proposed MTGBN method is built based on a few recent works that enhances the learning ability of valid Bayesian inference from observations. We improve with the proposed multitask learning technique and treat each subject in a class of observation as a task. The key ingredients in achieving this include the shared prior covariance among all tasks and the efficient MCEM algorithm. We perform extensive and systematic numerical experiments on simulation and real-life data and demonstrate the excellent performance of the proposed method. We further highlight that the proposed method is especially beneficial for learning multiple informative network structures from insufficient datasets with little shared prior knowledge.

\begin{appendix}
\section{Proof of Theorem \ref{Theorem1}}\label{Appendix a}
Here we provide the detailed proof of Theorem \ref{Theorem1}. 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)$ follows an inverse Wishart distribution $I W\left(\Sigma_{h}, \nu_{0}\right)$. Thus, the challenging of derivation here is that after we generate $\Sigma_{i}$ from $\Sigma_{h}$ (i.e., $\mathbb{P}\left(\Sigma_{i} \mid \Sigma_{h}\right)$), and how to generate a sparse Bayesian network $B_{i}$ from a covariance matrix $\Sigma_{i}$ (i.e., $\mathbb{P}\left(B_{i} \mid \Sigma_{i}\right)$).

Fortunately, \cite{dawid1993hyper} construct 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:
\begin{equation}
\begin{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}.\label{eq:Pbsigma}
\end{small}
\end{equation}

In fact, equation \eqref{eq:Pbsigma} 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 normalization 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:
\begin{equation}
\mathbb{P}\left(\widetilde{G}_{i} \mid \Sigma^{h}\right)=h( \Sigma_{h},\widetilde{G}_{i})=\frac{\prod_{C \in \mathcal{C}_{i}}\left|2^{-1} \Sigma_{h C C}\right|^{\left(\nu_{0}+|C|-1\right) / 2} \Gamma_{|C|}^{-1}\left(\frac{\nu_{0}+|C|-1}{2}\right)}{\prod_{S \in \mathcal{S}_{i}}\left|2^{-1} \Sigma_{h S S}\right|^{\left(\nu_{0}+|S|-1\right) / 2} \Gamma_{|S|}^{-1}\left(\frac{\nu_{0}+|S|-1}{2}\right)},
\label{eq:h}
\end{equation}
where $\Gamma(\cdot)$ denotes  Gamma function, $\mathcal{C}$ and $\mathcal{S}$ are the cliques and separators of the decomposable undirected graph  $\widetilde{G}$, $\Sigma_{h C C}$ and $\Sigma_{h S S}$ 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}).
\]

\section{Proof of Theorem \ref{Theorem2}} \label{Appendix b} Here we provide the detailed proof of Theorem \ref{Theorem2}. 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 \cite{geiger1994learning} shows that the conditional distribution of $D$ given the network structure $B$ can be factorized as
\begin{equation}
\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)},
\label{eq:Pdb}
\end{equation}
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 \eqref{eq:Pdb} follows a multivariate normal distribution.

Similar to equation \eqref{eq:Pdb}, 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)},
\]
where $\Pi_{k}$ is 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}$,  $D_{i}^{\Pi_{k}}$ is the dataset $D_{i}$ restricted to $\Pi_{k}$. Therefore, the main fact we need to do is to derive $\mathbb{P}(D_{i}|\Sigma_{h})$. From equation \eqref{eq:Psigma} and \eqref{eq:Pdsigmai}, we have
\begin{align}\nonumber
\mathbb{P}(D_{i},\Sigma_{i}|\Sigma_{h})\nonumber&=\mathbb{P}(D_{i}|\Sigma_{i})\mathbb{P}(\Sigma_{i}|\Sigma_{h})\nonumber\\
&=\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\} .
\label{eq:Pdsigma2}
\end{align}

Moreover, since $\mathbb{P}(D_{i}|\Sigma_{i}){\sim}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)}\label{eq:data_given_sigma_h}.
\]

\section{Proof of Theorem~\ref{Theorem3}} \label{Appendix c}
In this section we provide the detailed proof of Theorem \ref{Theorem3}. 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
\begin{align}
&\mathbb{P}(\mathcal{B},\mathcal{D}|\Sigma_{h})\nonumber\\
&  =\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}\nonumber\\
 & =\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}.
\label{eq:Pbdsigma}
\end{align}

The challenging of derivation in equation \eqref{eq:Pbdsigma} is how to integrate $\Sigma_{i}$ and generate $B_{i}$ and $D_{i}$ from $\Sigma_{h}$. Fortunately, Theorem \ref{Theorem1} 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 \eqref{eq:Pbdsigma} 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 theequation \eqref{eq:Pdsigma2}, 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
\begin{equation}
\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}.
\label{eq:integral}
\end{equation}
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 \eqref{eq:integral}. Therefore,
\begin{align}
&\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}\nonumber \\
= & \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}.\label{eq:integral2}
\end{align}

Comparing $\mathbb{P}(\Sigma_{i}|\Sigma_{h})$ in the equation \eqref{eq:Pbsigma} and $\mathbb{P}(\Sigma_{i}|\mathcal{D}_{i},\Sigma_{h})d\Sigma_{i}$ in equation \eqref{eq:integral2}, we find that they only differ in the parameters of the inverse Wishart distributions. 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 \eqref{eq:h}, 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,
\[
\begin{aligned}
&\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}\nonumber \\
= & \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},
\end{aligned}
\]
and
\begin{small}
\[
\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}}.
\]
\end{small}

\section{Details of HMC Algorithm \label{Appendix d}}
In this section we continue the discussion on the HMC algorithm in Step 5 of Algorithm \ref{mcem}. 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$ (\textbf{Transform}) and inverse-transform from $V$ to $\Sigma_{h}$ (\textbf{Inverse-Transform}).

\textbf{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.
\]

\textbf{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 \ref{sec:computing}, 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:
\begin{align}\nonumber
\log\{\mathbb{P}(V|\mathcal{B},\mathcal{D})\}& =\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}|\nonumber \\
 & \quad+\sum_{i=1}^{p}(m\nu_{0}+p-i+2)v_{ii}+C.
\label{eq:logV}
\end{align}

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
\[
\begin{aligned}
\mathrm{d}g & =\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)\}\\
 & =2\mathrm{tr}\{L^{\mathrm{T}}I_{C}^{\mathrm{T}}(\Sigma^{h}+A)_{CC}^{-1}I_{C}(\mathrm{d}L)\}.
\end{aligned}
\]

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
\[
\begin{aligned}
\mathrm{d}f & =\mathrm{tr}\{(\Sigma_{h}+A)^{-1}(L(\mathrm{d}L)^{\mathrm{T}}+(\mathrm{d}L)L^{\mathrm{T}})\}\\
 & =\mathrm{tr}\{(\mathrm{d}L)^{\mathrm{T}}(\Sigma_{h}+A)^{-1}L\}+\mathrm{tr}\{L^{\mathrm{T}}(\Sigma_{h}+A)^{-1}(\mathrm{d}L)\}\\
 & =2\mathrm{tr}\{L^{\mathrm{T}}(\Sigma_{h}+A)^{-1}(\mathrm{d}L)\}.
\end{aligned}
\]

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}
\]

According to the above definition, the specific solution process and iterative steps of HMC algorithm are shown in Algorithm \ref{alg:HMC}.
\begin{figure}[]
		\renewcommand{\algorithmicrequire}{\textbf{Input:}}
		\renewcommand{\algorithmicensure}{\textbf{Output:}}
		\let\@latex@error\@gobble
		\begin{algorithm}[H]
			\caption{HMC for Sampling $\Sigma_h$}
			\begin{algorithmic}[1]\label{alg:HMC}
				\REQUIRE   Observational datasets $\mathcal{D}$, BN structures $\mathcal{B}$ and initial $\Sigma_{h}^{(0)}$.
				\ENSURE $\{\Sigma_{h}^{(l)}\}_{l=1}^{N}$: $l$ samples of $\Sigma_{h}$.
				\STATE \textbf{Initialization:} The number of iterations $N$, the number of steps $K$, the step length $\epsilon$, and a specified covariance matrix to simulate momentum $M$.
				\STATE $V^{(0)}$ $\leftarrow $ Transform ($\Sigma_{h}^{(0)}$).
				\STATE $\log\{\mathbb{P}(V|\mathcal{B},\mathcal{D})\}$ = equation \eqref{eq:logV}.
				\STATE Calculate $\log\{\mathbb{P}(V^{(0)}|\mathcal{B},\mathcal{D})\}$.
				\FOR{$t=1, \ldots, N$}
					\STATE $\mathbf{p} \leftarrow \mathcal{N} (0, \mathbf{M})$,
					\STATE $V^{(t)} \leftarrow V^{(t-1)}, \tilde{V} \leftarrow V^{(t-1)}, \tilde{\mathbf{p}} \leftarrow \mathbf{p}$,
					\FOR{$i=1, \ldots, K$}
					   \STATE $\tilde{V}, \tilde{\mathbf{p}} \leftarrow$ Leapfrog $(\tilde{V}, \tilde{\mathbf{p}}, \epsilon, \mathbf{M})$.
					\ENDFOR
					\STATE $\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)$,
					\STATE With probability $\alpha, V^{(t)} \leftarrow \tilde{V} \text { and } \mathbf{p}^{(t)} \leftarrow-\tilde{\mathbf{p}}$,
					\STATE  $\Sigma_{h}^{(t)}$ $\leftarrow$ Inverse-Transform ($V^{(t)}$).
				\ENDFOR
			   	\RETURN $\Sigma_{h}^{(1)}, \ldots, \Sigma_{h}^{(N)}$.
				\STATE \textbf{function} Leapfrog $\left(V^{*}, \mathbf{p}^{*}, \epsilon, \mathbf{M}\right)$
				\STATE \quad $\tilde{\mathbf{p}} \leftarrow \mathbf{p}^{*}+(\epsilon / 2) \nabla_{V} \log f\left(V^{*}\right)$,
				\STATE \quad $\tilde{V} \leftarrow V^{*}+\epsilon \mathbf{M}^{-1} \tilde{\mathbf{p}}$,
				\STATE \quad $\tilde{\mathbf{p}} \leftarrow \tilde{\mathbf{p}}+(\epsilon / 2) \nabla_{V} \log f(\tilde{V})$,
				\STATE \quad \textbf{return} $\tilde{V}, \tilde{\mathbf{p}}$.
				\STATE \textbf{end function}
			\end{algorithmic}
		\end{algorithm}
\end{figure}

\section{ROI information of DMN networks}

\begin{table}[]
\centering
\caption{Parts of ROIs in the three DMN networks and their MNI coordinates }
\begin{tabular}{@{}cccccclllllll@{}}
\toprule
Region name & Abbreviation & Broadman Areas & x & y & z   \\ \midrule
\multicolumn{1}{l}{\textbf{aDMN (Andrews-Hanna Default Mode Networks)}}    \\
\multicolumn{1}{l}{\textit{\textbf{P (PCC-aMPFC Core)}}}   \\
Anterior medial prefrontal cortex & aDMN1 & 10, 32 & -6 & 52 & -2   \\
Posterior cingulate cortex & aDMN2 & 23, 31 & -8 & -56 & 26  \\
\multicolumn{1}{l}{\textit{\textbf{D(dMPFC Subsystem)}}} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} \\
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  \\
\multicolumn{1}{l}{\textit{\textbf{M(MTL Subsystem)}}} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{}  \\
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 &  \\
\multicolumn{1}{l}{\textbf{dDMN (Dorsal Default Mode Networks)}} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} \\
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   \\
\multicolumn{1}{l}{\textbf{vDMN (Ventral Default Mode Networks)}} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{} & \multicolumn{1}{l}{}  \\
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   \\ \bottomrule
\end{tabular}
\label{tab: roi}
\end{table}

\end{appendix}

\newpage

\bibliographystyle{imsart-nameyear}

redpajama_set_name: RedPajamaArXiv