4\sum_{j=1}^{p}\|\mathbf{x}_{j}-\mathbf{y}_{j}\|_{2}^{2}=4\|\mathbf{X}-\mathbf{Y}\|_{F}^{2},$

where the last inequality follows by Lemma H.2. ∎

######  Lemma H.2.

Let
$\widetilde{\mathcal{K}}(\tau_{1},\tau_{2})=\\{\mathbf{x}\in\mathbb{R}^{d}:\tau_{1}\leq\|\mathbf{x}\|_{2}\leq\tau_{2}\\}$
and $\widetilde{\mathcal{P}}_{\tau}$ be the projection operator that projects
$\mathbf{x}$ to $\widetilde{\mathcal{K}}(\tau_{1},\tau_{2})$. Then, for any
$\mathbf{y}\in\widetilde{\mathcal{K}}(\tau_{1},\tau_{2})$, we have

$\|\widetilde{\mathcal{P}}_{\tau}(\mathbf{x})-\mathbf{y}\|_{F}\leq
2\|\mathbf{x}-\mathbf{y}\|_{F}.$

###### Proof of Lemma H.2.

We consider the following three cases:
_Case 1_ : If $\|\mathbf{x}\|_{2}\geq\tau_{2}$, then projection to
$\widetilde{\mathcal{K}}(\tau_{1},\tau_{2})$ is equivalent as the projection
to a ball with radius $\tau_{2}$, which is a convex set. Then, following the
property that projection to the convex set is contraction, we have
$\|\widetilde{\mathcal{P}}_{\tau}(\mathbf{x})-\mathbf{y}\|_{F}\leq\|\mathbf{x}-\mathbf{y}\|_{F}$.
_Case 2_ : If $\tau_{1}\leq\|\mathbf{x}\|_{2}\leq\tau_{2}$, then it is clear
that $\widetilde{\mathcal{P}}_{\tau}(\mathbf{x})=\mathbf{x}$ and hence
$\|\widetilde{\mathcal{P}}_{\tau}(\mathbf{x})-\mathbf{y}\|_{F}=\|\mathbf{x}-\mathbf{y}\|_{F}$.
_Case 3_ : If $\|\mathbf{x}\|_{2}\leq\tau_{1}$, then define
$\widehat{\mathbf{x}}=\widetilde{\mathcal{P}}_{\tau}(\mathbf{x})$ we can apply
triangle inequality and obtain

$\displaystyle\|\widehat{\mathbf{x}}-\mathbf{y}\|_{2}\leq\|\mathbf{x}-\mathbf{y}\|_{2}+\|\widehat{\mathbf{x}}-\mathbf{x}\|_{2}\leq
2\|\mathbf{x}-\mathbf{y}\|_{2},$

where the last inequality follows by the fact that
$\|\widehat{\mathbf{x}}-\mathbf{x}\|_{2}\leq\|\mathbf{x}-\mathbf{y}\|_{2}$.
Therefore, we have
$\|\widetilde{\mathcal{P}}_{\tau}(\mathbf{x})-\mathbf{y}\|_{F}\leq
2\|\mathbf{x}-\mathbf{y}\|_{F}$.

Taking the maximum coefficient of three cases, we complete the proof. ∎

######  Lemma H.3.

[Lemma $B.3$ in (zhang2018unified)] Let
$\mathbf{B}^{\star}\in\mathbb{R}^{k\times k}$ and suppose that there are at
most $\alpha^{\star}$-fraction of non-zero entries per row and per column of
$\mathbf{B}^{\star}$. Then, for any $\mathbf{B}$ and
$1\geq\alpha\geq\alpha^{\star}$ and $\mathcal{H}_{\alpha}(\cdot)$ be the hard
thresholding operator defined in Section 4, we have

$\|\mathcal{H}_{\alpha}(\mathbf{B})-\mathbf{B}^{\star}\|_{F}\leq\left(1+\sqrt{\frac{2\alpha^{\star}}{\alpha-\alpha^{\star}}}\right)\|\mathbf{B}-\mathbf{B}^{\star}\|_{F}.$

######  Lemma H.4 (Group-sparse hard thresholding).

Let $s,s^{\star}$ be some integers such that $s>s^{\star}$,
$\mathbf{B},\mathbf{B}^{\star}\in\mathbb{R}^{k\times pk}$ and
$\mathbf{B}^{\star}$ is $s^{\star}$ group-sparse, i.e.,
$\|\mathbf{B}^{\star}\|_{r(k,0)}\leq s^{\star}$. Recall
${\mathcal{T}}_{s}(\cdot)$ defined in Section 4, we have

$\|{\mathcal{T}}_{s}(\mathbf{B})-\mathbf{B}^{\star}\|_{F}\leq\left(1+\sqrt{\frac{s^{\star}}{s-s^{\star}}}\right)\|\mathbf{B}-\mathbf{B}^{\star}\|_{F}.$

###### Proof.

First, let $\widehat{\mathbf{B}}={\mathcal{T}}_{s}(\mathbf{B})$,
$\text{supp}(\widehat{\mathbf{B}})=\widehat{S}$,
$\text{supp}(\mathbf{B}^{\star})=S^{\star}$ and $I=\widehat{S}\cup S^{\star}$.
We have

$\|\widehat{\mathbf{B}}-\mathbf{B}^{\star}\|_{F}=\|[\widehat{\mathbf{B}}-\mathbf{B}^{\star}]_{I}\|_{F}\leq\|[\widehat{\mathbf{B}}-\mathbf{B}]_{I}\|_{F}+\|[\mathbf{B}-\mathbf{B}^{\star}]_{I}\|_{F}.$

Let $U=S^{\star}/\widehat{S}$ and $V=\widehat{S}/S^{\star}$, then we have

$\|[\widehat{\mathbf{B}}-\mathbf{B}]_{I}\|_{F}^{2}=\|[\mathbf{B}]_{U}\|_{F}^{2}\leq\frac{s^{\star}}{s-s^{\star}}\|[\mathbf{B}]_{V}\|_{F}^{2}\leq\frac{s^{\star}}{s-s^{\star}}\|[\mathbf{B}-\mathbf{B}^{\star}]_{I}\|_{F}^{2}.$

Therefore, we have

$\|\widehat{\mathbf{B}}-\mathbf{B}^{\star}\|_{F}\leq\left(1+\sqrt{\frac{s^{\star}}{s-s^{\star}}}\right)\|[\mathbf{B}-\mathbf{B}^{\star}]_{I}\|_{F}\leq\left(1+\sqrt{\frac{s^{\star}}{s-s^{\star}}}\right)\|\mathbf{B}-\mathbf{B}^{\star}\|_{F}.$

Wwe complete the analysis. ∎

##  Appendix I Auxiliary Lemmas

This section introduces several useful properties that are used in analyses in
previous sections.

######  Lemma I.1.

Assume that
$(\mathbf{x}^{(1)},\mathbf{y}^{(1)}),\ldots,(\mathbf{x}^{(N)},\mathbf{y}^{(N)})$
are $N$ independent realizations of
$(\mathbf{x},\mathbf{y})\in(\mathbb{R}^{p},\mathbb{R}^{d})$, which are
distributed as $\mathbf{x}\sim\mathcal{N}({\bf 0},\bm{\Sigma}_{x})$ and
$\mathbf{y}\sim\mathcal{N}({\bf 0},\bm{\Sigma}_{y})$. Let
$\bar{\mathbf{x}}=N^{-1}\sum_{n=1}^{N}\mathbf{x}^{(n)}$,
$\bar{\mathbf{y}}=N^{-1}\sum_{n=1}^{N}\mathbf{y}^{(n)}$ and
$\widehat{\bm{\Sigma}}_{xy}=N^{-1}\sum_{n=1}^{N}(\mathbf{x}^{(n)}-\bar{\mathbf{x}})(\mathbf{y}^{(n)}-\bar{\mathbf{y}})^{\top}$
be the sample covariance. Given $\delta>0$, there exist an universal constant
$C_{0}>0$ such that

$\mathbb{P}\left\\{\frac{\|\widehat{\bm{\Sigma}}_{xy}-\mathbb{E}\widehat{\bm{\Sigma}}_{xy}\|_{2}}{\|\mathbb{E}\widehat{\bm{\Sigma}}_{xy}\|_{2}}\geq
C_{0}\left(\sqrt{\frac{p+d+\log(1/\delta)}{N}}+\frac{p+d+\log(1/\delta)}{N}\right)\right\\}\leq\delta.$

Suppose that there exist $\sigma_{p\wedge d},\sigma_{1}>0$ such that
$\sigma_{p\wedge
d}\leq\sigma_{min}(\mathbb{E}\bm{\Sigma}_{xy})\leq\sigma_{max}(\mathbb{E}\bm{\Sigma}_{xy})\leq\sigma_{1}$.
If $N=O(\kappa^{2}(p+d))$ with $\kappa=\sigma_{1}/\sigma_{p\wedge d}$, then

$\mathbb{P}\left\\{\frac{\sigma_{p\wedge
d}}{2}\leq\sigma_{\min}(\widehat{\bm{\Sigma}}_{xy})\leq\sigma_{\max}(\widehat{\bm{\Sigma}}_{xy})\leq\frac{3\sigma_{1}}{2}\right\\}\geq
1-\exp(-p-d).$

###### Proof.

The first part of the proof can be found in Theorem 6.1.1 in
tropp2015introduction, so here we only present the second part of the proof.
Note that when $N=O(\kappa^{2}(p+d))$, we have
$\|\widehat{\bm{\Sigma}}_{xy}-\mathbb{E}\widehat{\bm{\Sigma}}_{xy}\|_{2}\leq\frac{\sigma_{p\wedge
d}}{2}$ with probability at least $1-\exp\\{-(p+d)\\}$. By Weyl’s inequality,
we have

$\sigma_{\max}(\widehat{\bm{\Sigma}}_{xy})\leq\sigma_{\max}(\mathbb{E}\widehat{\bm{\Sigma}}_{xy})+\|\widehat{\bm{\Sigma}}_{xy}-\mathbb{E}\widehat{\bm{\Sigma}}_{xy}\|_{2}\leq\frac{3\sigma_{1}}{2},$

and

$\sigma_{\min}(\widehat{\bm{\Sigma}}_{xy})\geq\sigma_{\min}(\mathbb{E}\widehat{\bm{\Sigma}}_{xy})-\|\widehat{\bm{\Sigma}}_{xy}-\mathbb{E}\widehat{\bm{\Sigma}}_{xy}\|_{2}\geq\frac{\sigma_{p\wedge
d}}{2},$

and hence we complete the second part of the proof. ∎

######  Lemma I.2.

Let $\mathbf{Y}\in\mathbb{R}^{pd\times N}$ be a matrix whose columns are
i.i.d. random vectors drawn from a distribution. Define
$\widehat{\bm{\Sigma}}=\frac{1}{N}\mathbf{Y}(\mathbf{Y})^{\top}\in\mathbb{R}^{pd\times
pd}$, $\bm{\Sigma}=\mathbb{E}[\widehat{\bm{\Sigma}}]\in\mathbb{R}^{pd\times
pd}$. Let $\mathscr{M}_{i}\subset\\{1,\ldots,p\\}$, where
$\left|\mathscr{M}_{i}\right|\leq r$. Let
$\widehat{\bm{\Sigma}}_{\mathscr{M}_{i}\mathscr{M}_{i}}\in\mathbb{R}^{\left|\mathscr{M}_{i}\right|d\times\left|\mathscr{M}_{i}\right|d}$
be the sub-matrix of $\widehat{\bm{\Sigma}}$ whose rows and columns are
corresponding to nodes in $\mathscr{M}_{i}$. Define the event

$\displaystyle{\mathcal{E}_{i}}$ $\displaystyle(\mathscr{M}_{i})$
$\displaystyle=\left\\{\frac{\sigma_{\min}(\bm{\Sigma}_{\mathscr{M}_{i}\mathscr{M}_{i}})}{2}\leq\sigma_{\min}(\widehat{\bm{\Sigma}}_{\mathscr{M}_{i}\mathscr{M}_{i}})\leq\sigma_{\max}(\widehat{\bm{\Sigma}}_{\mathscr{M}_{i}\mathscr{M}_{i}})\leq\frac{3\sigma_{\max}(\bm{\Sigma}_{\mathscr{M}_{i}\mathscr{M}_{i}})}{2}\right\\}.$

Let $\kappa=\sigma_{\max}(\bm{\Sigma})/\sigma_{\min}(\bm{\Sigma})$. Under the
condition that $N=O(\kappa^{2}(dr+\log M+\log p))$, the event
${\mathcal{E}}=\cap_{i=1}^{p}{\mathcal{E}_{i}}(\mathscr{M}_{i})$ happens with
probability at least $1-M^{-1}\exp\\{-dr\\}$.

###### Proof.

Then, by Lemma I.1, under the condition that $N=O(\kappa^{2}(dr+\log M+\log
p))$, we have

$\mathbb{P}\left(|\\!|\\!|\widehat{\bm{\Sigma}}_{{\mathscr{M}_{i}}{\mathscr{M}_{i}}}-\bm{\Sigma}_{\mathscr{M}_{i}\mathscr{M}_{i}}|\\!|\\!|_{{2}}>\frac{\sigma_{\min}(\bm{\Sigma}_{\mathscr{M}_{i}\mathscr{M}_{i}})}{2}\right)\leq\frac{1}{pM}\exp(-dr),$

and we will be on the event ${\mathcal{E}_{i}}(\mathscr{M}_{i})$ with
probability at least $1-(pM)^{-1}\exp\left(-dr\right)$. Then, we define the
intersection of all events ${\mathcal{E}_{i}}(\mathscr{M}_{i})$ for $i\in V$
as

$\displaystyle{\mathcal{E}}=\cap_{i=1}^{p}{\mathcal{E}_{i}}(\mathscr{M}_{i}).$

Then, by De Morgan’s laws, event ${\mathcal{E}}$ happens with probability at
least $1-M^{-1}\exp(-dr)$.

∎

######  Lemma I.3 (Theorem $2.1.12$ in nesterov2003introductory).

For a $L$-smooth and $m$-strongly convex function $h$, we have

$\displaystyle\langle\nabla h_{i}(X)-\nabla
h_{i}(Y),X-Y\rangle\geq\frac{mL}{m+L}\|X-Y\|_{F}^{2}+\frac{1}{m+L}\|\nabla
h_{i}(X)-\nabla h_{i}(Y)\|_{F}^{2}.$

######  Lemma I.4 (Theorem $6.5$ in wainwright2019high).

Let $\mathbf{X}\in\mathbb{R}^{p\times d}$ and
$\mathcal{N}_{1}=\\{\mathbf{v}_{1},\ldots,\mathbf{v}_{n}\\}$ be the $1/8$-net
of $\mathbb{S}^{p-1}$ and
$\mathcal{N}_{2}=\\{\mathbf{u}_{1},\ldots,\mathbf{u}_{m}\\}$ be the $1/8$-net
of $\mathbb{S}^{d-1}$. Then,

$\sup_{\mathbf{v}\in\mathbb{S}^{p-1}}\sup_{\mathbf{u}\in\mathbb{S}^{d-1}}\left|\langle{\mathbf{v}},{\mathbf{A}\mathbf{u}}\rangle\right|\leq
2\max_{\mathbf{v}\in\mathcal{N}_{1}}\max_{\mathbf{u}\in\mathcal{N}_{2}}\left|\langle{\mathbf{v}},{\mathbf{A}\mathbf{u}}\rangle\right|.$

######  Lemma I.5.

Let $\frac{1}{p}+\frac{1}{q}=1$ and $\frac{1}{m}+\frac{1}{n}=1$. Then the dual
norm of $\|\cdot\|_{p,m}$ is $\|\cdot\|_{q,n}$.

###### Proof.

Assume that we have two matrices $\mathbf{X},\mathbf{Y}\in\mathbb{R}^{p\times
d}$ and recall the definition of the dual norm

$\|\mathbf{X}\|_{*}=\sup\\{\mathop{\mathrm{tr}}(\mathbf{Y}^{\top}\mathbf{X})|\|\mathbf{Y}\|\leq
1\\}.$

Then, it is easy to see that

$\mathop{\mathrm{tr}}(\mathbf{Y}^{\top}\mathbf{X})\leq\|\mathbf{Y}\|\|\mathbf{X}\|_{*},$

for every pair of $\mathbf{X}$, $\mathbf{Y}$.

Let $\mathbf{x}_{1},\ldots,\mathbf{x}_{d}$ be columns of $\mathbf{X}$ and
$\mathbf{y}_{1},\ldots,\mathbf{y}_{d}$ be columns of $\mathbf{Y}$. Apply
Hölder’s inequality, we have

$\displaystyle\mathop{\mathrm{tr}}(\mathbf{X}^{\top}\mathbf{Y})$
$\displaystyle=\sum_{j=1}^{d}\langle{\mathbf{x}_{j}},{\mathbf{y}_{j}}\rangle$
$\displaystyle\leq\sum_{j=1}^{d}\|\mathbf{x}_{j}\|_{p}\|\mathbf{y}_{j}\|_{q}.$
Apply the Hölder’s inequality again, we can obtain
$\displaystyle\leq\left(\sum_{j=1}^{d}\|\mathbf{x}_{j}\|_{p}^{m}\right)^{\frac{1}{m}}\left(\sum_{j=1}^{d}\|\mathbf{y}_{j}\|_{q}^{n}\right)^{\frac{1}{n}}$
$\displaystyle=\|\mathbf{X}\|_{p,m}\|\mathbf{Y}\|_{q,n}.$

The second step is to show that the upper bound is achievable. Let $y_{kj}$ be
the $kj$-th entry of $\mathbf{Y}$ such that

$y_{kj}=\frac{\mathop{\mathrm{sign}}(x_{kj})|x_{kj}|^{p-1}}{\|\mathbf{x}_{j}\|_{p}^{p-m}(\sum_{j^{\prime}=1}^{d}\|\mathbf{x}_{j^{\prime}}\|_{p}^{m})^{1/n}},\qquad
k=1,\ldots p,j=1,\ldots,d.$

We can easily verify that $\|\mathbf{Y}\|_{q,n}=1$. Note that we have

$\|\mathbf{y}_{j}\|_{q}=\frac{\|\mathbf{x}_{j}\|_{p}^{q/p}}{\|\mathbf{x}_{j}\|_{p}^{p-m}(\sum_{j^{\prime}=1}^{d}\|\mathbf{x}_{j^{\prime}}\|_{p}^{m})^{1/n}}=\frac{\|\mathbf{x}_{j}\|^{m-1}_{p}}{(\sum_{j^{\prime}=1}^{d}\|\mathbf{x}_{j^{\prime}}\|_{p}^{m})^{1/n}},\qquad
j=1,\ldots,d.$

Summing over $d$ columns, we then conclude that

$\|\mathbf{Y}\|_{q,n}^{n}=\sum_{j=1}^{d}\|\mathbf{y}_{j}\|_{q}^{n}=\frac{\sum_{j=1}^{d}\|\mathbf{x}_{j}\|_{p}^{n(m-1)}}{\sum_{j^{\prime}=1}^{d}\|\mathbf{x}_{j^{\prime}}\|_{p}^{m}}=1.$

Moreover, we can verify that

$\displaystyle\mathop{\mathrm{tr}}(\mathbf{X}^{\top}\mathbf{Y})=\sum_{j=1}^{d}\langle{\mathbf{x}_{j}},{\mathbf{y}_{j}}\rangle$
$\displaystyle=\sum_{j=1}^{d}\frac{\|\mathbf{x}_{j}\|_{p}^{p}}{\|\mathbf{x}_{j}\|_{p}^{p-m}(\sum_{j^{\prime}=1}^{d}\|\mathbf{x}_{j^{\prime}}\|_{p}^{m})^{1/n}}$
$\displaystyle=\frac{\sum_{j=1}^{d}\|\mathbf{x}_{j}\|_{p}^{m}}{(\sum_{j^{\prime}=1}^{d}\|\mathbf{x}_{j^{\prime}}\|_{p}^{m})^{1/n}}=\left(\sum_{j=1}^{d}\|\mathbf{x}_{j}\|_{p}^{m}\right)^{1/m}=\|\mathbf{X}\|_{p,m}.$

Here we complete the proof. ∎

######  Lemma I.6.

Let $\frac{1}{p}+\frac{1}{q}=1$ and recall the definition of
$\|\cdot\|_{r(k,p)}$ and $\|\cdot\|_{r(k,q)}$ in Definition 3.3. Then the dual
norm of $\|\cdot\|_{r(k,p)}$ is $\|\cdot\|_{r(k,q)}$.

###### Proof.

The proof is similar to proof of Lemma I.5. Assume that
$\mathbf{X},\mathbf{Y}\in\mathbb{R}^{d\times ks}$ and are represented as
$\mathbf{X}=(\mathbf{X}_{1}\cdots\mathbf{X}_{s})$ and
$\mathbf{Y}=(\mathbf{Y}_{1}\cdots\mathbf{Y}_{s})$ for
$\mathbf{X}_{i}\in\mathbb{R}^{d\times k}$ and
$\mathbf{Y}_{i}\in\mathbb{R}^{d\times k}$ for $i=2,\ldots,s$. Then, apply
Hölder’s inequality, we can write

$\displaystyle\mathop{\mathrm{tr}}(\mathbf{X}^{\top}\mathbf{Y})$
$\displaystyle=\sum_{i=1}^{m}\mathop{\mathrm{tr}}(\mathbf{X}_{i}^{\top}\mathbf{Y}_{i})$
$\displaystyle=\sum_{i=1}^{m}\operatorname{vec}\left(\mathbf{X}_{i}\right)^{\top}\operatorname{vec}\left(\mathbf{Y}_{i}\right)$
$\displaystyle\leq\sum_{i=1}^{m}\|\operatorname{vec}\left(\mathbf{X}_{i}\right)\|_{p}\|\operatorname{vec}\left(\mathbf{Y}_{i}\right)\|_{q}.$
Recall the definition of $\|\cdot\|_{r(k,p)}$ and $\|\cdot\|_{r(k,q)}$ in
Definition 3.3, where we can apply Lemma I.5 with $m=n=2$ and arrive at the
following conclusion
$\displaystyle\leq\|\mathbf{X}\|_{r(k,p)}\|\mathbf{Y}\|_{r(k,q)}.$

∎

######  Lemma I.7.

Let $\mathscr{A},\mathscr{B}\in\mathfrak{B}(\mathbb{H}_{1},\mathbb{H}_{2})$ be
two compact operators, where $\mathbb{H}_{1},\mathbb{H}_{2}$ are separable
Hilbert spaces. Let $\mathscr{P}_{A}$ be the projection operator to
$\overline{\operatorname{Im}(\mathscr{A})}$ and $\mathscr{P}_{B}$ be the
projection operator to $\overline{\operatorname{Im}(\mathscr{B})}$. If
$\operatorname{rank}(\mathscr{A})>\operatorname{rank}(\mathscr{B})$, then

$|\\!|\\!|\mathscr{P}_{B}^{\perp}\mathscr{P}_{A}|\\!|\\!|\geq|\\!|\\!|\mathscr{P}_{B}\mathscr{P}_{A}^{\perp}|\\!|\\!|.$

Lemma I.7 similarly holds for the Hilbert-Schimidt norm. Since we only use the
result for the operator norm, we do not provide the result for the Hilbert-
Schimidt norm.

###### Proof.

We start by decomposing the projection operator as
$\mathscr{P}_{A}=\mathscr{P}_{A_{1}}+\mathscr{P}_{A_{2}}$, where
$\operatorname{rank}(\mathscr{P}_{A_{1}})=\operatorname{rank}(\mathscr{P}_{B})$
and $\mathscr{P}_{A_{2}}\mathscr{P}_{B}=0$. Then

$\displaystyle|\\!|\\!|\mathscr{P}_{B}\mathscr{P}_{A}^{\perp}|\\!|\\!|$
$\displaystyle=|\\!|\\!|\mathscr{P}_{B}(I-\mathscr{P}_{A_{1}}-\mathscr{P}_{A_{2}})|\\!|\\!|=|\\!|\\!|\mathscr{P}_{B}(I-\mathscr{P}_{A_{1}})|\\!|\\!|=|\\!|\\!|\mathscr{P}_{B}\mathscr{P}_{A_{1}}^{\perp}|\\!|\\!|$
Since
$\operatorname{rank}(\mathscr{P}_{A_{1}})=\operatorname{rank}(\mathscr{P}_{B})$,
we have
$|\\!|\\!|\mathscr{P}_{B}\mathscr{P}_{A_{1}}^{\perp}|\\!|\\!|=|\\!|\\!|\mathscr{P}_{B}^{\perp}\mathscr{P}_{A_{1}}|\\!|\\!|$
(davis1970rotation). Then, the above display is equivalent to
$\displaystyle=|\\!|\\!|\mathscr{P}_{B}^{\perp}\mathscr{P}_{A_{1}}|\\!|\\!|\leq|\\!|\\!|\mathscr{P}_{B}^{\perp}\mathscr{P}_{A}|\\!|\\!|,$

which completes the proof. ∎

######  Lemma I.8.

Let $\mathscr{A},\mathscr{B}\in\mathfrak{B}(\mathbb{H}_{1},\mathbb{H}_{2})$ be
two compact operators, where $\mathbb{H}_{1},\mathbb{H}_{2}$ are separable
Hilbert spaces. Let $\mathscr{P}_{A}$ be the projection operator to
$\overline{\operatorname{Im}(\mathscr{A})}$ and $\mathscr{P}_{B}$ be the
projection operator to $\overline{\operatorname{Im}(\mathscr{B})}$. If
$\operatorname{rank}(\mathscr{A})\geq\operatorname{rank}(\mathscr{B})$, then

$|\\!|\\!|\mathscr{A}^{\dagger}-\mathscr{B}^{\dagger}|\\!|\\!|\leq
2\left(|\\!|\\!|\mathscr{A}^{\dagger}|\\!|\\!|^{2}\vee|\\!|\\!|\mathscr{B}^{\dagger}|\\!|\\!|^{2}\right)|\\!|\\!|\mathscr{A}-\mathscr{B}|\\!|\\!|.$

Lemma I.8 similarly holds for the Hilbert-Schmidt norm with constant 2
replaced by $2^{1/2}$. We only prove the result for the operator norm.

###### Proof.

We generalize Theorem $3.3$ in stewart1977perturbation to a Hilbert space. Let
$\mathscr{Q}_{A}$ be the projection operator to $\ker(\mathscr{A})$ and
$\mathscr{P}_{A}$ be the orthogonal projection to
$\overline{\operatorname{Im}(\mathscr{A})}$. Recalling the properties of a
pseudo-inverse, we have $\mathscr{P}_{A}=\mathscr{A}\mathscr{A}^{\dagger}$ and
define $\mathscr{R}_{A}=\mathscr{A}^{\dagger}\mathscr{A}=I-\mathscr{Q}_{A}$.
We define $\mathscr{P}_{B}$ and $\mathscr{R}_{B}$ in the same way. Then, we
can write

$\displaystyle\mathscr{B}^{\dagger}$
$\displaystyle=\mathscr{B}^{\dagger}\mathscr{B}\mathscr{B}^{\dagger}$
$\displaystyle=\mathscr{B}^{\dagger}\mathscr{P}_{B}$
$\displaystyle=\mathscr{B}^{\dagger}\mathscr{P}_{B}(\mathscr{A}\mathscr{A}^{\dagger}+I-\mathscr{A}\mathscr{A}^{\dagger})$
$\displaystyle=\mathscr{B}^{\dagger}\mathscr{P}_{B}(\mathscr{A}\mathscr{A}^{\dagger}\mathscr{A}\mathscr{A}^{\dagger}+I-\mathscr{A}\mathscr{A}^{\dagger})$
$\displaystyle=\mathscr{B}^{\dagger}\mathscr{P}_{B}(\mathscr{A}\mathscr{R}_{A}\mathscr{A}^{\dagger}+\mathscr{P}_{A}^{\perp}).$

Similarly, write

$\displaystyle\mathscr{A}^{\dagger}$
$\displaystyle=\mathscr{A}^{\dagger}\mathscr{A}\mathscr{A}^{\dagger}$
$\displaystyle=\mathscr{R}_{A}\mathscr{A}^{\dagger}$
$\displaystyle=(\mathscr{B}^{\dagger}\mathscr{B}+I-\mathscr{B}^{\dagger}\mathscr{B})\mathscr{R}_{A}\mathscr{A}^{\dagger}$
$\displaystyle=(\mathscr{B}^{\dagger}\mathscr{B}\mathscr{B}^{\dagger}\mathscr{B}+I-\mathscr{B}^{\dagger}\mathscr{B})\mathscr{R}_{A}\mathscr{A}^{\dagger}$
$\displaystyle=(\mathscr{B}^{\dagger}\mathscr{P}_{B}\mathscr{B}+\mathscr{R}_{B}^{\perp})\mathscr{R}_{A}\mathscr{A}^{\dagger}.$

Therefore, we have

$\displaystyle\mathscr{B}^{\dagger}-\mathscr{A}^{\dagger}=-\mathscr{B}^{\dagger}\mathscr{P}_{B}(\mathscr{B}-\mathscr{A})\mathscr{R}_{A}\mathscr{A}^{\dagger}+\mathscr{B}^{\dagger}\mathscr{P}_{B}\mathscr{P}_{A}^{\perp}-\mathscr{R}_{B}^{\perp}\mathscr{R}_{A}\mathscr{A}^{\dagger}=T_{1}+T_{2}+T_{3}.$

Note that $\overline{\operatorname{Im}(\mathscr{B}^{\dagger})}$ is orthogonal
to the null space $\ker(\mathscr{B})$ and we have
$\operatorname{Im}(T_{1}),\operatorname{Im}(T_{2})\in\overline{\operatorname{Im}(\mathscr{B}^{\dagger})}$
and $\operatorname{Im}(T_{3})\in\ker(\mathscr{B})$. We have

$|\\!|\\!|\mathscr{B}^{\dagger}-\mathscr{A}^{\dagger}|\\!|\\!|\leq|\\!|\\!|T_{1}+T_{2}|\\!|\\!|+|\\!|\\!|T_{3}|\\!|\\!|.$

Since
$T_{1}+T_{2}=\mathscr{B}^{\dagger}\left\\{-\mathscr{P}_{B}(\mathscr{B}-\mathscr{A})\mathscr{A}^{\dagger}\mathscr{P}_{A}+\mathscr{P}_{B}\mathscr{P}_{A}^{\perp}\right\\}$,
we have

$\displaystyle|\\!|\\!|T_{1}+T_{2}|\\!|\\!|$
$\displaystyle\leq|\\!|\\!|\mathscr{B}^{\dagger}|\\!|\\!|\left(|\\!|\\!|\mathscr{P}_{B}(\mathscr{B}-\mathscr{A})\mathscr{A}^{\dagger}|\\!|\\!|+|\\!|\\!|\mathscr{P}_{B}\mathscr{P}_{A}^{\perp}|\\!|\\!|\right)$
$\displaystyle\leq|\\!|\\!|\mathscr{B}^{\dagger}|\\!|\\!|\left(|\\!|\\!|\mathscr{P}_{B}(\mathscr{B}-\mathscr{A})\mathscr{A}^{\dagger}|\\!|\\!|+|\\!|\\!|\mathscr{P}_{B}^{\perp}\mathscr{P}_{A}|\\!|\\!|\right)$
$\displaystyle(\text{Lemma~{}\ref{lemma:projection_hs}})$
$\displaystyle=|\\!|\\!|\mathscr{B}^{\dagger}|\\!|\\!|\left\\{|\\!|\\!|\mathscr{P}_{B}(\mathscr{B}-\mathscr{A})\mathscr{A}^{\dagger}|\\!|\\!|+|\\!|\\!|\mathscr{P}_{B}^{\perp}(\mathscr{B}-\mathscr{A})\mathscr{A}^{\dagger}|\\!|\\!|\right\\}$
$\displaystyle(\mathscr{P}_{B}^{\perp}\mathscr{B}=0)$
$\displaystyle\leq|\\!|\\!|\mathscr{B}^{\dagger}|\\!|\\!||\\!|\\!|\mathscr{A}^{\dagger}|\\!|\\!|\left\\{|\\!|\\!|\mathscr{P}_{B}(\mathscr{B}-\mathscr{A})|\\!|\\!|+|\\!|\\!|\mathscr{P}_{B}^{\perp}(\mathscr{B}-\mathscr{A})|\\!|\\!|\right\\}$
$\displaystyle=|\\!|\\!|\mathscr{B}^{\dagger}|\\!|\\!||\\!|\\!|\mathscr{A}^{\dagger}|\\!|\\!||\\!|\\!|\mathscr{B}-\mathscr{A}|\\!|\\!|.$

For $T_{3}$, since
$|\\!|\\!|\mathscr{R}_{B}^{\perp}\mathscr{R}_{A}|\\!|\\!|=|\\!|\\!|(\mathscr{R}_{B}^{\perp}\mathscr{R}_{A})^{*}|\\!|\\!|=|\\!|\\!|\mathscr{R}_{A}^{*}(\mathscr{R}_{B}^{\perp})^{*}|\\!|\\!|=|\\!|\\!|\mathscr{R}_{A}\mathscr{R}_{B}^{\perp}|\\!|\\!|$,
we have

$\displaystyle|\\!|\\!|\mathscr{R}_{B}^{\perp}\mathscr{R}_{A}\mathscr{A}^{\dagger}|\\!|\\!|$
$\displaystyle\leq|\\!|\\!|\mathscr{A}^{\dagger}|\\!|\\!||\\!|\\!|\mathscr{R}_{A}\mathscr{R}_{B}^{\perp}|\\!|\\!|$
$\displaystyle=|\\!|\\!|\mathscr{A}^{\dagger}|\\!|\\!||\\!|\\!|\mathscr{A}^{\dagger}(\mathscr{A}-\mathscr{B})\mathscr{R}_{B}^{\perp}|\\!|\\!|$
$\displaystyle(\mathscr{B}\mathscr{R}_{B}^{\perp}=0)$
$\displaystyle\leq|\\!|\\!|\mathscr{A}^{\dagger}|\\!|\\!|^{2}|\\!|\\!|\mathscr{B}-\mathscr{A}|\\!|\\!|.$

Combining the last two displays, we have

$|\\!|\\!|\mathscr{A}^{\dagger}-\mathscr{B}^{\dagger}|\\!|\\!|\leq
2\left(|\\!|\\!|\mathscr{A}^{\dagger}|\\!|\\!|^{2}\vee|\\!|\\!|\mathscr{B}^{\dagger}|\\!|\\!|^{2}\right)|\\!|\\!|\mathscr{A}-\mathscr{B}|\\!|\\!|,$

which completes the proof. ∎

######  Lemma I.9.

[Adapted from Theorem 5.2.2 in hsing2015theoretical Let
$\mathscr{T},\widetilde{\mathscr{T}}\in\mathfrak{B}(\mathbb{H}_{1},\mathbb{H}_{2})$
be two compact operators associated with singular systems
$\\{\lambda_{j},f_{1j},f_{2j}\\}_{j\in\mathbb{N}}$ and
$\\{\widetilde{\lambda}_{j},\widetilde{f}_{1j},\widetilde{f}_{2j}\\}_{j\in\mathbb{N}}$,
respectively. Without the loss of generality, we assume that
$\\{f_{1j}\\}_{j\in\mathbb{N}}$ and $\\{f_{2j}\\}_{j\in\mathbb{N}}$ provide
CONS for $\mathbb{H}_{1}$ and $\mathbb{H}_{2}$, respectively. Define
$\Delta=\widetilde{\mathscr{T}}-\mathscr{T}$ and $\zeta_{j}=(2^{-1}\min_{k\neq
j}|\lambda_{k}^{2}-\lambda_{j}^{2}|)^{-1}(|\\!|\\!|\mathscr{T}|\\!|\\!|+|\\!|\\!|\Delta|\\!|\\!|)|\\!|\\!|\Delta|\\!|\\!|$.
Assume that $\zeta_{j}<1-\epsilon$ for any number $\epsilon\in(0,1)$. Then,
there exists a constant $C=C_{\epsilon}\in(0,\infty)$ such that

$\bigg{\|}(\widetilde{f}_{1j}-f_{1j})-\sum_{k\neq
j}\frac{\lambda_{k}\langle{f_{2k}},{\Delta
f_{1j}}\rangle_{\mathbb{H}_{2}}+\lambda_{j}\langle{f_{2j}},{\Delta
f_{1k}}\rangle_{\mathbb{H}_{2}}}{\lambda_{k}^{2}-\lambda_{j}^{2}}f_{1k}\bigg{\|}_{\mathbb{H}_{1}}\leq
C{\gamma_{j}},$

with
$\gamma_{j}=\zeta_{j}\\{\zeta_{j}+|\\!|\\!|\Delta|\\!|\\!|/(|\\!|\\!|\Delta|\\!|\\!|+|\\!|\\!|\mathscr{T}|\\!|\\!|)\\}$.

######  Lemma I.10.

Let $\mathbf{A}$ and $\mathbf{B}$ be $d\times d$ real symmetric matrices such
that $0\prec\mathbf{A},\mathbf{B}\preceq\lambda\mathbf{I}$ for some constant
$\lambda>0$, then we have

$\|\mathbf{A}^{-1/2}-\mathbf{B}^{-1/2}\|_{2}\leq\left(3\lambda^{1/2}\|\mathbf{A}^{-1/2}\|_{2}+1\right)\|\mathbf{B}^{-3/2}\|_{2}\|\mathbf{A}-\mathbf{B}\|_{2}.$

###### Proof.

First, we write

$\displaystyle\|\mathbf{A}^{-1/2}-\mathbf{B}^{-1/2}\|_{2}$
$\displaystyle=\|\mathbf{A}^{-1/2}(\mathbf{B}^{3/2}-\mathbf{A}^{3/2})\mathbf{B}^{-3/2}+(\mathbf{A}-\mathbf{B})\mathbf{B}^{-3/2}\|_{2}$
$\displaystyle\leq\|\mathbf{A}^{-1/2}\|_{2}\|\mathbf{B}^{-3/2}\|_{2}\|\mathbf{B}^{3/2}-\mathbf{A}^{3/2}\|_{2}+\|\mathbf{B}^{-3/2}\|_{2}\|\mathbf{A}-\mathbf{B}\|_{2}.$

Apply Lemma 8 in fukumizu2007statistical, we could bound
$\|\mathbf{B}^{3/2}-\mathbf{A}^{3/2}\|_{2}\leq
3\lambda^{1/2}\|\mathbf{A}-\mathbf{B}\|_{2}$. Then we complete the proof. ∎

##  Appendix J Discussion of the Transformation Operator

In this section, we discuss an alternative initialization approach as compared
to the approach proposed in Section 4.1. We aggregate the samples across nodes
$i\in V$ to initiate ${\mathbf{A}^{m}}$. The detail is documented in Appendix
J.1. In Appendix J.2, we discuss the construction of ${\mathbf{A}^{m}}$ varied
across nodes. In Appendix J.2, we discuss the generalization of the
initialization 4.3 to $M\geq 2$ modalities.

###  J.1 Aggregation of Samples Across Nodes

The initialization procedure introduced in Section 4.1 only use the sample
from node $1$. As an alternative, one could aggregate the samples across nodes
to conduct joint canonical correlation analysis. That is, we compute the
singular decomposition of

$\widehat{\mathbf{R}}_{12}^{\sharp}=\frac{1}{p}\sum_{i=1}^{p}(\widehat{\bm{\Sigma}}_{i,i}^{1})^{-1/2}\widehat{\bm{\Sigma}}_{i,i}^{12}(\widehat{\bm{\Sigma}}_{i,i}^{2})^{-1/2}.$

We denote $\widehat{\mathbf{V}}^{\sharp}_{1}$ as the matrix whose columns are
the top-k right singular vectors of $\widehat{\mathbf{R}}_{12}^{\sharp}$,
$\widehat{\mathbf{V}}^{\sharp}_{2}$ as the matrix whose columns are the top-k
left singular vectors of $\widehat{\mathbf{R}}_{12}^{\sharp}$ and
$\widehat{\bm{\Gamma}}^{\sharp}$ is a diagonal matrix whose diagonal entries
are top-k singular values of $\widehat{\mathbf{R}}_{12}^{\sharp}$. Then, we
initiate ${\mathbf{A}^{m}}$ as

${\mathbf{A}^{m}}^{\sharp}=(\widehat{\bm{\Gamma}}^{\sharp})^{-1/2}\widehat{\mathbf{V}}_{m}^{\sharp\top}\left(p^{-1}\sum_{i=1}^{p}\widehat{\bm{\Sigma}}^{m}_{ii}\right)^{-1/2}.$

The experimental setup is described in Section 7.2 and we test the experiments
for $20$ runs and take the average of the results. In the simulation, we vary
$N$ and measure the distance
$\sum_{m=1}^{M}\sigma_{A}^{-1}\|{\mathbf{A}^{m}}-\mathbf{A}^{m\star}\|_{F}^{2}$.
The result is displayed in Figure 4. The resulting figure indicates that the
normalized distance is consistent with respect to the sample size $N$.

Figure 4: We compare two approaches under $p=\\{50,100,150\\}$,
$\sigma_{A}=\max_{m=1,\ldots,M}\|\mathbf{A}^{m\star}\|_{2}$. The plots
indicate that both approaches have comparable performance, in particular in
the case when $N$ is large leading to small $(k\sum_{m=1}^{M}k_{m}\vee
k^{2}\log p)/N$.

###  J.2 Extension to Distinct Operators

Our model is easily extended to accommodate a distinct transformation operator
for each node $i$, i.e., $\mathbf{A}^{m}_{i}$. Then, the objective function
(3.7) becomes:

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

The initialization for $\\{\mathbf{A}_{i}^{m}\\}_{i=1,\ldots,p}$ can be
obtained by simply repeating (4.3) for each node. Similarly, we can tweak
Algorithm 1 to individually update $\mathbf{A}_{i}^{m}$ for this more flexible
model by alternately minimizing the objective (J.1).

###  J.3 Generalization to $M$ modalities

In this section, we discuss a simple extension that enables generalization to
$M\geq 2$ modalities. Note that Algorithm 1 can be applied to any $M$, and
hence we will focus on the generalization of the initialization method (4.3).
The initialization method (4.3) is the result of the maximum log-likelihood
estimator of the Gaussian linear latent model (bach2005probabilistic). It
follows that the generalization to $M$ modalities may be solved using the
maximum log-likelihood of the $M$-modal generative model. There are multiple
potential ways to carry out the estimation, e.g., the EM algorithm. Here, we
provide a simple extension from (4.3).
Given $m\in\\{1,\ldots,M\\}$, let $\backslash
m=\\{1,\ldots,M\\}\backslash\\{m\\}$. Note that we can generalize (4.2) to

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

where $\mathbf{y}^{\backslash
m}_{i}=(\mathbf{y}_{i}^{1\top},\ldots,\mathbf{y}_{i}^{m-1\top},\mathbf{y}_{i}^{m+1\top},\ldots,\mathbf{y}_{i}^{M\top})^{\top}$,
$\mathbf{L}^{\backslash
m\star}=(\mathbf{L}^{1\star\top}\cdots\mathbf{L}^{m-1\star\top}\mathbf{L}^{m+1\star\top}\cdots\mathbf{L}^{M\star\top})^{\top}$.
Here $\bm{\Sigma}_{i,i}^{\backslash m,\bf{q}}$ is a block-diagonal matrix
whose $j$-th block is $\bm{\Sigma}_{i,i}^{j,\bf{q}}$ for $j<m$ and
$\bm{\Sigma}_{i,i}^{j+1,\bf{q}}$ for $j\geq m$; $\bm{\Sigma}_{i,i}^{\backslash
m,\bm{\mu}}$ is a matrix of the covariance of
$(\bm{\mu}_{i}^{1\top},\ldots,\bm{\mu}_{i}^{m-1\top},\bm{\mu}_{i}^{m+1\top},\ldots,\bm{\mu}_{i}^{M\top})^{\top}$.
Then, we can view $(\mathbf{y}^{m}_{i},\mathbf{y}^{\backslash m}_{i})$ as two
modalities and apply Theorem 5.1 with $k_{2}=\sum_{i\in\backslash m}k_{i}$. We
generalize the result in the following paragraph. The idea is that since
Theorem 5.3 holds for any $M\geq 2$, if the generalized initialization method
can fulfill Assumption 8, then Theorem 5.4 can be generalized to any $M\geq
2$.
Define $\mathbf{R}^{m\backslash
m}=(\bm{\Sigma}^{m})^{-1/2}\bm{\Sigma}^{m\backslash m}(\bm{\Sigma}^{\backslash
m})^{-1/2}$. The following Assumption is a variant of Assumption 5.

######  Assumption 5’.

Suppose that there exists an $m\in\\{1,\ldots,M\\}$ such that
$\min(k_{m},\sum_{i\in\backslash m}k_{i})\geq k$ and the top-k singular values
of $\mathbf{R}^{m\backslash m}$ satisfy:
$\gamma_{1}>\gamma_{2}>\ldots>\gamma_{k}>0$.

######  Theorem 5.1’ (Generalization of Theorem 5.1).

Let
$\kappa_{m}=\sigma_{\max}(\bm{\Sigma}^{m})/\sigma_{\min}(\bm{\Sigma}^{m})$.
Suppose there exists a $m$ satisfies Assumption 5’ and Assumptions 6 holds.
Suppose that $N=O(\kappa_{m}^{2}k_{m}\vee\max_{i\in\backslash
m}\kappa_{i}^{2}\sum_{i\in\backslash m}k_{i})$. 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 for $m^{\prime}=1,\ldots,M$

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

with probability at least $1-5\exp(-\sum_{i=1}^{M}k_{i})$, where
$C_{\gamma_{k},\nu,\rho}>0$ is a constant.

###### Proof of Theorem 5.1’.

Note that if Assumption 5’ is satisfied, this implies that there exists a $m$
such that $\mathbf{R}^{m\backslash m}$ satisfies Assumption 5. Then we can
apply Theorem 5.1 and obtain the result. ∎

##  Appendix K Additional Simulation Results

In this section, we present additional simulation results. Section K.1
introduces the graph and data generation processes. Section K.2 introduces the
simulation process of the edge sets used in Graph 4. Section K.3–K.4 present
additional simulation results with different sample size and noise model.

###  K.1 Data Generation Procedures

We first introduce the precision structures followed by the data generation
processes.

Construction of inverse covariance ${\bm{\Omega}}$: Instead of generating
graphs from $\beta^{\star}_{ij}$, we directly construct sparse inverse
covariance operators. In simulations, we require the ranks of the operators to
be finite. Assume that the true latent space is $r$-dimensional. Note that $r$
and $k$ are different: $r$ is the true dimension that generates the process
and $k$ is the estimated (low-rank) dimension via Section 4.2. We follow graph
generation processes introduced in Zapata2019functional and
qiao2019functional. Let $\Psi\in\mathbb{R}^{r\times r}$ be a tri-diagonal
matrix such that $\Psi_{i,i}=1$ for $i=1,\ldots,r$ and
$\Psi_{i,i+1}=\Psi_{i+1,i}=0.5$ for $i=1,\ldots,r-1$. Let
${\bm{\Omega}}\in\mathbb{R}^{pr\times pr}=({\bm{\Omega}}_{i,j})$ be the
precision matrix with $r=9$ and ${\bm{\Omega}}_{i,j}\in\mathbb{R}^{r\times r}$
for $i,j=1,\ldots,p$. Since ${\bm{\Omega}}_{i,j}={\bf 0}$ if and only if
$\mathbf{B}_{ij}^{r\star}={\bf 0}$ for $i\neq j$, we consider consider the
following structures on ${\bm{\Omega}}$:

  * •

Graph 1: This model is similar to Graph 1 in Section $5.1$ in
(qiao2019functional). The diagonal blocks have and
${\bm{\Omega}}_{i,i}=\mathbf{I}$ for $i\in V$. For any $i=1,\ldots,p-1$, the
off-diagonal blocks have
${\bm{\Omega}}_{i,i+1}={\bm{\Omega}}_{i+1,i}=0.4\Psi$. For any
$i=1,\ldots,p-2$, the off-diagonal blocks have
${\bm{\Omega}}_{i,i+2}={\bm{\Omega}}_{i+2,i}=0.2\Psi$. For all other off-
diagonal blocks, we have ${\bm{\Omega}}_{i,j}={\bf 0}$.

  * •

Graph 2: This model has similar structure as the Graph 2 in Section $5.1$ in
(qiao2019functional) with the assumption that $p$ must be a constant multiple
of $10$. For $t=1,\ldots,p/10$ and let
${\bm{\Omega}}_{t}=({\bm{\Omega}}_{i,j})_{i=10(t-1)+1,\ldots,10t,j=10(t-1)+1,\ldots,10t}$
be a $10r\times 10r$ sub-matrix of ${\bm{\Omega}}$. If $t$ is an odd number,
${\bm{\Omega}}_{t}$ comes from Graph 1 with $p=10$. If $t$ is an even number,
then ${\bm{\Omega}}_{t}=\mathbf{I}$.

  * •

Graph 3: This model is the same as Graph 1 except that we have an additional
structure: ${\bm{\Omega}}_{i,i+3}={\bm{\Omega}}_{i+3,i}=0.1\Psi$ for
$i=1,\ldots,p-3$.

  * •

Graph 4: The structure is similar to example used in Zapata2019functional,
which violates the partial separability structure (Zapata2019functional). We
adopt the modified graph structure, Model D in zhao2021high as the graph
candidate for Graph 4. First, each node has the number of neighbors following
a power law distribution $f(y)=y^{-2}$ and the candidates of neighbors are
selected uniformly. Then, we partition the edge set $F$ into $r$ edge sets,
$F_{1},\ldots,F_{r}$, where the construction is described in Appendix A.4 of
(Zapata2019functional) and we restate the simulation procedure in Appendix K.2
for completeness. Given $F_{l}$, we construct a
$\widetilde{\bm{\Omega}}_{l}\in\mathbb{R}^{p\times p}$ precision matrix for
$l=1,\ldots,r$. The $ij$-th entry of $\widetilde{\bm{\Omega}}_{l}$ is
constructed as follows

$\widetilde{\bm{\Omega}}_{l,ij}=\left\\{\begin{array}[]{ll}1,&i=j;\\\
0,&(i,j)\not\in F_{l}\text{ or }i<j;\\\
\operatorname{Unif}\left(\left[-\frac{2}{3},-\frac{1}{3}\right]\cup\left[\frac{1}{3},\frac{2}{3}\right]\right),&(i,j)\in
F_{l}.\end{array}\right.$

We then normalize $\widetilde{\bm{\Omega}}_{l}$ such that each row of
$\widetilde{\bm{\Omega}}_{l}$ has unit norm. Then, we symmetrize
$\widetilde{\bm{\Omega}}_{l}$ by computing
$(\widetilde{\bm{\Omega}}_{l}+\widetilde{\bm{\Omega}}_{l}^{\top})/2$ and
setting the diagonal entries to be $1$. We define
$\bm{\Sigma}_{\text{ps}}={\rm diag}(\bm{\Sigma}_{1},\ldots,\bm{\Sigma}_{r})$
where $\bm{\Sigma}_{l}=3l^{-1.8}\widetilde{\bm{\Omega}}_{l}^{-1}$ for
$l=1,\ldots,r$. Define $\bar{\bm{\Omega}}\in\mathbb{R}^{pr\times pr}$ be a
precision matrix whose $l$-th $p\times p$ block diagonal is
$\bar{\bm{\Omega}}_{l,l}=\widetilde{\bm{\Omega}}_{l}$ and off-diagonal blocks
are
$\bar{\bm{\Omega}}_{l,l+1}=\bar{\bm{\Omega}}_{l+1,l}=\\{\widetilde{\bm{\Omega}}_{l}-{\rm
diag}(\widetilde{\bm{\Omega}}_{l})+\widetilde{\bm{\Omega}}_{l+1}-{\rm
diag}(\widetilde{\bm{\Omega}}_{l+1})\\}/2$. We then obtain

$\bm{\Omega}={\rm diag}(\bm{\Sigma}_{\text{ps}})^{-1/2}{{\rm
diag}(\bar{\bm{\Omega}})^{-1/2}\bar{\bm{\Omega}}{\rm
diag}(\bar{\bm{\Omega}})^{-1/2}}{\rm diag}(\bm{\Sigma}_{\text{ps}})^{-1/2}.$

Construction of ${\mathscr{A}^{m}}$: Let
$\mathbf{A}^{m,r}=(\mathbf{L}^{m,r})^{\dagger}\in\mathbb{R}^{r\times r_{m}}$
be the matrix representation of ${\mathscr{A}^{m}}$: for each entry
$a_{\ell\ell^{\prime}}^{m,r}$ in $\mathbf{A}^{m,r}$, we have
$a_{\ell\ell^{\prime}}^{m,r}=\langle{{\mathscr{A}^{m}}\phi_{\ell^{\prime}}^{m}},{\phi_{\ell}}\rangle$.
Noting that $\mathbf{A}^{m\star}$ is a sub-matrix of $\mathbf{A}^{m,r}$. To
construct $\mathbf{A}^{m,r}$ for $m=1,\ldots,M$, we first generate sparse
orthonormal rows, where the ratio of the non-zero entries are $1/3$. Then we
scale the magnitude of row $i$ in $\mathbf{A}^{m,r}$ with $0.2(i+1)+1$ for
$i=1,\ldots,r_{m}$. The covariance matrix of $\mathbf{x}^{m}$ for $m=1,2$ is
hence
$\\{(\mathbf{I}_{p}\otimes\mathbf{L}^{m,r}){\bm{\Omega}}^{-1}(\mathbf{I}_{p}\otimes(\mathbf{L}^{m,r})^{\top})\\}$.

Construction of Noise models: We consider two covariance structures of noise
models.

  * •

Noise Model 1: $\bm{\Sigma}^{m,\mathbf{q}}=\sigma\mathbf{I}$, where
$\mathbf{I}\in\mathbb{R}^{pr\times pr}$ for $n=1,\ldots,N$ and $m=1,\ldots,M$.
In the simulation, we set $\sigma=0.05$

  * •

Noise Model 2: $\Sigma^{m,\mathbf{q}}$ is a block-diagonal matrix with $p/10$
blocks. First, we generate a block-diagonal matrix ${\mathbf{F}^{m}}$ with
$p/10$ blocks and each row in a block is orthonormal to the other rows in the
same block. For $m=2$, we rotate ${\mathbf{F}^{m}}$ $90$ degree clockwise. Let
$\lambda_{m,1}\geq\ldots\geq\lambda_{m,pr}$ be the eigenvalues of
$\widetilde{\mathbf{F}}^{m}=({\mathbf{F}^{m}}+{\mathbf{F}^{m}})/2$. Then, we
make $\widetilde{\mathbf{F}}^{m}$ positive definite by taking the low-rank
parts of $\widetilde{\mathbf{F}}^{m}$ such that the corresponding eigenvalues
are greater than zero. Then, we normalize the remaining eigenvalues
$\lambda_{m,i}$ by
$0.01\lambda_{m,i}\\{\max_{j}(\lambda_{m,j})\\}^{-1}\sigma_{\max}(\bm{\Sigma}^{m,\mathbf{q}})$.

The major difference between Noise Model 1 and Noise Model 2 is that data from
all modalities corrupted with Noise Model 1 have identical graph structures in
the observational space. Moreover, when magnitude of the noise is small, i.e.,
$\sigma$ is small, the graph structures of the observed graphs and the latent
graph will look almost identical. In contrast, data corrupted with Noise Model
2 has very different observed graph structures between modalities. In
addition, the latent graph has distinct graph structure from the observed
graphs. We construct Noise Model 2 to mimic the real world situation that data
from different modalities are corrupted with distinct structured noise and to
test the robustness of the proposed model.
Data Simulation Process: We generate $N$ samples from the latent space,
$\mathbf{z}^{(1)},\ldots\mathbf{z}^{(N)}\in\mathbb{R}^{pr}$ under the
distribution $\mathcal{N}({\bf 0},{\bm{\Omega}}^{-1})$, where $\bm{\Omega}$ is
constructed by one of the Graph 1–4. Then, we compute $\mathbf{x}^{m,(n)}$ as
$\mathbf{x}^{m,(n)}=(\mathbf{I}_{p}\otimes\mathbf{L}^{m,r})\mathbf{z}^{(n)}$
for $n=1,\ldots,N$ and $m=1,\ldots,M$. The observed samples are

$\mathbf{y}^{m,(n)}=\mathbf{x}^{m,(n)}+\mathbf{q}^{m,(n)},\quad
n=1,\ldots,N,\;m=1,\ldots,M,$

where $\mathbf{q}^{m,(n)}\sim\mathcal{N}({\bf 0},\bm{\Sigma}^{m,\mathbf{q}})$
is generated independently from either of Noise Model 1 or Noise Model 2.

###  K.2 Simulation Details of Graph 4

Algorithm 3 Partition of the edge sets $F_{1},\ldots,F_{r}$

Input: $F$, $\tau\in[0,1]$

$F_{c}\leftarrow$ randomly select $\tau\left|F\right|$ of the edges from $F$

for $l=1,\ldots,r$ do

$F_{l}\leftarrow F_{c}$

$l\leftarrow 1$, $c\leftarrow 1$

for $e\in E\backslash E_{c}$ do

$F_{l}\leftarrow F_{l}\bigcup e$

$l\leftarrow l+1$

if $l>c$ then

$l\leftarrow 1$

$c\leftarrow(c+1)\text{mod }r$

This section discusses the details of the simulation process of the edge set
$F_{1},\ldots,F_{r}$ of Graph 4. We follow the same simulation process
introduced in Section A4.1 in Zapata2019functional and restate it here for
clarity. First, we generate a graph with the edge set $F$ whose edges follow
the power law distribution $f(y)=y^{-2}$. Then, we partition the edge set $F$
into $r$ edge sets $F_{1},\ldots,F_{r}$ such that $F=\bigcup_{l=1}^{r}F_{i}$.
The partition procedure is described in Algorithm 3.

###  K.3 Additional Results for Graph Simulations

Figure 5: The ROC curves of Graph 1–4. The additive noise is generated from
noise model 1.The AUC is discussed in Table 2. The proposed method has
consistent performance across four graphs and $p=\\{50,100,150\\}$.

We run the experiments for noise model 1. Figure 5 shows the ROC curves under
noise model 1. Table 3–2 indicate that the proposed method has gained in
performance in both the AUC and AUC15 under different graph settings and
dimension $p$. It is worth noting that the neighborhood regression approach on
single data modality proposed by zhao2021high has smaller AUC and AUC15
compared to our method, suggesting that integrating data modality might
improve the performance. Graph 2 has a simpler and sparser graph structure
compared to the other three graphs, and hence most estimators perform well in
this setting. In the case of Graph 1 and Graph 3 which have much more
complicated graph structures, qiao2019functional and moysidis2021joint achieve
much lower AUC. Finally, while most methods fail in Graph 4, our method still
retains good performance under $p=50,100,150$.

|  | AUC | AUC15
---|---|---|---
|  | Dimension (p)
Graph | Method | $50$ | $100$ | $150$ | $50$ | $100$ | $150$
Graph 1 | Proposed | $0.88(0.02)$ | $0.93(0.01)$ | $0.90(0.01)$ | $0.64(0.03)$ | $0.77(0.03)$ | $0.70(0.02)$
FGGM | $0.55(0.05)$ | $0.63(0.13)$ | $0.57(0.07)$ | $0.12(0.05)$ | $0.26(0.19)$ | $0.11(0.04)$
PSFGGM | $0.83(0.01)$ | $0.92(0.01)$ | $0.90(0.00)$ | $0.47(0.01)$ | $0.64(0.02)$ | $0.56(0.01)$
FPCA | $0.84(0.29)$ | $0.89(0.12)$ | $0.87(0.15)$ | $0.51(0.03)$ | $0.70(0.01)$ | $0.68(0.02)$
| JFGGM | $0.47(0.02)$ | $0.64(0.01)$ | – | $0.07(0.00)$ | $0.15(0.02)$ | –
Graph 2 | Proposed | $0.97(0.01)$ | $0.98(0.01)$ | $0.98(0.00)$ | $0.82(0.03)$ | $0.90(0.03)$ | $0.88(0.01)$
FGGM | $0.96(0.00)$ | $0.96(0.01)$ | $0.95(0.00)$ | $0.74(0.03)$ | $0.81(0.02)$ | $0.75(0.02)$
PSFGGM | $0.89(0.01)$ | $0.89(0.02)$ | $0.89(0.01)$ | $0.63(0.03)$ | $0.65(0.03)$ | $0.65(0.02)$
FPCA | $0.89(0.01)$ | $0.89(0.01)$ | $0.91(0.01)$ | $0.74(0.02)$ | $0.74(0.04)$ | $0.77(0.01)$
| JFGGM | $0.83(0.03)$ | $0.91(0.01)$ | – | $0.28(0.07)$ | $0.62(0.03)$ | –
Graph 3 | Proposed | $0.79(0.03)$ | $0.80(0.03)$ | $0.88(0.02)$ | $0.39(0.03)$ | $0.48(0.04)$ | $0.56(0.03)$
FGGM | $0.68(0.02)$ | $0.52(0.00)$ | $0.73(0.01)$ | $0.36(0.02)$ | $0.08(0.00)$ | $0.40(0.02)$
PSFGGM | $0.65(0.01)$ | $0.67(0.01)$ | $0.70(0.01)$ | $0.35(0.00)$ | $0.34(0.00)$ | $0.36(0.00)$
FPCA | $0.76(0.01)$ | $0.68(0.02)$ | $0.75(0.01)$ | $0.44(0.02)$ | $0.41(0.02)$ | $0.44(0.01)$
| JFGGM | $0.63(0.04)$ | $0.68(0.01)$ | – | $0.12(0.03)$ | $0.16(0.01)$ | –
Graph 4 | Proposed | $0.83(0.00)$ | $0.85(0.00)$ | $0.86(0.00)$ | $0.47(0.00)$ | $0.46(0.00)$ | $0.46(0.00)$
FGGM | $0.72(0.00)$ | $0.63(0.00)$ | $0.67(0.00)$ | $0.15(0.00)$ | $0.11(0.00)$ | $0.27(0.00)$
PSFGGM | $0.54(0.00)$ | $0.56(0.00)$ | $0.54(0.00)$ | $0.09(0.00)$ | $0.17(0.00)$ | $0.23(0.00)$
FPCA | $0.62(0.00)$ | $0.49(0.00)$ | $0.49(0.00)$ | $0.35(0.01)$ | $0.08(0.00)$ | $0.07(0.00)$
| JFGGM | $0.68(0.02)$ | $0.69(0.00)$ | – | $0.20(0.01)$ | $0.29(0.00)$ | –

Table 2: The average AUC of Graph 1–4 over 10 runs. The value inside the
parentheses denotes the standard deviation. The additive noise is generated
from noise model 1, where the AUC plot is shown in Figure 5. The proposed
method consistently achieves the largest AUC and AUC15 across four different
graphs and $p=\\{50,100,150\\}$.

|  | AUC | AUC15
---|---|---|---
|  | Dimension (p)
Graph | Method | $50$ | $100$ | $150$ | $50$ | $100$ | $150$
Graph 1 | Proposed | $0.87(0.02)$ | $0.92(0.02)$ | $0.87(0.01)$ | $0.59(0.02)$ | $0.72(0.04)$ | $0.62(0.02)$
FGGM | $0.60(0.02)$ | $0.76(0.01)$ | $0.64(0.01)$ | $0.17(0.02)$ | $0.45(0.02)$ | $0.15(0.02)$
PSFGGM | $0.83(0.01)$ | $0.91(0.01)$ | $0.90(0.01)$ | $0.46(0.01)$ | $0.63(0.02)$ | $0.56(0.01)$
FPCA | $0.84(0.01)$ | $0.89(0.00)$ | $0.85(0.02)$ | $0.52(0.03)$ | $0.70(0.01)$ | $0.65(0.02)$
| JFGGM | $0.48(0.01)$ | $0.66(0.01)$ | – | $0.07(0.00)$ | $0.18(0.02)$ | –
Graph 2 | Proposed | $0.96(0.01)$ | $0.98(0.01)$ | $0.98(0.00)$ | $0.82(0.02)$ | $0.88(0.02)$ | $0.82(0.02)$
FGGM | $0.95(0.01)$ | $0.96(0.01)$ | $0.96(0.01)$ | $0.71(0.02)$ | $0.78(0.03)$ | $0.76(0.03)$
PSFGGM | $0.88(0.02)$ | $0.89(0.00)$ | $0.89(0.01)$ | $0.63(0.03)$ | $0.64(0.02)$ | $0.65(0.02)$
FPCA | $0.89(0.02)$ | $0.90(0.01)$ | $0.90(0.02)$ | $0.70(0.02)$ | $0.76(0.03)$ | $0.75(0.02)$
| JFGGM | $0.85(0.02)$ | $0.89(0.03)$ | – | $0.33(0.07)$ | $0.64(0.04)$ | –
Graph 3 | Proposed | $0.77(0.05)$ | $0.80(0.48)$ | $0.83(0.03)$ | $0.38(0.04)$ | $0.48(0.05)$ | $0.52(0.05)$
FGGM | $0.71(0.02)$ | $0.52(0.00)$ | $0.73(0.01)$ | $0.36(0.02)$ | $0.09(0.01)$ | $0.40(0.02)$
PSFGGM | $0.66(0.01)$ | $0.67(0.01)$ | $0.70(0.01)$ | $0.35(0.00)$ | $0.33(0.01)$ | $0.37(0.01)$
FPCA | $0.76(0.02)$ | $0.69(0.01)$ | $0.74(0.00)$ | $0.45(0.01)$ | $0.41(0.01)$ | $0.44(0.01)$
| JFGGM | $0.66(0.03)$ | $0.69(0.02)$ | – | $0.15(0.02)$ | $0.18(0.01)$ | –
Graph 4 | Proposed | $0.84(0.00)$ | $0.84(0.00)$ | $0.86(0.00)$ | $0.47(0.00)$ | $0.45(0.00)$ | $0.48(0.00)$
FGGM | $0.72(0.00)$ | $0.63(0.00)$ | $0.67(0.00)$ | $0.15(0.00)$ | $0.11(0.00)$ | $0.27(0.00)$
PSFGGM | $0.53(0.00)$ | $0.56(0.00)$ | $0.51(0.00)$ | $0.10(0.00)$ | $0.14(0.00)$ | $0.12(0.00)$
FPCA | $0.61(0.00)$ | $0.48(0.00)$ | $0.51(0.00)$ | $0.34(0.00)$ | $0.07(0.00)$ | $0.08(0.00)$
| JFGGM | $0.75(0.04)$ | $0.63(0.00)$ | – | $0.23(0.03)$ | $0.24(0.00)$ | –

Table 3: The average AUC of Model 1–4 over 10 runs. The value inside the
parentheses denotes the standard deviation. The additive noise is generated
from noise model 2, where the AUC plot is shown in Figure 2. The rightmost
three columns denote the average AUC for FPR between $[0,0.15]$, normalized to
have a maximum area $1$. The proposed method achieves the largest AUC and
AUC15 across four different graphs and $p=\\{50,100,150\\}$. Figure 6: The AUC
of proposed model where data is corrupted by noise model 2. Left to right:
Graph 1–Graph 4

###  K.4 Experiment 1: p v.s. N

We plot the ROC curve of various graph with $p\in\\{50,100,150\\}$ and
$\in\\{50,100,200,500\\}$. The results are displayed in Figure 6–7. The
corresponding AUC is documented in Table 4 and Table 5, respectively. For
Graph 1–3, the AUC consistently increases as the sample size increases.
However, the sample size has mild effect on Graph 4 in both cases.

|  | AUC | AUC15
---|---|---|---
|  | Dimension (p)
Graph | Method | $50$ | $100$ | $150$ | $50$ | $100$ | $150$
Graph 1 | $N=50$ | $0.84(0.01)$ | $0.87(0.01)$ | $0.83(0.02)$ | $0.55(0.04)$ | $0.64(0.03)$ | $0.53(0.04)$
$N=100$ | $0.87(0.02)$ | $0.92(0.02)$ | $0.85(0.02)$ | $0.59(0.02)$ | $0.72(0.04)$ | $0.53(0.06)$
$N=200$ | $0.88(0.02)$ | $0.95(0.01)$ | $0.86(0.01)$ | $0.59(0.02)$ | $0.78(0.05)$ | $0.53(0.03)$
$N=500$ | $0.89(0.01)$ | $0.96(0.00)$ | $0.89(0.02)$ | $0.59(0.02)$ | $0.79(0.01)$ | $0.56(0.03)$
Graph 2 | $N=50$ | $0.94(0.02)$ | $0.96(0.01)$ | $0.96(0.01)$ | $0.75(0.03)$ | $0.81(0.02)$ | $0.81(0.03)$
$N=100$ | $0.96(0.01)$ | $0.98(0.01)$ | $0.98(0.00)$ | $0.82(0.02)$ | $0.88(0.02)$ | $0.88(0.02)$
$N=200$ | $0.98(0.00)$ | $0.99(0.00)$ | $0.99(0.00)$ | $0.87(0.01)$ | $0.93(0.01)$ | $0.94(0.01)$
$N=500$ | $0.98(0.00)$ | $0.99(0.00)$ | $1.00(0.00)$ | $0.90(0.01)$ | $0.96(0.00)$ | $0.97(0.00)$
Graph 3 | $N=50$ | $0.76(0.04)$ | $0.71(0.02)$ | $0.76(0.02)$ | $0.43(0.03)$ | $0.35(0.02)$ | $0.43(0.03)$
$N=100$ | $0.77(0.05)$ | $0.80(0.03)$ | $0.83(0.03)$ | $0.38(0.04)$ | $0.48(0.05)$ | $0.52(0.05)$
$N=200$ | $0.80(0.05)$ | $0.84(0.03)$ | $0.87(0.04)$ | $0.40(0.04)$ | $0.53(0.05)$ | $0.54(0.09)$
$N=500$ | $0.81(0.05)$ | $0.92(0.01)$ | $0.92(0.01)$ | $0.39(0.06)$ | $0.62(0.04)$ | $0.64(0.05)$
Graph 4 | $N=50$ | $0.86(0.00)$ | $0.84(0.00)$ | $0.84(0.00)$ | $0.44(0.00)$ | $0.44(0.00)$ | $0.45(0.00)$
$N=100$ | $0.84(0.00)$ | $0.84(0.00)$ | $0.86(0.00)$ | $0.47(0.00)$ | $0.45(0.00)$ | $0.48(0.00)$
$N=200$ | $0.76(0.00)$ | $0.83(0.00)$ | $0.86(0.00)$ | $0.47(0.00)$ | $0.46(0.00)$ | $0.46(0.00)$
$N=500$ | $0.76(0.00)$ | $0.83(0.00)$ | $0.86(0.00)$ | $0.47(0.00)$ | $0.46(0.00)$ | $0.46(0.00)$

Table 4: The average AUC of Graph 1–4 over 10 runs with different sample size.
The corresponding plots are displayed in Figure 6. The value inside the
parentheses denotes the standard deviation. The additive noise is generated
from noise model 2. Figure 7: The AUC of proposed model where data is
corrupted by noise model 1. Left to right: Graph 1–Graph 4

|  | AUC | AUC15
---|---|---|---
|  | Dimension (p)
Graph | Method | $50$ | $100$ | $150$ | $50$ | $100$ | $150$
Graph 1 | $N=50$ | $0.84(0.02)$ | $0.88(0.02)$ | $0.82(0.03)$ | $0.57(0.03)$ | $0.66(0.04)$ | $0.57(0.03)$
$N=100$ | $0.85(0.01)$ | $0.92(0.01)$ | $0.84(0.01)$ | $0.58(0.02)$ | $0.74(0.02)$ | $0.48(0.02)$
$N=200$ | $0.89(0.01)$ | $0.95(0.01)$ | $0.86(0.01)$ | $0.60(0.02)$ | $0.78(0.04)$ | $0.51(0.03)$
$N=500$ | $0.88(0.01)$ | $0.97(0.01)$ | $0.87(0.01)$ | $0.59(0.01)$ | $0.83(0.03)$ | $0.52(0.03)$
Graph 2 | $N=50$ | $0.95(0.01)$ | $0.96(0.01)$ | $0.95(0.01)$ | $0.76(0.03)$ | $0.80(0.02)$ | $0.79(0.02)$
$N=100$ | $0.97(0.01)$ | $0.98(0.01)$ | $0.98(0.00)$ | $0.82(0.03)$ | $0.90(0.03)$ | $0.88(0.01)$
$N=200$ | $0.98(0.00)$ | $0.99(0.00)$ | $0.99(0.00)$ | $0.88(0.01)$ | $0.93(0.01)$ | $0.94(0.01)$
$N=500$ | $0.98(0.00)$ | $0.99(0.00)$ | $1.00(0.00)$ | $0.90(0.01)$ | $0.96(0.00)$ | $0.97(0.00)$
Graph 3 | $N=50$ | $0.72(0.04)$ | $0.75(0.03)$ | $0.77(0.02)$ | $0.38(0.06)$ | $0.42(0.04)$ | $0.44(0.04)$
$N=100$ | $0.77(0.05)$ | $0.80(0.03)$ | $0.85(0.02)$ | $0.38(0.04)$ | $0.48(0.04)$ | $0.55(0.04)$
$N=200$ | $0.81(0.02)$ | $0.87(0.02)$ | $0.90(0.01)$ | $0.38(0.03)$ | $0.52(0.06)$ | $0.60(0.03)$
$N=500$ | $0.84(0.03)$ | $0.92(0.02)$ | $0.94(0.01)$ | $0.40(0.04)$ | $0.61(0.04)$ | $0.69(0.03)$
Graph 4 | $N=50$ | $0.86(0.00)$ | $0.84(0.00)$ | $0.85(0.00)$ | $0.45(0.00)$ | $0.45(0.00)$ | $0.45(0.00)$
$N=100$ | $0.83(0.00)$ | $0.85(0.00)$ | $0.86(0.00)$ | $0.47(0.00)$ | $0.46(0.00)$ | $0.46(0.00)$
$N=200$ | $0.88(0.00)$ | $0.86(0.00)$ | $0.86(0.00)$ | $0.48(0.00)$ | $0.46(0.00)$ | $0.46(0.00)$
$N=500$ | $0.87(0.00)$ | $0.82(0.02)$ | $0.86(0.00)$ | $0.49(0.00)$ | $0.48(0.00)$ | $0.45(0.00)$

Table 5: The average AUC of Graph 1–4 over 10 runs with different sample size.
The corresponding plots are displayed in Figure 7. The value inside the
parentheses denotes the standard deviation. The additive noise is generated
from noise model 1. Figure 8: Distance v.s. statistical error. We fix $N=324$
for $p=50,100$ and $N=289$ for $p=150$ and vary $k={3,5,7,9}$. When the value
on the $x$-axis is small, both the average error and maximum error are small
in all graphs. The errors increase as the value on the $x$-axis increases.

###  K.5 Distance v.s. $k$

In addition to the sample complexity experiment discussed in Section 7.2, we
verify Theorem 5.3 with varying $k$ as well. We choose $k=\\{3,5,7,9\\}$,
where we set $k=k_{m}$ and $N=324$ here. We run the simulation for $20$
independent simulated datasets and the average result is shown in Figure 8.
Note that the lines are nearly linear for three different graphs and
$p=50,100,150$, supporting the result from Theorem 5.3.

###  K.6 Sensitivity of the Tuning Parameters

Figure 9: Top row: The average residuals of signals from data modality 1 after
projecting to $k_{1}$ number of basis functions, where the ground truth is
$k_{1}=9$. Bottom row: The average residuals of signals from data modality 2
after projecting to $k_{2}$ number of basis functions, where the ground truth
is $k_{2}=9$. Each line is plotted as the average result of $N=100$ subjects
and over $10$ independent runs. The standard deviation of each point is at the
scale of $0.1$. Both figures indicate that there is a turning point in the
residual when $k_{1}$, $k_{2}$ exceed the $9$, suggesting that the elbow
method could give us relatively accurate estimates.  Figure 10:  The top-$k$
canonical correlations averaged over $10$ independent runs, where the true
$k=9$. The canonical correlation plot has an abrupt drop to $0$ at $k=10$,
when $k$ exceeds the ground truth value.

Figure 11: The BIC score plot with $5$-fold cross-validation. In this
simulation, we choose $N=100$, $p=100$, $s=4$, and $\alpha=0.44$. Both the
plots from the top row and bottom row suggest that the score is sensitive to
the values of $\tau_{1}$ but indifferent to the choices of $\tau_{2}$. Figure
12: The sample complexity plot with varying $\alpha$. We select other tunning
parameters $\\{s,\tau_{1},\tau_{2}\\}$ using the procedure discussed in
Section 4.2 and apply $(4.3)$ to initialize $\mathbf{A}^{m}$. The true
$\alpha^{\star}$ is $0.33$ for Graph 1–2, and $0.44$ for Graph 3. As the plots
indicate, under-selection of $\alpha$ leads to a larger error when $x$
approaches zero, where in this case we have a large sample size $N$. In
contrast, large $\alpha$ has a weaker influence on the plot. We find that the
correct choice of $\alpha=\alpha^{\star}$ will lead to optimal results. Figure
13: The ROC plot of varying $k$ (ground truth $k=9$ and $k_{1}=k_{2}=9$). We
use $N=100$ and each plot is the average result over $10$ independent runs.
The underlying AUCs are smaller when $k$ and $k_{m}$ are greatly under-
selected. We found that when $k$ is slightly under-selected, i.e. $k=7$, the
AUCs are close to the AUCs of $k=9$ in most settings.  Figure 14: The ROC plot
of varying $k$ and $k_{m}$ with $k_{m}\geq k$. The ground truth values of $k$
and $k_{m}$ are $9$. We use $N=100$, $p=50$, and noise model 1. Each plot is
the average result over $10$ independent runs. Top row: when both $k_{1}$ and
$k_{2}$ are significantly under-selected, the choice of $k$ is insensitive to
the ROC. Centered two rows: when both $k_{1}$ and $k_{2}$ are slightly under-
selected or equal to the true value, larger $k$ has a larger AUC. Bottom row:
when both $k_{1}$ and $k_{2}$ are over-selected, the ROC curves of $k=9$ and
$k=11$ almost overlap. It only has a suboptimal curve in Graph 3.

In this section, we analyze the variable selection method introduced in
Section 4.2. First, we discuss the practice of using elbow method to select
$k$ and $k_{m}$. We then discuss how mis-specification of $k$ and $k_{m}$
affects the results. Finally, we discuss the sensitivity of choosing
$\tau_{1},\tau_{2}$ and $\alpha$ in the following.
First, we run the elbow algorithm discussed in Section 4.2 to select $k$ and
$k_{m}$. The result is displayed in Figure 9–10, where there are clear turning
points in all cases. Furthermore, the elbow points match the true values.
Figure 13 shows the ROC plot when $k$ is under-selected, and we select $k_{m}$
based on the elbow method. It shows that the AUCs are smaller when $k$ is
significantly smaller than the true $k=9$. In the case when $k=7$, the AUCs
are close to the case when $k=9$, as Table 6 shows. In Figure 14, we vary
$k_{m}=\\{5,7,9,11\\}$ and $k=\\{3,5,7,9,11\\}$. Furthermore, we restrict
$k_{m}\geq k$ due to the technical constraints of CCA. The corresponding AUC
is documented in Table 7. The result indicates that if the difference between
$k_{m}$ and $k$ is small, i.e., $k_{m}-k=2$, the underlying AUC is close to
the case when $k_{m}=k$. However, if $k$ is much smaller than $k_{m}$, we see
the decay in the underlying AUC. When we over-select $k_{m}=11$, the
corresponding AUCs of $k=9$ and $k=11$ are close to the case when $k_{m}=k=9$.
However, under-selection of $k$ results in smaller AUCs.
We select the candidates of $\tau_{1}$ to be $\\{1/n:n=1,\ldots,5\\}$ and the
candidates for $\tau_{2}$ to be $\\{n:n=1,\ldots,5\\}$. Then, we apply 5-fold
cross-validation with BIC metric to select the optimal values. In the
following simulation, we fix $s=4$ and $\alpha=4/9$ and plot the BIC score
with respect to $\tau_{1}$ and $\tau_{2}$, shown in Figure 11. The plot
indicates that the score is insensitive to the upper bound, the choice of
$\tau_{2}$, and more sensitive to the lower bound, the choice of $\tau_{1}$.
Finally, we find varying $\alpha$ is less sensitive to the result so long as
$\alpha\geq\alpha^{\star}$. We present the simulation result in Figure 12. The
simulation details are discussed in Section 7.2.

|  | AUC | AUC15
---|---|---|---
|  | Dimension (p)
Graph | Method | $50$ | $100$ | $150$ | $50$ | $100$ | $150$
Graph 1 | $k=3$ | $0.79(0.02)$ | $0.87(0.04)$ | $0.80(0.02)$ | $0.44(0.02)$ | $0.62(0.04)$ | $0.47(0.04)$
$k=5$ | $0.84(0.01)$ | $0.91(0.02)$ | $0.80(0.02)$ | $0.50(0.03)$ | $0.69(0.06)$ | $0.45(0.05)$
$k=7$ | $0.85(0.03)$ | $0.92(0.01)$ | $0.85(0.02)$ | $0.55(0.03)$ | $0.72(0.02)$ | $0.52(0.06)$
$k=9$ | $0.85(0.02)$ | $0.92(0.01)$ | $0.85(0.02)$ | $0.55(0.03)$ | $0.74(0.03)$ | $0.52(0.06)$
Graph 2 | $k=3$ | $0.96(0.01)$ | $0.97(0.01)$ | $0.97(0.01)$ | $0.81(0.04)$ | $0.86(0.02)$ | $0.88(0.02)$
$k=5$ | $0.96(0.01)$ | $0.98(0.00)$ | $0.97(0.01)$ | $0.81(0.02)$ | $0.89(0.01)$ | $0.88(0.02)$
$k=7$ | $0.96(0.01)$ | $0.98(0.01)$ | $0.98(0.00)$ | $0.83(0.03)$ | $0.88(0.03)$ | $0.89(0.02)$
$k=9$ | $0.97(0.00)$ | $0.98(0.00)$ | $0.98(0.00)$ | $0.82(0.02)$ | $0.90(0.01)$ | $0.88(0.01)$
Graph 3 | $k=3$ | $0.72(0.04)$ | $0.78(0.02)$ | $0.75(0.04)$ | $0.44(0.02)$ | $0.37(0.04)$ | $0.43(0.03)$
$k=5$ | $0.71(0.05)$ | $0.78(0.01)$ | $0.80(0.04)$ | $0.37(0.05)$ | $0.36(0.01)$ | $0.44(0.05)$
$k=7$ | $0.81(0.04)$ | $0.79(0.01)$ | $0.81(0.02)$ | $0.42(0.04)$ | $0.40(0.03)$ | $0.45(0.03)$
$k=9$ | $0.78(0.06)$ | $0.79(0.02)$ | $0.84(0.04)$ | $0.41(0.05)$ | $0.43(0.03)$ | $0.51(0.06)$
Graph 4 | $k=3$ | $0.80(0.00)$ | $0.81(0.00)$ | $0.84(0.00)$ | $0.38(0.00)$ | $0.41(0.00)$ | $0.41(0.00)$
$k=5$ | $0.83(0.00)$ | $0.84(0.00)$ | $0.85(0.00)$ | $0.43(0.00)$ | $0.45(0.00)$ | $0.44(0.00)$
$k=7$ | $0.85(0.00)$ | $0.85(0.00)$ | $0.86(0.00)$ | $0.45(0.00)$ | $0.45(0.00)$ | $0.46(0.00)$
$k=9$ | $0.86(0.00)$ | $0.85(0.00)$ | $0.86(0.00)$ | $0.46(0.00)$ | $0.46(0.00)$ | $0.46(0.00)$

Table 6: The average AUC of Model 1–4 over 10 runs. The ground truth
$k=k_{m}=9$. We vary $k=\\{3,5,7,9\\}$. The value inside the parentheses
denotes the standard deviation. The additive noise is generated from noise
model 1, where the AUC plot is shown in Figure 13. The rightmost three columns
denote the average AUC for FPR between $[0,0.15]$ divided by $0.15$. The
division is a normalization such that the maximum area will be $1$. When $k$
is significant under selected, $k=3$, the AUC and AUC15 are smaller. In
contrast, the AUC and AUC15 increase as $k$ increases close to the true value
$9$.

|  | AUC | AUC15
---|---|---|---
|  | $k_{m}$
Graph | $k$ | $5$ | $7$ | $9$ | $11$ | $5$ | $7$ | $9$ | $11$
Graph 1 | $3$ | $0.87(0.02)$ | $0.80(0.03)$ | $0.79(0.02)$ | $0.79(0.02)$ | $0.60(0.02)$ | $0.48(0.05)$ | $0.44(0.02)$ | $0.45(0.05)$
$5$ | $0.87(0.02)$ | $0.82(0.01)$ | $0.84(0.01)$ | $0.81(0.02)$ | $0.61(0.04)$ | $0.51(0.01)$ | $0.50(0.03)$ | $0.47(0.04)$
$7$ | – | $0.84(0.02)$ | $0.85(0.03)$ | $0.85(0.02)$ | – | $0.54(0.03)$ | $0.55(0.03)$ | $0.56(0.05)$
$9$ | – | – | $0.85(0.02)$ | $0.86(0.02)$ | – | – | $0.55(0.03)$ | $0.57(0.05)$
$11$ | – | – | – | $0.86(0.02)$ | – | – | – | $0.57(0.05)$
Graph 2 | $3$ | $0.96(0.01)$ | $0.95(0.01)$ | $0.96(0.01)$ | $0.95(0.01)$ | $0.81(0.03)$ | $0.80(0.05)$ | $0.81(0.04)$ | $0.80(0.03)$
$5$ | $0.97(0.01)$ | $0.97(0.01)$ | $0.96(0.01)$ | $0.97(0.00)$ | $0.82(0.03)$ | $0.81(0.03)$ | $0.81(0.02)$ | $0.82(0.02)$
$7$ | – | $0.96(0.01)$ | $0.96(0.01)$ | $0.97(0.01)$ | – | $0.80(0.05)$ | $0.83(0.03)$ | $0.82(0.02)$
$9$ | – | – | $0.97(0.00)$ | $0.97(0.01)$ | – | – | $0.82(0.02)$ | $0.82(0.03)$
$11$ | – | – | – | $0.97(0.00)$ | – | – | – | $0.82(0.03)$
Graph 3 | $3$ | $0.78(0.03)$ | $0.76(0.02)$ | $0.72(0.04)$ | $0.75(0.03)$ | $0.45(0.04)$ | $0.42(0.03)$ | $0.44(0.03)$ | $0.45(0.06)$
$5$ | $0.76(0.04)$ | $0.74(0.04)$ | $0.71(0.05)$ | $0.73(0.04)$ | $0.39(0.04)$ | $0.39(0.03)$ | $0.37(0.05)$ | $0.38(0.03)$
$7$ | – | $0.74(0.03)$ | $0.81(0.04)$ | $0.76(0.05)$ | – | $0.38(0.04)$ | $0.42(0.04)$ | $0.38(0.04)$
$9$ | – | – | $0.78(0.06)$ | $0.79(0.04)$ | – | – | $0.41(0.05)$ | $0.41(0.04)$
$11$ | – | – | – | $0.77(0.04)$ | – | – | – | $0.42(0.03)$
Graph 4 | $3$ | $0.81(0.00)$ | $0.80(0.00)$ | $0.80(0.00)$ | $0.80(0.00)$ | $0.40(0.00)$ | $0.40(0.00)$ | $0.38(0.00)$ | $0.38(0.00)$
$5$ | $0.84(0.00)$ | $0.84(0.00)$ | $0.83(0.00)$ | $0.83(0.00)$ | $0.42(0.00)$ | $0.43(0.00)$ | $0.43(0.00)$ | $0.43(0.00)$
$7$ | – | $0.86(0.00)$ | $0.85(0.00)$ | $0.85(0.00)$ | – | $0.44(0.00)$ | $0.45(0.00)$ | $0.45(0.00)$
$9$ | – | – | $0.86(0.00)$ | $0.86(0.00)$ | – | – | $0.46(0.00)$ | $0.46(0.00)$
$11$ | – | – | – | $0.86(0.00)$ | – | – | – | $0.46(0.00)$

Table 7: The average AUC of Graph 1–4 over 10 runs. The ground truth
$k=k_{m}=9$. We vary $k=\\{3,5,7,9,11\\}$ and $k_{m}=\\{5,7,9,11\\}$. Since
$k$ must be smaller or equal to $k_{m}$ due to the initialization method, we
only show the results for $k\leq k_{m}$. The value inside the parentheses
denotes the standard deviation. The additive noise is generated from noise
model 1, where the AUC plot is shown in Figure 14. The rightmost three columns
denote the average AUC for FPR between $[0,0.15]$, normalized to have a
maximum area $1$. Fix $k_{m}$, both AUC and AUC15 are optimal when $k$ is
close to $k_{m}$.

##  Appendix L Treatment to Discrete Observed Data

Throughout the paper we have assumed that observations across modalities are
continuous functions. As a result, this article focuses mainly on the
construction of the latent model, and we assume continuous observations for
simplicity. In practice, data from different modalities are expected to be
recorded with different temporal resolutions and can be viewed as discrete
data. When the observations are discrete, we can still compute the function
score by projecting the discrete samples to the discretized basis functions.
Recall the basis $\\{\phi_{1}^{m},\ldots,\phi^{m}_{k_{m}}\\}$ and let
$t_{i1}^{m,(n)},\ldots,t_{iq}^{m,(n)}$ be the sampling time points of subject
$n$ of the modality $m$ at node $i$. We can obtain the function score by
solving a least-squares problem (cf. Section 3.1 of zhao2022fudge) and obtain

$\widetilde{\mathbf{y}}^{m,(n)}_{i}=(\mathbf{C}_{i}^{m,(n)\top}\mathbf{C}_{i}^{m,(n)})^{-1}\mathbf{C}_{i}^{m,(n)\top}\mathbf{h}_{i}^{m,(n)},$

where

$\mathbf{C}_{i}^{m,(n)}=\begin{bmatrix}\phi_{1}^{m}(t_{i1}^{m,(n)})&\ldots&\phi_{k_{m}}^{m}(t_{i1}^{m,(n)})\\\
\vdots&\ddots&\vdots\\\
\phi_{1}^{m}(t_{iq}^{m,(n)})&\ldots&\phi_{k_{m}}^{m}(t_{iq}^{m,(n)})\end{bmatrix},\quad\mathbf{h}_{i}^{m,(n)}=\begin{bmatrix}\mathcal{Y}_{i}^{m,(n)}(t_{i1}^{m,(n)})\\\
\vdots\\\ \mathcal{Y}_{i}^{m,(n)}(t_{iq}^{m,(n)}).\end{bmatrix}$

Then, we can replace $\mathbf{Y}_{i}^{m}$ with
$\widetilde{\mathbf{Y}}_{i}^{m}=(\widetilde{\mathbf{y}}_{i}^{m,(1)},\ldots,\widetilde{\mathbf{y}}_{i}^{m,(N)})$
in (3.7). The question is then how well can we estimate the parameters with
the new objective function? zhao2022fudge analyzed the conditions when the
covariance of $\widetilde{\mathbf{y}}_{i}^{m}$ is close to the covariance of
${\mathbf{y}}_{i}^{m}$. This provides insight that we might be able to
quantify the error due to discretization under some regularity conditions,
i.e., when the basis functions have smooth structures and the samples are
evenly spaced.

The family of basis functions is often unknown in practice. Hence, one
approach uses functional PCA (fPCA) to estimate the basis functions. Existing
literature has studied this setting from both theoretical and methodological
perspectives (yao2005functional; li2010uniform; cai2011optimal;
amini2012sampled; zhang2016sparse). Motivated by qiao2020doubly who studied
functional graphical model under the setting of the discrete sample, we
briefly outline an extension of the methodology from yao2005functional to our
model and propose a simple treatment.In this first stage, we apply the
algorithm proposed by yao2005functional to estimate the functional score
$\mathbf{y}_{i}^{m,(n)}$, denoted as
$\widehat{\mathbf{y}}^{m,(n)}_{i}=(\widehat{y}_{i,1}^{m,(n)},\ldots,\widehat{y}_{i,k_{m}}^{m,(n)})^{\top}$
for $i\in V$ and the estimated basis functions
$\\{\widehat{\phi}_{1}^{m},\ldots,\widehat{\phi}_{k_{m}}^{m}\\}$ individually
for each modality $m=1,\ldots,M$. Assume we observe $J$ discrete samples
randomly sampled at $T_{1}^{m,(n)},\ldots,T_{J}^{m,(n)}$:

$\mathcal{Y}_{i}\left(T_{j}^{m,(n)}\right)=\chi_{i}\left(T_{j}^{m,(n)}\right)+\xi_{i}\left(T_{j}^{m,(n)}\right),$

where $\mathcal{Y}_{i}({T_{j}^{m,(n)}})$ is a shorthand for
$\mathcal{Y}_{i}^{m,(n)}({T_{j}^{m,(n)}})$,
$\chi_{i}({T_{j}^{m,(n)}})=\chi_{i}^{m,(n)}({T_{j}^{m,(n)}})$ and
$\xi_{i}({T_{j}^{m,(n)}})=\xi_{i}^{m,(n)}({T_{j}^{m,(n)}})$. Define
$\mathscr{K}^{m}_{i}(u,v)$ as the covariance estimator of
$\mathcal{Y}^{m}_{i}(u)$ and $\mathcal{Y}^{m}_{i}(v)$. The first step is to
estimate $\mathscr{K}^{m}_{i}(u,v)$ with discrete observations. Let $h>0$ be a
smoothing constant and denote $K_{h}(\cdot)=h^{-1}K(\cdot/h)$ be a smoothing
kernel function. Given $u,v$, we consider the minimization of the following
function with respect to $(\beta_{0},\beta_{1},\beta_{2})$:

$\sum_{n=1}^{N}\sum_{T_{j}^{m,(n)}\neq
T_{j^{\prime}}^{m,(n)}}\left\\{\mathcal{Y}_{i}\left(T_{j}^{m,(n)}\right)\mathcal{Y}_{i}\left(T_{j^{\prime}}^{m,(n)}\right)-\beta_{0}-\beta_{1}\left(T_{j}^{m,(n)}-u\right)-\beta_{2}\left(T_{j^{\prime}}^{m,(n)}-v\right)\right\\}^{2}\\\
\times
K_{h}\left(T_{j}^{m,(n)}-u\right)K_{h}\left(T_{j^{\prime}}^{m,(n)}-v\right).$
(L.1)

Define $\widehat{\beta}_{0}$ be the optimal solution of $\beta_{0}$ in (L.1).
The estimator is obtained as
$\widehat{\mathscr{K}}_{i}^{m}(u,v)=\widehat{\beta}_{0}$. Let
$\\{(\widehat{\lambda}_{i,\ell},\widehat{\phi}^{m}_{i,\ell})\\}_{\ell=1}^{k_{m,i}^{\prime}}$
be the eigen-pairs of $\widehat{\mathscr{K}}_{i}^{m}(\cdot,\cdot)$. Define
$\\{\widehat{\phi}_{1}^{m},\ldots,\widehat{\phi}_{k_{m}}^{m}\\}$ as the set of
orthonormal basis functions that spans
$\\{\widehat{\phi}_{1,1}^{m},\ldots,\widehat{\phi}_{1,k_{m,1}^{\prime}}^{m},\ldots,\widehat{\phi}_{p,1}^{m},\ldots,\widehat{\phi}_{p,k_{m,p}^{\prime}}^{m}\\}$.
Let

$\widetilde{\bm{\Sigma}}^{m,(n)}_{i}=[\widehat{\mathscr{K}}_{i}^{m}(u,v)+\widehat{\sigma}^{2}I(u=v)]_{u,v\in\\{T_{1}^{m,(n)},\ldots,T_{J}^{m,(n)}\\}}\in\mathbb{R}^{J\times
J},$

where the selection of $\widehat{\sigma}$ has been discussed in
yao2005functional,

$\widetilde{\mathbf{y}}_{i}^{m,(n)}=(\mathcal{Y}_{i}({T_{1}^{m,(n)}}),\ldots,\mathcal{Y}_{i}({T_{J}^{m,(n)}}))^{\top}\in\mathbb{R}^{J}.$

Then, the estimator (yao2005functional) for ${y}_{i,\ell}^{m,(n)}$ under the
discrete setting is

$\widehat{y}_{i,\ell}^{m,(n)}=\widehat{\mathbf{d}}_{i,\ell}^{m,(n)\top}(\widetilde{\bm{\Sigma}}^{m,(n)}_{i})^{-1}\widetilde{\mathbf{y}}_{i}^{m,(n)},\quad\ell=1,\ldots,k_{m},$

where
$\widehat{\mathbf{d}}_{i,\ell}^{m,(n)}=(\widehat{\mathbf{d}}_{i,\ell}^{m}(T_{1}^{m,(n)}),\ldots,\widehat{\mathbf{d}}_{i,\ell}^{m}(T_{J}^{m,(n)}))^{\top}\in\mathbb{R}^{J}$
and

$\widehat{\mathbf{d}}_{i,\ell}^{m}(T_{j}^{m,(n)})=\int_{u}\widehat{\mathscr{K}}_{i}^{m}(T_{j}^{m,(n)},u)\widehat{\phi}^{m}_{\ell}(u)du,\quad
j=T_{1}^{m,(n)},\ldots,T_{J}^{m,(n)}.$

After obtaining the estimates
$\widehat{\mathbf{Y}}_{i}^{m}=(\widehat{\mathbf{y}}_{i}^{m,(1)},\ldots,\widehat{\mathbf{y}}_{i}^{m,(N)})\in\mathbb{R}^{k_{m}\times
N}$ from Stage 1, we can replace $\mathbf{Y}_{i}^{m}$ with
$\widehat{\mathbf{Y}}_{i}^{m}$ in (3.7). Then, the rest of the estimation
procedures follow the proposal.

##  Appendix M Experiments on Real Data

This section discusses the implementation details of the concurrent EEG-fMRI
recordings (sadaghiani2010intrinsic). Section M.1 introduces the data
preprocessing pipeline. Section M.2 discusses the details of the regression
procedure introduced in Section 8. Section M.3 shows the visualization of the
precision matrices.

###  M.1 Preprocessing Pipeline

First, we conduct the z-transform of the time-series for both EEG and fMRI
data, which is a standard preprocessing step. In the second step, we project
the data to different bases as outlined in the following. We first regress out
the global signal and standardize each time-series. The truncated fMRI time-
series has $143$ time points and the original EEG time-series has $71680$ time
points, after removing a few time points in the beginning
(poldrack2011handbook). We span the fMRI data using the Fourier basis
functions. Then, we down-sample the EEG data evenly to $1024$ time points so
that we can project it to wavelet family basis functions, a common basis
family used to decompose EEG signals (gandhi2011comparative). The candidates
of the wavelet bases are Daubechies’ extremal phase wavelets, Daubechies’
“least-asymmetric” wavelets, and Coiflets wavelets, provided by the
‘wavethresh’ R-package. To select the best wavelet basis, we use the Shannon
entropy type function as the evaluation metric. The steps are as the
following. For each subject $n$ and each region $i$, we compute the wavelet
coefficients $w^{m,(n)}_{i,\ell}$ for $\ell=1,2,\ldots,1024$. We normalize
each coefficient as
$\bar{w}_{\ell}^{m}=\sum_{i,n}w^{m,(n)}_{i,\ell}/\max_{\ell}(\sum_{i,n}w^{m,(n)}_{i,\ell})$.
We then select the basis family that has the smallest entropy
$-\sum_{\ell}(\bar{w}_{i}^{m})^{2}\log(\bar{w}_{\ell}^{m})^{2}$. The rest of
the tuning parameters are selected using the procedure discussed in Section
4.2.

###  M.2 Details of the Neighborhood Regression Procedure

The functional neighborhood regression method for data of single modality
follows closely from (3.7), except we do not need to estimate the
transformation operator ${\mathbf{A}^{m}}$ here. Hence, given
$\mathbf{y}_{i}^{m,(1)},\ldots,\mathbf{y}_{i}^{m,(N)}$ independent samples and
let
${\mathbf{Y}_{i}^{m}}=(\mathbf{y}_{i}^{m,(1)},\ldots,\mathbf{y}_{i}^{m,(N)})\in\mathbb{R}^{k_{m}\times
N}$ for $i\in V$, we define

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

Then, we optimize the objective function:

$\displaystyle\mathop{\mathrm{argmin}}_{\mathbf{B}}g_{\mathbf{Y}}(\mathbf{B})$
$\displaystyle\text{s.t.}\quad\mathbf{B}_{i}\in\mathcal{K}_{B}(s,\alpha),\quad
i\in V.$ (M.1)

The optimization problem (M.2) is carried out by projected gradient descent,
as shown in Algorithm 4. Next, we use the estimated $\widehat{\mathbf{B}}_{i}$
for $i\in V$ to construct the edge set via AND operation (4.1). Finally, we
estimate the inverse covariance matrix by solving (8.1).

Algorithm 4 Sparse Neighborhood Regression

Input: $\\{\mathbf{Y}^{m}\\}_{m=1,\ldots,M}$

for $i\in V$ do

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

while Not converged do

for $i\in V$ do

$\mathbf{B}_{i}^{(t+.25)}\leftarrow\mathbf{B}_{i}^{(t)}-\eta_{B}\nabla_{\mathbf{B}_{i}}g_{\mathbf{Y}}$

$\mathbf{B}_{i}^{(t+.5)}\leftarrow{\mathcal{T}}_{\vartheta_{1}s^{\star}/2}(\mathbf{B}_{i}^{(t)})$

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

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

Figure 15:  Estimated graph structure from three edge set candidates. Left to
right: the graphs of data modality 1 and the graphs of data modality 2. Top to
bottom: the individual edge set, the latent edge set, and the fused edge set.
The block diagonal box presents the partition using Yeo 7 network
(yeo2011organization).

###  M.3 Visualization of Precision Matrices

As the number of bases used to span two data modalities is different, the
precision matrices $\widehat{\Omega}^{1}\in\mathbb{R}^{pk_{1}\times pk_{1}}$
and $\widehat{\Omega}^{2}\in\mathbb{R}^{pk_{2}\times pk_{2}}$ have different
dimensions. This makes it hard to compare the precision matrices of two
graphs. As an alternative, given estimated $\widehat{\Omega}^{m}$ estimated
from (8.1) with
$E\in\\{\widehat{E}^{m},\widehat{E}^{z},\widehat{E}^{m}\cup\widehat{E}^{z}\\}$
for $m=1,2$, we compute the magnitude of the $(i,j)$th submatrix
$\widehat{\omega}_{i,j}^{m}$ of $\widehat{\Omega}^{m}$, in Frobenious norm.
The resulting plots are shown Figure 15, where the $(i,j)$th is the magnitude
for $\omega_{i,j}^{m}$ and the label is documented in Table
LABEL:tab:DK_label. The figure indicates that the latent dependency graph
includes some edge connections that were not originally detected by any
modality individually.

###  M.4 Computation Complexity

In this section, we analyze the computation time for running the experiments.
The sample size, dimension, and number of time points are documented in
Section 8. The computing time of the model selection of
$\\{s,\alpha,\tau_{1},\tau_{2}\\}$ introduced in Section 4.2 is shown in the
step 5 of Table 8. We run grid search on
$s=\\{3,6,9,12,15,18,21,24,27,30,33\\}$, $\alpha=\\{3,7,10,13\\}$,
$\tau_{1}=\\{0.5,0.25\\}$ and $\tau_{2}=\\{2,4\\}$. The computing time of edge
selection (proposed algorithm) is shown in Step 7. In the experiment, we use
learning rate $10^{-4}$ to initialize $\mathbf{B}$, the learning rate
$10^{-3}$ for updating $\mathbf{B}$ and $10^{-4}$ for updating $\mathbf{A}$ in
Algorithm 1. Then, we select
$\lambda=\\{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0\\}$ in (8.1) using 5-fold
cross-validation with the BIC metric. Then, we use $skggm$ package to
implement the graph estimation with selected $\lambda$. The computing time of
model selection of $\lambda$ along with the graph estimation is shown in Step
8–9 of Table 8. Noting that Step 5, 8–9 take much longer time to finish as the
running time depends on the size of the tuning parameter candidates, other
steps are executed only once and could be done within a few minutes as shown
in Table 8.

Task |  | Running Time (s)
---|---|---
Step 1: | standard preprocessing | $202.60$
Step 2: | model Selection: basis family, $k$, $k_{m}$ | $201.64$
Step 3: | project data to selected basis | $16.91$
Step 4: | initialization of $\mathbf{A}$ | $2.374$
Step 5: | model selection: $s,\alpha,\tau_{1},\tau_{2}$ | $23610.33^{*}$
Step 6: | initialization of $\mathbf{B}$ (Algorithm 2) | $251.17$
Step 7: | Algorithm 1 | $277.19$
Step 8: | estimate graph of fMRI | $2777.08^{*}$
Step 9: | estimate graph of EEG | $4289.77^{*}$

Table 8: Running time of the fMRI-EEG experiment. We use Intel Xeon Processor E5-2620 v3 @ 2.40 GHz to run all the steps. The * implies that we parallel the experiment with 24 CPUs and the remaining experiments are executed with single CPU.  Table 9:  The index table of Desikan-Killiany atlas (desikan2006automated) reordered to Yeo 7 network (yeo2011organization). VIS = visual; SM = somatomotor; DA = dorsal attention; VA = ventral attention; L = limbic; FP = fronto-parietal; DMN = default mode network.  Matrix Index |  Desikan-Killiany Atlas Index |  Abbreviation |  Yeo 7 Network |  Name
---|---|---|---|---
1 |  4 |  lCUN |  VIS |  L cuneus
2 |  6 |  lFUS |  VIS |  L fusiform
3 |  10 |  lLOG |  VIS |  L lateral occipital gyrus
4 |  12 |  lLING |  VIS |  L lingual
5 |  20 |  lperiCAL |  VIS |  L pericalcarine
6 |  38 |  rCUN |  VIS |  R cuneus
7 |  40 |  rFUS |  VIS |  R fusiform
8 |  44 |  rLOG |  VIS |  R lateral occipital gyrus
9 |  46 |  rLING |  VIS |  R lingual
10 |  54 |  rperiCAL |  VIS |  R pericalcarine
11 |  16 |  lparaC |  SM |  L paracentral
12 |  21 |  lpostC |  SM |  L postcentral
13 |  23 |  lpreC |  SM |  L precentral
14 |  29 |  lSTG |  SM |  L superior temporal gyrus
15 |  33 |  lTT |  SM |  L transverse temporal
16 |  50 |  rparaC |  SM |  R paracentral
17 |  55 |  rpostC |  SM |  R postcentral
18 |  56 |  rPCC |  SM |  R posterior cingulate cortex
19 |  57 |  rpreC |  SM |  R precentral
20 |  63 |  rSTG |  SM |  R superior temporal gyrus
21 |  67 |  rTT |  SM |  R transverse temporal
22 |  28 |  lSPL |  DA |  L superior parietal lobule
23 |  62 |  rSPL |  DA |  R superior parietal lobule
24 |  2 |  lcACC |  VA |  L caudal anterior cingulate
25 |  17 |  lpOPER |  VA |  L pars opercularis
26 |  30 |  lSMAR |  VA |  L supramarginal gyrus
27 |  34 |  lINS |  VA |  L insula
28 |  36 |  rcACC |  VA |  R caudal anterior cingulate
29 |  64 |  rSMAR |  VA |  R supramarginal gyrus
30 |  68 |  rINS |  VA |  R insula
31 |  5 |  lENT |  L |  L entorhinal
32 |  8 |  lITG |  L |  L inferior temporal gyrus
33 |  11 |  lLOF |  L |  L lateral orbitofrontal
34 |  13 |  lMOF |  L |  L medial orbitofrontal
35 |  31 |  lFP |  L |  L frontal pole
36 |  32 |  lTP |  L |  L temporal pole
37 |  39 |  rENT |  L |  R entorhinal
38 |  42 |  rITG |  L |  R inferior temporal gyrus
39 |  45 |  rLOF |  L |  R lateral orbitofrontal
40 |  47 |  rMOF |  L |  R medial orbitofrontal
41 |  65 |  rFP |  L |  R frontal pole
42 |  66 |  rTP |  L |  R temporal pole
43 |  26 |  lrMFG |  FP |  L rostral middle frontal gyrus
44 |  37 |  rcMFG |  FP |  R caudal middle frontal gyrus
45 |  51 |  rpOPER |  FP |  R pars opercularis
46 |  53 |  rpTRI |  FP |  R pars triangularis
47 |  60 |  rrMFG |  FP |  R rostral middle frontal gyrus
48 |  1 |  lBSTS |  DMN |  L bank of the superior temporal sulcus
49 |  3 |  lcMFG |  DMN |  L caudal middle frontal gyrus
50 |  7 |  lIPL |  DMN |  L inferior parietal lobule
51 |  9 |  liCC |  DMN |  L isthmus cingulate cortex
52 |  14 |  lMTG |  DMN |  L middle temporal gyrus
53 |  18 |  lpORB |  DMN |  L pars orbitalis
54 |  19 |  lpTRI |  DMN |  L pars triangularis
55 |  22 |  lPCC |  DMN |  L posterior cingulate cortex
56 |  24 |  lPCUN |  DMN |  L precuneus
57 |  25 |  lrACC |  DMN |  L rostral anterior cingulate cortex
58 |  27 |  lSFG |  DMN |  L superior frontal gyrus
59 |  35 |  rBSTS |  DMN |  R bank of the superior temporal sulcus
60 |  41 |  rIPL |  DMN |  R inferior parietal lobule
61 |  43 |  riCC |  DMN |  R isthmus cingulate cortex
62 |  48 |  rMTG |  DMN |  R middle temporal gyrus
63 |  52 |  rpORB |  DMN |  R pars orbitalis
64 |  58 |  rPCUN |  DMN |  R precuneus
65 |  59 |  rrACC |  DMN |  R rostral anterior cingulate cortex
66 |  61 |  rSFG |  DMN |  R superior frontal gyrus
67 |  15 |  lPARH |  DMN |  L parahippocampal
68 |  49 |  rPARH |  DMN |  R parahippocampal

## References