# Latent Multimodal Functional Graphical Model Estimation

Katherine Tsai  Department of Electrical and Computer Engineering, University
of Illinois Urbana-Champaign  Boxin Zhao  Booth School of Business, University
of Chicago  Sanmi Koyejo  Department of Computer Science, Stanford University
Mladen Kolar  Booth School of Business, University of Chicago

###### Abstract

Joint multimodal functional data acquisition, where functional data from
multiple modes are measured simultaneously from the same subject, has emerged
as an exciting modern approach enabled by recent engineering breakthroughs in
the neurological and biological sciences. One prominent motivation to acquire
such data is to enable new discoveries of the underlying connectivity by
combining multimodal signals. Despite the scientific interest, there remains a
gap in principled statistical methods for estimating the graph underlying
multimodal functional data. To this end, we propose a new integrative
framework that models the data generation process and identifies operators
mapping from the observation space to the latent space. We then develop an
estimator that simultaneously estimates the transformation operators and the
latent graph. This estimator is based on the partial correlation operator,
which we rigorously extend from the multivariate to the functional setting.
Our procedure is provably efficient, with the estimator converging to a
stationary point with quantifiable statistical error. Furthermore, we show
recovery of the latent graph under mild conditions. Our work is applied to
analyze simultaneously acquired multimodal brain imaging data where the graph
indicates functional connectivity of the brain. We present simulation and
empirical results that support the benefits of joint estimation.

Keywords: integrative analysis; multimodal data; functional Gaussian graphical
model; Neighborhood regression

##  1 Introduction

Recent engineering breakthroughs have enabled new ways to acquire rich
multimodal data from individual subjects. For example, high-throughput
sequencing enables the acquisition of genotype, and gene expression, among
other signals, from the same set of subjects (hao2021integrated). Our work is
motivated by emerging technology that simultaneously acquires data from
electroencephalogram (EEG) and functional magnetic resonance imaging (fMRI)
(morillon2010neurophysiological). FMRI and EEG data are multivariate time
series, where, after standard preprocessing (wirsich2020concurrent), each
dimension represents a region of the brain. We study the case where both the
EEG and fMRI data are parcellated into the same atlas, resulting in the same
number of dimensions. We view the time series of each region as a function of
time, namely the functional data (ramsay2005fitting), given the continuous
underlying brain signals and the high sampling rate of the measurements. Our
goal is to estimate the functional connectivity of the brain by solving a
graph estimation problem (qiu2016joint; qiao2019functional).

There are several challenges in integrating multimodal functional data for our
problem. First, estimation is often performed in a high-dimensional setting,
where the ambient dimension is much larger than the sample size. Furthermore,
multimodal data are often highly correlated across modes, since they measure
related features of the same subject. Thus, many high-dimensional methods,
such as graphical lasso (yuan2007model), are not easily applied, as they
require restrictions on correlations to achieve the desired statistical
properties. Furthermore, in our setting of interest, the EEG and fMRI data are
noisy and confounded by non-neural activity (murphy2013resting; goto2015head).
Therefore, the graphs estimated by either modality alone are inaccurate. Due
to the observed structured noise and confounders, the graphs estimated
separately from the two modalities may be dissimilar, although they contain
partial information about the same underlying brain network
(wirsich2020multi). To address these shortcomings, we present a novel
generative framework for observed processes (EEG-fMRI measurements) and latent
processes (the underlying brain networks). Then, we provide a framework to
jointly estimate the inverse operator mapping from the observed to the latent
space, along with the latent graph, which encodes conditional independence.

Our framework uses functional graphical models (qiao2019functional;
zhao2021high) as a building block. We study the setting in which multimodal
temporal data are viewed as functions belonging to different Hilbert spaces or
subspaces of the same Hilbert space. This allows us to encode temporal
characteristics using functional scores, a vector of real numbers obtained by
projection to specific bases, effectively circumventing the need to address
temporal discrepancies between modalities; for example, the bottom right of
Figure 1 indicates that fMRI and EEG recordings have distinct frequency
characteristics, belonging to different subspaces of the space of continuous
functions. Our estimation approach handles such a challenging setting by first
constructing linear operators that transform the functional data from the
observation spaces to a shared latent space. Then, we estimate the latent
graph using functional neighborhood regression. Our algorithm jointly
estimates the graphical model and the linear operators through an alternating
iterative procedure.

Our work offers several contributions. From the modeling perspective, we
propose an integrative function-on-function regression framework to estimate a
latent functional Gaussian graphical model. We develop a novel initialization
method inspired by the equivalence of the maximum log-likelihood estimator of
the latent model and canonical correlation analysis (bach2005probabilistic),
when there are two data modalities. While our focus is on two data modalities,
the proposed initialization method can be generalized to $M\geq 2$ by merging
the data modalities, as discussed in Appendix J.3. Subsequently, we devise an
efficient iterative method that boasts a linear rate of convergence.
Specifically, our optimization analysis draws inspiration from recent advances
in nonconvex optimization (zhang2018unified). Regarding the theoretical
contributions, we demonstrate that, under mild conditions, we can recover the
latent graph with high probability. We showcase the effectiveness of the model
through simulations and the analysis of concurrent EEG-fMRI. Empirical results
indicate improved prediction scores when using the estimated latent graph.

Figure 1: An illustration of the generative processes. Top left: The latent
graph of interest. Bottom left: The diagram of the generative model. Top
right: Observed processes from multiple sources. Simultaneous fMRI recordings,
$\mathcal{Y}^{1}_{i}(t)$, and EEG recordings, $\mathcal{Y}_{i}^{2}(t)$, from
morillon2010neurophysiological. Bottom right: The spectrum of fMRI dominates
in the lower frequency domain while that of EEG spans to a higher frequency
domain.

##  2 Methodology

This section introduces a generative model for multimodal data. As our goal is
to estimate the latent graphical model using regression, we establish the
relation between function-on-function regression and conditional independence
in the Gaussian setting. We rigorously define the partial covariance operator
and prove its equivalence to conditional independence.

###  2.1 Notation

Given a separable Hilbert space $\mathbb{H}$ of continuous functions with a
Complete OrthoNormal System (CONS): $\\{\phi_{\ell}:\ell\in\mathbb{N}\\}$, the
inner product is defined as
$\langle{f},{g}\rangle_{\mathbb{H}}=\sum_{i=1}^{\infty}\langle{f},{\phi_{i}}\rangle\langle{g},{\phi_{i}}\rangle$,
where
$\langle{f},{\phi_{i}}\rangle=\int_{t\in{\mathcal{T}}}f(t)\phi_{i}(t)dt$. The
induced norm is defined as
$\|f\|_{\mathbb{H}}=(\langle{f},{f}\rangle_{\mathbb{H}})^{1/2}$ for any
$f\in\mathbb{H}$. We let $\mathbb{H}^{p}$ be the Cartesian product of
$\mathbb{H}^{(1)}\times\cdots\times\mathbb{H}^{(p)}$ where
$\mathbb{H}^{(j)}\equiv\mathbb{H}$ for all $j=1,\ldots,p$ and its inner
product is defined as
$\langle{f},{g}\rangle_{\mathbb{H}^{p}}=\sum_{i=1}^{p}\langle{f_{i}},{g_{i}}\rangle_{\mathbb{H}}$
for any $f,g\in\mathbb{H}^{p}$. Let
$\mathfrak{B}(\mathbb{H}_{1},\mathbb{H}_{2})$ be the class of linear bounded
operators from $\mathbb{H}_{1}$ to $\mathbb{H}_{2}$. A linear operator
$\mathscr{K}\in\mathfrak{B}(\mathbb{H}_{1},\mathbb{H}_{2})$ is in the
equivalence class of a zero operator if
$\mathscr{K}\in\\{\mathscr{A}\in\mathfrak{B}(\mathbb{H}_{1},\mathbb{H}_{2}):\langle{f},{\mathscr{K}g}\rangle=0,\;\forall
g\in\mathbb{H}_{1},f\in\mathbb{H}_{2}\\}$. For any
$\mathscr{T}\in\mathfrak{B}(\mathbb{H}_{1},\mathbb{H}_{2})$, we define the
Hilbert-Schmidt norm as
$\left\|\mathscr{T}\right\|_{\text{HS}}^{2}=\sum_{\ell=1}^{\infty}\|\mathscr{T}\phi_{1,\ell}\|_{\mathbb{H}_{2}}^{2}$,
where $\\{\phi_{1,\ell}:\ell\in\mathbb{N}\\}$ is a CONS for $\mathbb{H}_{1}$
and $|\\!|\\!|\mathscr{T}|\\!|\\!|=\sup_{f\in\mathbb{H}_{1},f\neq
0}\|\mathscr{T}f\|_{\mathbb{H}_{2}}/\|f\|_{\mathbb{H}_{1}}$. $\mathscr{T}$ is
a Hilbert-Schmidt operator if $\left\|\mathscr{T}\right\|_{\text{HS}}<\infty$.
Given an operator $\mathscr{T}$, $\operatorname{Im}(\mathscr{T})$ is the image
space of $\mathscr{T}$. Let $\chi$ be a random element taking values in a
measurable space $(\mathbb{H},\mathscr{B}(\mathbb{H}))$, where
$\mathscr{B}(\cdot)$ denotes the Borel $\sigma$-algebra. A random element
$\chi$ in $L^{2}(\Omega,\mathscr{F},\mu)$ satisfies
$\mathbb{E}\|\chi\|_{\mathbb{H}}^{2}=\int_{\Omega}\|\chi\|_{\mathbb{H}}^{2}d\mu<\infty$.
Let $x_{1},y\in\mathbb{H}_{1}$ and $x_{2}\in\mathbb{H}_{2}$, we define the
tensor product $(x_{2}\otimes x_{1}):\mathbb{H}_{1}\rightarrow\mathbb{H}_{2}$
as $(x_{2}\otimes x_{1})y=\langle{x_{1}},{y}\rangle_{\mathbb{H}_{1}}x_{2}$.
The mean element is defined as $\mathfrak{m}=\int_{\Omega}\chi d\mu$ and the
covariance operator $\mathscr{K}\in\mathfrak{B}({\mathbb{H}},{\mathbb{H}})$
for $\chi$ is defined as
$\mathscr{K}=\mathbb{E}[(\chi-\mathfrak{m})\otimes(\chi-\mathfrak{m})]=\int_{\Omega}(\chi-\mathfrak{m})\otimes(\chi-\mathfrak{m})d\mu$.
Given a matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$ and for $p,q\geq 1$, we
define
$\|\mathbf{A}\|_{p,q}=(\sum_{j=1}^{n}(\sum_{i=1}^{m}|a_{ij}|^{p})^{q/p})^{1/q}$,
$\|\mathbf{A}\|_{p}=\sup_{\mathbf{x}\neq
0}\|\mathbf{A}\mathbf{x}\|_{p}/\|\mathbf{x}\|_{p}$ and $\|\mathbf{A}\|_{F}$ as
the Frobenius norm of $\mathbf{A}$. We use $\mathbb{N}$ to denote positive
integers.

###  2.2 The Latent and Observed Processes

Suppose that there are $M$ modalities of functional data, each living in a
Hilbert space $\mathbb{H}_{m}$, $m=1,\ldots,M$. For modality $m$, we observe a
$p$-dimensional random function
$\mathcal{Y}^{m}=(\mathcal{Y}_{1}^{m},\ldots,\mathcal{Y}_{p}^{m})^{\top}\in\mathbb{H}_{m}^{p}$,
where $\mathbb{H}_{m}$ is a separable Hilbert space of continuous functions
defined on a closed interval $\mathcal{T}\subseteq\mathbb{R}$. Additionally,
denote by $\mathbb{H}$ the latent separable Hilbert space of continuous
functions defined on the closed interval $\mathcal{T}$. Assume that
$\mathcal{Y}^{m}$ is driven by some $p$-dimensional latent processes
$\mathcal{Z}=(\mathcal{Z}_{1},\ldots,\mathcal{Z}_{p})\in\mathbb{H}^{p}$ where
$\mathcal{Z}_{i}:L^{2}(\Omega,\mathcal{F},\mu)\rightarrow(\mathbb{H},\mathscr{B}(\mathbb{H}))$
is a centered Gaussian random element. In neuroscience applications, we have
$M=2$ modalities, ${\mathcal{Y}^{m}}$, $m=1,2$, which represent fMRI and EEG
measurements, and $\mathcal{Z}$ is the latent functional brain process.

Motivated by recent findings that EEG and fMRI measurements can be modeled as
linear transformations of brain signals (calhoun2012multisubject;
chen2013dynamic), we define
${\mathscr{L}^{m}}\in\mathfrak{B}({\mathbb{H}},{\mathbb{H}_{m}})$,
$m=1,\ldots,M$, as data generation operators that transform data from latent
space to observed spaces. We assume that the transformation operator is the
same across nodes; hence ${\mathcal{Y}_{i}^{m}}$ can be decomposed as:

$\displaystyle\mathcal{Y}_{i}^{m}=\chi_{i}^{m}+\xi_{i}^{m},\qquad\chi_{i}^{m}={\mathscr{L}^{m}}\mathcal{Z}_{i},\qquad
i\in V,\;m=1,\ldots,M,$ (2.1)

where $V=\\{1,\ldots,p\\}$ is the set of vertices, $\chi_{i}^{m}$ is obtained
by a deterministic transformation of $\mathcal{Z}_{i}\in\mathbb{H}$, the noise
$\xi_{i}^{m}$, with $\mathbb{E}[\xi_{i}^{m}]=0$ and
$\mathbb{E}[\|\xi_{i}^{m}\|_{\mathbb{H}_{m}}^{2}]<\infty$, is the nuisance
independent of $(\mathcal{Z}_{i})_{i\in V}$,
$(\xi_{j}^{m})_{j\in{V\backslash\\{i\\}}}$, and $(\xi_{j}^{m^{\prime}})_{m\neq
m^{\prime},j=1,\ldots,p}$. See Figure 1 for an illustration of the model.

To avoid identifiability issues, we assume the following regularity
conditions.

######  Assumption 1.

Let ${\mathscr{L}^{m}}$, $m=1,\ldots,M$, be compact operators. The set of
eigenfunctions associated with nonzero distinct eigenvalues, i.e.,
$\lambda_{1}>\lambda_{2}>\ldots$, of ${\mathscr{L}^{m}}^{*}{\mathscr{L}^{m}}$
is denoted as $\\{\phi_{\ell}\in\mathbb{H}:\ell\in\mathbb{N}\\}$, where
${\mathscr{L}^{m}}^{*}$ is the adjoint operator of ${\mathscr{L}^{m}}$. We
assume that

$z_{i,\ell}:=\langle{\mathcal{Z}_{i}},{\phi_{\ell}}\rangle_{\mathbb{H}}\sim\mathcal{N}(0,1),\qquad\ell\in\mathbb{N},$
(2.2)

and $z_{i,\ell}$ is not correlated with $z_{i,\ell^{\prime}}$ for
$\ell\neq\ell^{\prime}$.

This assumption is common in factor model analysis (hsing2015theoretical,
Chapter 10.4). Note that for $i\neq j$, $z_{i,\ell}$ and $z_{j,\ell^{\prime}}$
may be correlated. The distinct eigenvalues assumption allows the set of
eigenfunctions to be unique.

Let $V=\\{1,\ldots,p\\}$ be the set of vertices, and let $E=\\{(i,j):i,j\in
V,i\neq j\\}$ be the set of edges. Our goal is to estimate the undirected
latent graph $G=(V,E)$ underlying $\mathcal{Z}$ from the observed processes
$\mathcal{Y}^{m}$, $m=1,\ldots,M$. The conditional independence for functional
graphical models is defined as in qiao2019functional.

######  Definition 2.1.

A centered Gaussian random vector $\mathcal{Z}\in\mathbb{H}^{p}$ follows a
functional graphical model with respect to an undirected graph $G=(V,E)$ if we
have

$\displaystyle\mathcal{Z}_{i}\perp\\!\\!\\!\perp\mathcal{Z}_{j}\mid\mathcal{Z}_{-(i,j)}\text{
if and only if }(i,j)\not\in E,$

where $\mathcal{Z}_{-(i,j)}=\\{\mathcal{Z}_{j^{\prime}}:j^{\prime}\in
V\backslash\\{i,j\\}\\}$ denotes the components of $\mathcal{Z}$ indexed by
$V\backslash\\{i,j\\}$.

qiao2019functional estimate $E$ by studying the inverse covariance operator of
$\mathcal{Z}$, which generally does not exist in infinite dimensional space
and can only be well-approximated under some restrictions on the eigenvalues
of the covariance operator. In contrast, in the multivariate Gaussian setting,
it is well known that conditional covariance is zero if and only if the
partial covariance is zero, which motivates neighborhood regression as an
estimator of conditional independence (meinshausen2006high; peng2009partial).
We develop a neighborhood regression estimator for the functional data
setting, bypassing the need to compute the inverse covariance operator
directly.

###  2.3 Functional Partial Covariance Operator

Although the relation between the partial and conditional correlation has been
implicitly stated in the RKHS setting (fukumizu2009kernel), to our knowledge,
the formal notion of partial covariance operators has not been clearly
established in the literature on functional Gaussian graphical models. The
form of the regression model that establishes the equivalence of conditional
independence is unclear: when does a zero _regression element_ imply
conditional independence? To answer this, we introduce the partial covariance
operator and establish the corresponding equivalence with respect to
conditional independence.

The regression of $\mathcal{Z}_{i}$ on $\mathcal{Z}_{-(i,j)}$, is defined as

$\displaystyle\\{\widetilde{\beta}_{ij^{\prime}}\\}_{j^{\prime}\in
V\backslash\\{i,j\\}}=\mathop{\mathrm{argmin}}_{\begin{subarray}{c}\beta_{ij^{\prime}}\in\mathfrak{B}(\mathbb{H},\mathbb{H})\\\
j^{\prime}\in
V\backslash\\{i,j\\}\end{subarray}}\mathbb{E}\left[\bigg{\|}\mathcal{Z}_{i}-\sum_{j^{\prime}\in
V\backslash\\{i,j\\}}\beta_{ij^{\prime}}\mathcal{Z}_{j^{\prime}}\bigg{\|}_{\mathbb{H}}^{2}\right],$
(2.3)

where $\beta_{ij}\in\mathfrak{B}(\mathbb{H},\mathbb{H})$ is a bounded linear
operator. For any $f,g\in\mathbb{H}$, we define the partial cross covariance
operator $\mathscr{K}_{ij\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}}$
as:

$\left\langle{\mathscr{K}_{ij\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}}f},{g}\right\rangle_{\mathbb{H}}=\mathbb{E}\left[\left\langle{\mathcal{Z}_{i}-\sum_{j^{\prime}\in
V\backslash\\{i,j\\}}\widetilde{\beta}_{ij^{\prime}}\mathcal{Z}_{j^{\prime}}},{g}\right\rangle_{\mathbb{H}}\left\langle{\mathcal{Z}_{j}-\sum_{j^{\prime\prime}\in
V\backslash\\{i,j\\}}\widetilde{\beta}_{jj^{\prime\prime}}\mathcal{Z}_{j^{\prime\prime}}},{f}\right\rangle_{\mathbb{H}}\right],$
(2.4)

that is, the cross covariance operator of $\mathcal{Z}_{i}$ and
$\mathcal{Z}_{j}$ after removing the effect of $\mathcal{Z}_{-(i,j)}$. The
partial cross covariance operator
$\mathscr{K}_{ij\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}}$ is a
bounded linear operator, which can be shown using the same technique as in
Lemma A.7. We show the following result.

######  Theorem 2.2.

The partial cross covariance operators
$\mathscr{K}_{ij\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}}$ and
$\mathscr{K}_{ji\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}}$ are in the
equivalence class of the zero operator if and only if
$\mathcal{Z}_{i}\perp\\!\\!\\!\perp\mathcal{Z}_{j}\mid\mathcal{Z}_{-(i,j)}$
or, equivalently, if $(i,j)\not\in E$.

This result allows us to estimate the partial covariance operator to measure
the conditional independence of two nodes given the remaining nodes. An
immediate corollary to Theorem 2.2 is that we can write $\mathcal{Z}_{i}$ as a
linear combination of $\mathcal{Z}_{j}$, $j\in{V\backslash\\{i\\}}$.

######  Corollary 2.3.

There exists $\beta_{ij}^{\star}\in\mathfrak{B}(\mathbb{H},\mathbb{H})$,
$j\in{V\backslash\\{i\\}}$, such that

$\displaystyle\mathcal{Z}_{i}=\sum_{j\in{V\backslash\\{i\\}}}\beta_{ij}^{\star}\mathcal{Z}_{j}+\mathcal{W}_{i},$
(2.5)

and $\mathcal{W}_{i}\in\mathbb{H}$ is a Gaussian random function independent
of $\mathcal{Z}_{j}$, $j\in{V\backslash\\{i\\}}$. Furthermore,
$\beta_{ij}^{\star}$ is in the equivalence class of the zero operator if and
only if
$\mathcal{Z}_{i}\perp\\!\\!\\!\perp\mathcal{Z}_{j}\mid\mathcal{Z}_{-(i,j)}$.

Based on Corollary 2.3, specifically the relationship (2.5), we can apply
functional regression to measure conditional independence. However, since the
space of bounded linear operators $\mathfrak{B}(\mathbb{H},\mathbb{H})$ is
infinite dimensional, it is computationally intractable to estimate the
parameters in (2.5). Thus, it is important to study a subclass of models that
can be well-approximated by finite-rank operators. As a consequence, we study
the case where the true regression operators $\beta_{ij}^{\star}$ are compact
and $\|\beta_{ij}^{\star}\|_{\text{HS}}<\infty$.

######  Assumption 2.

The operators $\beta^{\star}_{ij}$, $i,j\in V$ and $i\neq j$, are Hilbert-
Schmidt operators.

Our formulation of functional neighborhood regression closely follows the
formulation in zhao2021high. Since
$\beta_{ij}^{\star}\in\mathfrak{B}_{\text{HS}}(\mathbb{H},\mathbb{H})$ is a
Hilbert-Schmidt operator, it admits a singular system (hsing2015theoretical,
Theorem 4.3.1). Hence, we are able to represent the singular system in terms
of $\\{\phi_{\ell},\ell\in\mathbb{N}\\}$ and obtain

$\beta_{ij}^{\star}=\sum_{\ell,\ell^{\prime}=1}^{\infty}b_{ij,\ell\ell^{\prime}}^{\star}\phi_{\ell}\otimes\phi_{\ell^{\prime}},\quad
i,j\in V,i\neq j,$ (2.6)

where $b_{ij,\ell\ell^{\prime}}^{\star}\in\mathbb{R}$. Under Theorem 2.2, one
can verify that
$\mathcal{Z}_{i}\perp\\!\\!\\!\perp\mathcal{Z}_{j}\mid\mathcal{Z}_{-(i,j)}$ if
and only if $\|\beta_{ij}^{\star}\|_{HS}=0$. Therefore, the set of neighbors
of $\mathcal{Z}_{i}$ is defined as

$\displaystyle\mathscr{N}_{i}^{\star}=\\{j\in{V\backslash\\{i\\}}:\|\beta_{ij}^{\star}\|_{HS}>0\\}.$
(2.7)

In the high-dimensional setting, where $p$ is larger than the sample size $N$,
we often consider that the coefficients are sparse (meinshausen2006high),
rendering a sparse network structure. Such structure is often observed in
biological and clinical experiments (tsai2022joint). In the functional data
scenario, a sparse functional network is equivalent to assuming that many of
$\beta_{ij}^{\star}$ are zero operators:

######  Assumption 3.

The set $\\{\beta_{ij}^{\star}\\}_{j\in{V\backslash\\{i\\}}}$ has at most
$s^{\star}$ nonzero operators.

Assumption 3 implies that the neighborhood $\mathscr{N}_{i}^{\star}$, defined
in (2.7), satisfies $\left|\mathscr{N}_{i}^{\star}\right|\leq s^{\star}$.

For $\beta_{ij}^{\star}$ not a zero operator, we expect that
$\beta_{ij}^{\star}$ can be well approximated by some finite-rank operator:
consider a $k$-dimensional subspace of the original space, denoted as
$\\{\phi_{\ell}:\ell=1,\ldots,k\\}$, we want
$\beta_{ij}^{k\star}:=\sum_{\ell,\ell^{\prime}=1}^{k}b_{ij,\ell\ell^{\prime}}^{k}\phi_{\ell}\otimes\phi_{\ell^{\prime}}\approx\beta_{ij}^{\star}$.
Observe that
$\|\beta_{ij}^{\star}\|_{HS}=\|\beta_{ij}^{\star}-\beta_{ij}^{\star{k}}+\beta_{ij}^{\star{k}}\|_{HS}\leq\|\beta_{ij}^{\star}-\beta_{ij}^{\star{k}}\|_{HS}+\|\beta_{ij}^{\star{k}}\|_{HS}$.
Hence, if there exists a positive $\epsilon_{0}$ such that $\min_{i\in
V}\min_{j\in\mathscr{N}_{i}^{\star}}\|\beta_{ij}^{\star}\|_{HS}\geq
2\epsilon_{0}$, then we can choose a ${k}$ such that the truncated signal
still has half the magnitude of the original signal:
$\|\beta_{ij}^{{k\star}}\|_{HS}\geq\epsilon_{0}$. We define

$\Lambda(k):=\min_{i\in
V}\min_{j\in\mathscr{N}_{i}^{\star}}\|\beta_{ij}^{k\star}\|_{HS}.$

In addition, for $\beta_{ij}^{\star}$ not a zero operator, we expect that most
of $b_{ij,\ell\ell^{\prime}}^{\star}$, $\ell,\ell^{\prime}\in\mathbb{N}$, are
zero or have a small magnitude. We make the following assumption on
$\beta^{\star}_{ij}$.

######  Assumption 4.

The rank $k$ is selected such that $\Lambda(k)\geq(1/2)\min_{i\in
V}\min_{j\in\mathscr{N}_{i}^{\star}}\|\beta_{ij}^{\star}\|_{HS}=:\epsilon_{0}$.
There exists a constant $\alpha^{\star}\in(0,1]$ such that the sets
$\\{b_{ij,\ell 1}^{\star},\ldots,b_{ij,\ell k}^{\star}\\}$ and
$\\{b_{ij,1\ell^{\prime}}^{\star},\ldots,b_{ij,k\ell^{\prime}}^{\star}\\}$
each have at most $\alpha^{\star}k$ non-zero coefficients, for
$\ell,\ell^{\prime}=1,\ldots,k$.

The first assumption implies that (2.7) is equivalent as

$\displaystyle\mathscr{N}_{i}^{\star}=\\{j\in{V\backslash\\{i\\}}:\|\beta_{ij}^{k\star}\|_{HS}>\epsilon_{0}\\}.$
(2.8)

The second structural assumption implies that each basis function is only
correlated with at most $\alpha k-1$ other basis functions but we allow the
rank to grow with $k$. From the neuroscience perspective
(olsen2012hippocampus; vidaurre2017brain), the first level of sparsity,
Assumption 3, corresponds to the sparse connectivity between sub-networks of
the whole connectome. We consider $\phi_{\ell}$ to be the $\ell$-th cognitive
process (posner1988localization), e.g., visual imaging, word reading, shifting
visual attention, and etc. Assumption 4 implies two connectivity structures,
as outlined in the following. First, each node $i$ only has a few ongoing
cognitive processes. Secondly, each connectivity process $\ell$ in node $i$ is
sparsely correlated with the other cognitive process $\ell^{\prime}$ in node
$j$ (park2013structural; vidaurre2017brain). This assumption further implies
that the operator $\beta_{ij}^{\star}$ is approximately low-rank, because many
of the $\\{b_{ij,\ell\ell^{\prime}}^{\star}\\}$ are zero. Additionally, the
variability of sparsity patterns imposed by Assumption 3–4 are subject to the
cognitive tasks performed. For example, cohen2016segregation have found that
during the motor cognitive process, within sub-network connectivity is highly
activated, meaning that $\alpha^{\star}$ may be high. In contrast, during the
memory task, between-network communication is more active, meaning that
$s^{\star}$ might be high.  While the former assumption is well recognized in
functional graphical model literature (qiao2019functional), we argue that
adding the second level sparsity assumption provides more informative
interpretability of the underlying functional processes.

##  3 Estimation

We describe the estimation of the parameters for the generative model
described in Section 2. We first introduce an estimation procedure for the
transformation operator, which we formulate as an inverse problem. Then, the
objective function in infinite dimensional space is outlined in Section 3.2.
Finally, the objective function in a finite dimensional space along with the
structured assumptions are discussed in Section 3.3.

###  3.1 Transformations to the Latent Space

Our problem is to estimate the graphical model $G=(V,E)$ in the latent space
through data from multiple modalities
$\mathcal{Y}^{1},\ldots,\mathcal{Y}^{M}$. While the mathematical formulation
in (2.1) outlines the generative model, from the estimation point of view, we
focus on its inverse, that is, the linear operators that transform the
functional processes from the observation spaces to the latent space. In
general, the inverse operator of ${\mathscr{L}^{m}}$ does not exist; even if
it does exist, it might not be continuous. However, the inverse operator is
well defined when the operator domain is restricted to
${\operatorname{Im}({\mathscr{L}^{m}})}$, giving us hope to recover the
inverse of ${\mathscr{L}^{m}}$ on the restricted domain. Following Definition
3.5.7 in hsing2015theoretical, we define the Moore-Penrose (generalized)
inverse of ${\mathscr{L}^{m}}$.

######  Definition 3.1.

Let $\widetilde{\mathscr{L}}^{m}$ be the operator ${\mathscr{L}^{m}}$
restricted to $\ker({\mathscr{L}^{m}})^{\perp}$, $m=1,\ldots,M$. The domain of
${\mathscr{L}^{m}}^{\dagger}$ is defined as
$\operatorname{Dom}({\mathscr{L}^{m}}^{\dagger})=\operatorname{Im}({\mathscr{L}^{m}})\oplus\operatorname{Im}({\mathscr{L}^{m}})^{\perp}$.
Then, for any $y\in\operatorname{Dom}({\mathscr{L}^{m}}^{\dagger})$,
${\mathscr{L}^{m}}^{\dagger}y=(\widetilde{\mathscr{L}}^{m})^{-1}y$ if
$y\in\operatorname{Im}({\mathscr{L}^{m}})$ and
${\mathscr{L}^{m}}^{\dagger}y=0$ if
$y\in\operatorname{Im}({\mathscr{L}^{m}})^{\perp}$.

We use ${\mathscr{A}^{m}}$ to denote the pseudoinverse operator
${\mathscr{L}^{m}}^{\dagger}$ for simplicity of notation. Since
${\chi_{i}^{m}}\in\operatorname{Im}({\mathscr{L}^{m}})$, we have the following
relationship:

$\mathcal{Z}_{i}={\mathscr{A}^{m}}(\chi_{i}^{m}),\qquad m=1,\ldots,M,\;i\in
V.$ (3.1)

The operator ${\mathscr{A}^{m}}$ maps the random elements from the observation
spaces back to the latent space, where we can jointly estimate the latent
graph from multiple views $\mathcal{Y}^{1},\ldots,\mathcal{Y}^{M}$.

The following representation of ${\mathscr{A}^{m}}$ will be useful for
estimation from the observed data. Let
$\\{\phi_{\ell}^{m}\in\mathbb{H}_{m}:\ell\in\mathbb{N}\\}$ be the set of
eigenfunctions of $\mathscr{L}^{m}\mathscr{L}^{m*}$. Since
${\mathscr{L}^{m}}\in\mathfrak{B}({\mathbb{H}},{\mathbb{H}_{m}})$ is a compact
operator, when the domain is restricted to
$\operatorname{Im}({\mathscr{L}^{m}})\oplus\operatorname{Im}({\mathscr{L}^{m}})^{\perp}$
we have:

${\mathscr{A}^{m}}=\sum_{\ell=1}^{\infty}\sum_{\ell^{\prime}=1}^{\infty}a_{\ell\ell^{\prime}}^{m\star}\phi_{\ell}\otimes\phi_{\ell^{\prime}}^{m},\quad{\mathscr{A}^{m}}\chi=\sum_{\ell=1}^{\infty}\sum_{\ell^{\prime}=1}^{\infty}a_{\ell\ell^{\prime}}^{m\star}\langle{\chi},{\phi_{\ell^{\prime}}^{m}}\rangle_{\mathbb{H}_{m}}\phi_{\ell},$
(3.2)

where $a_{\ell\ell^{\prime}}^{m\star}\in\mathbb{R}$ are coefficients and
$\chi\in\operatorname{Dom}({\mathscr{A}^{m}})$.

Although ${\mathscr{A}^{m}}$ is well defined on the restricted domain, it is
not necessarily a bounded operator. Since $\mathscr{L}^{m}$ is a compact
operator, to have a bounded inverse ${\mathscr{A}^{m}}$, $\mathscr{L}^{m}$
must be finite dimensional (groetsch1984theory, Corollary 1.3.3). The
assumption is mild in our application because the rank of ${\mathscr{L}^{m}}$
can be arbitrarily large, so long as it is finite.

###  3.2 Approximation of Infinite Dimensional Estimator

When the noise ${\xi_{i}^{m}}$ in (2.1) is small in magnitude compared to
${\chi_{i}^{m}}$, then
$\mathcal{Z}_{i}={\mathscr{A}^{m}}{\chi_{i}^{m}}\approx{\mathscr{A}^{m}}{\mathcal{Y}_{i}^{m}}$.
We combine (2.3) and the estimation of ${\mathscr{A}^{m}}$ and optimize the
objective:

$\min_{\\{{\mathscr{A}^{m}}\\},\\{\beta_{ij}\\}}\qquad\sum_{i=1}^{p}\sum_{m=1}^{M}\mathbb{E}\left[\bigg{\|}{\mathscr{A}^{m}}{\mathcal{Y}_{i}^{m}}-\sum_{j\in{V\backslash\\{i\\}}}\beta_{ij}{\mathscr{A}^{m}}\mathcal{Y}_{m,j}\bigg{\|}_{\mathbb{H}}^{2}\right],$
(3.3)

where ${\mathscr{A}^{m}}$ is of the form (3.2) and $\beta_{ij}$ is of the form
(2.6).

The temporal realization of random elements $\mathcal{Z}_{i}(\omega,\cdot)$
and ${\mathcal{Y}_{i}^{m}}(\omega,\cdot)$, $\omega\in\Omega$, are bivariate
functions $\mathcal{Z}_{i}(\omega,t)$ and ${\mathcal{Y}_{i}^{m}}(\omega,t)$
that we assume are jointly measurable with respect to the product
$\sigma$-field $\mathcal{F}\times\mathscr{B}({\mathcal{T}})$. Since both
$\mathcal{Z}_{i}$ and ${\mathcal{Y}_{i}^{m}}$ are mean-squared integrable
random elements taking values in spaces of continuous functions, by the
Karhunen–Loève theorem (karhunen1947uber; loeve1945fonctions), we can uniquely
represent them as

$\mathcal{Z}_{i}(\omega,t)=\sum_{\ell=1}^{\infty}z_{i,\ell}(\omega)\phi_{\ell}(t),\qquad{\mathcal{Y}_{i}^{m}}(\omega,t)=\sum_{\ell=1}^{\infty}y_{i,\ell}^{m}(\omega)\phi_{\ell}^{m}(t),\qquad\omega\in\Omega,\;t\in{\mathcal{T}},$
(3.4)

where we recall that $z_{i,\ell}(\omega)$ is defined in (2.2) and
$y_{i,\ell}^{m}(\omega)=\langle{{\mathcal{Y}_{i}^{m}}(\omega,\cdot)},{\phi_{\ell}^{m}(\cdot)}\rangle_{\mathbb{H}_{m}}$.
For conciseness, we omit $\omega$ from the notation on
$y_{i,\ell}^{m}(\omega)$ and $z_{i,\ell}(\omega)$, and $t$ on $\phi_{\ell}(t)$
and $\phi_{\ell}^{m}(t)$.

######  Proposition 3.2.

The optimization problem in (3.3) is equivalent to

$\min_{\\{a_{\ell\ell^{\prime}}^{m}\\},\\{b_{ij,\ell\ell^{\prime}}\\}}\qquad\sum_{m=1}^{M}\sum_{i=1}^{p}\mathbb{E}\left[\sum_{\ell\in\mathbb{N}}\left(\sum_{\ell^{\prime\prime}\in\mathbb{N}}a_{\ell\ell^{\prime\prime}}^{m}y_{i,\ell^{\prime\prime}}^{m}-\sum_{j\in{V\backslash\\{i\\}}}\sum_{\ell^{\prime},\ell^{\prime\prime}\in\mathbb{N}}b_{ij,\ell\ell^{\prime}}a_{\ell^{\prime}\ell^{\prime\prime}}^{m}y_{j,\ell^{\prime\prime}}^{m}\right)^{2}\right],$
(3.5)

where
${\mathscr{A}^{m}}=\sum_{\ell,\ell^{\prime}}a_{\ell\ell^{\prime}}^{m}\phi_{\ell}\otimes\phi_{\ell^{\prime}}^{m}$
and
$\beta_{ij}=\sum_{\ell,\ell^{\prime\prime}}b_{ij,\ell\ell^{\prime\prime}}\phi_{\ell}\otimes\phi_{\ell^{\prime\prime}}$.

Based on Proposition 3.2, the optimization over the operators
${\mathscr{A}^{m}}$ and $\beta_{ij}$ in (3.3) is transformed to optimization
over sequences of real numbers
$\\{a_{\ell\ell^{\prime}}^{m}\\}_{\ell,\ell^{\prime}\in\mathbb{N}}$ for
$m=1,\ldots,M$ and
$\\{b_{ij,\ell\ell^{\prime}}\\}_{\ell,\ell^{\prime}\in\mathbb{N}}$ for $i\in
V$ and $j\in{V\backslash\\{i\\}}$. While functional data live in infinite
dimensional spaces, in practice, one is faced with finite computational
resources and is forced to truncate both the latent functions and observed
functions to perform estimation. In addition, Assumption 3–4 indicates that
most of $b_{ij,\ell\ell^{\prime}}^{\star}$ are zero. Hence, using the result
from (3.4), the ${k}$-dimensional approximation of $\mathcal{Z}_{i}$ via
projection to $\\{\phi_{1},\ldots,\phi_{k}\\}$ is
$\sum_{\ell=1}^{k}z_{i,\ell}\phi_{\ell}$. For the observed space, we assume
that the signals are smooth enough such that the $k_{m}$-dimensional
projection $\sum_{\ell=1}^{k_{m}}y_{i,\ell}^{m}\phi_{\ell}^{m}$ can well
approximate ${\mathcal{Y}_{i}^{m}}$. These observations imply that the
following optimization problem can recover the model parameters sufficiently
well:

$\mathbf{A}^{\star},\mathbf{B}^{\star}:=\arg\min_{\mathbf{A},\mathbf{B}}\sum_{m=1}^{M}\sum_{i=1}^{p}\mathbb{E}\left[\bigg{\|}{\mathbf{A}^{m}}\mathbf{y}_{i}^{m}-\sum_{j\in{V\backslash\\{i\\}}}\mathbf{B}_{ij}{\mathbf{A}^{m}}\mathbf{y}_{j}^{m}\bigg{\|}_{F}^{2}\right],$
(3.6)

where $\mathbf{A}^{\star}=\\{\mathbf{A}^{m\star}\\}_{m=1}^{M}$,
$\mathbf{B}^{\star}=\\{\mathbf{B}_{i}^{\star}\\}_{i=1}^{p}$,
$\mathbf{y}_{i}^{m}=(y_{i,1}^{m},\ldots,y_{i,k_{m}}^{m})^{\top}\in\mathbb{R}^{k_{m}}$,
${\mathbf{A}^{m}}=(a_{\ell\ell^{\prime}}^{m})_{\ell=1,\ell^{\prime}=1}^{k,k_{m}}$,
and
$\mathbf{B}_{ij}=(b_{ij,\ell\ell^{\prime}})_{\ell=1,\ell^{\prime}=1}^{k,k}\in\mathbb{R}^{k\times
k}$. The parameters $k$ and $k_{m}$ are suitably chosen – we discuss the
selection of $k$ and $k_{m}$ in more detail in Section 4.2.

###  3.3 Empirical Objective Function and Structured Assumptions

Throughout the section, we assume the sets of basis
$\\{\phi_{1}^{m},\ldots,\phi_{k_{m}}^{m}\\}$, $m=1,\ldots,M$ are provided. We
are given independent observations from $N$ individuals, where for each
individual $n$ and each vertex $i$, we have data from $M$ modalities,
$\\{\mathcal{Y}_{i}^{m,(n)}\\}_{m=1}^{M}$. The vector $\mathbf{y}_{i}^{m,(n)}$
is obtained by projecting $\mathcal{Y}_{i}^{m,(n)}$ onto the set of basis
$\\{\phi_{1}^{m},\ldots,\phi_{k_{m}}^{m}\\}$, and we let
${\mathbf{Y}_{i}^{m}}=(\mathbf{y}_{i}^{m,(1)},\ldots,\mathbf{y}_{i}^{m,(N)})\in\mathbb{R}^{k_{m}\times
N}$. The sample version of the objective in (3.6) is:

$\displaystyle f(\mathbf{A},\mathbf{B})$
$\displaystyle=\sum_{i=1}^{p}\sum_{m=1}^{M}\frac{1}{2N}\bigg{\|}{\mathbf{A}^{m}}{\mathbf{Y}_{i}^{m}}-\sum_{j\in{V\backslash\\{i\\}}}\mathbf{B}_{ij}{\mathbf{A}^{m}}{\mathbf{Y}_{j}^{m}}\bigg{\|}_{F}^{2}.$
(3.7)

To enforce the sparsity of $\mathbf{B}$, we define the following norm.

######  Definition 3.3.

Let $\mathbf{U}=(\mathbf{U}_{1}\cdots\mathbf{U}_{p})\in\mathbb{R}^{d\times
kp}$, where $\mathbf{U}_{i}\in\mathbb{R}^{d\times k}$. We define the norm

$\|\mathbf{U}\|_{r(k,\ell)}=\left\|\begin{pmatrix}\operatorname{vec}\left(\mathbf{U}_{1}\right)\cdots\operatorname{vec}\left(\mathbf{U}_{p}\right)\end{pmatrix}\right\|_{2,\ell},\quad\ell\geq
1,$

and let
$\|\mathbf{U}\|_{r(k,0)}=\sum_{i=1}^{p}\mathbbm{1}(\|\operatorname{vec}\left(\mathbf{U}_{i}\right)\|_{2}>0)$,
where $\mathbbm{1}(\cdot)$ denotes the indicator function.

Let
$\mathbf{B}_{i}^{\star}=(\mathbf{B}_{i1}^{\star}\cdots\mathbf{B}_{i(i-1)}^{\star}\mathbf{B}_{i(i+1)}^{\star}\cdots\mathbf{B}_{ip}^{\star})\in\mathbb{R}^{k\times
k(p-1)}$. Under Assumption 3, we have that
$\|\mathbf{B}_{i}^{\star}\|_{r(k,0)}\leq s^{\star}$. Assumption 4 implies that
each nonzero submatrix $\mathbf{B}_{ij}^{\star}$ of $\mathbf{B}_{i}^{\star}$
has columns and rows with at most $\alpha$-fraction of nonzero entries. We
define the following constraint set

$\mathcal{K}_{B}(s,\alpha)=\\{\mathbf{U}:\|\mathbf{U}\|_{r(k,0)}\leq
s,\|\mathbf{U}_{j,m\cdot}\|_{0},\|\mathbf{U}_{j,\cdot n}\|_{0}\leq\alpha
k,m,n=1,\ldots,k,j=1,\ldots,p-1\\},$

where $\mathbf{U}_{j}\in\mathbb{R}^{k\times k}$ denotes the sub-matrix of
$\mathbf{U}\in\mathbb{R}^{k\times k(p-1)}$. We assume that
${\mathbf{A}^{m}}\in\mathbb{R}^{k\times k_{m}}$ satisfies the following
constraints. Let $\mathbf{W}_{j\cdot}$ be the $j$-th row of a matrix
$\mathbf{W}$, we define

$\mathcal{K}_{m}(\tau_{1},\tau_{2})=\left\\{\mathbf{W}:\tau_{1}^{1/2}\sigma_{\min}(\mathbf{A}^{m(0)})\leq\|\mathbf{W}_{j\cdot}\|_{2}\leq(\tau_{2}/k_{m})^{1/2}\sigma_{\max}(\mathbf{A}^{m(0)}),\;j=1,\ldots
k\right\\},$ (3.8)

where $\mathbf{A}^{m(0)}$ is the initial guess. We will later verify that
there exist some $\tau_{1},\tau_{2}>0$ such that $\mathbf{A}^{m\star}$ lies in
$\mathcal{K}_{m}(\tau_{1},\tau_{2})$. The lower bound $\tau_{1}$ ensures that
any row of ${\mathbf{A}^{m}}$ is always bounded away from zero.111We need to
make sure that ${\mathbf{A}^{m}}$ is nonzero so that the solution is
nontrivial, since ${\mathbf{A}^{m}}={\bf 0}$ and $\mathbf{B}_{i}={\bf 0}$ for
$m=1,\ldots,M$ and $i\in V$ is always an optimal solution in the unconstrained
optimization problem. If $\tau_{2}$ is proportional to $k^{1/2}$, then the
upper bound constraint can be viewed as the incoherence condition
(candes2011robust). We want to optimize the following objective function:

$\displaystyle\mathop{\mathrm{argmin}}_{\mathbf{A},\mathbf{B}}$
$\displaystyle\quad f(\mathbf{A},\mathbf{B}),$ (3.9)
$\displaystyle\quad\text{s.t.}\quad\mathbf{B}_{i}\in\mathcal{K}_{B}(s,\alpha),\;i\in
V;$
$\displaystyle\quad\quad\quad\;{\mathbf{A}^{m}}\in\mathcal{K}_{m}(\tau_{1},\tau_{2}),\;m=1,\ldots,M.$

##  4 Algorithm

Algorithm 1 Latent Graph Estimation

Input:
$\\{\mathbf{Y}_{m}\\}_{m=1,\ldots,M}$,$\\{\mathbf{A}^{m(0)}\\}_{m=1,\ldots,M}$,
$\\{\mathbf{B}_{i}^{(0)}\\}_{i\in V}$,$\tau=(\tau_{1},\tau_{2})$, $s$,
$\alpha$, $\eta_{A}$, $\eta_{B}$;

Output: $\\{\widehat{\mathbf{A}}^{m}\\}_{m=1,\ldots,M}$,
$\\{\widehat{\mathbf{B}}_{i}\\}_{i\in V}$

while Not converged do

for m=1,…,M do

$\mathbf{A}^{m(t+1)}\leftarrow\mathcal{P}_{m,\tau}(\mathbf{A}^{m(t)}-\eta_{A}\nabla_{{\mathbf{A}^{m}}}f$)

for $i\in V$ do

$\mathbf{B}_{i}^{(t+.5)}\leftarrow{\mathcal{T}}_{s}(\mathbf{B}_{i}^{(t)}-\eta_{B}\nabla_{\mathbf{B}_{i}}f$)

for $j\in{V\backslash\\{i\\}}$ do

$\mathbf{B}_{ij}^{(t+1)}\leftarrow\mathcal{H}_{\alpha}(\mathbf{B}_{ij}^{(t+.5)})$

We propose a two-stage algorithm to minimize (3.9) efficiently. We start by
introducing the second stage. Given a suitably chosen initial tuple
$(\mathbf{A}^{(0)},\mathbf{B}^{(0)})$ from the first stage, we use the
alternating projected gradient descent with the group-sparse hard-thresholding
operator. We first define the $s$-group sparse truncation operator as

$\displaystyle[{\mathcal{T}}_{s}(\mathbf{B}_{i})]_{\cdot,k(j-1)+1:kj}=\left\\{\begin{array}[]{cc}\mathbf{B}_{ij}\in\mathbb{R}^{k\times
k},&\text{the magnitude $\|\mathbf{B}_{ij}\|_{F}$ is among the top-s of
$\\{\mathbf{B}_{ij}\\}$};\\\ {\bf 0}\in\mathbb{R}^{k\times
k},&\text{otherwise}.\end{array}\right.$

The operator ${\mathcal{T}}_{s}$ keeps the largest $s$ submatrices of
$\mathbf{B}_{i}\in\mathbb{R}^{k\times k(p-1)}$. Next, we define the
$\alpha$-truncation operator as

$[\mathcal{H}_{\alpha}(\mathbf{B}_{ij})]_{u,v}=\left\\{\begin{array}[]{ll}[\mathbf{B}_{ij}]_{u,v},&\text{$[\mathbf{B}_{ij}]_{u,v}$
is one of the $\alpha k$ largest elements in}\\\ &\text{magnitude of
$[\mathbf{B}_{ij}]_{u,\cdot}$ and $[\mathbf{B}_{ij}]_{\cdot,v}$.}\\\
0,&\text{otherwise}.\end{array}\right.$

The operator $\mathcal{H}_{\alpha}$ keeps the largest $\alpha$-fraction of
entries in each row and column of $\mathbf{B}_{ij}\in\mathbb{R}^{k\times k}$.
Then, the projection to the set $\mathcal{K}_{B}(s,\alpha)$ can be implemented
by the composition of ${\mathcal{T}}_{s}$ and $\mathcal{H}_{\alpha}$, as shown
in Algorithm 1. The operator $\mathcal{P}_{m,\tau}$ as the projection operator
to the set $\mathcal{K}_{m}(\tau_{1},\tau_{2})$ is defined as

$\displaystyle\mathcal{P}_{m,\tau}({\mathbf{A}^{m}})=\mathop{\mathrm{argmin}}_{\mathbf{W}\in\mathcal{K}_{m}(\tau_{1},\tau_{2})}\|\mathbf{W}-{\mathbf{A}^{m}}\|_{F}.$

After Algorithm 1 converges, we select the edges of the graph either using the
AND or OR operation with a threshold $\epsilon_{0}>0$:

$\widehat{\mathscr{N}}_{i}(\epsilon_{0})=\\{(i,j):\|\widehat{\mathbf{B}}_{ij}\|_{F}\geq\epsilon_{0}\text{
and (or) }\|\widehat{\mathbf{B}}_{ji}\|_{F}\geq\epsilon_{0}\\}.$ (4.1)

###  4.1 Initialization Procedure

We describe a procedure to find a good initial tuple
$(\\{\mathbf{A}^{m(0)}\\}_{m=1}^{M},\\{\mathbf{B}_{i}^{(0)}\\}_{i=1}^{p})$. We
focus on the case when $M=2$, as this is the case for the motivating example.
To find initial estimates, we connect the problem with probabilistic canonical
correlation analysis. We first construct ${\mathbf{A}^{m}}^{(0)}$ by finding
the canonical correlation between two views and then compute
$\mathbf{B}_{i}^{(0)}$ with an iterative method by fixing ${\mathbf{A}^{m}}$.

Let $\mathbf{L}^{m\star}\in\mathbb{R}^{k_{m}\times k}$ be such that the
$\ell\ell^{\prime}$-th entry of $\mathbf{L}^{m\star}$ is
$L_{\ell\ell^{\prime}}^{m\star}=\langle{{\mathscr{L}^{m}}\phi_{\ell^{\prime}}^{m}},{\phi_{\ell}^{m}}\rangle$.
Let
${\mathbf{q}}_{i}^{m}=(q_{i,1}^{m},\ldots,q_{i,k_{m}}^{m})^{\top}\in\mathbb{R}^{k_{m}}$,
where $q_{i,j}^{m}=\langle{\xi_{i}^{m}},{\phi_{j}^{m}}\rangle$, be the
truncation of $\xi_{i}^{m}$;
${\mu}_{i}^{m}=(\mu_{i,1}^{m},\ldots,\mu_{i,k_{m}}^{m})^{\top}\in\mathbb{R}^{k_{m}}$
be the bias induced by finite-dimensional truncation and
$\mu_{i,\ell^{\prime}}^{m}=\sum_{\ell=k+1}^{\infty}\langle{{\mathscr{L}^{m}}(z_{i,\ell}\phi_{\ell})},{\phi_{\ell^{\prime}}^{m}}\rangle$
for $\ell^{\prime}=1,\ldots,k_{m}$. Recalling that
$\mathbf{z}_{i}\sim\mathcal{N}({\bf 0},\mathbf{I}_{k})$, it follows that

$\displaystyle\mathbf{y}_{i}^{m}|\mathbf{z}_{i}$
$\displaystyle\sim\mathcal{N}(\mathbf{L}^{m\star}\mathbf{z}_{i},\bm{\Sigma}^{m,\mathbf{q}}_{i,i}+\bm{\Sigma}^{m,{\bm{\mu}}}_{i,i}),$
(4.2)

where $\bm{\Sigma}^{m,{\bm{\mu}}}_{i,i}$ is the covariance of $\mu_{i}^{m}$
and $\bm{\Sigma}^{m,{\mathbf{q}}}_{i}$ is the covariance of
$\mathbf{q}_{i}^{m}$. From (2.2) we have that $\bm{\mu}_{i}^{m}$ is
uncorrelated with $\mathbf{z}_{i}$. Recall that we assume that
$\mathbf{L}^{m\star}$ is the same across $i\in V$. We briefly discuss how
extend to the setting where $\mathbf{L}_{i}^{m\star}$ differ across the nodes
in Appendix J.2. To estimate $\mathbf{L}^{m\star}$ from data, we pick the
first node $i=1$ for convenience. We discuss an alternative to this strategy
in Appendix J.1 The log-likelihood is defined as

$\displaystyle\ell({\bf S})=(k_{1}+k_{2})\log 2\pi$
$\displaystyle+\log\left|{\bf S}\right|+\mathop{\mathrm{tr}}({\bf
S}^{-1}\widehat{\bm{\Sigma}}),$ $\displaystyle{\bf
S}=\begin{pmatrix}\mathbf{L}^{1}{\mathbf{L}^{1}}^{\top}+\bm{\Sigma}^{1,\mathbf{q}}_{11}&\mathbf{L}^{1}{\mathbf{L}^{2}}^{\top}\\\
\mathbf{L}^{2}{\mathbf{L}^{1}}^{\top}&\mathbf{L}^{2}{\mathbf{L}^{2}}^{\top}+\bm{\Sigma}^{2,\mathbf{q}}_{11}\end{pmatrix},$
$\displaystyle\quad\widehat{\bm{\Sigma}}=\frac{1}{N}\sum_{n=1}^{N}\begin{pmatrix}\mathbf{y}_{1}^{1,(n)}\mathbf{y}_{1}^{1,(n)\top}&\mathbf{y}_{1}^{1,(n)}\mathbf{y}_{1}^{2,(n)\top}\\\
\mathbf{y}_{1}^{2,(n)}\mathbf{y}_{1}^{1,(n)\top}&\mathbf{y}_{1}^{2,(n)}\mathbf{y}_{1}^{2,(n)\top}\end{pmatrix}.$

Let $\bm{\Sigma}^{m}_{11}$, $\widehat{\bm{\Sigma}}^{m}_{11}$ be the covariance
matrix and sample covariance of node $1$ for modality $m$, respectively. Let
$\bm{\Sigma}^{12}_{11}$, $\widehat{\bm{\Sigma}}^{12}_{11}$ be the cross-
covariance matrix and sample cross-covariance matrix for two modalities at
node $1$, respectively. By Theorem 2 in bach2005probabilistic, the maximum
log-likelihood estimator for $\mathbf{L}^{m\star}$ is
$\widehat{\mathbf{L}}^{m}=\widehat{\bm{\Sigma}}^{m}_{11}\widehat{\mathbf{V}}^{m}\widehat{\bm{\Gamma}}^{1/2}$,
where $\widehat{\bm{\Gamma}}$ is a diagonal matrix whose diagonal entries are
the top-k singular values of
$\widehat{\mathbf{R}}^{12}=(\widehat{\bm{\Sigma}}^{1}_{11})^{-1/2}\widehat{\bm{\Sigma}}^{12}_{11}(\widehat{\bm{\Sigma}}^{2}_{11})^{-1/2}$,
columns of $\widehat{\mathbf{V}}^{1}$ are the corresponding top-k left
singular vectors and columns of $\widehat{\mathbf{V}}^{2}$ are the
corresponding right singular vectors. The initial estimate ${\mathbf{A}^{m}}$
is obtained by taking the pseudo-inverse of $\widehat{\mathbf{L}}^{m}$,
denoted as

$\displaystyle{\mathbf{A}^{m}}^{(0)}=\bm{\Gamma}^{-1/2}\widehat{\mathbf{V}}^{m\top}(\widehat{\bm{\Sigma}}^{m}_{11})^{-1/2}.$
(4.3)

After obtaining ${\mathbf{A}^{m}}^{(0)}$, we compute the initial estimate for
$\mathbf{B}:=\\{\mathbf{B}_{i}\\}_{i=1}^{p}$. We solve the following
constrained optimization problem:

$\min_{\begin{subarray}{c}\mathbf{B}_{i}\in\mathcal{K}_{B}(s,\alpha)\end{subarray}}h_{i}(\mathbf{B}_{i})=\min_{\begin{subarray}{c}\mathbf{B}_{i}\in\mathcal{K}_{B}(s,\alpha)\end{subarray}}\sum_{m=1}^{M}\frac{1}{2MN}\bigg{\|}{\mathbf{A}^{m}}^{(0)}{\mathbf{Y}_{i}^{m}}-\mathbf{B}_{i}\\{\mathbf{I}_{p-1}\otimes{({\mathbf{A}^{m}}^{(0)})}\\}{\mathbf{Y}_{\backslash\\{i\\}}^{m}}\bigg{\|}_{F}^{2}$

using the projected gradient descent described in Algorithm 2.

Algorithm 2 Initialization of $\mathbf{B}$

Input:
$\\{\mathbf{Y}_{m}\\}_{m=1,\ldots,M}$,$\\{\mathbf{A}^{m(0)}\\}_{m=1,\ldots,M}$,
$\\{\mathbf{B}_{0,i}={\bf 0}_{i\in V}\\}$, $s$, $\alpha$, $\eta_{B_{0}}$;

for $i\in V$, $t=1,\ldots,L-1$ do

$\mathbf{B}_{0,i}^{(t+.5)}\leftarrow{\mathcal{T}}_{s}(\mathbf{B}_{0,i}^{(t)}-\eta_{B_{0}}\nabla_{\mathbf{B}_{0,i}}h)$

for $j\in{V\backslash\\{i\\}}$ do

$\mathbf{B}_{0,ij}^{(t+1)}\leftarrow\mathcal{H}_{\alpha}(\mathbf{B}_{0,ij}^{(t+.5)})$

$\mathbf{B}_{i}^{(0)}\leftarrow\mathbf{B}_{0,i}^{(L)}$

###  4.2 Selection of Basis Parameters

We select $k_{m}$ based on the projection score to the basis, where we employ
the elbow method to decide $k_{m}$. Specifically, we compute the mean-squared
error of the projected signals with the original signals and pick the elbow
point. Then, we compute the canonical correlation of two views and select $k$
based on the canonical correlation score using the elbow method. Although
Theorem 5.4 requires knowledge of $\varepsilon_{0}$, in practice we find
setting $\varepsilon_{0}$ to be a small constant $10^{-3}$ suffices to give
stable and reproducible results in various simulations and real data. The
parameters $s$, $\alpha$, $\tau_{1}$, $\tau_{2}$ are selected based on the
$5$-fold cross-validation with the BIC score function used in the functional
regression (zhao2021high):

$\displaystyle\operatorname{BIC}(s,\alpha,\tau_{1},\tau_{2})=\sum_{i=1}^{p}\left\\{\sum_{m=1}^{M}N\log\det\left(\frac{1}{2N}\mathbf{G}_{i}^{m}\mathbf{G}_{i}^{m\top}\right)+|\widehat{\mathscr{N}}_{i}(s,\alpha,\tau_{1},\tau_{2})|\log
N\right\\},$

where
$\mathbf{G}_{i}^{m}={\mathbf{A}^{m}}{\mathbf{Y}_{i}^{m}}-\sum_{j\in{V\backslash\\{i\\}}}\mathbf{B}_{ij}{\mathbf{A}^{m}}{\mathbf{Y}_{j}^{m}}$
for $i\in V$ and $m=1,\ldots,M$.

##  5 Theory

We show that Algorithm 1–2 can recover the underlying latent graphs with high
probability under mild conditions. The quality of the graph estimate depends
on how well we can estimate $\mathbf{A}^{\star}$ and $\mathbf{B}^{\star}$. We
start by showing the convergence guarantee for both parameters. The
convergence analysis provides the rate of the algorithmic convergence, along
with the statistical error at the stationary points.

In Section 5.1, we assess the quality of the initial estimates. The theory
suggests the consistency of the estimator with a small truncation error. The
convergence guarantee of Algorithm 1 is provided in Section 5.2 along with the
quantification of the statistical error at the convergence points. Section 5.3
presents the guarantee of latent graph recovery.

###  5.1 Analysis of the Initialization

We quantify the initial error of ${\mathbf{A}^{m}}^{(0)}$ to
$\mathbf{A}^{m\star}$ and $\mathbf{B}_{i}^{(0)}$ to $\mathbf{B}_{i}^{\star}$
in terms of the sample size $N$. When the truncation error is negligible, the
result suggests that ${\mathbf{A}^{m}}^{(0)}$ converges in Frobenius norm at a
rate $O(\sqrt{k(k_{1}+k_{2})/N})$. Furthermore, we show that
$\mathbf{B}_{i}^{(0)}$ converges in Frobenius norm at a rate
$O(\sqrt{k(k_{1}+k_{2})/N}+\sqrt{s^{\star}k^{2}/N})$, where the first term is
induced by the estimation error of ${\mathbf{A}^{m}}^{(0)}$ and the second
term is is the statistical error of $\mathbf{B}_{i}^{(0)}$.

We make the following assumption on the matrix
$\mathbf{R}^{12}=(\bm{\Sigma}^{1}_{1,1})^{-1/2}\bm{\Sigma}^{12}_{1,1}(\bm{\Sigma}^{2}_{1,1})^{-1/2}$.

######  Assumption 5.

Suppose that $\min\\{k_{1},k_{2}\\}\geq k$ and the top-k singular values of
$\mathbf{R}^{12}$ satisfy: $\gamma_{1}>\gamma_{2}>\ldots>\gamma_{k}>0$.

By definition, the singular values of $\mathbf{R}^{12}$ are the canonical
correlations of two views. Assumption 5 is for identifiability purpose. It
implies that the canonical vectors are unique, up to sign changes. Let the
sets of left and right singular vectors of $\widehat{\mathbf{R}}^{12}$ and
$\mathbf{R}^{12}$ be
$\\{(\widehat{\mathbf{v}}_{i}^{1},\widehat{\mathbf{v}}_{i}^{2})\\}_{i=1,\ldots,k}$
and
$\\{({\mathbf{v}}_{i}^{1\star},{\mathbf{v}}_{i}^{2\star})\\}_{i=1,\ldots,k}$,
respectively. Hence we define $\mathbf{Q}$ as a diagonal matrix whose $i$-th
diagonal entry $q_{i}\in\\{-1,1\\}$ satisfies the condition that
$\langle{q_{i}\widehat{\mathbf{v}}_{i}^{1}},{\mathbf{v}_{i}^{1\star}}\rangle\geq
0$.

######  Assumption 6.

The covariance
$\bm{\Sigma}^{m}=\frac{1}{N}\mathbb{E}[\mathbf{Y}^{m}\mathbf{Y}^{m\top}]\in\mathbb{R}^{pk_{m}\times
pk_{m}}$ satisfies

$0<\nu_{x}\leq\sigma_{\min}(\bm{\Sigma}^{m})<\sigma_{\max}(\bm{\Sigma}^{m})\leq\rho_{x}<\infty,\quad
m=1,\ldots,M.$

We now establish the theory of the distance of ${\mathbf{A}^{m}}^{(0)}$ and
$\mathbf{A}^{m\star}$ under sign matrix $\mathbf{Q}$.

######  Theorem 5.1.

Let
$\kappa_{m}=\sigma_{\max}(\bm{\Sigma}^{m})/\sigma_{\min}(\bm{\Sigma}^{m})$.
Suppose that $N=O(\max_{m=1,2}\kappa_{m}^{2}k_{m})$ and Assumption 5–6 hold.
Let $\mathscr{L}^{m,r}=\sum_{(\ell,\ell^{\prime})\in
I}\langle{\mathscr{L}^{m}\phi_{\ell^{\prime}}},{\phi_{\ell}^{m}}\rangle_{\mathbb{H}_{m}}\phi_{\ell}^{m}\otimes\phi_{\ell^{\prime}}$,
with
$I=\\{(\ell,\ell^{\prime})\in\mathbb{N}\times\mathbb{N}:(\ell,\ell^{\prime})\neq(a,b),a=1,\ldots,k_{m},b=1,\ldots,k\\}$,
denote the truncation term. Then

$\|{\mathbf{A}^{m}}^{(0)}-\mathbf{Q}\mathbf{A}^{m\star}\|_{F}\leq
C_{\gamma_{k},\nu,\rho}\left(1+\max_{j\neq
i}\frac{1}{|\gamma_{j}-\gamma_{i}|}\right)\sqrt{k\frac{k_{1}+k_{2}}{N}}+\sqrt{8k}|\\!|\\!|\mathscr{A}^{m}|\\!|\\!|^{2}|\\!|\\!|\mathscr{L}^{m,r}|\\!|\\!|,$
(5.1)

with probability at least $1-5\exp(-\sum_{m=1}^{2}k_{m})$, where
$C_{\gamma_{k},\nu,\rho}>0$ is a constant depending on $\gamma_{k}$,
$\nu_{x}$, and $\rho_{x}$ defined in Appendix C.

######  Remark 1.

The second term in (5.1) is the truncation error, which becomes small as $k$
and $k_{m}$ grow. To obtain a non-trivial upper bound, we must have
$|\\!|\\!|{\mathscr{A}^{m}}|\\!|\\!|<\infty$. Therefore, the smallest non-zero
singular value of ${\mathscr{L}^{m}}$ should be bounded away from zero,
suggesting that assuming ${\mathscr{L}^{m}}$ to be a finite-rank operator
(hsing2015theoretical, Theorem 4.2.3) is a necessary condition.  If the
truncation error scales the same as the statistical error, the first term of
(5.1), the theorem tells us that ${\mathbf{A}^{m}}^{(0)}$ is reasonably close
to $\mathbf{A}^{m\star}$, up to a sign change, at the rate of
$O(\sqrt{k(k_{1}+k_{2})/N})$. We can achieve this by choosing the appropriate
$k_{m}$ because $|\\!|\\!|\mathscr{L}^{m,r}|\\!|\\!|$ decreases with
increasing $k_{m}$. Therefore, there exists a $k_{m}$ such that
$|\\!|\\!|\mathscr{L}^{m,r}|\\!|\\!|=O(\sqrt{(k_{1}+k_{2})/N})$. Additionally,
the first term of (5.1) is inversely proportional to the difference between
two consecutive canonical correlations. Hence, the canonical correlations must
be distinct, as stated in Assumption 5, and a gap between two consecutive
canonical correlations is required to obtain a non-trivial upper bound.

Note that multiplications with a signed diagonal matrix do not break
structural assumptions; that is,
$\mathbf{Q}\mathbf{A}^{m\star}\in\mathcal{K}_{m}(\tau_{1},\tau_{2})$ and
$(\mathbf{Q}\mathbf{B}_{i1}^{\star}\mathbf{Q}\cdots\mathbf{Q}\mathbf{B}_{ip}^{\star}\mathbf{Q}){\top}\in\mathcal{K}_{B}(s,\alpha)$.
Furthermore, we have
$f(\mathbf{A},\mathbf{B})=f(\\{\mathbf{Q}{\mathbf{A}^{m}}\\},\\{\mathbf{Q}\mathbf{B}_{i}(\mathbf{I}_{p-1}\otimes\mathbf{Q}^{\top})\\})$
and
$h_{i}(\mathbf{B}_{i})=h_{i}(\mathbf{Q}\mathbf{B}_{i}(\mathbf{I}_{p-1}\otimes\mathbf{Q}^{\top}))$
for any $\mathbf{Q}$. Therefore, it would be cumbersome to write out
$\mathbf{Q}$ explicitly for the rest of the analysis. We abuse the notation
here by writing $\mathbf{Q}\mathbf{A}^{m\star}$ as $\mathbf{A}^{m\star}$ and
$\mathbf{Q}\mathbf{B}_{i}^{\star}\mathbf{Q}^{\top}$ as
$\mathbf{B}_{i}^{\star}$. Define
$\widetilde{\mathbf{B}}_{i}^{\star}=(-\mathbf{B}_{i1}^{\star}\cdots-\mathbf{B}_{ip}^{\star})\in\mathbb{R}^{k\times
pk}$, where $\mathbf{B}_{ii}^{\star}=\mathbf{I}_{k}$.

######  Assumption 7.

For every $i\in V$, $\widetilde{\mathbf{B}}_{i}^{\star}$ satisfies

$0<\nu_{b}^{1/2}\leq\sigma_{\min}(\widetilde{\mathbf{B}}_{i}^{\star})<\sigma_{\max}(\widetilde{\mathbf{B}}_{i}^{\star})\leq\rho_{b}^{1/2}<\infty.$

For every $m=1,\ldots,M$, $\mathbf{A}^{m\star}$ satisfies

$0<\nu_{a}^{1/2}\leq\sigma_{\min}(\mathbf{A}^{m\star})<\sigma_{\max}(\mathbf{A}^{m\star})\leq\rho_{a}^{1/2}<\infty.$

Assumptions 6–7 guarantee that $f(\mathbf{A},\mathbf{B}^{\star})$ is strongly
convex and smooth with respect to ${\mathbf{A}^{m}}$, $m=1,\ldots,M$, and
$f(\mathbf{A}^{\star},\mathbf{B})$ is strongly convex and smooth with respect
to $\mathbf{B}_{i}$, $i\in V$. Furthermore, combined with the result of
Theorem 5.1, it follows that $h_{i}(\mathbf{B})$ is strongly convex and smooth
with respect to $\mathbf{B}_{i}$ for $i\in V$ with high probability.

Let $s=\vartheta_{1}s^{\star}/2$, $\alpha=\vartheta_{2}\alpha^{\star}/2$ for
some constants $\vartheta_{1}$ and $\vartheta_{2}$. Our theory requires more
stringent conditions: $\vartheta_{1}$,$\vartheta_{2}$ are some constants
greater than $2$, which implies that
$\mathbf{B}_{i}^{\star}\in\mathcal{K}_{B}(s^{\star},\alpha^{\star})\subset\mathcal{K}_{B}(s,\alpha)$
for $i\in V$. We state the convergence rate of Algorithm 2.

######  Lemma 5.2.

Suppose that Assumptions 3–4 and 6–7 hold. Let
$m^{\prime}=\mathop{\mathrm{argmin}}\sigma_{\min}^{2}({\mathbf{A}^{m}}^{(0)})$,
$\eta_{B_{0}}=\\{\max_{m=1,\ldots,M}\sigma_{\max}^{2}(\mathbf{A}^{m(0)})\|\widehat{\bm{\Sigma}}^{m}\|_{2}\\}^{-1}$,
$N=O(\max_{m=1,\ldots,M}\kappa_{m}^{2}(k_{m}s^{\star}+\log M+\log p))$ and
$\pi_{0}=C_{\vartheta}\left\\{1-{\sigma_{\min}^{2}(\mathbf{A}^{m^{\prime}(0)})\eta_{B_{0}}\nu_{x}}/2\right\\}$,
where $C_{\vartheta}$, $C_{\gamma}$ and $C_{\delta,i}$ are constants whose
explicit values are given in Appendix C. Then, for each $1\leq i\leq p$, after
$L$ iterations of Algorithm 2 we have

$\displaystyle\|\mathbf{B}_{i}^{(L)}-\mathbf{B}_{i}^{\star}\|_{F}\leq\pi_{0}^{L}\|\mathbf{B}_{i}^{\star}\|_{F}+\frac{C_{\gamma}C_{\vartheta}}{1-\pi_{0}}\frac{\eta_{B_{0}}}{M}\sum_{m=1}^{M}\|{\mathbf{A}^{m}}^{(0)}-\mathbf{A}^{m\star}\|_{F}+\frac{C_{\delta,i}C_{\vartheta}\eta_{B_{0}}}{1-\pi_{0}}\sqrt{\frac{\alpha^{\star}s^{\star}k^{2}}{N}},$

with probability at least $1-3\max_{m=1,\ldots,M}\exp(-k_{m}s^{\star})$.

######  Remark 2.

The quality of the initial guess of $\mathbf{B}_{i}^{\star}$ depends on
$M^{-1}\sum_{m=1}^{M}\|{\mathbf{A}^{m}}^{(0)}-\mathbf{A}^{m\star}\|_{F}$.
While $s^{\star}k^{2}$ in the second term seems large at first glance, it is
worth pointing out that the maximum degree of a node in functional graphical
model scales with $s^{\star}k^{2}\ll k^{2}p$.

###  5.2 Convergence Analysis of Algorithm 1

We show the minimum number of iterations required for Algorithm 1 to produce a
useful result. Additionally, Algorithm 1 converges in a linear rate under mild
assumptions.

The following assumption constrains the space of the local region where
Algorithm 1 takes place. The condition is further fulfilled by Theorem 5.1 and
Lemma 5.2.

######  Assumption 8.

There exist constants $C_{A}$ and $C_{B}$ such that the initial guesses
satisfy $\|\mathbf{A}^{m(0)}-\mathbf{A}^{m\star}\|_{2}\leq
C_{A}\|\mathbf{A}^{m\star}\|_{2}$,
$\|\mathbf{B}_{i}^{(0)}-\mathbf{B}_{i}^{\star}\|\leq
C_{B}\|\mathbf{B}_{i}^{\star}\|$ for any unitarily invariant norm and
$\|\mathbf{B}_{i}^{(0)\top}-\mathbf{B}_{i}^{\star\top}\|_{1}\leq
C_{B}\|\mathbf{B}_{i}^{\star\top}\|_{1}$ for $m=1,\ldots,M$ and $i\in V$.

Let $C_{\alpha},C_{\beta}$, $C_{\eta_{A}},C_{\eta_{B}}$ be constants depending
on the data, as detailed in Appendix C. We introduce the statistical error
$\Xi:=\sqrt{N^{-1}(\Xi_{1}+\Xi_{2})}$ with

$\displaystyle\Xi_{1}:=C_{\alpha}\max_{m=1,\ldots,M}{k(k_{m}+k+\log
M)},\quad\Xi_{2}:=C_{\beta}{s^{\star}(k^{2}+\log p)},$

where $\Xi_{1}$ is induced by the statistical error of estimating
${\mathbf{A}^{m}}$ and $\Xi_{2}$ is by the statistical error of estimating
$\mathbf{B}_{i}$. Under Assumption 8, we define $R_{0}$ as the smallest
positive real number satisfying
$\max_{m=1,\ldots,M}\|\mathbf{A}^{m(0)}-{\mathbf{A}^{m}}\|_{F}^{2}+\max_{i\in
V}\|\mathbf{B}_{i}^{(0)}-\mathbf{B}_{i}^{\star}\|_{F}^{2}\leq R_{0}^{2},$ and

$\displaystyle\delta_{0}$
$\displaystyle:=\exp\left[-\\{(s^{\star}+1)\min_{m=1,\ldots,M}k_{m}\wedge(\alpha
k)^{2}\\}\right],$ $\displaystyle\pi$
$\displaystyle:=\left[4\left\\{1-p(s^{\star}+1)C_{\eta_{A}}\eta_{A}\right\\}\vee
C_{\vartheta}\left(1-MC_{\eta_{B}}\eta_{B}\right)\right].$

We now state the convergence guarantee for Algorithm 1, with proof in Appendix
D.4.

######  Theorem 5.3.

Suppose that Assumptions 3–4,6–8 hold and that, $\Xi^{2}\leq
C_{1}(C_{A}^{2}\rho_{a}+C_{B}^{2}\rho_{b})$, and
$N=O(\max_{m=1,\ldots,M}\kappa_{m}^{2}\\{k_{m}(s^{\star}+1)+\log p+\log
M\\})$. The step sizes $\eta_{A}$ and $\eta_{B}$ are selected to satisfy (i)
$\eta_{A}\leq{C_{2}}\\{p(s^{\star}+1)\\}^{-1}$ and
$\eta_{B}\leq{C_{3}}{M}^{-1}$, (ii) $p\eta_{A}=C_{\vartheta}M\eta_{B}$, and
(iii) $\pi<1$, where $\\{C_{i}\\}_{i=1}^{3}$ are positive constants. After $L$
iterations of Algorithm 1, we have

$\displaystyle\max_{m=1,\ldots,M}\|\mathbf{A}^{m(L)}-{\mathbf{A}^{m}}^{\star}\|_{F}^{2}+\max_{i\in
V}\|\mathbf{B}_{i}^{(L)}-\mathbf{B}_{i}^{\star}\|_{F}^{2}\leq\pi^{L}R_{0}^{2}+\frac{1}{1-\pi}\Xi^{2},$

with probability at least $1-10\delta_{0}$.

######  Remark 3.

From the theorem it follows that Algorithm 1 converges to a local optimum that
is close to the population parameters with distance of up to
$2(1-\pi)^{-1}\Xi^{2}$ if the number of iterations $L$ exceeds
$o((-\log(1-\pi)+2\log\Xi-2\log R_{0})/\log\pi)$.

###  5.3 Graph Recovery

The final step is to combine Theorem 5.1, 5.3, and Lemma 5.2 together and
choose the appropriate parameters and the number of iterations. We have the
following main result.

######  Theorem 5.4.

Suppose that Assumptions 1–7 hold and that $N$, $\eta_{A}$, $\eta_{B}$,
$\eta_{B_{0}}$ satisfy the conditions stated in Theorem 5.3 and Lemma 5.2.
Additionally, $N$ is large enough so that
$\Xi^{2}\leq(1-\pi)2^{-3}\Lambda^{2}(k)$ and there exists a $k_{m}$ such that
$|\\!|\\!|\mathscr{L}^{m,r}|\\!|\\!|=o(1/\sqrt{k_{m}})$. We set
${\mathbf{A}^{m}}^{(0)}$ as (4.3), run Algorithm 2 with the number of
iterations $-C_{1}\log\\{(k^{2}+\log s^{\star})/{N}\\}$ and then run Algorithm
1 with $(-\log(1-\pi)+2\log\Xi-2\log R_{0})/\log\pi$ iterations. Let
$\widehat{\mathscr{N}}_{i}$ be the estimated edges based on the outcomes from
Algorithm 1 with the edge selection threshold $\epsilon_{0}=2^{-1}\Lambda(k)$
defined in (4.1). Then, we have
$\mathbb{P}(\widehat{\mathscr{N}}_{i}=\mathscr{N}_{i}^{\star};1\leq i\leq
p)\geq 1-14\delta_{0}.$

######  Remark 4.

The theorem states that we can recover the edge set $\mathscr{N}_{i}^{\star}$
with probability at least $1-14\delta_{0}$. Note that to establish graph
recovery with a reasonable sample size, the truncated latent signal
$\Lambda(k)$ cannot be too small. As the minimum sample size scales with
$\Lambda^{-2}(k)$, a sufficiently large sample size is required to meet the
condition $\Xi^{2}\leq(1-\pi)2^{-3}\Lambda^{2}(k)$.

##  6 Related Work

We review relevant research on multimodal estimation and functional graphical
models.

The low-rank latent space assumption has been widely explored in multimodal
estimation (zhou2015group; yang2016non). JIVE (lock2013joint) jointly
estimates shared low-rank components and individual low-rank components and
has been successfully applied to biomedical data (o2016r). Our work differs in
the way that we find the low-rank model of the original data generation
process, while previous work finds a low-dimensional representation of the
data, providing a better interpretation of the underlying mechanism.
Multimodal integration framework has been widely used in joint prediction
tasks. li2018integrative constructed linear additive low-rank predictors for
joint multitask regression. li2021integrative proposed a statistical inference
procedure to select significant modalities for integrative linear regression
models. Another methodology is to stack multimodal data in a tensor and
perform low-rank tensor regression (zhou2013tensor). The predictive power of
various models was improved by using rich multimodal data in applications,
including clinical diagnoses (wolfers2015estimating) and biomarker detection
(mimitou2021scalable).

There is increasing interest in functional graphical models.
qiao2019functional established the penalized maximum log-likelihood framework
for Gaussian functional data. qiao2020doubly successively extended it to the
discrete sampling setting and further considered estimation of time-varying
graphs. Zapata2019functional developed a separability condition for the
multivariate covariance operator and applied it to learning functional
graphical models. zhao2021high developed a function-on-function regression
model, as the functional data version of the neighborhood regression method.
While the aforementioned work focuses on the Gaussian distribution, another
active line of work develops functional graphical model estimators under the
non-Gaussian setting. li2018nonparametric and lee2021nonparametric used
additive conditional independence (li2014additive) and proposed nonparametric
estimators. moysidis2021joint studied the joint estimation of multiple
functional Gaussian graphical models by building a hierarchical structure on
inverse covariance. Although this procedure is most closely related to our
work, there are considerable differences in both methodology and theory.
moysidis2021joint extended the approach of guo2011joint to the functional data
setting. They imposed a hierarchical structure on the inverse covariance, and
jointly estimated multiple graphs. Our method, in contrast, takes a generative
perspective and studies a “single” latent graph shared over multiple views.

##  7 Simulations

The simulations focus on comparisons of sparse precision matrices, where we
test four types of sparse graphs. First, we compare our models with three
existing methods by computing the Area under the ROC Curve (AUC). Second, we
compute the convergence distance with respect to the sample size, verifying
the results from Theorem 5.3. Details of data generation processes are
presented in Appendix K.1, where we synthesize four types of different graph
structures and two noise models.

###  7.1 Comparison

Figure 2: The ROC curves of Graph 1–4. The additive noise is generated from
noise model 2. The AUC is discussed in Table 3. The proposed method has
consistent performance across four graphs and $p=\\{50,100,150\\}$.  Figure 3:
Distance v.s. statistical error. The solid lines denote the proposed
algorithm. The dashed lines denote the case where ${\mathbf{A}^{m}}$ is first
initiated with true $\mathbf{A}^{m\star}$ for $m=1,\ldots,M$ and then refined
by Algorithm 1–2 subsequently.

We compare the model with 3 other methods: FGGM (qiao2019functional), PSFGGM
(Zapata2019functional), FPCA (zhao2021high), and JFGGM (moysidis2021joint).
Since those candidates were not originally designed to learn latent graphs, we
apply them to learn the graph of each modality separately and then take the
intersection of the graphs as the surrogate of the latent graph. Let
$\widehat{E}^{1},\widehat{E}^{2}$ be the edge set estimated by one of the
algorithms, then we construct the latent graph as
$\widehat{E}^{z}=\\{(i,j):(i,j)\in\widehat{E}^{1},(i,j)\in\widehat{E}^{2}\\}$.
We fix $k$, $k_{m}$ to be the true values and set $N=100$,
${\mathcal{T}}=[0,1]$ and $p=\\{50,100,150\\}$. Since JFGGM requires computing
an eigen-decomposition at each iterate, which is computationally expensive
when $p$ is large, we omit testing JFGGM under the setting of $p=150$.

To evaluate the performance of competing methods, we vary the sparsity
parameters and plot the Receiver Operating Characteristic (ROC) curves, where
the x-axis denotes the False Positive Rate (FPR) of the estimated edges and
the y-axis denotes the True Positive Rate (TPR) of the estimated edges. Under
the “AND” operation of edge selection and $i\neq j$, the TPR and FPR are
defined respectively as:

$\displaystyle\text{TPR}=\frac{|\\{(i,j):(j,i),(i,j)\in
E^{z}\cap\widehat{E}^{z}\\}|}{|\\{(i,j):(i,j)\in
E^{z}\\}|},\quad\text{FPR}=\frac{|\\{(i,j):(i,j),(j,i)\in\widehat{E}^{z}\backslash
E^{z}\\}|}{|\\{(i,j):(i,j)\not\in E^{z}\\}|}.$

To plot the ROC curve, we vary the sparsity level $s$ for the proposed method
and vary the sparse regularization coefficient $\lambda$ for the competing
methods. Each point in the ROC plot is $(\text{FPR}(s),\text{TPR}(s))$ or
$(\text{FPR}(\lambda),\text{TPR}(\lambda))$ for $1\leq s\leq p$ and $\lambda$
is chosen from a grid of values $[0,\lambda_{\max}]$, where
$\lambda_{\max}=\max_{i\neq
j}N^{-1}\|{\mathbf{Y}_{i}^{m}}({\mathbf{Y}_{j}^{m}})^{\top}\|_{F}$ is
sufficient as discussed in Proposition 1 in zhao2021high. Furthermore, we test
the proposed method with step size $\eta_{A}$, $\eta_{B}$, and $\eta_{B_{0}}$
chosen from $\\{0.1,0.01,0.001,0.0001\\}$ and select the combination of step
size $(\eta_{A},\eta_{B},\eta_{B_{0}})$ that has the largest AUC. We run the
simulation for $10$ runs and take the average of the results. The simulation
result is plotted in Figure 2 and the corresponding AUC is presented in Table
3 in Appendix K.3. The proposed method achieves the best or comparable
performance for Graph 1–4. Appendix K includes more simulation results along
with discussion on the comparing methods.

###  7.2 Distance v.s. sample size

Empirically, we demonstrate that the metric
$\max_{m=1,\ldots,M}\|{\mathbf{A}^{m}}-\mathbf{A}^{m\star}\|_{F}^{2}+\max_{i\in
V}\|\mathbf{B}_{i}-\mathbf{B}_{i}^{\star}\|_{F}^{2}$ at convergence scales
almost linearly with respect to $(k\sum_{m=1}^{M}k_{m}\vee k^{2}\log p)/N$,
which matches the result of Theorem 5.3.222We assume that $s^{\star}$ is a
constant here. We test the graph model 1–3 with the magnitude of the off-
diagonal entries further scaled down by half. For additive noise, we use noise
model 1. For each simulation, we fix $k$, $p$, and $k_{m}$ for $m=1,\ldots,M$
and vary $N$. We take the average of the results for $20$ runs of simulations.
From Figure 3, we see that when $(k\sum_{m=1}^{M}k_{m}\vee k^{2}\log p)/N$ is
smaller than $0.45$,
$\sum_{m=1}^{M}\|{\mathbf{A}^{m}}-\mathbf{A}^{m\star}\|_{F}^{2}+\sum_{i=1}^{p}\|\mathbf{B}_{i}-\mathbf{B}_{i}^{\star}\|_{F}^{2}$
and
$\max_{m=1,\ldots,M}\|{\mathbf{A}^{m}}-\mathbf{A}^{m\star}\|_{F}^{2}+\max_{i\in
V}\|\mathbf{B}_{i}-\mathbf{B}_{i}^{\star}\|_{F}^{2}$ scale almost linearly
with $(k\sum_{m=1}^{M}k_{m}\vee k^{2}\log p)/N$ for most tasks. Note that the
solid lines converge to the dashed lines as the x-axis approaches zero. This
implies that the CCA initialization (4.3) is a consistent estimator, matching
the result of Theorem 5.1.

##  8 Brain Network Estimation with EEG-fMRI

We apply our proposed model to estimate the connectivity patterns of brain
networks using concurrent measurements of EEG-fMRI
(morillon2010neurophysiological; sadaghiani2010intrinsic). The dataset
includes $23$ test subjects in each of the two sessions. Both EEG and fMRI are
first source-localized (wirsich2020concurrent) and then parcellated according
to the Desikan-Killiany cortical atlas (desikan2006automated), which has $68$
parcels. Each recorded session is $600$ seconds. In the first $300$ seconds,
all subjects are synced to watch a short movie clip; in the following $300$
seconds, all the subjects are resting. It is believed that brain networks
display similar patterns across subjects during movie-watching as compared to
the resting-state (vanderwal2019movies). Hence, we only use the first $300$
seconds for estimation and testing. Our goal is to predict the brain networks
of session 2 using the data from session 1. Specifically, we first learn the
precision matrices using the data from session 1. Then, we apply the results
to predict the brain network of session 2, where the log-likelihood is used
for evaluation. The preprocessing pipeline is presented in Appendix M.1.

Since it is difficult to directly evaluate the prediction error of the latent
graph, we design a surrogate task to demonstrate the effectiveness of the
proposed model. We used the estimated edge set $\widehat{E}^{z}$ as auxiliary
information to learn the graph of EEG and fMRI. After obtaining the estimated
edge set $\widehat{E}_{z}$ through the proposed algorithm, we use it to
estimate the precision matrices of fMRI and EEG. Let $m=1$ denote the index
for the fMRI modality and $m=2$ denote the index for the EEG modality. Given
an edge set $E$, we denote ${\mathcal{S}}^{+}_{E}(k_{m})$ be the set of
$pk_{m}\times pk_{m}$ positive definite matrices with support associated to
$E$:
${\mathcal{S}}^{+}_{E}(k_{m}):=\\{{\bm{\Omega}}\in{\mathcal{S}}^{+}:\omega_{i,j}=0,\;(i,j)\not\in
E\\}$, where $\omega_{i,j}\in\mathbb{R}^{k_{m}\times k_{m}}$ is the $(i,j)$-th
submatrix of ${\bm{\Omega}}$. Let $\lambda\geq 0$. We estimate the graph by
solving the problem:

$\widehat{\bm{\Omega}}^{m}=\mathop{\mathrm{argmin}}_{{\bm{\Omega}}^{m}\in{\mathcal{S}}^{+}_{E}(k_{m})}\mathop{\mathrm{tr}}(\widehat{\bm{\Sigma}}^{m}{\bm{\Omega}}^{m})-\log\det({\bm{\Omega}}^{m})+\lambda\sum_{i\neq
j,(i,j)\in E}\|\omega_{i,j}^{m}\|_{1,1}.$ (8.1)

When $\lambda=0$, this is a standard procedure to estimate the graph after
obtaining an edge set through neighborhood regression (ma2016joint). When
$\lambda>0$, this means that the submatrix $\omega_{i,j}$ is sparse, following
a similarly structured assumption of the second condition of
$\mathcal{K}_{B}(s,\alpha)$. We then use the skggm package
(laska_narayan_2017_830033) to optimize (8.1). The parameter $\lambda$ is
selected using 5-fold cross-validation with the BIC metric:
$-2\widehat{L}_{N}+s\log N$, where $\widehat{L}_{N}$ is the log-likelihood and
$s$ is the number of parameters. Let $\widehat{E}_{1}$, $\widehat{E}^{2}$ be
the edge sets obtained by running the regression type algorithm, as detailed
in Appendix M.2, separately for fMRI and EEG data. We validate the
effectiveness of $\widehat{E}^{z}$ by testing three types of edge sets used in
(8.1): individual edge set, latent edge set, and fused edge set. The
individual edge set is $E=\widehat{E}^{1}$ for the fMRI modality or
$E=\widehat{E}^{2}$ for the EEG modality; the latent edge set is
$E=\widehat{E}^{z}$ for both modalities; and the fused edge set is
$E=\widehat{E}^{1}\cup\widehat{E}^{z}$ for the fMRI modality or
$E=\widehat{E}^{2}\cup\widehat{E}^{z}$ for the EEG modality. Then, for the
edge candidate and each data modality, we refit (8.1) to estimate the
precision matrix.

To address the small sample size, we employ sampling with replacement and
repeat the experiments for 5 runs. The result is shown in Table 1, where we
use the log-likelihood as the prediction score. The result indicates that,
with the auxiliary edge information, the fused edge set has maximum in-sample
log-likelihood and out-of-sample log-likelihood for both data modalities. Note
that, with the latent edge set alone, we might neglect the individual graph
structure for each modality, and hence the result is suboptimal.

Data modality | Edge candidate | In-sample Log-likelihood | Out-of-sample Log-likelihood
---|---|---|---
Modality 1 (fMRI) | Individual edge set | $-118411.88(2543.47)$ | $-124217.55(315.57)$
Latent edge set | $-115916.02(3515.06)$ | $-122497.21(377.08)$
| Fused edge set | $-104131.66(1928.12)$ | $-117812.87(248.49)$
Modality 2 (EEG) | Individual edge set | $-22417.26(1189.99)$ | $-34798.19(401.08)$
Latent edge set | $-28523.50(2059.87)$ | $-36820.28(321.84)$
Fused edge set | $-21023.02(701.56)$ | $-32842.76(447.48)$

Table 1: The in-sample and out-of-sample log-likelihood with $3$ different
edge set candidates. The result is the average over $5$ runs and the value in
each parenthesis denotes the standard deviation.

##  9 Discussion

We have developed a new procedure for integrating multimodal functional data
to estimate the underlying latent graph, providing a new statistical solution
to answer scientific questions. The functions are assumed to be continuous
whereas we might only have access to the discrete samples in practice. We
included a discussion to address this issue in Appendix L.  We studied the
setting where data across modalities share the same set of nodes but have
heterogeneous temporal characteristics. In practice, many applications might
have a different number of nodes across modalities. For example, the spatial
resolution of the fMRI data is much higher than that of the EEG data.
Therefore, fMRI is known to capture richer spatial information. On the
contrary, the temporal resolution of the EEG data is much higher than that of
the fMRI data, offering more expressive temporal information. It is left for
future work to develop new statistical models and estimation procedures that
accommodate more complex spatial and temporal discrepancies.

##  10 Acknowledgement

We thank Thomas Alderson and Sepideh Sadaghiani for their help with data pre-
processing and insightful discussion. We thank Benjamin Morillon, Katia
Lehongre, and Anne-Lise Giraud for their generosity of sharing the original
EEG-fMRI raw data. We thank Jonathan Wirsich for providing source-localized
pre-processed EEG-fMRI data. We thank Percy Zhai for sharing the code of the
implementation of the functional Gaussian graphical model (qiao2019functional)
and functional neighborhood regression (zhao2021high). The research project is
funded by the National Science Foundation (NSF) Graduate Research Fellowships
Program, NSF 2046795, 1909577, 1934986, 2216912, and NIFA award
2020-67021-32799.

## References

##  Appendix A Conditional Covariance Operator and Partial Covariance Operator

This section shows the analysis of Section 2.3. Appendix A.1 introduces the
notation. The main result is presented in in Appendix A.2. We start by
discussing the univariate case and the extension to multivariate case follows
similarly. Appendix A.3 shows the proof of Theorem 2.2; Appendix A.5 shows the
proof of Corollary 2.3; Appendix A.6 presents the proofs of auxiliary lemmas.

###  A.1 Notation

Let $\mathbb{H}$ be a separable Hilbert space and $\mathscr{B}(\cdot)$ be the
Borel set. Define $L^{2}(\Omega,\mathcal{F},\mu)$ as the space of square-
integrable random element endowed with inner product
$\mathbb{E}[\langle{\cdot},{\cdot}\rangle_{\mathbb{H}}]$. Let
$\mathfrak{B}(\mathbb{H},\mathbb{H})$ be the space of bounded linear operator
from $\mathbb{H}$ to $\mathbb{H}$. Since $\mathbb{H}$ is a complete normed
space, it follows that $\mathfrak{B}(\mathbb{H},\mathbb{H})$ is a Banach
space. Given a linear operator $\mathscr{T}$, we say that $\mathscr{T}$ is in
the equivalence class of a zero operator if $|\\!|\\!|\mathscr{T}|\\!|\\!|=0$.

###  A.2 Main Results

In this section, we connect the relation of the partial covariance operator,
defined later in (A.3) with the covariance operator defined in
fukumizu2007kernel. We first start with the simple case that discusses
conditioning on the single random element. Consider $\mathcal{Z}_{1}$,
$\mathcal{Z}_{2}$ be mean-squared integrable random elements in a separable
Hilbert space. Let $\mathscr{K}_{12}$ be the covariance operator for
$\mathcal{Z}_{1}$ $\mathcal{Z}_{2}$ and $\mathscr{K}_{1}$, $\mathscr{K}_{1}$
be the variance operator of $\mathcal{Z}_{1}$ and $\mathcal{Z}_{2}$
respectively. It is known (baker1973joint) that there exists a unique bounded
linear operator $\mathscr{V}_{12}$ with norm
$|\\!|\\!|\mathscr{V}_{12}|\\!|\\!|\leq 1$ such that
$\mathscr{K}_{12}=\mathscr{K}_{1}^{1/2}\mathscr{V}_{12}\mathscr{K}_{2}^{1/2}$,
$\operatorname{Im}(\mathscr{V}_{12})\subset\overline{\operatorname{Im}(\mathscr{K}_{1})}$
and
$\ker(\mathscr{V}_{12})^{\perp}\subset\overline{\operatorname{Im}(\mathscr{K}_{2})}$.
Then the conditional covariance operator is defined as

$\displaystyle\mathscr{V}_{12\mid
3}=\mathscr{V}_{12}-\mathscr{V}_{13}\mathscr{V}_{32},$ (A.1)

which could be applied to measure conditional inference.
fukumizu2004dimensionality shows that $\mathscr{V}_{12\mid 3}=0$ if and only
if $\mathcal{Z}_{1}\perp\\!\\!\\!\perp\mathcal{Z}_{2}\mid\mathcal{Z}_{3}$ when
$\mathbb{H}$ is a RKHS. We will show that this property also holds under the
Gaussian assumption in any separable Hilbert space. Before establishing the
theoretical results, we introduce the definition of the Gaussian random
element. First, we introduce the definition of Gaussian linear space.

######  Definition A.1 (Gaussian linear space, adopted from Definition 1.2 in
janson1997gaussian).

A Gaussian linear space is a real linear space of random variables, defined on
$(\Omega,\mathcal{F},\mu)$, such that each variable in the space is centered
Gaussian.

Furthermore, by Theorem 1.5 in janson1997gaussian, we know that any set of
random variables in linear Gaussian space is jointly normal. Next, we
introduce the definition of Gaussian random element.

######  Definition A.2.

Let $\mathcal{Z}\in L^{2}(\Omega,\mathcal{F},\mu)$ be a centered random
element taking values in $(\mathbb{H},\mathscr{B}(\mathbb{H}))$, we say
$\mathcal{Z}$ is a Gaussian random element if
$\\{\langle{\mathcal{Z}},{f}\rangle_{\mathbb{H}}:f\in\mathbb{H}\\}$ forms a
Gaussian linear space.

Consider three centered Gaussian random elements $\mathcal{Z}_{i}\in
L^{2}(\Omega,\mathcal{F},\mu)$ for $i=1,2,3$ taking values in $\mathbb{H}$.
Our goal is to measure the conditional independence of $\mathcal{Z}_{1}$ and
$\mathcal{Z}_{2}$ given the random element $\mathcal{Z}_{3}$. Let
$L^{2}(\Omega,\sigma(\mathcal{Z}_{3}),\mu)$ be a closed linear subspace of
$L^{2}(\Omega,\mathcal{F},\mu)$, where $\sigma(\mathcal{Z}_{3})$ denote the
smallest $\sigma$-algebra generated by $\mathcal{Z}_{3}$. Let
$\widetilde{\beta}_{i3}\in\mathfrak{B}(\mathbb{H},\mathbb{H})$ be a bounded
linear operator, where $\mathfrak{B}(\mathbb{H},\mathbb{H})$ is a Banach
space. Since $\widetilde{\beta}_{i3}$ is continuous, it is measurable and
hence $\widetilde{\beta}_{i3}\mathcal{Z}_{3}\in
L^{2}(\Omega,\sigma(\mathcal{Z}_{3}),\mu)$. Additionally,
$\mathcal{Z}_{i}-\widetilde{\beta}_{i3}\mathcal{Z}_{3}$ is the residual of
$\mathcal{Z}_{i}$ after regressing on $\mathcal{Z}_{3}$ for $i=1,2$. The
regression procedure is defined as:

$\widetilde{\beta}_{i3}=\mathop{\mathrm{argmin}}_{\begin{subarray}{c}\beta\in\mathfrak{B}(\mathbb{H},\mathbb{H})\end{subarray}}\mathbb{E}\left[\|\mathcal{Z}_{i}-\beta\mathcal{Z}_{3}\|^{2}_{\mathbb{H}}\right],\quad
i=1,2.$ (A.2)

Since $\mathcal{Z}_{i}$ for $i=1,2,3$ are centered Gaussian random elements,
the mean element of $\mathcal{Z}_{i}-\widetilde{\beta}_{i3}\mathcal{Z}_{3}$ is
a zero element of $\mathbb{H}$ for $i=1,2$. The partial cross covariance
operator $\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$
is defined as

$\langle{\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}f},{g}\rangle=\mathbb{E}\left[\langle{\mathcal{Z}_{1}-\widetilde{\beta}_{13}\mathcal{Z}_{3}},{g}\rangle_{\mathbb{H}}\langle{\mathcal{Z}_{2}-\widetilde{\beta}_{23}\mathcal{Z}_{3}},{f}\rangle_{\mathbb{H}}\right],\quad
f,g\in\mathbb{H},$ (A.3)

which are bounded linear operators, as shown in Lemma A.7. We define the
partial cross covariance operator
$\mathscr{K}_{21\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ as the
adjoint operator of
$\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$. We
establish the following properties.

######  Theorem A.3.

Let $\mathcal{Z}_{i}\in L^{2}(\Omega,\mathcal{F},\mu)$ for $i=1,2,3$ be
centered Gaussian random elements taking values in
$(\mathbb{H},\mathscr{B}(\mathbb{H}))$.
$\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ and
$\mathscr{K}_{21\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ are
defined in (A.3). The following properties hold.

  1. 1.

$\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ and
$\mathscr{K}_{21\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ are in
the equivalence class of the zero operator if and only if $\mathscr{V}_{12\mid
3}$ and $\mathscr{V}_{21\mid 3}$ are in the equivalence class of the zero
operator.

  2. 2.

$\mathcal{Z}_{1}$ and $\mathcal{Z}_{2}$ are conditionally independent given
$\mathcal{Z}_{3}$, denoted as
$\mathcal{Z}_{1}\perp\\!\\!\\!\perp\mathcal{Z}_{2}\mid\mathcal{Z}_{3}$, if and
only if $\mathscr{V}_{12\mid 3}$ and $\mathscr{V}_{21\mid 3}$ are in the
equivalence class of zero operator.

To evaluate the conditional independence of two Gaussian random elements in
the Hilbert space, one could check the partial covariance operator, realized
by computing (A.2). Theorem A.3 shows the conditional independence properties
when $\mathcal{Z}_{1}$, $\mathcal{Z}_{2}$ and $\mathcal{Z}_{3}$ are random
elements. Oftentimes, we encounter the situation when $\mathcal{Z}_{3}$ is not
a single element, but a set of elements. For instance, this is the case for
estimating Gaussian graphical models. Let
$\mathcal{Z}=(\mathcal{Z}_{1},\ldots,\mathcal{Z}_{p})^{\top}\in\mathbb{H}^{p}$
be centered Gaussian random vector with $\mathcal{Z}_{i}\in
L^{2}(\Omega,\mathcal{F},\mu)$ for $i=1,\ldots,p$. We could easily generalize
Theorem A.3 to a $p$-dimensional centered random vector by studying the linear
subspace $L^{2}(\Omega,\sigma(\mathcal{Z}_{i}:i=1,\ldots,p),\mu)$. The result
is stated in Theorem 2.2.

###  A.3 Proof of Theorem A.3

It suffices to verify the properties of
$\mathscr{K}_{12,\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$, and
similar properties hold for
$\mathscr{K}_{21,\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$. We
begin with showing the first statement. Notice that since
$\operatorname{Im}(\mathscr{V}_{ij})\subset\overline{\operatorname{Im}(\mathscr{K}_{i})}$
and
$\ker(\mathscr{V}_{ij})^{\perp}\subset\overline{\operatorname{Im}(\mathscr{K}_{j})}$
for $i,j\in\\{1,2,3\\}$, $\mathscr{V}_{12\mid 3}=0$ if and only if
$\mathscr{K}_{1}^{1/2}\mathscr{V}_{12\mid 3}\mathscr{K}_{2}^{1/2}=0$.
Therefore, from the result of Lemma A.4, we could establish the following
equivalence: for all $f,g\in\mathbb{H}$,

$\mathscr{V}_{12\mid 3}(f,g)=0,\;\mathscr{V}_{21\mid
3}(f,g)=0\Longleftrightarrow\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}(f,g)=0,\;\mathscr{K}_{21\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}(f,g)=0,$

which proves the first statement. It remains to show that, under the Gaussian
setting, the following equivalence holds

$\mathcal{Z}_{1}\perp\\!\\!\\!\perp\mathcal{Z}_{2}\mid\mathcal{Z}_{3}\Longleftrightarrow\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3},\mathscr{K}_{21\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}\text{
are in the equivalent class of the zero operator}.$

Note that under the Gaussian assumption, by Lemma A.6, the conditional
estimator is the linear estimator. Then, one could look at residual element
$\mathcal{Z}_{i\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ as the
residual as
$\mathcal{Z}_{i\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}=\mathcal{Z}_{i}-\mathbb{E}[\mathcal{Z}_{i}\mid\mathcal{Z}_{3}]$
for $i=1,2$. In addition, since $\mathcal{Z}_{i}$ for $i=1,2,3$ are zero mean,
the mean elements
$\mathfrak{m}_{i\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ for
$i=1,2$ are zero elements. Consequently, for any $f,g\in\mathbb{H}$, we could
write

$\displaystyle\langle{f},{\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}g}\rangle_{\mathbb{H}}$
$\displaystyle=\mathbb{E}\left[\langle{\mathcal{Z}_{1}-\mathbb{E}[\mathcal{Z}_{1}\mid\mathcal{Z}_{3}]},{f}\rangle_{\mathbb{H}}\langle{\mathcal{Z}_{2}-\mathbb{E}[\mathcal{Z}_{2}\mid\mathcal{Z}_{3}]},{g}\rangle_{\mathbb{H}}\right]$
$\displaystyle=\mathbb{E}_{\mathcal{Z}_{3}}\big{[}\underbrace{\mathbb{E}_{\mathcal{Z}_{1},\mathcal{Z}_{2}\mid\mathcal{Z}_{3}}\left[\langle{\mathcal{Z}_{1}-\mathbb{E}[\mathcal{Z}_{1}\mid\mathcal{Z}_{3}]},{f}\rangle_{\mathbb{H}}\langle{\mathcal{Z}_{2}-\mathbb{E}[\mathcal{Z}_{2}\mid\mathcal{Z}_{3}]},{g}\rangle_{\mathbb{H}}\mid\mathcal{Z}_{3}\right]}_{=:\mathscr{T}(f,g)}\big{]}.$
(A.4)

By definition, if
$\mathcal{Z}_{1}\perp\\!\\!\\!\perp\mathcal{Z}_{2}\mid\mathcal{Z}_{3}$, for
any $f,g\in\mathbb{H}$, we have $\mathscr{T}(f,g)=0$, which implies that
$\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}=0$.
Conversely, under the Gaussian setting, the conditional covariance does not
depend on the value of $\mathcal{Z}_{3}$. Therefore, we can ignore the outer
expectation of (A.4) and see that
$\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}=0$ implies
$\mathcal{Z}_{1}\perp\\!\\!\\!\perp\mathcal{Z}_{2}\mid\mathcal{Z}_{3}$.

###  A.4 Proof of Theorem 2.2

The proof step follows similarly to the proof step of the second statement of
Theorem A.3. Noting that, we could generalize the argument of Lemma A.5 to a
set of random element $\\{\mathcal{Z}_{1},\ldots,\mathcal{Z}_{d}\\}$ and show
that the space

$\left\\{\sum_{j^{\prime}\in\\{1,\ldots,p\\}\backslash\\{i,j\\}}\beta_{j^{\prime}}\mathcal{Z}_{j^{\prime}}:\beta_{j^{\prime}}\in\mathfrak{B}(\mathbb{H},\mathbb{H}),j^{\prime}\in\\{1,\ldots,p\\}\backslash\\{i,j\\}\right\\}\\\
\subseteq
L^{2}(\Omega,\sigma(\mathcal{Z}_{j^{\prime}};j^{\prime}\in\\{1,\ldots,p\\}\backslash\\{i,j\\}),\mu),$

is linear and closed, following a similar proof step of Lemma A.5. By the
orthogonality principle of Hilbert space (see alma99954915694905899 for an
exposition), we have

$\left\\{\mathcal{Z}_{i}-\sum_{j^{\prime}\in\\{1,\ldots,p\\}\backslash\\{i,j\\}}\widetilde{\beta}_{ij^{\prime}}\mathcal{Z}_{j^{\prime}}\right\\}\perp\mathcal{Z}_{j^{\prime}},\quad
j\in\\{1,\ldots,p\\}\backslash\\{i,j\\}.$

Then, under the jointly Gaussian assumption, we could write

$\displaystyle\mathcal{Z}_{i}-\sum_{j^{\prime}\in\\{1,\ldots,p\\}\backslash\\{i,j\\}}\widetilde{\beta}_{ij^{\prime}}\mathcal{Z}_{j^{\prime}}=\mathcal{Z}_{i}-\mathbb{E}[\mathcal{Z}_{i}\mid\mathcal{Z}_{-(i,j)}];$
$\displaystyle\mathcal{Z}_{j}-\sum_{j^{\prime}\in\\{1,\ldots,p\\}\backslash\\{i,j\\}}\widetilde{\beta}_{jj^{\prime}}\mathcal{Z}_{j^{\prime}}=\mathcal{Z}_{i}-\mathbb{E}[\mathcal{Z}_{j}\mid\mathcal{Z}_{-(i,j)}],$

following similar argument of Lemma A.6. Then, the remaining step is the same
as the second part of the proof of Theorem A.3, where it is left to verify the
condition: for any $f,g\in\mathbb{H}$

$\mathbb{E}\left[\langle{\mathcal{Z}_{i}-\mathbb{E}[\mathcal{Z}_{i}\mid\mathcal{Z}_{-(i,j)}]},{f}\rangle_{\mathbb{H}}\langle{\mathcal{Z}_{i}-\mathbb{E}[\mathcal{Z}_{j}\mid\mathcal{Z}_{-(i,j)}]},{g}\rangle_{\mathbb{H}}\mid
Z_{j^{\prime}},j^{\prime}\in\\{1,\ldots,p\\}\backslash\\{i,j\\}\right]=0,$

similar as (A.4). Therefore, following Definition 2.1, we can show that
$\mathscr{K}_{ij{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}}}$ is in
the equivalence class of a zero operator if and only if
$\mathcal{Z}_{i}\perp\\!\\!\\!\perp\mathcal{Z}_{j}\mid\mathcal{Z}_{-(i,j)}$
for $(i,j)\not\in E$.

###  A.5 Proof of Corollary 2.3

The first statement is easy to verify as we can apply the orthogonality
principle of the Hilbert space and obtain

$\mathcal{Z}_{i}=\sum_{j^{\prime}\in{V\backslash\\{i\\}}}\beta_{ij^{\prime}}^{\star}\mathcal{Z}_{j^{\prime}}+\mathscr{W}_{i},\quad\mathcal{Z}_{j}=\sum_{j^{\prime}\in\\{1,\ldots,p\\}\backslash\\{j\\}}\beta_{jj^{\prime}}^{\star}\mathcal{Z}_{j^{\prime}}+\mathscr{W}_{j},$

where $\mathscr{W}_{i}$ is independent of $\mathcal{Z}_{j^{\prime}}$ for
$j^{\prime}\in{V\backslash\\{i\\}}$ and $\mathscr{W}_{j}$ is independent of
$\mathcal{Z}_{j^{\prime}}$ for
$j^{\prime}\in\\{1,\ldots,p\\}\backslash\\{j\\}$. Furthermore,
$\mathscr{W}_{i}$ and $\mathscr{W}_{j}$ are centered and independent of each
other.

We then prove the second statement. Here we prove an equivalent statement:
$\beta_{ij}^{\star}$ and $\beta_{ji}^{\star}$ are in the equivalent class of
zero operator if and only if
$\mathcal{Z}_{i}\perp\\!\\!\\!\perp\mathcal{Z}_{j}\mid\mathcal{Z}_{-(i,j)}$.
We first show the forward direction: if $\beta_{ij}^{\star}$ and
$\beta_{ji}^{\star}$ are in the equivalent class of zero operator, then
$\mathcal{Z}_{i}\perp\\!\\!\\!\perp\mathcal{Z}_{j}\mid\mathcal{Z}_{-(i,j)}$.
By Theorem 2.2, we only need to show that
$\mathscr{K}_{ij\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}}$ is in the
equivalence class of zero operator. For any $f,g\in\mathbb{H}$, we can write

$\displaystyle\langle{\mathscr{K}_{ij{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}}}f},{g}\rangle_{\mathbb{H}}$
$\displaystyle=\mathbb{E}\left[\langle{\beta_{ij}^{\star}\mathcal{Z}_{j}+\mathscr{W}_{i}},{f}\rangle_{\mathbb{H}}\langle{\beta_{ji}^{\star}\mathcal{Z}_{i}+\mathscr{W}_{j}},{g}\rangle_{\mathbb{H}}\right]$
$\displaystyle=\mathbb{E}\left[\langle{\mathscr{W}_{i}},{f}\rangle_{\mathbb{H}}\langle{\mathscr{W}_{j}},{g}\rangle_{\mathbb{H}}\right]=0.$

To show the other direction, it suffices to show that
$\mathscr{K}_{ij\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}}$ is in the
equivalence class of a zero operator implies that $\beta_{ij}^{\star}$ in the
equivalence class of a zero operator, which is equivalent to show that

$\langle{\mathbb{E}[\mathcal{Z}_{i}\mid\mathcal{Z}_{-(i,j)}]},{f}\rangle_{\mathbb{H}}=\langle{\mathbb{E}[\mathcal{Z}_{i}\mid\mathcal{Z}_{j\in{V\backslash\\{i\\}}}]},{f}\rangle_{\mathbb{H}}.$

First, we denote the residuals of $\mathcal{Z}_{i}$ and $\mathcal{Z}_{j}$ on
$\mathcal{Z}_{j^{\prime}}$ for
$j^{\prime}\in\\{1,\ldots,p\\}\backslash\\{i,j\\}$ as $r_{i}$ and $r_{j}$,
respectively. Then for any $f\in\mathbb{H}$, we have

$\displaystyle\langle{\mathbb{E}[\mathcal{Z}_{i}\mid\mathcal{Z}_{j\in{V\backslash\\{i\\}}}]},{f}\rangle_{\mathbb{H}}$
$\displaystyle=\langle{\mathbb{E}[\mathcal{Z}_{i}\mid\mathbb{E}[Z_{j}\mid\mathcal{Z}_{-(i,j)}]+r_{j},\mathcal{Z}_{-(i,j)}]},{f}\rangle_{\mathbb{H}}$
$\displaystyle=\langle{\mathbb{E}[\mathcal{Z}_{i}\mid
r_{j},\mathcal{Z}_{-(i,j)}]},{f}\rangle_{\mathbb{H}}$
$\displaystyle=\langle{\mathbb{E}[\mathbb{E}[\mathcal{Z}_{i}\mid\mathcal{Z}_{-(i,j)}]\mid
r_{j},\mathcal{Z}_{-(i,j)}]},{f}\rangle_{\mathbb{H}}+\langle{\mathbb{E}[r_{i}\mid
r_{j},\mathcal{Z}_{-(i,j)}]},{f}\rangle_{\mathbb{H}}.$ Apply the fact that
$\mathbb{E}[\mathcal{Z}_{i}\mid\mathcal{Z}_{-(i,j)}]$ is independent of
$r_{j}$ and $r_{i}$ is independent of $\mathcal{Z}_{-(i,j)}$, the above
display is equivalent as
$\displaystyle=\langle{\mathbb{E}[\mathcal{Z}_{i}\mid\mathcal{Z}_{-(i,j)}]},{f}\rangle_{\mathbb{H}}+\langle{\mathbb{E}[r_{i}\mid
r_{j}]},{f}\rangle_{\mathbb{H}}.$ Since
$\mathscr{K}_{ij\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}}$ is in the
equivalence class of zero operator, we know that the residual $r_{i}$ and
$r_{j}$ are independent. Therefore, the above display is equivalent as
$\displaystyle=\langle{\mathbb{E}[\mathcal{Z}_{i}\mid\mathcal{Z}_{-(i,j)}]},{f}\rangle_{\mathbb{H}}+\langle{\mathbb{E}[r_{i}]},{f}\rangle_{\mathbb{H}}$
$\displaystyle=\langle{\mathbb{E}[\mathcal{Z}_{i}\mid\mathcal{Z}_{-(i,j)}]},{f}\rangle_{\mathbb{H}}.$

Hence, we complete the proof.

###  A.6 Auxiliary Lemmas

######  Lemma A.4.

Under the conditions of Theorem A.3, for any $f,g\in\mathbb{H}$, the
properties hold

$\langle{\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}g},{f}\rangle_{\mathbb{H}}=\langle{\mathscr{K}_{1}^{1/2}\mathscr{V}_{12\mid
3}\mathscr{K}_{2}^{1/2}g},{f}\rangle_{\mathbb{H}},\quad\langle{\mathscr{K}_{21\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}f},{g}\rangle_{\mathbb{H}}=\langle{\mathscr{K}_{2}^{1/2}\mathscr{V}_{21\mid
3}\mathscr{K}_{1}^{1/2}f},{g}\rangle_{\mathbb{H}}.$

###### Proof of Lemma A.4.

Recall the definition of the residual element
$\mathcal{Z}_{i\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ in (A.2)
and based on Lemma A.8, we could write
$\mathcal{Z}_{i\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ as
$\mathcal{Z}_{i}-\beta_{i}\mathcal{Z}_{i}$ for some $\beta_{i}$ that satisfies
$\mathscr{K}_{i3}-\beta_{i}\mathscr{K}_{3}=0$ for $i=1,2$. Since
$\mathcal{Z}_{i}$ for $i=1,2,3$ are zero-mean, it follows that the residual
mean elements $\mathfrak{m}_{i\cdot 3}$ for $i=1,2$ are zero elements. Then,
for any $f,g\in\mathbb{H}$, write

$\displaystyle\langle{f},{\mathscr{K}_{12,\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}g}\rangle_{\mathbb{H}}$
$\displaystyle=\mathbb{E}[\langle{\mathcal{Z}_{1\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}},{f}\rangle_{\mathbb{H}}\langle{\mathcal{Z}_{2\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}},{g}\rangle_{\mathbb{H}}]$
$\displaystyle=\mathbb{E}[\langle{\mathcal{Z}_{1}-\beta_{1}\mathcal{Z}_{3}},{f}\rangle_{\mathbb{H}}\langle{\mathcal{Z}_{2}-\beta_{2}\mathcal{Z}_{3}},{g}\rangle_{\mathbb{H}}]$
$\displaystyle=\langle{f},{(\mathscr{K}_{12}+\beta_{1}\mathscr{K}_{3}\beta_{2}^{*}-\beta_{1}\mathscr{K}_{32}-\mathscr{K}_{13}\beta_{2}^{*})g}\rangle_{\mathbb{H}},$
where $\beta_{2}^{*}$ is the adjoint operator of $\beta_{2}$. Apply the fact
that $\beta_{1}\mathscr{K}_{3}=\mathscr{K}_{13}$. The above display is
equivalent as
$\displaystyle=\langle{f},{(\mathscr{K}_{12}-\beta_{1}\mathscr{K}_{32})g}\rangle_{\mathbb{H}}.$

Next, write $\mathscr{K}_{1}^{1/2}\mathscr{V}_{12\mid
3}\mathscr{K}_{2}^{1/2}=\mathscr{K}_{12}-\mathscr{K}_{1}^{1/2}\mathscr{V}_{13}\mathscr{V}_{32}\mathscr{K}_{2}^{1/2}$.
Then, for any $f,g\in\mathbb{H}$, we have

$\displaystyle\langle{f},{\mathscr{K}_{1}^{1/2}\mathscr{V}_{12\mid
3}\mathscr{K}_{2}^{1/2}g}\rangle_{\mathbb{H}}=\langle{f},{\mathscr{K}_{12}g}\rangle_{\mathbb{H}}-\langle{f},{\mathscr{K}_{1}^{1/2}\mathscr{V}_{13}\mathscr{V}_{32}\mathscr{K}_{2}^{1/2}g}\rangle_{\mathbb{H}}.$

Therefore, to show the first statement, it suffices to show that
$\langle{f},{\mathscr{K}_{1}^{1/2}\mathscr{V}_{13}\mathscr{V}_{32}\mathscr{K}_{2}^{1/2}g}\rangle_{\mathbb{H}}=\langle{f},{\beta_{1}\mathscr{K}_{32}g}\rangle_{\mathbb{H}}$
for any $f,g\in\mathbb{H}$.

Recall that
$\operatorname{Im}(\mathscr{V}_{32})\subset\overline{\operatorname{Im}(\mathscr{K}_{3})}=\overline{\operatorname{Im}(\mathscr{K}_{3}^{1/2})}$
and define the projection operator to
$\overline{\operatorname{Im}(\mathscr{K}_{3}^{1/2})}$ as
$\mathscr{P}_{3}=\mathscr{K}_{3}^{1/2}(\mathscr{K}_{3}^{1/2})^{\dagger}$.
Then, we have $\mathscr{P}_{3}\mathscr{V}_{32}=\mathscr{V}_{32}$. Therefore,
we can write

$\displaystyle\langle{f},{\mathscr{K}_{1}^{1/2}\mathscr{V}_{13}\mathscr{P}_{3}\mathscr{V}_{32}\mathscr{K}_{2}^{1/2}g}\rangle_{\mathbb{H}}$
$\displaystyle=\langle{f},{\mathscr{K}_{1}^{1/2}\mathscr{V}_{13}\mathscr{K}_{3}^{1/2}(\mathscr{K}_{3}^{1/2})^{\dagger}\mathscr{V}_{32}\mathscr{K}_{2}^{1/2}g}\rangle_{\mathbb{H}}$
$\displaystyle=\langle{f},{\mathscr{K}_{13}(\mathscr{K}_{3}^{1/2})^{\dagger}\mathscr{V}_{32}\mathscr{K}_{2}^{1/2}g}\rangle_{\mathbb{H}}.$
Applying the fact that $\mathscr{K}_{13}=\beta_{1}\mathscr{K}_{3}$, the above
display yields
$\displaystyle=\langle{f},{\beta_{1}\mathscr{K}_{3}^{1/2}\mathscr{K}_{3}^{1/2}(\mathscr{K}_{3}^{1/2})^{\dagger}\mathscr{V}_{32}\mathscr{K}_{2}^{1/2}g}\rangle_{\mathbb{H}}$
$\displaystyle=\langle{f},{\beta_{1}\mathscr{K}_{3}^{1/2}\mathscr{P}_{3}\mathscr{V}_{32}\mathscr{K}_{2}^{1/2}g}\rangle_{\mathbb{H}}$
$\displaystyle=\langle{f},{\beta_{1}\mathscr{K}_{3}^{1/2}\mathscr{V}_{32}\mathscr{K}_{2}^{1/2}g}\rangle_{\mathbb{H}}$
$\displaystyle=\langle{f},{\beta_{1}\mathscr{K}_{32}g}\rangle_{\mathbb{H}}.$

Hence, we complete the proof. ∎

######  Lemma A.5.

Let $\mathcal{Z}\in L^{2}(\Omega,\mathcal{F},\mu)$ be a centered random
element taking values in the measurable Hilbert space
$(\mathbb{H},\mathscr{B}(\mathbb{H}))$, where $\mathbb{H}$ is a separable
Hilbert space. Define the space

$V=\\{\beta\mathcal{Z}:\beta\in\mathfrak{B}(\mathbb{H},\mathbb{H})\\},$

where $\mathfrak{B}(\mathbb{H},\mathbb{H})$ is the space of all linear bounded
operator mapping to its own Hilbert space $\mathbb{H}$. Then, $V$ is a closed
linear subspace of $L^{2}(\Omega,\sigma(\mathcal{Z}),\mu)$.

###### Proof of Lemma A.5.

First, we verify that $V\subseteq L^{2}(\Omega,\sigma(\mathcal{Z}),\mu)$.
Since $\beta$ is a bounded linear operator from a complete normed space to
itself, it follows that $\beta$ is a continuous operator and hence measurable.
Therefore, the composite $\beta\mathcal{Z}$ is a measurable operator from
$(\Omega,\sigma(\mathcal{Z}))$ to $(\mathbb{H},\mathscr{B}(\mathbb{H}))$.
Moreover, we have
$\mathbb{E}[\|\beta\mathcal{Z}\|_{\mathbb{H}}^{2}]\leq|\\!|\\!|\beta|\\!|\\!|^{2}\mathbb{E}[\|\mathcal{Z}\|_{\mathbb{H}}^{2}]<\infty$.
Therefore, we have $V\subseteq L^{2}(\Omega,\sigma(\mathcal{Z}),\mu)\subset
L^{2}(\Omega,\mathcal{F},\mu)$.

It is easy to see that $V$ forms a linear subspace and hence what is left is
to verify that $V$ is closed. Let $\\{\beta_{n}\\}_{n\in\mathbb{N}}$ be a
sequence of indexed bounded linear operators
$\beta_{n}\in\mathfrak{B}(\mathbb{H},\mathbb{H})$ for $n\in\mathbb{N}$
converging to $\beta$. Since $\mathfrak{B}(\mathbb{H},\mathbb{H})$ is a Banach
space, it follows that $\beta\in\mathfrak{B}(\mathbb{H},\mathbb{H})$.
Therefore, we have $\\{\beta_{n}\mathcal{Z}\\}_{n\geq 1}$ converging to
$\beta\mathcal{Z}\in V$. Hence $V$ is closed.

∎

######  Lemma A.6.

Let $\mathcal{Z}_{1},\mathcal{Z}_{2}\in L^{2}(\Omega,\mathcal{F},\mu)$ be
centered Gaussian random elements taking values in
$(\mathbb{H},\mathscr{B}(\mathbb{H}))$, where $\mathbb{H}$ is a separable
Hilbert space. Let

$\widehat{\beta}=\mathop{\mathrm{argmin}}_{\beta\in\mathfrak{B}(\mathbb{H},\mathbb{H})}\mathbb{E}\left[\|\mathcal{Z}_{1}-\beta\mathcal{Z}_{2}\|_{\mathbb{H}}^{2}\right].$

Then, we have

$\mathbb{E}[\mathcal{Z}_{1}\mid\mathcal{Z}_{2}]=\widehat{\beta}\mathcal{Z}_{2}.$

This is a classical result when $\mathcal{Z}_{1},\mathcal{Z}_{2}$ are Gaussian
random variables taking values in the real line. The proof follows by the fact
that uncorrelatedness implies independence under the jointly Gaussian setting.
The extension to Gaussian random element in the Hilbert space is
straightforward and yet we provide the statement for clarity.

###### Proof of Lemma A.6.

By Lemma A.8, we can write

$\mathcal{Z}_{1}=e+\widehat{\beta}\mathcal{Z}_{2},$

where we have
$\mathbb{E}[\langle{e},{\widehat{\beta}\mathcal{Z}_{2}}\rangle_{\mathbb{H}}]=0$
by the orthogonality principle of the Hilbert space. Since ${e}$ and
${\widehat{\beta}\mathcal{Z}_{2}}$ are uncorrelated, they are independent
under the joint Gaussian assumption. Therefore, given $\mathcal{Z}_{2}=f$,
$\widehat{\beta}\mathcal{Z}_{2}$ is deterministic and the randomness
$\mathcal{Z}_{2}\mid\mathcal{Z}_{1}$ only comes from $e$. As a result,
$\mathcal{Z}_{1}\mid\mathcal{Z}_{2}\sim\mathcal{N}(\widehat{\beta}\mathcal{Z}_{2},\mathscr{K}_{e})$,
where $\mathscr{K}_{e}$ is the covariance operator of $e$. Note that
$\widehat{\beta}\mathcal{Z}_{2}$ is the mean of
$\mathcal{Z}_{1}\mid\mathcal{Z}_{2}$ and hence we have
$\mathbb{E}[\mathcal{Z}_{1}\mid\mathcal{Z}_{2}]=\widehat{\beta}\mathcal{Z}_{2}$.
∎

######  Lemma A.7.

Let $\mathcal{Z}_{i}\in L^{2}(\Omega,\mathcal{F},\mu)$ for $i=1,2,3$ random
elements taking values in $(\mathbb{H},\mathscr{B}(\mathbb{H}))$. We define
$\mathcal{Z}_{1\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$,
$\mathcal{Z}_{2\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ be the
residuals of $\mathcal{Z}_{1}$ and $\mathcal{Z}_{2}$ after regressing on
$\mathcal{Z}_{3}$, respectively, where recall the definition of regression in
(A.2). The mean elements of
$\mathcal{Z}_{1\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$,
$\mathcal{Z}_{2\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ are
defined as
$\mathfrak{m}_{1\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ and
$\mathfrak{m}_{2\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$
respectively. Define the partial cross covariance operator

$\langle{\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}f},{g}\rangle=\mathbb{E}\left[\langle{(\mathcal{Z}_{1\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}-\mathfrak{m}_{1\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3})},{g}\rangle_{\mathbb{H}}\langle{\mathcal{Z}_{2\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}-\mathfrak{m}_{2\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}},{f}\rangle_{\mathbb{H}}\right],\quad
f,g\in\mathbb{H},$ (A.5)

which is linear and bounded. Let
$\mathscr{K}_{21\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ be the
adjoint operator of
$\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$. Then,
$\mathscr{K}_{21\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ is linear
and bounded.

If $\mathcal{Z}_{i}\in L^{2}(\Omega,\mathcal{F},\mu)$ for $i=1,2,3$ are zero-
mean. Then the mean elements
$\mathfrak{m}_{1\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ and
$\mathfrak{m}_{2\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ are zero
elements. Hence (A.5) is equivalent as (A.3).

###### Proof of Lemma A.7.

For notation simplicity, we define
$\widetilde{\mathcal{Z}}_{i}=\mathcal{Z}_{i\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}-\mathfrak{m}_{i\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$,
the marginal probability measure of $\widetilde{\mathcal{Z}}_{i}$ as
$\widetilde{\mu}_{i}$ for $i=1,2$ and the joint probability measure of
$\widetilde{\mathcal{Z}}_{1}$ and $\widetilde{\mathcal{Z}}_{2}$ as
$\widetilde{\mu}_{12}$. First, for any $f,g\in\mathbb{H}$, we define a
bilinear functional
$\widetilde{G}:\mathbb{H}\otimes\mathbb{H}\rightarrow\mathbb{R}$:

$\widetilde{G}(f,g)=\int_{\Omega\times\Omega}\langle{\widetilde{\mathcal{Z}}_{1}},{f}\rangle_{\mathbb{H}}\langle{\widetilde{\mathcal{Z}}_{2}},{g}\rangle_{\mathbb{H}}d\widetilde{\mu}_{12}.$

We show that $\widetilde{G}(f,g)$ is bounded for any bounded
$f,g\in\mathbb{H}$ by applying Cauchy-Schwarz inequality:

$\displaystyle\left|\widetilde{G}(f,g)\right|^{2}$
$\displaystyle=\left(\int_{\Omega\times\Omega}\langle{\widetilde{\mathcal{Z}}_{1}},{f}\rangle_{\mathbb{H}}\langle{\widetilde{\mathcal{Z}}_{2}},{g}\rangle_{\mathbb{H}}d\widetilde{\mu}_{12}\right)^{2}$
$\displaystyle\leq\left\\{\int_{\Omega}(\langle{\widetilde{\mathcal{Z}}_{1}},{f}\rangle_{\mathbb{H}})^{2}d\widetilde{\mu}_{1}\right\\}\left\\{\int_{\Omega}(\langle{\widetilde{\mathcal{Z}}_{2}},{g}\rangle_{\mathbb{H}})^{2}d\widetilde{\mu}_{2}\right\\}.$
Using the fact that
$\mathcal{Z}_{i}=\widetilde{\mathcal{Z}}_{i}+(\mathcal{Z}_{i}-\widetilde{\mathcal{Z}}_{i})$
and
$E[\langle{\widetilde{\mathcal{Z}}_{i}},{\mathcal{Z}_{i}-\widetilde{\mathcal{Z}}_{i}}\rangle_{\mathbb{H}}]=0$
for $i=1,2$, we could upper bound the above display as
$\displaystyle\leq\left\\{\int_{\Omega}(\langle{{\mathcal{Z}}_{1}},{f}\rangle_{\mathbb{H}})^{2}d{\mu}_{1}\right\\}\left\\{\int_{\Omega}(\langle{{\mathcal{Z}}_{2}},{g}\rangle_{\mathbb{H}})^{2}d{\mu}_{2}\right\\}$
$\displaystyle=\|\mathscr{K}_{1}^{1/2}f\|_{\mathbb{H}}^{2}\|\mathscr{K}_{2}^{1/2}g\|_{\mathbb{H}}^{2}<\infty.$
(A.6)

Therefore, for any fixed $g\in\mathbb{H}$, by Riesz’ representation theorem,
there exists a unique element $q_{g}$ in $\mathbb{H}$ such that
$G(f,g)=\langle{f},{q_{g}}\rangle$. Then, we define
$\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}:\mathbb{H}\rightarrow\mathbb{H}$
by
$\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}g=q_{g}$.
It is clear that
$\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ is linear
and bounded by (A.6):

$\displaystyle\|\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}g\|_{\mathbb{H}}^{2}=\sup_{f\in\mathbb{H}}\frac{|\widetilde{G}(f,g)|^{2}}{\|f\|_{\mathbb{H}}^{2}}\leq\sup_{f\in\mathbb{H}}\frac{\|\mathscr{K}_{1}^{1/2}f\|_{\mathbb{H}}^{2}}{\|f\|_{\mathbb{H}}^{2}}\|\mathscr{K}_{2}^{1/2}g\|_{\mathbb{H}}^{2}\leq\|\mathscr{K}_{1}\|\|\mathscr{K}_{2}\|\|g\|_{\mathbb{H}}^{2},$

which implies that
$\|\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}\|\leq\|\mathscr{K}_{1}\|\|\mathscr{K}_{2}\|<\infty$.
$\mathscr{K}_{12\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}3}$ is
defined everywhere in $\mathbb{H}$. ∎

######  Lemma A.8.

Under the same condition of Lemma A.7, define covariance operator as
$\mathscr{K}_{12}$. Define

$\widehat{\beta}=\mathop{\mathrm{argmin}}_{\begin{subarray}{c}\beta\in\mathfrak{B}(\mathbb{H},\mathbb{H}).\end{subarray}}\mathbb{E}\left[\|\mathcal{Z}_{1}-\beta\mathcal{Z}_{2}\|_{\mathbb{H}}^{2}\right].$

Then, $\widehat{\beta}$ satisfies
$\mathscr{K}_{12}=\widehat{\beta}\mathscr{K}_{2}$.

since
$\operatorname{Im}(\mathscr{K}_{21})\subseteq\operatorname{Im}(\mathscr{K}_{2})$,
it follows that
$\widehat{\beta}=(\mathscr{K}_{2}^{\dagger}\mathscr{K}_{21})^{*}$.

###### Proof of Lemma A.8.

Define the space
$V=\\{\beta\mathcal{Z}_{2}:\beta\in\mathfrak{B}(\mathbb{H},\mathbb{H})\\}$. By
Lemma A.5, we see that $V$ is a closed subspace of
$L^{2}(\Omega,\mathcal{F},\mu)$.

Then, by the orthogonality principle of the Hilbert space (see
alma99954915694905899 for an exposition), we have
$\mathbb{E}[(\mathcal{Z}_{1}-\widehat{\beta}\mathcal{Z}_{2})\otimes\mathcal{Z}_{2}]=0$.
This implies that for any $f,g\in\mathbb{H}$, we have

$\displaystyle\mathbb{E}[\langle{\mathcal{Z}_{1}-\widehat{\beta}\mathcal{Z}_{2}},{f}\rangle_{\mathbb{H}}\langle{\mathcal{Z}_{2}},{g}\rangle_{\mathbb{H}}]$
$\displaystyle=\langle{(\mathscr{K}_{12}-\widehat{\beta}\mathscr{K}_{2})g},{f}\rangle_{\mathbb{H}}=0.$

Therefore, we can conclude that
$\mathscr{K}_{12}-\widehat{\beta}\mathscr{K}_{2}=0$. ∎

##  Appendix B Truncation in Finite Dimensional Space

In this section, we first show the proof of Proposition 3.2, which transforms
the optimization problem (3.3) in the Hilbert space into an optimization
problem of sequences of real numbers. Then, in Appendix B.2, we derive and
propose a finite dimensional optimization problem that computing resources can
realize.

###  B.1 Proof of Proposition 3.2

Recall that for $m=1,\ldots,M$ and $i\in V$ we have

${\mathscr{A}^{m}}=\sum_{\ell,\ell^{\prime\prime}\in\mathbb{N}}a_{\ell\ell^{\prime\prime}}^{m}\phi^{m}_{\ell}\otimes\phi_{\ell^{\prime\prime}}^{m},\quad\beta_{ij}=\sum_{\ell,\ell^{\prime}\in\mathbb{N}}b_{ij,\ell\ell^{\prime}}\phi_{\ell}\otimes\phi_{\ell^{\prime}},\quad{\mathcal{Y}_{i}^{m}}=\sum_{\ell^{\prime\prime}\in\mathbb{N}}y_{i,\ell^{\prime\prime}}^{m}\phi_{\ell^{\prime\prime}}^{m}.$

Therefore, we can write

$\displaystyle\mathbb{E}\left[\bigg{\|}{\mathscr{A}^{m}}{\mathcal{Y}_{i}^{m}}-\sum_{j\in{V\backslash\\{i\\}}}\beta_{ij}{\mathscr{A}^{m}}\mathcal{Y}_{j}^{m}\bigg{\|}_{\mathbb{H}}^{2}\right]$
$\displaystyle\quad=\mathbb{E}\Bigg{[}\Bigg{\|}\left(\sum_{\ell,\ell^{\prime\prime}\in\mathbb{N}}a_{\ell\ell^{\prime\prime}}^{m}\phi^{m}_{\ell}\otimes\phi_{\ell^{\prime\prime}}^{m}\right)\left(\sum_{\ell^{\prime\prime}\in\mathbb{N}}y_{i,\ell^{\prime\prime}}^{m}\phi_{\ell^{\prime\prime}}^{m}\right)$
$\displaystyle\quad\quad\quad-\sum_{j\in{V\backslash\\{i\\}}}\left(\sum_{\ell,\ell^{\prime}\in\mathbb{N}}b_{ij,\ell\ell^{\prime}}\phi_{\ell}\otimes\phi_{\ell^{\prime}}\right)\left(\sum_{\ell^{\prime},\ell^{\prime\prime}\in\mathbb{N}}a_{\ell^{\prime}\ell^{\prime\prime}}^{m}\phi^{m}_{\ell^{\prime}}\otimes\phi_{\ell^{\prime\prime}}^{m}\right)\left(\sum_{\ell^{\prime\prime}\in\mathbb{N}}y_{j,\ell^{\prime\prime}}^{m}\phi_{\ell^{\prime\prime}}^{m}\right)\Bigg{\|}_{\mathbb{H}}^{2}\Bigg{]}$
$\displaystyle\quad=\mathbb{E}\Bigg{[}\Bigg{\|}\left(\sum_{\ell,\ell^{\prime\prime}\in\mathbb{N}}a_{\ell\ell^{\prime\prime}}^{m}y_{i,\ell^{\prime\prime}}^{m}\phi_{\ell}\right)-\sum_{j\in{V\backslash\\{i\\}}}\left(\sum_{\ell,\ell^{\prime},\ell^{\prime\prime}\in\mathbb{N}}b_{ij,\ell\ell^{\prime}}a_{\ell^{\prime}\ell^{\prime\prime}}^{m}y_{j,\ell^{\prime\prime}}^{m}\phi_{\ell}\right)\Bigg{\|}_{\mathbb{H}}^{2}\Bigg{]}$
$\displaystyle\quad=\mathbb{E}\left[\sum_{\ell\in\mathbb{N}}\left(\sum_{\ell^{\prime\prime}\in\mathbb{N}}a_{\ell\ell^{\prime\prime}}^{m}y_{i,\ell^{\prime\prime}}^{m}-\sum_{j\in{V\backslash\\{i\\}}}\sum_{\ell^{\prime},\ell^{\prime\prime}\in\mathbb{N}}b_{ij,\ell\ell^{\prime}}a_{\ell^{\prime}\ell^{\prime\prime}}^{m}y_{j,\ell^{\prime\prime}}^{m}\right)^{2}\right].$

Hence we complete the proof.

###  B.2 Derivation of Equation (3.6)

Similar to the latent space, we define the $k_{m}$-dimensional approximation
of $\chi_{i}^{m}$ as

$\chi_{i}^{k_{m}}:=\sum_{\ell=1}^{k_{m}}x_{i,\ell}^{m}\phi_{\ell}^{m},\quad
i\in V,\;m=1,\ldots,M,$

where $x_{i,\ell}^{m}=\langle{\chi_{i}^{m}},{\phi_{\ell}^{m}}\rangle$ for
$\ell=1,\ldots,k_{m}$. Then, we can represent the truncated
$\mathbf{z}_{i}=(z_{i,1},\ldots,z_{i,k})^{\top}\in\mathbb{R}^{k}$ from the
infinite dimensional vector $(z_{i,1},z_{i,2},\ldots)$ in the following form

$\displaystyle\mathbf{z}_{i}$ $\displaystyle=\begin{pmatrix}z_{i,1}\\\
\vdots\\\
z_{i,k}\end{pmatrix}=\sum_{j\in[p]\backslash\\{i\\}}\begin{pmatrix}b_{ij,11}^{\star}&\ldots&b_{ij,1k}^{\star}\\\
&\ddots&\\\
b_{ij,k1}^{\star}&\ldots&b_{ij,kk}^{\star}\end{pmatrix}\begin{pmatrix}z_{j,1}\\\
\vdots\\\
z_{j,k}\end{pmatrix}+\mathbf{r}_{i}+\mathbf{w}_{i}=\sum_{j\in[p]\backslash\\{i\\}}\mathbf{B}_{ij}^{\star}\mathbf{z}_{j}+\mathbf{r}_{i}+\mathbf{w}_{i},$
(B.1)

where
$r_{i,\ell}=\sum_{j\in{V\backslash\\{i\\}}}\sum_{\ell^{\prime}=k+1}^{\infty}\langle{b_{ij,\ell\ell^{\prime}}^{\star}\phi_{\ell^{\prime}}},{\mathcal{Z}_{j}}\rangle$
, $w_{i,\ell}=\langle{\mathcal{W}_{i}},{\phi_{\ell}}\rangle$ for
$\ell=1,\ldots,k$ and
$\mathbf{r}_{i}=(r_{i,1},\ldots,r_{i,k})^{\top}\in\mathbb{R}^{k}$ denotes the
residual vector. Let
$\mathbb{E}(\mathbf{r}_{i}\mathbf{r}_{i}^{\top})=\bm{\Sigma}_{i}^{\mathbf{r}}$
and
$\mathbb{E}(\mathbf{w}_{i}\mathbf{w}_{i}^{\top})=\bm{\Sigma}_{i}^{\mathbf{w}}$
for $i\in V$.

Following (3.1), we can also express (B.1) as

$\displaystyle\mathcal{Z}^{k}_{i}=\sum_{\ell=1}^{k}\langle{{\mathscr{A}^{m}}(\chi_{i}^{m})},{\phi_{\ell}}\rangle\phi_{\ell}$
$\displaystyle=\sum_{\ell^{\prime}=1}^{\infty}\sum_{\ell=1}^{k}\langle{a_{\ell\ell^{\prime}}^{m\star}\phi_{\ell^{\prime}}^{m}},{\chi_{i}^{m}}\rangle\phi_{\ell}$
$\displaystyle=\sum_{\ell=1}^{k}\left(\sum_{\ell^{\prime}=1}^{k_{m}}\langle{a_{\ell\ell^{\prime}}^{m\star}\phi_{\ell^{\prime}}^{m}},{\chi_{i}^{k_{m}}}\rangle\phi_{\ell}+\sum_{\ell^{\prime}=k_{m}+1}^{\infty}\langle{a_{\ell\ell^{\prime}}^{m\star}\phi_{\ell^{\prime}}^{m}},{\chi_{i}^{m}}\rangle\phi_{\ell}\right).$
(B.2)

To simplify the analysis, we assume that the second term is $0$. Then,
Equation (B.2) has the following equivalent expression in matrix form

$\mathbf{z}_{i}=\begin{pmatrix}z_{i,1}\\\ \vdots\\\
z_{i,k}\end{pmatrix}=\begin{pmatrix}a_{11}^{m\star}&\ldots&a_{1k_{m}}^{m\star}\\\
&\ddots&\\\
a_{k_{m}1}^{m\star}&\ldots&a_{kk_{m}}^{m\star}\end{pmatrix}\begin{pmatrix}x_{i,1}^{m}\\\
\vdots\\\ x_{i,k_{m}}^{m}\end{pmatrix}+\begin{pmatrix}u_{i,1}^{m}\\\ \vdots\\\
u_{i,k}^{m}\end{pmatrix}=\mathbf{A}^{m\star}\mathbf{x}_{i}^{m}+\mathbf{u}_{i}^{m},$
(B.3)

where
$u_{i,\ell}^{m}=\sum_{\ell^{\prime}=k+1}^{\infty}\langle{a_{\ell\ell^{\prime}}^{m\star}\phi_{\ell^{\prime}}^{m}},{\chi_{i}^{m}}\rangle$
for $\ell=1,\ldots,k$ and
$\mathbf{u}_{i}^{m}=(u_{i,1}^{m},\ldots,u_{i,k}^{m})^{\top}\in\mathbb{R}^{k}$
is the truncation error in the observation space. Let
$\mathbf{u}^{m}=(\mathbf{u}_{1}^{m\top},\ldots,\mathbf{u}_{p}^{m\top})^{\top}\in\mathbb{R}^{kp}$
and $\mathbb{E}(\mathbf{u}^{m}\mathbf{u}^{m\top})=\bm{\Sigma}^{m,\mathbf{u}}$.
Combining (B.1) and (B.3) together, we then have the relationship

$\mathbf{A}^{m\star}\mathbf{x}_{i}^{m}+\mathbf{u}_{i}^{m}-\left\\{\sum_{j\in{V\backslash\\{i\\}}}\mathbf{B}_{ij}^{\star}\left(\mathbf{A}^{m\star}\mathbf{x}_{j}^{m}+\mathbf{u}_{j}^{m}\right)+\mathbf{r}_{i}+\mathbf{w}_{i}\right\\}={\bf
0},$ (B.4)

where $\mathbf{w}_{i}$ follows independently from $\mathcal{N}({\bf
0},\bm{\Sigma}^{\mathbf{w}}_{i})$. Let
${\mathbf{q}}_{i}^{m}=(q_{i,1}^{m},\ldots,q_{i,k_{m}}^{m})^{\top}\in\mathbb{R}^{k_{m}}$
where $q_{i,j}^{m}=\langle{\xi_{i}^{m}},{\phi_{j}^{m}}\rangle$ for $j=1,\ldots
k_{m}$, we define the finite truncation of the observed random vector as

$\mathbf{y}_{i}^{m}=\mathbf{x}_{i}^{m}+\mathbf{q}_{i}^{m},\qquad
m=1,\ldots,M,i\in V,$ (B.5)

where $\mathbf{q}_{i}^{m}$ is the noise in the observation space and follows
the distribution $\mathcal{N}({\bf 0},\bm{\Sigma}^{m,\mathbf{q}}_{i})$.
Therefore, under the assumption that $k_{m}$ is large enough such that
$\mathbf{q}_{i}^{m}$ are small in magnitude, we can approximate (B.4) as

$\mathbf{A}^{m\star}\mathbf{y}_{i}^{m}+\mathbf{u}_{i}^{m}-\left\\{\sum_{j\in{V\backslash\\{i\\}}}\mathbf{B}_{ij}^{\star}\left(\mathbf{A}^{m\star}\mathbf{y}_{j}^{m}+\mathbf{u}_{j}^{m}\right)+\mathbf{r}_{i}+\mathbf{w}_{i}\right\\}\approx{\bf
0}.$

Therefore, a reasonable choice of the finite dimensional realization of the
original programming (3.5) is

$\min_{\\{{\mathbf{A}^{m}}\\},\\{\mathbf{B}_{ij}\\}}\sum_{m=1}^{M}\sum_{i=1}^{p}\mathbb{E}\left[\bigg{\|}{\mathbf{A}^{m}}\mathbf{y}_{i}^{m}-\sum_{j\in{V\backslash\\{i\\}}}\mathbf{B}_{ij}{\mathbf{A}^{m}}\mathbf{y}_{j}^{m}\bigg{\|}_{F}^{2}\right].$

##  Appendix C Notations

We define the following quantities

$\displaystyle\Upsilon_{1}=\frac{3\nu_{x}\nu_{b}\rho_{x}\rho_{b}}{16},\quad\Upsilon_{2}=\frac{\nu_{x}\nu_{b}+3\rho_{x}\rho_{b}}{4};$
$\displaystyle\Upsilon_{3}=\|\widehat{\bm{\Sigma}}^{m}\|_{2}(C_{B}+1)\max_{m=1,\ldots,M}\|\mathbf{A}^{m(0)}\|_{2}\max_{i\in
V}(\|\mathbf{B}_{i}^{\star}\|_{F}\vee\|\widetilde{\mathbf{B}}_{i}^{\star}\|_{1});$
$\displaystyle\Upsilon_{4}=2\sqrt{(1+\vartheta_{1})(1+\vartheta_{2})\tau\alpha^{\star}},\quad\Upsilon_{5}=4(1+C_{A})^{2}(2+C_{B})^{2}\rho_{b}\rho_{a}\max_{m=1,\ldots,M}\|\widehat{\bm{\Sigma}}^{m}\|_{2}^{2};$
$\displaystyle\Upsilon_{6}=\frac{\nu_{x}\nu_{a}}{2},\quad\Upsilon_{7}=\rho_{a}^{2}\max_{m=1,\ldots,M}\|\widehat{\bm{\Sigma}}^{m}\|_{2}^{2};$
$\displaystyle\Upsilon_{8}=C_{0}\rho_{x}\rho_{b}\left\\{\rho_{b}\left(\max_{m=1,\ldots,M}\|\bm{\Sigma}^{m,\mathbf{u}}\|_{2}+\rho_{a}\|\bm{\Sigma}^{m,\mathbf{q}}\|\right)+\max_{i\in
V}\left(\|\bm{\Sigma}_{i}^{\mathbf{r}}\|_{2}+\|\bm{\Sigma}_{i}^{\mathbf{w}}\|_{2}\right)\right\\};$
$\displaystyle\Upsilon_{9}=C_{1}C_{\vartheta}(1+\vartheta_{1})\rho_{x}\rho_{a}\left\\{\rho_{b}\left(\max_{m=1,\ldots,M}\|\bm{\Sigma}^{m,\mathbf{u}}\|_{2}+\rho_{a}\|\bm{\Sigma}^{m,\mathbf{q}}\|\right)+\max_{i\in
V}\left(\|\bm{\Sigma}_{i}^{\mathbf{r}}\|_{2}+\|\bm{\Sigma}_{i}^{\mathbf{w}}\|_{2}\right)\right\\};$
$\displaystyle
C_{\vartheta}=\left(1+\sqrt{\frac{2}{\vartheta_{1}-2}}\right)^{2}\left(1+\sqrt{\frac{4}{\vartheta_{2}-2}}\right)^{2},\quad
C_{\alpha}=\frac{\Upsilon_{8}}{\Upsilon_{1}},\quad
C_{\beta}=\frac{\Upsilon_{9}}{\Upsilon_{7}};$ $\displaystyle
C_{\gamma}=\max_{m=1,\ldots,M}\rho_{b}^{1/2}\|\widehat{\bm{\Sigma}}^{m}\|_{2}\left\\{\sigma_{\max}({\mathbf{A}^{m}}^{(0)})+\rho_{a}^{1/2}\right\\}\sqrt{(1+\vartheta_{1})s^{\star}};$
$\displaystyle
C_{\delta,i}=C_{2}\sqrt{(1+\vartheta_{2})(1+\vartheta_{1})}\max_{m=1,\ldots,M}\rho_{x}^{1/2}\rho_{a}^{1/2}\\{\rho_{b}^{1/2}\|\bm{\Sigma}^{m,\mathbf{u}}\|_{2}^{1/2}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+(\rho_{a}\rho_{b})^{1/2}\|\bm{\Sigma}^{m,\mathbf{q}}\|_{2}^{1/2}+\|\bm{\Sigma}_{i}^{\mathbf{r}}\|_{2}^{1/2}+\|\bm{\Sigma}_{i}^{\mathbf{w}}\|_{2}^{1/2}\\};$
$\displaystyle
C_{\gamma_{k},\nu,\rho}=C_{3}\bigg{\\{}\gamma_{k}^{-1/2}(\rho_{x}^{1/2}\nu_{x}^{-2}+\nu_{x}^{-3/2})+\rho_{x}^{2}\nu_{x}^{-5/2}(\gamma_{k}^{-1/2}+\gamma_{k}^{-2})(\rho_{x}^{1/2}\nu_{x}^{-1/2}+1)\bigg{\\}}$
$\displaystyle C_{\eta_{A}}=\frac{3\Upsilon_{1}}{4\Upsilon_{2}},\quad
C_{\eta_{B}}=\frac{\Upsilon_{6}}{4},$

where $C_{0},C_{1},C_{2},C_{3},C_{4},C_{5}$ are absolute constants, which are
independent of data parameters: $M,N,p,s^{\star},\alpha$, upper and lower
bounds of the singular values:
$\nu_{x},\rho_{x},\nu_{b},\rho_{b},\nu_{a},\rho_{a}$, model parameters:
$k,k_{m},C_{A},C_{B}$, and tuning parameters:
$\gamma_{1},\gamma_{2},\tau_{1},\tau_{2}$.

##  Appendix D Analysis of Edge Estimation

The section discusses the analysis of edge estimation, where the results are
presented in Theorem 5.3–5.4. We begin with introducing some notations and
math formulations that will assist the analysis in Section D.1. In Section
D.2, we discuss the shrinkage of the distance metric
$\max_{m=1,\ldots,M}\|{\mathbf{A}^{m}}-\mathbf{A}^{m\star}\|_{F}^{2}$ and
$\max_{i\in V}\|\mathbf{B}_{i}-\mathbf{B}_{i}^{\star}\|_{F}^{2}$ per iterate
of Algorithm 1 along with the statistical error induced by the finite sample
size setting. In section D.3, we combine the results from Section D.2 to state
the contraction of the distance metric
$\max_{m=1,\ldots,M}\|{\mathbf{A}^{m}}-\mathbf{A}^{m\star}\|_{F}^{2}+\max_{i\in
V}\|\mathbf{B}_{i}-\mathbf{B}_{i}^{\star}\|_{F}^{2}$ per iterate of Algorithm
1 and the underlying conditions. With results from Section D.3 in hand, we are
able to apply the telescoping technique to show convergence of Algorithm 1.
The analysis is presented in Section D.4. Section D.5–D.6 present the proof of
lemmas used in Section D.2–D.4.

###  D.1 Notations

For analysis simplicity, we can write the objective function (3.7) as

$\displaystyle f(\mathbf{A},\mathbf{B})$
$\displaystyle=\sum_{i=1}^{p}\sum_{m=1}^{M}\frac{1}{2N}\bigg{\|}{\mathbf{A}^{m}}{\mathbf{Y}_{i}^{m}}-\sum_{j\in{V\backslash\\{i\\}}}\mathbf{B}_{ij}{\mathbf{A}^{m}}{\mathbf{Y}_{j}^{m}}\bigg{\|}_{F}^{2}$
$\displaystyle=\sum_{i=1}^{p}\sum_{m=1}^{M}\frac{1}{2N}\left\|{\mathbf{A}^{m}}{\mathbf{Y}_{i}^{m}}-\mathbf{B}_{i}(\mathbf{I}_{p-1}\otimes{\mathbf{A}^{m}}){\mathbf{Y}_{\backslash\\{i\\}}^{m}}\right\|_{F}^{2}$
$\displaystyle=\sum_{i=1}^{p}\sum_{m=1}^{M}\frac{1}{2N}\|\widetilde{\mathbf{B}}_{i}(\mathbf{I}_{p}\otimes{\mathbf{A}^{m}})\mathbf{Y}^{m}\|_{F}^{2},$
(D.1)

where $\otimes$ denote the Kronecker product. Define

$\widehat{\bm{\Sigma}}^{m}=N^{-1}\mathbf{Y}^{m}(\mathbf{Y}^{m})^{\top},\quad\widehat{\bm{\Sigma}}^{m}_{\cdot\backslash\\{i\\}}=N^{-1}\mathbf{Y}^{m}({\mathbf{Y}_{\backslash\\{i\\}}^{m}})^{\top}.$

Then, for any matrix $\Delta\in\mathbb{R}^{k\times k_{m}}$, we can write

$\langle{\nabla_{{\mathbf{A}^{m}}}f(\mathbf{A},\mathbf{B})},{\Delta}\rangle=\sum_{i=1}^{p}\langle{\widetilde{\mathbf{B}}_{i}^{\top}\widetilde{\mathbf{B}}_{i}(\mathbf{I}_{p}\otimes{\mathbf{A}^{m}})\widehat{\bm{\Sigma}}^{m}},{(\mathbf{I}_{p}\otimes\Delta)}\rangle$
(D.2)

for $m=1,\ldots,M$. Moreover, the gradient of $f(\mathbf{A},\mathbf{B})$ with
respect to ${\mathbf{A}^{m}}$ for $m=1,\ldots,M$ is

$\nabla_{{\mathbf{A}^{m}}}f(\mathbf{A},\mathbf{B})=\sum_{j=1}^{p}(\mathbf{e}_{j}^{\top}\otimes\mathbf{I}_{k})\left\\{\sum_{i=1}^{p}\widetilde{\mathbf{B}}_{i}^{\top}\widetilde{\mathbf{B}}_{i}(\mathbf{I}_{p}\otimes{\mathbf{A}^{m}})\widehat{\bm{\Sigma}}^{m}\right\\}(\mathbf{e}_{j}\otimes\mathbf{I}_{k_{m}}).$
(D.3)

Similarly, for any matrix $\Delta\in\mathbb{R}^{k\times k(p-1)}$, we can write

$\langle{\nabla_{\mathbf{B}_{i}}f(\mathbf{A},\mathbf{B})},{\Delta}\rangle=\sum_{m=1}^{M}\langle{-\widetilde{\mathbf{B}}_{i}(\mathbf{I}_{p}\otimes{\mathbf{A}^{m}})\widehat{\bm{\Sigma}}^{m}_{\cdot\backslash\\{i\\}}(\mathbf{I}_{p-1}\otimes{\mathbf{A}^{m}})^{\top}},{\Delta}\rangle$
(D.4)

for every $i\in V$ and $j\in{V\backslash\\{i\\}}$.

###  D.2 Preliminaries: One-Step Update and statistical quantities

We show Theorem 5.3 in several steps. We begin with showing the contraction of
distance metric per iterate of Algorithm 1. First, given the current estimate
${\mathbf{A}^{m}}$ for $m=1,\ldots,M$, we define an update of
${\mathbf{A}^{m}}$ from Algorithm 1 as

$\displaystyle\mathbf{A}^{m+}=\mathcal{P}_{\tau}({\mathbf{A}^{m}}-\eta_{A}\nabla_{{\mathbf{A}^{m}}}f(\mathbf{A},\mathbf{B})),\quad
m=1,\ldots,M;$

Similarly, for every $i\in V$, we define the update of $\mathbf{B}_{i}$ as
$\mathbf{B}_{i}^{+}$. Let

$\displaystyle{\mathcal{S}}_{i}^{\star}=\text{supp}(\mathbf{B}_{i}^{\star}),\quad{\mathcal{S}}_{i}=\text{supp}(\mathbf{B}_{i}),\quad{\mathcal{S}}^{+}_{i}=\text{supp}(\mathbf{B}_{i}^{+}),\quad\bar{\mathcal{S}}_{i}={\mathcal{S}}^{\star}_{i}\cup{\mathcal{S}}_{i}\cup{\mathcal{S}}_{i}^{+}.$
(D.5)

By definition of $\mathbf{B}_{i}^{\star}$, we have
$\left|{\mathcal{S}}_{i}^{\star}\right|\leq s^{\star}(\alpha^{\star}k)^{2}$
and $\left|{\mathcal{S}}_{i}\right|,\left|{\mathcal{S}}_{i}^{+}\right|\leq
s(\alpha k)^{2}=(\vartheta_{1}s^{\star})(\vartheta_{2}\alpha^{\star}k)^{2}/8$.
Taking the union bound, it follows that
$\left|\bar{\mathcal{S}}_{i}\right|\leq(1+\vartheta_{1})s^{\star}(1+\vartheta_{2})^{2}(\alpha^{\star}k)^{2}$.

Similarly, we define the corresponding sets of node indices as

$\displaystyle\mathscr{N}_{i}^{\star}=\\{j:\|\mathbf{B}_{ij}^{\star}\|_{F}>\epsilon_{0}\\},\quad\mathscr{N}_{i}=\\{j:\|\mathbf{B}_{ij}\|_{F}>\epsilon_{0}\\};\;$
$\displaystyle\mathscr{N}_{i}^{+}=\\{j:\|\mathbf{B}_{i,j}^{+}\|_{F}>\epsilon_{0}\\},\quad\bar{\mathscr{N}}_{i}=\mathscr{N}_{i}^{\star}\cup\mathscr{N}_{i}\cup\mathscr{N}_{i}^{+},$
(D.6)

with $|\bar{\mathscr{N}}_{i}|\leq(1+\vartheta_{1})s^{\star}$ for $i\in V$.

By definition, we have

$\displaystyle\mathbf{B}_{i}^{+}$
$\displaystyle=\mathcal{H}_{\alpha}\circ{\mathcal{T}}_{s}\left(\mathbf{B}_{i}-\eta_{B}\nabla_{\mathbf{B}_{i}}f(\mathbf{A},\mathbf{B})\right)$
$\displaystyle=\mathcal{H}_{\alpha}\circ{\mathcal{T}}_{s}\left(\mathbf{B}_{i}-\eta_{B}[\nabla_{\mathbf{B}_{i}}f(\mathbf{A},\mathbf{B})]_{\bar{\mathcal{S}}_{i}}\right),$
(D.7)

for $i\in V$.

With definitions of ${\mathbf{A}^{m}}^{+}$ for $m=1,\ldots,M$ and
$\mathbf{B}_{i}^{+}$ for $i\in V$, Lemma D.1 states the contraction of the
distance $\|{\mathbf{A}^{m}}^{+}-\mathbf{A}^{m\star}\|_{F}^{2}$ with respect
to $\|{\mathbf{A}^{m}}-\mathbf{A}^{m\star}\|_{F}^{2}$ for $m=1,\ldots,M$ while
Lemma D.2 states the contraction of
$\|\mathbf{B}_{i}^{+}-\mathbf{B}_{i}^{\star}\|_{F}^{2}$ with respect to
$\|\mathbf{B}_{i}-\mathbf{B}_{i}^{\star}\|_{F}^{2}$ for $i\in V$ per iterate
of Algorithm 1. We leave the details of the proof in Section D.5.

######  Lemma D.1 (One-Step Update of ${\mathbf{A}^{m}}$).

Recall that $\Upsilon_{j}$ for $j=1,\ldots,6$ are constants defined in Section
C. Suppose that $k\leq\min_{m=1,\ldots,M}k_{m}$ and Assumption 6–7 hold. There
exists constants $C_{A}$ and $C_{B}$ such that
$\|{\mathbf{A}^{m}}-\mathbf{A}^{m\star}\|_{2}\leq
C_{A}\|\mathbf{A}^{m\star}\|_{2}$,
$\|\mathbf{B}_{i}-\mathbf{B}_{i}^{\star}\|\leq
C_{B}\|\mathbf{B}_{i}^{\star}\|$ for any unitarily invariant norm and
$\|\mathbf{B}_{i}^{\top}-\mathbf{B}_{i}^{\star\top}\|_{1}\leq
C_{B}\|\mathbf{B}_{i}^{\star\top}\|_{1}$ for $m=1,\ldots,M,\;i\in V$. Let
$\kappa_{m}=\sigma_{\max}(\bm{\Sigma}^{m})/\sigma_{\min}(\bm{\Sigma}^{m})$,
$N=O(\kappa_{m}^{2}(\max_{m=1,\ldots,M}k_{m}s^{\star}+\log p+\log M))$ and
$\eta_{A}\leq{2}/\\{3p(s^{\star}+1)\Upsilon_{2}\\}$, then

$\displaystyle\max_{m=1,\ldots,M}\|\mathbf{A}^{m+}-\mathbf{A}^{m\star}\|_{F}^{2}$
$\displaystyle\leq
4\bigg{[}\left\\{1-\eta_{A}\frac{3p(s^{\star}+1)\Upsilon_{1}}{2\Upsilon_{2}}+\eta_{A}\frac{2ps^{\star}\Upsilon_{4}^{2}\Upsilon_{3}^{2}}{\Upsilon_{6}}\right\\}\max_{m=1,\ldots,M}\|{\mathbf{A}^{m}}-\mathbf{A}^{m\star}\|_{F}^{2}$
$\displaystyle\quad+k\eta_{A}\left\\{\frac{4\Upsilon_{2}}{p(s^{\star}+1)\Upsilon_{1}}+3\eta_{A}\right\\}\max_{m=1,\ldots,M}\|\nabla_{{\mathbf{A}^{m}}}f(\mathbf{A}^{\star},\mathbf{B}^{\star})\|_{2}^{2}$
$\displaystyle\quad+p\eta_{A}\left\\{\frac{\Upsilon_{6}}{2}+3p(s^{\star}+1)\Upsilon_{5}\eta_{A}\right\\}\max_{i\in
V}\|\mathbf{B}_{i}-\mathbf{B}_{i}^{\star}\|_{F}^{2}\bigg{]},$ (D.8)

with probability at least
$1-\max_{m=1,\ldots,M}6\exp\\{-k_{m}(s^{\star}+1)\\}$.

######  Lemma D.2 (One-Step Update of $\mathbf{B}_{ij}$).

Recall that $\Upsilon_{j}$ for $j=3,\ldots,7$ are constants defined in Section
C. Suppose that $k\leq\min_{m=1,\ldots,M}k_{m}$ and Assumption 3–4, 6–7 hold.
There exists constants $C_{A}$ and $C_{B}$ such that
$\|{\mathbf{A}^{m}}-\mathbf{A}^{m\star}\|_{2}\leq
C_{A}\|\mathbf{A}^{m\star}\|_{2}$,
$\|\mathbf{B}_{i}-\mathbf{B}_{i}^{\star}\|\leq
C_{B}\|\mathbf{B}_{i}^{\star}\|$ for any unitarily invariant norm and
$\|\mathbf{B}_{i}^{\top}-\mathbf{B}_{i}^{\star\top}\|_{1}\leq
C_{B}\|\mathbf{B}_{i}^{\star\top}\|_{1}$ for $m=1,\ldots,M,\;i\in V$. Let
$\kappa_{m}=\sigma_{\max}(\bm{\Sigma}^{m})/\sigma_{\min}(\bm{\Sigma}^{m})$ and
$N=O(\kappa_{m}^{2}(\max_{m=1,\ldots,M}k_{m}s^{\star}+\log p+\log M))$. Then,
we have

$\displaystyle\max_{i\in
V}\|\mathbf{B}_{i}^{+}-\mathbf{B}_{i}^{\star}\|_{F}^{2}$ $\displaystyle\leq
C_{\vartheta}\bigg{\\{}\left(1-\eta_{B}M\Upsilon_{6}+3\eta_{B}^{2}M^{2}\Upsilon_{7}\right)\max_{i\in
V}\|\mathbf{B}_{i}-\mathbf{B}_{i}^{\star}\|_{F}^{2}$
$\displaystyle\quad+Ms^{\star}\eta_{B}\left(\frac{2\Upsilon_{4}^{2}\Upsilon_{3}^{2}}{\Upsilon_{6}}+3\Upsilon_{5}M\eta_{B}\right)\max_{m=1,\ldots,M}\|{\mathbf{A}^{m}}-\mathbf{A}^{m\star}\|_{F}^{2}$
$\displaystyle\quad+(1+\vartheta_{1})s^{\star}\eta_{B}\left(\frac{2}{M\Upsilon_{6}}+3\eta_{B}\right)\max_{i\in
V}\|[\nabla_{\mathbf{B}_{i}}f(\mathbf{A}^{\star},\mathbf{B}^{\star})]_{\bar{\mathcal{S}}_{i}}\|^{2}_{r(k,\infty)}\bigg{\\}},$
(D.9)

with probability at least $1-\max_{m=1,\ldots,M}\exp(-k_{m}s^{\star})$.

Observe that the right hand sides of (D.8) and (D.9) has the term of
$\|\nabla_{{\mathbf{A}^{m}}}f(\mathbf{A}^{\star},\mathbf{B}^{\star})\|_{2}^{2}$
and
$\|[\nabla_{\mathbf{B}_{i}}f(\mathbf{A}^{\star},\mathbf{B}^{\star})]_{\bar{\mathcal{S}}_{i}}\|^{2}_{r(k,\infty)}$,
respectively. By law of large number, ideally, these quantities would be
negligible as the number of sample size goes to infinity. Specifically,
$(\mathbf{A}^{\star},\mathbf{B}^{\star})$ is a stationary point of
$\mathbb{E}f(\mathbf{A}^{\star},\mathbf{B}^{\star})$. Hence, we expect that
$\|\nabla_{{\mathbf{A}^{m}}}f(\mathbf{A}^{\star},\mathbf{B}^{\star})\|_{2}^{2}$
and
$\|[\nabla_{\mathbf{B}_{i}}f(\mathbf{A}^{\star},\mathbf{B}^{\star})]_{\bar{\mathcal{S}}_{i}}\|^{2}_{r(k,\infty)}$
to be small when we have some reasonable number of sample size. Lemma D.3
states the upper bound of the
$\|\nabla_{{\mathbf{A}^{m}}}f(\mathbf{A}^{\star},\mathbf{B}^{\star})\|_{2}$
for $m=1,\ldots,M$ while Lemma D.4 provides the upper bound of
$\|[\nabla_{\mathbf{B}_{i}}f(\mathbf{A}^{\star},\mathbf{B}^{\star})]_{\bar{\mathcal{S}}_{i}}\|^{2}_{r(k,\infty)}$
for $i\in V$.

######  Lemma D.3 (statistical error bound of
$\nabla_{{\mathbf{A}^{m}}}f(\mathbf{A}^{\star},\mathbf{B}^{\star})$).

Under conditions of Lemma D.1 and assume that
$N=O(C_{\epsilon}^{2}(k+k_{m}+\log M))$. Let $C_{\epsilon}=C_{3}\max_{i\in
V}\rho_{x}^{1/2}\rho_{b}^{1/2}\big{\\{}\rho_{b}^{1/2}\|\bm{\Sigma}^{m,\mathbf{u}}\|_{2}^{1/2}+(\rho_{b}\rho_{a})^{1/2}\|\bm{\Sigma}^{m,\mathbf{q}}\|_{2}^{1/2}+\|\bm{\Sigma}_{i}^{\mathbf{r}}\|_{2}^{1/2}+\|\bm{\Sigma}_{i}^{\mathbf{w}}\|_{2}^{1/2}\big{\\}}$
for some absolute constant $C_{3}>0$.

Then, we have

$\frac{1}{p(s^{\star}+1)}\|\nabla_{{\mathbf{A}^{m}}}f(\mathbf{A}^{\star},\mathbf{B}^{\star})\|_{2}>C_{\epsilon}\sqrt{\frac{k_{m}+k+\log
M}{N}},$

with probability smaller than $M^{-1}\exp\\{-(k_{m}+k)\\}$.

######  Lemma D.4 (statistical error bound of
$\nabla_{\mathbf{B}_{i}}f(\mathbf{A}^{\star},\mathbf{B}^{\star})$).

Let $C_{4},C_{5}>0$ be some absolute constants and