(C.84) $\displaystyle=\left\langle\left(\int_{\mathscr{X}}K_{x}\otimes
K_{x}\,\mathrm{d}\left(\rho^{(k+1)}_{\mathscr{X}}-\rho^{(k)}_{\mathscr{X}}\right)\right)f,\left(\int_{\mathscr{X}}K_{x}\otimes
K_{x}\,\mathrm{d}\left(\rho^{(k+1)}_{\mathscr{X}}-\rho^{(k)}_{\mathscr{X}}\right)\right)f\right\rangle_{K}$
(C.85) $\displaystyle=\left\langle\left(\int_{\mathscr{X}}K_{x}\otimes
K_{x}\,\mathrm{d}\left(\rho^{(k+1)}_{\mathscr{X}}-\rho^{(k)}_{\mathscr{X}}\right)\right)\left(\int_{\mathscr{X}}K_{x}\otimes
K_{x}\,\mathrm{d}\left(\rho^{(k+1)}_{\mathscr{X}}-\rho^{(k)}_{\mathscr{X}}\right)\right)f,f\right\rangle_{K}$
(C.86)
$\displaystyle=\left\langle\int_{\mathscr{X}}K_{y}\left(\int_{\mathscr{X}}f(x)K(y,x)\,\mathrm{d}\Delta_{k}(x)\right)\,\mathrm{d}\Delta_{k}(y),f\right\rangle_{K}$
(C.87) $\displaystyle=\int_{\mathscr{X}}\left\langle
K_{y}\left(\int_{\mathscr{X}}f(x)K(y,x)\,\mathrm{d}\Delta_{k}(x)\right),f\right\rangle_{K}\,\mathrm{d}\Delta_{k}(y)$
(C.88)
$\displaystyle=\int_{\mathscr{X}}\left(\int_{\mathscr{X}}f(x)K(y,x)\,\mathrm{d}\Delta_{k}(x)\right)\left\langle
K_{y},f\right\rangle_{K}\,\mathrm{d}\Delta_{k}(y)$ (C.89)
$\displaystyle=\int_{\mathscr{X}}f(y)\left(\int_{\mathscr{X}}f(x)K(y,x)\,\mathrm{d}\Delta_{k}(x)\right)\,\mathrm{d}\Delta_{k}(y),~{}\forall~{}k\in\mathbb{N},$
(C.90)

where
$\Delta_{k}=\rho^{(k+1)}_{\mathscr{X}}-\rho^{(k)}_{\mathscr{X}}\in\mathcal{M}(\mathscr{X})$.
Since $C^{s}(\mathscr{X})\subseteq C(\mathscr{X})$ and
$(C(\mathscr{X}))^{*}=\mathcal{M}(\mathscr{X})$, then
$\mathcal{M}(\mathscr{X})\subseteq(C^{s}(\mathscr{X}))^{*}$, from which we
have $\Delta_{k}\in(C^{s}(\mathscr{X}))^{*}$. Denote

$\displaystyle
g_{k}(\cdot)=f(\cdot)\left(\int_{\mathscr{X}}f(x)K(\cdot,x)\,\mathrm{d}\Delta_{k}(x)\right),~{}\forall~{}k\in\mathbb{N}.$
(C.91)

Noting that $g\in\mathscr{H}_{K}\subseteq C^{s}(\mathscr{X})$, it follows from
$\Delta_{k}\in(C^{s}(\mathscr{X}))^{*}$ and the definition of
$(C^{s}(\mathscr{X}))^{*}$ that

$\displaystyle\int_{\mathscr{X}}g_{k}(y)\,\mathrm{d}\Delta_{k}(y)\leq\left\|\Delta_{k}\right\|_{(C^{s}(\mathscr{X}))^{*}}\left\|g_{k}\right\|_{C^{s}(\mathscr{X})}=\left\|\Delta_{k}\right\|_{(C^{s}(\mathscr{X}))^{*}}\left(\left\|g_{k}\right\|_{\infty}+\left|g_{k}\right|_{C^{s}(\mathscr{X})}\right).$
(C.92)

We now estimate $\|g_{k}\|_{\infty}$ and $|g_{k}|_{C^{s}(\mathscr{X})}$,
respectively.

It follows from Assumption 3 that $K\in
C^{s}(\mathscr{X}\times\mathscr{X})\subseteq C(\mathscr{X}\times\mathscr{X})$,
which shows that there exists a constant $\kappa_{1}<\infty$, such that
$\kappa_{1}=\sup_{x\in\mathscr{X}}\sqrt{K(x,x)}$. Noting that
$\|g\|_{\infty}\leq\kappa_{1}\|g\|_{K}$ and
$\|g\|_{C^{s}(\mathscr{X})}\leq(\kappa_{1}+\tau_{s})\|g\|_{K}$,
$\forall~{}g\in\mathscr{H}_{K}$, by Lemma D.4 and the reproducing property of
RKHS, we get

$\displaystyle~{}~{}~{}\left\|fK_{y}\right\|_{C^{s}(\mathscr{X})}$
$\displaystyle=\left\|fK_{y}\right\|_{\infty}+\left|fK_{y}\right|_{C^{s}(\mathscr{X})}$
$\displaystyle\leq\left\|f\right\|_{\infty}\left\|K_{y}\right\|_{\infty}+\left|f\right|_{C^{s}(\mathscr{X})}\left\|K_{y}\right\|_{\infty}+\left\|f\right\|_{\infty}\left|K_{y}\right|_{C^{s}(\mathscr{X})}$
$\displaystyle\leq\kappa_{1}\left\|f\right\|_{\infty}\left\|K_{y}\right\|_{K}+\kappa_{1}\left|f\right|_{C^{s}(\mathscr{X})}\left\|K_{y}\right\|_{K}+\left\|f\right\|_{C^{s}(\mathscr{X})}\left\|K_{y}\right\|_{C^{s}(\mathscr{X})}$
$\displaystyle\leq\kappa_{1}\sup_{y\in\mathscr{X}}\sqrt{K(y,y)}\left\|f\right\|_{\infty}+\kappa_{1}\sup_{y\in\mathscr{X}}\sqrt{K(y,y)}\left|f\right|_{C^{s}(\mathscr{X})}+\left(\kappa_{1}+\tau_{s}\right)\left\|f\right\|_{C^{s}(\mathscr{X})}\left\|K_{y}\right\|_{K}$
$\displaystyle\leq\kappa^{2}_{1}\left(\left\|f\right\|_{\infty}+\left|f\right|_{C^{s}(\mathscr{X})}\right)+\left(\kappa_{1}+\tau_{s}\right)\sup_{y\in\mathscr{X}}\sqrt{K(y,y)}\left\|f\right\|_{C^{s}(\mathscr{X})}$
$\displaystyle=\left(2\kappa^{2}_{1}+\kappa_{1}\tau_{s}\right)\left\|f\right\|_{C^{s}(\mathscr{X})},~{}\forall~{}y\in\mathscr{X},$

which shows that

$\displaystyle~{}~{}~{}\left|\int_{\mathscr{X}}f(x)K(y,x)\,\mathrm{d}\Delta_{k}(x)\right|$
(C.93)
$\displaystyle\leq\left\|\Delta_{k}\right\|_{(C^{s}(\mathscr{X}))^{*}}\left\|fK_{y}\right\|_{C^{s}(\mathscr{X})}$
(C.94)
$\displaystyle\leq\left(2\kappa^{2}_{1}+\kappa_{1}\tau_{s}\right)\left\|\Delta_{k}\right\|_{(C^{s}(\mathscr{X}))^{*}}\left\|f\right\|_{C^{s}(\mathscr{X})},~{}\forall~{}y\in\mathscr{X},~{}\forall~{}k\in\mathbb{N}.$
(C.95)

Thus, it follows from (C.91) and (C.93) that

$\displaystyle\left\|g_{k}\right\|_{\infty}\leq\left\|f\right\|_{\infty}\sup_{y\in\mathscr{X}}\left|\int_{\mathscr{X}}f(x)K(y,x)\,\mathrm{d}\Delta_{k}(x)\right|\leq\left(2\kappa^{2}_{1}+\kappa_{1}\tau_{s}\right)\left\|f\right\|^{2}_{C^{s}(\mathscr{X})}\left\|\Delta_{k}\right\|_{(C^{s}(\mathscr{X}))^{*}}.$
(C.96)

By Lemma D.4 and (C.93), we obtain

$\displaystyle~{}~{}~{}\left|g_{k}\right|_{C^{s}(\mathscr{X})}$ (C.97)
$\displaystyle\leq\left|f\right|_{C^{s}(\mathscr{X})}\left\|\int_{\mathscr{X}}f(x)K_{x}\,\mathrm{d}\Delta_{k}(x)\right\|_{\infty}+\left\|f\right\|_{\infty}\left|\int_{\mathscr{X}}f(x)K_{x}\,\mathrm{d}\Delta_{k}(x)\right|_{C^{s}(\mathscr{X})}$
(C.98)
$\displaystyle=\left|f\right|_{C^{s}(\mathscr{X})}\sup_{y\in\mathscr{X}}\left|\int_{\mathscr{X}}f(x)K(y,x)\,\mathrm{d}\Delta_{k}(x)\right|+\left\|f\right\|_{\infty}\left|\int_{\mathscr{X}}f(x)K_{x}\,\mathrm{d}\Delta_{k}(x)\right|_{C^{s}(\mathscr{X})}$
(C.99)
$\displaystyle\leq\left(2\kappa^{2}_{1}+\kappa_{1}\tau_{s}\right)\left\|\Delta_{k}\right\|_{(C^{s}(\mathscr{X}))^{*}}\left\|f\right\|^{2}_{C^{s}(\mathscr{X})}+\left\|f\right\|_{C^{s}(\mathscr{X})}\left|\int_{\mathscr{X}}f(x)K_{x}\,\mathrm{d}\Delta_{k}(x)\right|_{C^{s}(\mathscr{X})}.$
(C.100)

By the definition of $(C^{s}(\mathscr{X}))^{*}$, we get

$\displaystyle~{}~{}~{}\left|\int_{\mathscr{X}}f(x)K_{x}\,\mathrm{d}\Delta_{k}(x)\right|_{C^{s}(\mathscr{X})}$
(C.101) $\displaystyle=\sup_{z_{1}\neq
z_{2}\in\mathscr{X}}\left|\int_{\mathscr{X}}f(x)\frac{K(z_{1},x)-K(z_{2},x)}{\left\|z_{1}-z_{2}\right\|^{s}}\,\mathrm{d}\Delta_{k}(x)\right|$
(C.102)
$\displaystyle\leq\left\|\Delta_{k}\right\|_{(C^{s}(\mathscr{X}))^{*}}\sup_{z_{1}\neq
z_{2}\in\mathscr{X}}\left\|f\frac{K_{z_{1}}-K_{z_{2}}}{\left\|z_{1}-z_{2}\right\|^{s}}\right\|_{C^{s}(\mathscr{X})},~{}\forall~{}k\in\mathbb{N}.$
(C.103)

By the definition of $\|\cdot\|_{C^{s}(\mathscr{X})}$ and Assumption 3, we
have

$\displaystyle~{}~{}~{}\left\|f\frac{K_{z_{1}}-K_{z_{2}}}{\left\|z_{1}-z_{2}\right\|^{s}}\right\|_{C^{s}(\mathscr{X})}$
(C.104)
$\displaystyle=\left\|f\frac{K_{z_{1}}-K_{z_{2}}}{\left\|z_{1}-z_{2}\right\|^{s}}\right\|_{\infty}+\left|f\frac{K_{z_{1}}-K_{z_{2}}}{\left\|z_{1}-z_{2}\right\|^{s}}\right|_{C^{s}(\mathscr{X})}$
(C.105)
$\displaystyle\leq\left\|f\right\|_{\infty}\left\|\frac{K_{z_{1}}-K_{z_{2}}}{\left\|z_{1}-z_{2}\right\|^{s}}\right\|_{\infty}+\left|f\frac{K_{z_{1}}-K_{z_{2}}}{\left\|z_{1}-z_{2}\right\|^{s}}\right|_{C^{s}(\mathscr{X})}$
(C.106)
$\displaystyle\leq\left\|f\right\|_{C^{s}(\mathscr{X})}\sup_{(z_{1},x)\neq(z_{2},x)\in\mathscr{X}\times\mathscr{X}}\frac{\left|K(z_{1},x)-K(z_{2},x)\right|}{\left\|z_{1}-z_{2}\right\|^{s}}+\left|f\frac{K_{z_{1}}-K_{z_{2}}}{\left\|z_{1}-z_{2}\right\|^{s}}\right|_{C^{s}(\mathscr{X})}$
(C.107)
$\displaystyle\leq\left|K\right|_{C^{s}(\mathscr{X}\times\mathscr{X})}\left\|f\right\|_{C^{s}(\mathscr{X})}+\left|f\frac{K_{z_{1}}-K_{z_{2}}}{\left\|z_{1}-z_{2}\right\|^{s}}\right|_{C^{s}(\mathscr{X})},~{}\forall~{}z_{1}\neq
z_{2}\in\mathscr{X}.$ (C.108)

It follows from Lemma D.4 and Assumption 3 that

$\displaystyle~{}~{}~{}\left|f\frac{K_{z_{1}}-K_{z_{2}}}{\left\|z_{1}-z_{2}\right\|^{s}}\right|_{C^{s}(\mathscr{X})}$
$\displaystyle\leq\left|f\right|_{C^{s}(\mathscr{X})}\left\|\frac{K_{z_{1}}-K_{z_{2}}}{\left\|z_{1}-z_{2}\right\|^{s}}\right\|_{\infty}+\left\|f\right\|_{\infty}\left|\frac{K_{z_{1}}-K_{z_{2}}}{\left\|z_{1}-z_{2}\right\|^{s}}\right|_{C^{s}(\mathscr{X})}$
$\displaystyle\leq\left|f\right|_{C^{s}(\mathscr{X})}\sup_{(z_{1},x)\neq(z_{2},x)\in\mathscr{X}\times\mathscr{X}}\frac{\left|K(z_{1},x)-K(z_{2},x)\right|}{\left\|z_{1}-z_{2}\right\|^{s}}+\left\|f\right\|_{C^{s}(\mathscr{X})}\left|\frac{K_{z_{1}}-K_{z_{2}}}{\left\|z_{1}-z_{2}\right\|^{s}}\right|_{C^{s}(\mathscr{X})}$
$\displaystyle\leq\left\|f\right\|_{C^{s}(\mathscr{X})}\left(\left|K\right|_{C^{s}(\mathscr{X}\times\mathscr{X})}+\sup_{w_{1}\neq
w_{2}\in\mathscr{X}}\frac{\left|K(z_{1},w_{1})-K(z_{2},w_{1})-K(z_{1},w_{2})+K(z_{2},w_{2})\right|}{\left\|z_{1}-z_{2}\right\|^{s}\left\|w_{1}-w_{2}\right\|^{s}}\right)$
$\displaystyle\leq\left\|f\right\|_{C^{s}(\mathscr{X})}\left(\left|K\right|_{C^{s}(\mathscr{X}\times\mathscr{X})}+\tau_{s}\right),~{}\forall~{}z_{1}\neq
z_{2}\in\mathscr{X},$

Thus, by the above and (C.104), we get

$\displaystyle\left\|f\frac{K_{z_{1}}-K_{z_{2}}}{\left\|z_{1}-z_{2}\right\|^{s}}\right\|_{C^{s}(\mathscr{X})}\leq\left\|f\right\|_{C^{s}(\mathscr{X})}\left(2\left|K\right|_{C^{s}(\mathscr{X}\times\mathscr{X})}+\tau_{s}\right).$
(C.109)

Putting (C.104)-(C.109) into (C.101) leads to

$\displaystyle\left|\int_{\mathscr{X}}f(x)K_{x}\,\mathrm{d}\Delta_{k}(x)\right|_{C^{s}(\mathscr{X})}\leq\left\|\Delta_{k}\right\|_{(C^{s}(\mathscr{X}))^{*}}\left\|f\right\|_{C^{s}(\mathscr{X})}\left(2\left|K\right|_{C^{s}(\mathscr{X}\times\mathscr{X})}+\tau_{s}\right),~{}\forall~{}k\in\mathbb{N},$

which together with (C.97) gives

$\displaystyle\left|g_{k}\right|_{C^{s}(\mathscr{X})}\leq\left\|\Delta_{k}\right\|_{(C^{s}(\mathscr{X}))^{*}}\left\|f\right\|^{2}_{C^{s}(\mathscr{X})}\left(2\kappa^{2}_{1}+\kappa_{1}\tau_{s}+2\left|K\right|_{C^{s}(\mathscr{X}\times\mathscr{X})}+\tau_{s}\right),~{}\forall~{}k\in\mathbb{N}.$
(C.110)

Putting (C.91)-(C.92), (C.96) and (C.110) into (C.82) shows

$\displaystyle\left\|\left(\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right]-\mathbb{E}\left[K_{x_{k}}\otimes
K_{x_{k}}|{\mathcal{F}}_{k-1}\right]\right)f\right\|^{2}_{K}$
$\displaystyle\leq\left\|\Delta_{k}\right\|^{2}_{(C^{s}(\mathscr{X}))^{*}}\left\|f\right\|^{2}_{C^{s}(\mathscr{X})}\left(4\kappa^{2}_{1}+2\kappa_{1}\tau_{s}+2\left|K\right|_{C^{s}(\mathscr{X}\times\mathscr{X})}+\tau_{s}\right)$
$\displaystyle\leq\left\|\Delta_{k}\right\|^{2}_{(C^{s}(\mathscr{X}))^{*}}\left\|f\right\|^{2}_{K}\left(\kappa_{1}+\tau_{s}\right)^{2}\left(4\kappa^{2}_{1}+2\kappa_{1}\tau_{s}+2\left|K\right|_{C^{s}(\mathscr{X}\times\mathscr{X})}+\tau_{s}\right),~{}\forall~{}f\in\mathscr{H}_{K},~{}\forall~{}k\in\mathbb{N},$

from which we have

$\displaystyle~{}~{}~{}\left\|\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right]-\mathbb{E}\left[K_{x_{k}}\otimes
K_{x_{k}}|{\mathcal{F}}_{k-1}\right]\right\|_{\mathscr{L}(\mathscr{H}_{K})}$
$\displaystyle\leq\left\|\Delta_{k}\right\|_{(C^{s}(\mathscr{X}))^{*}}\left(\kappa_{1}+\tau_{s}\right)\sqrt{4\kappa^{2}_{1}+2\kappa_{1}\tau_{s}+2\left|K\right|_{C^{s}(\mathscr{X}\times\mathscr{X})}+\tau_{s}}$
$\displaystyle=\left\|\rho^{(k+1)}_{\mathscr{X}}-\rho^{(k)}_{\mathscr{X}}\right\|_{(C^{s}(\mathscr{X}))^{*}}\left(\kappa_{1}+\tau_{s}\right)\sqrt{4\kappa^{2}_{1}+2\kappa_{1}\tau_{s}+2\left|K\right|_{C^{s}(\mathscr{X}\times\mathscr{X})}+\tau_{s}},~{}\forall~{}k\in\mathbb{N},$

by the above and (35), we know that

$\displaystyle\left\|\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right]-\mathbb{E}\left[K_{x_{k}}\otimes
K_{x_{k}}|{\mathcal{F}}_{k-1}\right]\right\|_{\mathscr{L}(\mathscr{H}_{K})}=O\left(a_{k}\lambda^{2}_{k}\right).$
(C.111)

Noting that $\lim_{x\to 0}\frac{1-(1-x)^{a}}{x}=a$,
$\forall~{}a\in\mathbb{R}$, by Condition 1, we obtain

$\displaystyle\lim_{k\to\infty}\frac{\lambda_{k}-\lambda_{k+1}}{a_{k}\lambda^{2}_{k}}=\lim_{k\to\infty}\left(\frac{(k+1)^{\tau_{1}+\tau_{2}}}{k+2}\times\frac{1-\left(1-\frac{1}{k+2}\right)^{\tau_{2}}}{\frac{1}{k+2}}\right)=0.$
(C.112)

It follows from Assumption 3 that $K\in C(\mathscr{X}\times\mathscr{X})$,
which shows that Assumption 1 holds. Hence, by Lemma D.5, (C.111)-(C.112), we
get

$\displaystyle\left\|f_{\lambda}(k+1)-f_{\lambda}(k)\right\|_{K}=O\left(a_{k}\lambda_{k}\left\|f_{\lambda}(k)-f^{\star}\right\|_{K}\right).$
(C.113)

It follows from the definition of the random Tikhonov regularization path
$f_{\lambda}(k)$ of $f^{\star}$ that
$\|f_{\lambda}(k)\|_{K}\leq\|f^{\star}\|_{K}$, then we have
$\|f_{\lambda}(k)-f^{\star}\|_{K}\leq 2\|f^{\star}\|_{K}<\infty$, which
together with (C.113) and Condition 1 gives

$\displaystyle\left\|f_{\lambda}(k+1)-f_{\lambda}(k)\right\|_{L^{2}\left(\Omega;\mathscr{H}_{K}\right)}=o\left(\lambda_{k}\right).$
(C.114)

It follows from Assumptions 2-3, Condition 1, (34) and Proposition 3 that the
online data streams $\\{(x_{k},y_{k}),k\in\mathbb{N}\\}$ satisfy the RKHS
persistence of excitation condition, then by (C.114) and Lemma C.2, we get

$\displaystyle\lim_{k\to\infty}\left\|f_{\lambda}(k)-f^{\star}\right\|_{K}=\lim_{k\to\infty}\left\|f_{\lambda}(k)-f^{\star}\right\|_{L^{2}\left(\Omega;\mathscr{H}_{K}\right)}=0.$
(C.115)

Combining (C.113) with (C.115) leads to

$\displaystyle\left\|f_{\lambda}(k+1)-f_{\lambda}(k)\right\|_{L^{2}\left(\Omega;\mathscr{H}_{K}\right)}=\left\|f_{\lambda}(k+1)-f_{\lambda}(k)\right\|_{K}=o\left(a_{k}\lambda_{k}\right).$
(C.116)

Noting that the online data streams $\\{(x_{k},y_{k}),k\in\mathbb{N}\\}$
satisfy the RKHS persistence of excitation condition, by (C.116) and Theorem
1, we have

$\displaystyle\lim_{k\to\infty}\left\|f_{k}-f^{\star}\right\|_{L^{2}\left(\Omega;\mathscr{H}_{K}\right)}=0.$

∎

##  Appendix D Key lemmas

######  Lemma D.1.

If the sequences $\\{a_{k},k\in\mathbb{N}\\}$ and
$\\{\lambda_{k},k\in\mathbb{N}\\}$ satisfy

$a_{k}=\frac{1}{(k+1)^{\tau_{1}}},\quad\lambda_{k}=\frac{1}{(k+1)^{\tau_{2}}},~{}\forall~{}k\in\mathbb{N},$

where
$0.1<\tau_{2}<0.5<\tau_{1}<1,~{}\tau_{1}+\tau_{2}<1,~{}3\tau_{2}<\tau_{1}$,
then

  * (A).

$\displaystyle\sum_{i=1}^{k}a^{2}_{i}\prod_{j=i+1}^{k}(1-a_{j}\lambda_{j})=O\left((k+1)^{\tau_{2}-\tau_{1}}\mbox{ln}(k+1)\right)$.

  * (B).

$\displaystyle\sum_{i=1}^{k}a^{2}_{i}\prod_{j=i+1}^{k}(1-a_{j}\lambda_{j})\sqrt{k-i+1}=O\left((k+1)^{\frac{3\tau_{2}-\tau_{1}}{2}}\mbox{ln}^{\frac{3}{2}}(k+1)\right)$.

###### Proof.

Noting that $1-x\leq e^{-x}$, $\forall~{}x\geq 0$, then we have

$\displaystyle\prod_{j=i+1}^{k}\left(1-\frac{1}{(j+1)^{\tau_{1}+\tau_{2}}}\right)\leq\exp\left(-\sum_{j=i+1}^{k}\frac{1}{(j+1)^{\tau_{1}+\tau_{2}}}\right).$
(D.1)

By directly computing, we get

$\displaystyle\sum_{j=i+1}^{k}\frac{1}{(k+1)^{\tau_{1}+\tau_{2}}}\geq\int_{i+1}^{k}\frac{1}{(k+1)^{\tau_{1}+\tau_{2}}}\,\mathrm{d}x=\frac{1}{1-\tau_{1}-\tau_{2}}\left((k+1)^{1-\tau_{1}-\tau_{2}}-(i+2)^{1-\tau_{1}-\tau_{2}}\right),$

which together with (D.1) leads to

$\displaystyle\prod_{j=i+1}^{k}\left(1-\frac{1}{(j+1)^{\tau_{1}+\tau_{2}}}\right)\leq\exp\left(-\frac{1}{1-\tau_{1}-\tau_{2}}\left((k+1)^{1-\tau_{1}-\tau_{2}}-(i+2)^{1-\tau_{1}-\tau_{2}}\right)\right).$
(D.2)

Denote

$\displaystyle\epsilon_{k}=\left\lceil
2(k+1)^{\tau_{1}+\tau_{2}}\mbox{ln}(k+1)\right\rceil,~{}\forall~{}k\in\mathbb{N}.$
(D.3)

Noting that $\epsilon_{k}=o(k)$ and $\epsilon^{-1}_{k}=o(1)$, there exists a
positive integer $k_{0}$, such that

$k_{0}\leq\epsilon_{k}\leq 2\epsilon_{k}<k.$

On one hand, when $k_{0}\leq i\leq k-1-\epsilon_{k}$,

$\displaystyle i+2\leq k+1-\epsilon_{k}.$ (D.4)

Noting that $(1-x)^{\alpha}\leq 1-\alpha x,~{}\forall~{}\alpha,x\in[0,1]$,
then we obtain

$\displaystyle\left(\frac{k+1-\epsilon_{k}}{k+1}\right)^{1-\tau_{1}-\tau_{2}}=\left(1-\frac{\epsilon_{k}}{k+1}\right)^{1-\tau_{1}-\tau_{2}}\leq
1-\frac{\epsilon_{k}(1-\tau_{1}-\tau_{2})}{k+1},$

which shows that

$\displaystyle(k+1)^{1-\tau_{1}-\tau_{2}}-(k+1-\epsilon_{k})^{1-\tau_{1}-\tau_{2}}$
$\displaystyle\geq$
$\displaystyle(k+1)^{-\tau_{1}-\tau_{2}}\epsilon_{k}(1-\tau_{1}-\tau_{2})$
$\displaystyle\geq$ $\displaystyle 2(1-\tau_{1}-\tau_{2})\mbox{ln}(k+1).$

Combing the above with (D.4) gives

$\displaystyle~{}~{}~{}\frac{1}{1-\tau_{1}-\tau_{2}}\left((k+1)^{1-\tau_{1}-\tau_{2}}-(i+2)^{1-\tau_{1}-\tau_{2}}\right)$
(D.5)
$\displaystyle\geq\frac{1}{1-\tau_{1}-\tau_{2}}\left((k+1)^{1-\tau_{1}-\tau_{2}}-(k+1-\epsilon_{k})^{1-\tau_{1}-\tau_{2}}\right)$
(D.6) $\displaystyle\geq 2\mbox{ln}(k+1).$ (D.7)

By putting (D.5) into (D.2), we get

$\displaystyle\prod_{j=i+1}^{k}\left(1-\frac{1}{(j+1)^{\tau_{1}+\tau_{2}}}\right)\leq\exp\left(-2\mbox{ln}(k+1)\right)=\frac{1}{(k+1)^{2}},~{}k_{0}\leq
i\leq k-1-\epsilon_{k},$

which shows that

$\displaystyle~{}~{}~{}\sum_{i=1}^{k-1-\epsilon_{k}}a^{2}_{i}\prod_{j=i+1}^{k}(1-a_{j}\lambda_{j})$
(D.8)
$\displaystyle=\sum_{i=1}^{k-1-\epsilon_{k}}\frac{1}{(i+1)^{2\tau_{1}}}\prod_{j=i+1}^{k}\left(1-\frac{1}{(j+1)^{\tau_{1}+\tau_{2}}}\right)$
(D.9)
$\displaystyle=\left(\sum_{i=1}^{k_{0}-1}+\sum_{i=k_{0}}^{k-1-\epsilon_{k}}\right)\frac{1}{(i+1)^{2\tau_{1}}}\prod_{j=i+1}^{k}\left(1-\frac{1}{(j+1)^{\tau_{1}+\tau_{2}}}\right)$
(D.10)
$\displaystyle\leq\sum_{i=1}^{k_{0}-1}\exp\left(-\sum_{j=i+1}^{k}\frac{1}{(j+1)^{\tau_{1}+\tau_{2}}}\right)+\frac{k-\epsilon_{k}}{(k+1)^{2}}$
(D.11)
$\displaystyle=O\left(\frac{1}{(k+1)^{\tau_{1}+\tau_{2}}}\right)+O\left(\frac{1}{k+1}\right)$
(D.12) $\displaystyle=O\left(\frac{1}{(k+1)^{\tau_{1}+\tau_{2}}}\right),$
(D.13)

similarly, we have

$\displaystyle~{}~{}~{}\sum_{i=1}^{k-1-\epsilon_{k}}a^{2}_{i}\prod_{j=i+1}^{k}(1-a_{j}\lambda_{j})\sqrt{k-i+1}$
(D.14)
$\displaystyle=\sum_{i=1}^{k-1-\epsilon_{k}}\frac{1}{(i+1)^{2\tau_{1}}}\prod_{j=i+1}^{k}\left(1-\frac{1}{(j+1)^{\tau_{1}+\tau_{2}}}\right)\sqrt{k-i+1}$
(D.15)
$\displaystyle=\left(\sum_{i=1}^{k_{0}-1}+\sum_{i=k_{0}}^{k-1-\epsilon_{k}}\right)\frac{1}{(i+1)^{2\tau_{1}}}\prod_{j=i+1}^{k}\left(1-\frac{1}{(j+1)^{\tau_{1}+\tau_{2}}}\right)\sqrt{k-i+1}$
(D.16)
$\displaystyle\leq\sum_{i=1}^{k_{0}-1}\exp\left(-\sum_{j=i+1}^{k}\frac{1}{(j+1)^{\tau_{1}+\tau_{2}}}\right)\sqrt{k}+\frac{(k-\epsilon_{k})\sqrt{k}}{(k+1)^{2}}$
(D.17)
$\displaystyle=O\left(\frac{1}{(k+1)^{\tau_{1}+\tau_{2}-0.5}}\right)+O\left(\frac{1}{\sqrt{k+1}}\right)$
(D.18) $\displaystyle=O\left(\frac{1}{(k+1)^{\tau_{1}+\tau_{2}-0.5}}\right).$
(D.19)

On the other hand, when $k-\epsilon_{k}\leq i\leq k$, we have $k\leq
2k-2\epsilon_{k}\leq 2i$, from which we get

$\displaystyle\frac{1}{(i+1)^{2\tau_{1}}}\leq\frac{4^{\tau_{1}}}{(k+2)^{2\tau_{1}}},~{}k-\epsilon_{k}\leq
i\leq k.$ (D.20)

Noting that

$\displaystyle\prod_{j=i+1}^{k}\left(1-\frac{1}{(j+1)^{\tau_{1}+\tau_{2}}}\right)\leq
1,~{}\forall i\in\mathbb{N},$

then by (D.3) and (D.20), we obtain

$\displaystyle~{}~{}~{}\sum_{i=k-\epsilon_{k}}^{k}a^{2}_{i}\prod_{j=i+1}^{k}(1-a_{j}\lambda_{j})$
(D.21)
$\displaystyle=\sum_{i=k-\epsilon_{k}}^{k}\frac{1}{(i+1)^{2\tau_{1}}}\prod_{j=i+1}^{k}\left(1-\frac{1}{(j+1)^{\tau_{1}+\tau_{2}}}\right)$
(D.22)
$\displaystyle\leq\frac{4^{\tau_{1}}(\epsilon_{k}+1)}{(k+2)^{2\tau_{1}}}$
(D.23)
$\displaystyle\leq\frac{4^{\tau_{1}}\left(2(k+1)^{\tau_{1}+\tau_{2}}\mbox{ln}(k+1)+2\right)}{(k+2)^{2\tau_{1}}}$
(D.24) $\displaystyle=O\left((k+1)^{\tau_{2}-\tau_{1}}\mbox{ln}(k+1)\right),$
(D.25)

following the same way, we also get

$\displaystyle~{}~{}~{}\sum_{i=k-\epsilon_{k}}^{k}a^{2}_{i}\prod_{j=i+1}^{k}(1-a_{j}\lambda_{j})\sqrt{k-i+1}$
(D.26)
$\displaystyle=\sum_{i=k-\epsilon_{k}}^{k}\frac{1}{(i+1)^{2\tau_{1}}}\prod_{j=i+1}^{k}\left(1-\frac{1}{(j+1)^{\tau_{1}+\tau_{2}}}\right)\sqrt{k-i+1}$
(D.27)
$\displaystyle\leq\frac{4^{\tau_{1}}(\epsilon_{k}+1)\sup_{k-\epsilon_{k}\leq
i\leq k}\sqrt{k-i+1}}{(k+2)^{2\tau_{1}}}$ (D.28)
$\displaystyle\leq\frac{4^{\tau_{1}}(\epsilon_{k}+1)\sqrt{\epsilon_{k}+1}}{(k+2)^{2\tau_{1}}}$
(D.29)
$\displaystyle\leq\frac{4^{\tau_{1}}\left(2(k+1)^{\tau_{1}+\tau_{2}}\mbox{ln}(k+1)+2\right)^{\frac{3}{2}}}{(k+2)^{2\tau_{1}}}$
(D.30)
$\displaystyle=O\left((k+1)^{\frac{3\tau_{2}-\tau_{1}}{2}}\mbox{ln}^{\frac{3}{2}}(k+1)\right).$
(D.31)

Combining (D.8) with (D.21) leads to

$\displaystyle~{}~{}~{}\sum_{i=1}^{k}a^{2}_{i}\prod_{j=i+1}^{k}(1-a_{j}\lambda_{j})$
$\displaystyle=\left(\sum_{i=1}^{k-1-\epsilon_{k}}+\sum_{i=k-\epsilon_{k}}^{k}\right)a^{2}_{i}\prod_{j=i+1}^{k}(1-a_{j}\lambda_{j})$
$\displaystyle=O\left(\frac{1}{(k+1)^{\tau_{1}+\tau_{2}}}\right)+O\left((k+1)^{\tau_{2}-\tau_{1}}\mbox{ln}(k+1)\right)$
$\displaystyle=O\left((k+1)^{\tau_{2}-\tau_{1}}\mbox{ln}(k+1)\right).$

By (D.14) and (D.26), we conclude that

$\displaystyle~{}~{}~{}\sum_{i=1}^{k}a^{2}_{i}\prod_{j=i+1}^{k}(1-a_{j}\lambda_{j})\sqrt{k-i+1}$
$\displaystyle=\left(\sum_{i=1}^{k-1-\epsilon_{k}}+\sum_{i=k-\epsilon_{k}}^{k}\right)a^{2}_{i}\prod_{j=i+1}^{k}(1-a_{j}\lambda_{j})\sqrt{k-i+1}$
$\displaystyle=O\left(\frac{1}{(k+1)^{\tau_{1}+\tau_{2}-0.5}}\right)+O\left((k+1)^{\frac{3\tau_{2}-\tau_{1}}{2}}\mbox{ln}^{\frac{3}{2}}(k+1)\right)$
$\displaystyle=O\left((k+1)^{\frac{3\tau_{2}-\tau_{1}}{2}}\mbox{ln}^{\frac{3}{2}}(k+1)\right).$

∎

######  Lemma D.2.

If Assumptions 1-2 and Condition 1 hold, and

$\displaystyle\left\|f_{\lambda}(k+1)-f_{\lambda}(k)\right\|_{L^{2}\left(\Omega;\mathscr{H}_{K}\right)}=o\left(\lambda_{k}\right),$
(D.32)

then

$\displaystyle\lim_{k\to\infty}\left\|f_{\lambda}(k,h)-f_{\lambda}(k)\right\|_{L^{2}\left(\Omega;\mathscr{H}_{K}\right)}=0.$
(D.33)

###### Proof.

It follows from the definitions of $f_{\lambda}(k)$ and $f_{\lambda}(k,h)$
that

$\displaystyle f_{\lambda}(k,h)-f_{\lambda}(k)$ (D.34)
$\displaystyle=\left(\sum_{i=k+1}^{k+h}\mathbb{E}\left[H_{i}|{\mathcal{F}}_{k}\right]+\left(\sum_{i=k+1}^{k+h}\lambda_{i}\right)I\right)^{-1}\left(\sum_{i=k+1}^{k+h}\mathbb{E}\left[H_{i}|{\mathcal{F}}_{k}\right]\right)f^{\star}$
(D.35)
$\displaystyle-\left(\sum_{i=k+1}^{k+h}\mathbb{E}\left[H_{i}|{\mathcal{F}}_{k}\right]+\left(\sum_{i=k+1}^{k+h}\lambda_{i}\right)I\right)^{-1}\left(\sum_{i=k+1}^{k+h}\mathbb{E}\left[H_{i}|{\mathcal{F}}_{k}\right]+\left(\sum_{i=k}^{k+h}\lambda_{i}\right)I\right)f_{\lambda}(k)$
(D.36)
$\displaystyle=\left(\sum_{i=k+1}^{k+h}\mathbb{E}\left[H_{i}|{\mathcal{F}}_{k}\right]+\left(\sum_{i=k+1}^{k+h}\lambda_{i}\right)I\right)^{-1}\left(\left(\sum_{i=k+1}^{k+h}\mathbb{E}\left[H_{i}|{\mathcal{F}}_{k}\right]\right)f^{\star}\right.$
(D.37)
$\displaystyle\left.-\left(\sum_{i=k+1}^{k+h}\mathbb{E}\left[H_{i}|{\mathcal{F}}_{k}\right]+\left(\sum_{i=k+1}^{k+h}\lambda_{i}\right)I\right)f_{\lambda}(k)\right)$
(D.38)
$\displaystyle=\left(\sum_{i=k+1}^{k+h}\mathbb{E}\left[H_{i}|{\mathcal{F}}_{k}\right]+\left(\sum_{i=k+1}^{k+h}\lambda_{i}\right)I\right)^{-1}\left(\left(\sum_{i=k+2}^{k+h}\mathbb{E}\left[H_{i}|{\mathcal{F}}_{k}\right]\right)f^{\star}\right.$
(D.39)
$\displaystyle\left.-\left(\sum_{i=k+2}^{k+h}\mathbb{E}\left[H_{i}|{\mathcal{F}}_{k}\right]+\left(\sum_{i=k+2}^{k+h}\lambda_{i}\right)I\right)f_{\lambda}(k)\right)$
(D.40)
$\displaystyle=\left(\sum_{i=k+1}^{k+h}\mathbb{E}\left[H_{i}|{\mathcal{F}}_{k}\right]+\left(\sum_{i=k+1}^{k+h}\lambda_{i}\right)I\right)^{-1}\left(\sum_{i=k+2}^{k+h}\left(\mathbb{E}\left[H_{i}|{\mathcal{F}}_{k}\right]+\lambda_{i}I\right)(f_{\lambda}(i)-f_{\lambda}(k))\right).$
(D.41)

Noting that

$\displaystyle\left\|\left(\sum_{i=k+1}^{k+h}\mathbb{E}\left[H_{i}|{\mathcal{F}}_{k}\right]+\left(\sum_{i=k+1}^{k+h}\lambda_{i}\right)I\right)^{-1}\right\|_{K}\leq\left(\sum_{i=k+1}^{k+h}\lambda_{i}\right)^{-1}~{}\mathrm{a.s.},~{}\forall~{}k\in\mathbb{N},$

then by Assumption 1, Condition 1, Minkowski inequality and (D.34), we get

$\displaystyle\left\|f_{\lambda}(k,h)-f_{\lambda}(k)\right\|_{L^{2}\left(\Omega;\mathscr{H}_{K}\right)}$
$\displaystyle\leq$
$\displaystyle\frac{\kappa+1}{h}(k+h+1)^{\tau_{2}}\sum_{i=k+1}^{k+h}\left\|f_{\lambda}(i)-f_{\lambda}(k)\right\|_{L^{2}\left(\Omega;\mathscr{H}_{K}\right)}$
(D.42) $\displaystyle=$ $\displaystyle
O\left((k+1)^{\tau_{2}}\sum_{i=k+1}^{k+h}\|f_{\lambda}(i+1)-f_{\lambda}(i)\|_{L^{2}\left(\Omega;\mathscr{H}_{K}\right)}\right).$
(D.43)

By the algorithm (11), we have

$\displaystyle
f_{k}-f_{k+1}=a_{k}H_{k}f_{k}-a_{k}y_{k}K_{x_{k}}+a_{k}\lambda_{k}f_{k},~{}\forall~{}k\in\mathbb{N}.$

It follows from Lemma 4 that
$\sup_{k\in\mathbb{N}}\|f_{k}\|_{L^{2}(\Omega;\mathscr{H}_{K})}<\infty$, then
by the statistical learning model (1), Assumptions 1-2, Condition 1 and
Minkowski inequality, we obtain

$\displaystyle\sum_{i=k+1}^{k+h}\|f_{\lambda}(i+1)-f_{\lambda}(i)\|_{L^{2}\left(\Omega;\mathscr{H}_{K}\right)}=O\left((k+1)^{-\tau_{1}}\right).$
(D.44)

It follows from Condition 1 that $\tau_{1}>\tau_{2}$, by putting (D.44) into
(D.42), we have (D.33).

∎

######  Lemma D.3.

If Assumption 1 and Condition 1 hold, and the online data streams
$\\{(x_{k},y_{k}),k\in\mathbb{N}\\}$ satisfy the RKHS persistence of
excitation condition, then

$\lim_{k\to\infty}\|f_{\lambda}(k,h)-f^{\star}\|_{L^{2}\left(\Omega;\mathscr{H}_{K}\right)}=0.$

###### Proof.

By Assumption 1 and Proposition A.3, without loss of generality, we assume the
eigensystem of $\sum_{i=k+1}^{k+h}\mathbb{E}[H_{i}|{\mathcal{F}}_{k}]$ is
$\\{\Lambda_{k}(i),e_{k}(i),i\in\mathbb{N}\\}$. It follows from the spectral
theorem of the compact operator that

$\displaystyle f^{\star}=\sum_{i=0}^{\infty}\langle
f^{\star},e_{k}(i)\rangle_{K}e_{k}(i)~{}\mathrm{a.s.},~{}\forall~{}k\in\mathbb{N}.$
(D.45)

Noting that

$\displaystyle\left(\sum_{i=k+1}^{k+h}\mathbb{E}\left[H_{i}|{\mathcal{F}}_{k}\right]\right)f^{\star}$
$\displaystyle=$
$\displaystyle\left(\sum_{i=k+1}^{k+h}\mathbb{E}\left[H_{i}|{\mathcal{F}}_{k}\right]\right)\left(\sum_{i=0}^{\infty}\langle
f^{\star},e_{k}(i)\rangle_{K}e_{k}(i)\right)$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{\infty}\Lambda_{k}(i)\langle
f^{\star},e_{k}(i)\rangle_{K}e_{k}(i)~{}\mathrm{a.s.},~{}\forall~{}k\in\mathbb{N},$

which together with the definition of $f_{\lambda}(k,h)$ gives

$\displaystyle f_{\lambda}(k,h)$ $\displaystyle=$
$\displaystyle\left(\sum_{i=k+1}^{k+h}\mathbb{E}\left[H_{i}|{\mathcal{F}}_{k}\right]+\left(\sum_{i=k+1}^{k+h}\lambda_{i}\right)I\right)^{-1}\left(\sum_{i=k+1}^{k+h}\mathbb{E}\left[H_{i}|{\mathcal{F}}_{k}\right]\right)f^{\star}$
(D.46) $\displaystyle=$
$\displaystyle\sum_{i=0}^{\infty}\frac{\Lambda_{k}(i)}{\Lambda_{k}(i)+\displaystyle\sum_{i=k+1}^{k+h}\lambda_{i}}\langle
f^{\star},e_{k}(i)\rangle_{K}e_{k}(i)~{}\mathrm{a.s.},~{}\forall~{}k\in\mathbb{N}.$
(D.47)

Since the online data streams $\\{(x_{k},y_{k}),k\in\mathbb{N}\\}$ satisfy the
RKHS persistence of excitation condition, then we have

$\displaystyle\inf_{k\in\mathbb{N}}\Lambda_{k}(i)>0~{}\mathrm{a.s.},~{}\forall~{}i\in\mathbb{N}.$
(D.48)

Noting that

$\displaystyle\sum_{i=0}^{\infty}\left(\frac{\displaystyle\sum_{i=k+1}^{k+h}\lambda_{i}}{\Lambda_{k}(i)+\displaystyle\sum_{i=k+1}^{k+h}\lambda_{i}}\right)^{2}\left|\langle
f^{\star},e_{k}(i)\rangle_{K}\right|^{2}\leq\sum_{i=0}^{\infty}\left|\langle
f^{\star},e_{k}(i)\rangle_{K}\right|^{2}=\|f^{\star}\|^{2}_{K},~{}\forall~{}k\in\mathbb{N},$

by (D.45)-(D.48) and the dominated convergence theorem, we obtain

$\displaystyle\|f_{\lambda}(k,h)-f^{\star}\|^{2}_{L^{2}\left(\Omega;\mathscr{H}_{K}\right)}$
$\displaystyle=$
$\displaystyle\left\|\sum_{i=0}^{\infty}\frac{\displaystyle\sum_{i=k+1}^{k+h}\lambda_{i}}{\Lambda_{k}(i)+\displaystyle\sum_{i=k+1}^{k+h}\lambda_{i}}\langle
f^{\star},e_{k}(i)\rangle_{K}e_{k}(i)\right\|^{2}_{L^{2}\left(\Omega;\mathscr{H}_{K}\right)}$
(D.49) $\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{i=0}^{\infty}\left(\frac{\displaystyle\sum_{i=k+1}^{k+h}\lambda_{i}}{\Lambda_{k}(i)+\displaystyle\sum_{i=k+1}^{k+h}\lambda_{i}}\right)^{2}\left|\langle
f^{\star},e_{k}(i)\rangle_{K}\right|^{2}\right]$ (D.50) $\displaystyle=$
$\displaystyle\sum_{i=0}^{\infty}\mathbb{E}\left[\left(\frac{\displaystyle\sum_{i=k+1}^{k+h}\lambda_{i}}{\Lambda_{k}(i)+\displaystyle\sum_{i=k+1}^{k+h}\lambda_{i}}\right)^{2}\left|\langle
f^{\star},e_{k}(i)\rangle_{K}\right|^{2}\right]$ (D.51) $\displaystyle\leq$
$\displaystyle\sum_{i=0}^{\infty}\mathbb{E}\left[\left(\frac{\displaystyle\sum_{i=k+1}^{k+h}\lambda_{i}}{\displaystyle\inf_{k\in\mathbb{N}}\Lambda_{k}(i)+\displaystyle\sum_{i=k+1}^{k+h}\lambda_{i}}\right)^{2}\left|\langle
f^{\star},e_{k}(i)\rangle_{K}\right|^{2}\right].$ (D.52)

Noting that $\inf_{k\in\mathbb{N}}\Lambda_{k}(i)>0$ a.s. together with
Condition 1 implies that

$\displaystyle\lim_{k\to\infty}\left(\frac{\displaystyle\sum_{i=k+1}^{k+h}\lambda_{i}}{\displaystyle\inf_{k\in\mathbb{N}}\Lambda_{k}(i)+\displaystyle\sum_{i=k+1}^{k+h}\lambda_{i}}\right)^{2}\left|\langle
f^{\star},e_{k}(i)\rangle_{K}\right|^{2}=0~{}\mathrm{a.s.},$ (D.53)

thus, it follows from (D.49)-(D.53) and the dominated convergence theorem that

$\lim_{k\to\infty}\|f_{\lambda}(k,h)-f^{\star}\|_{L^{2}\left(\Omega;\mathscr{H}_{K}\right)}=0.$

∎

######  Lemma D.4.

If $(X,d)$ is a compact metric space, $0\leq s\leq 1$, then

$\displaystyle|fg|_{C^{s}(X)}\leq|f|_{C^{s}(X)}\|g\|_{\infty}+\|f\|_{\infty}|g|_{C^{s}(X)},~{}\forall~{}f,g\in
C^{s}(X).$

###### Proof.

It follows from the definitions of $|\cdot|_{C^{s}(X)}$ and
$\|\cdot\|_{C^{s}(X)}$ that

$\displaystyle~{}~{}~{}|fg|_{C^{s}(X)}$ $\displaystyle=\sup_{x\neq y\in
X}\frac{|f(x)g(x)-f(y)g(y)|}{(d(x,y))^{s}}$ $\displaystyle=\sup_{x\neq y\in
X}\frac{|f(x)g(x)-f(y)g(x)+f(y)g(x)-f(y)g(y)|}{(d(x,y))^{s}}$
$\displaystyle\leq\sup_{x\neq y\in
X}\frac{|f(x)g(x)-f(y)g(x)|+|f(y)g(x)-f(y)g(y)|}{(d(x,y))^{s}}$
$\displaystyle\leq\sup_{x\neq y\in
X}\frac{|f(x)g(x)-f(y)g(x)|}{(d(x,y))^{s}}+\sup_{x\neq y\in
X}\frac{|f(y)g(x)-f(y)g(y)|}{(d(x,y))^{s}}$
$\displaystyle\leq\left(\sup_{x\neq y\in
X}\frac{|f(x)-f(y)|}{(d(x,y))^{s}}\right)\left(\sup_{x\in
X}|g(x)|\right)+\left(\sup_{y\in X}|g(y)|\right)\left(\sup_{x\neq y\in
X}\frac{|g(x)-g(y)|}{(d(x,y))^{s}}\right)$
$\displaystyle=|f|_{C^{s}(X)}\|g\|_{\infty}+\|f\|_{\infty}|g|_{C^{s}(X)},~{}\forall~{}f,g\in
C^{s}(X).$

∎

######  Lemma D.5.

If Assumption 1 holds, and $\\{\lambda_{k},k\in\mathbb{N}\\}$ is a sequence of
positive real numbers, then

$\displaystyle\left\|f_{\lambda}(k+1)-f_{\lambda}(k)\right\|_{K}$
$\displaystyle\leq$
$\displaystyle\left(\frac{\left\|\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right]-\mathbb{E}\left[K_{x_{k}}\otimes
K_{x_{k}}|{\mathcal{F}}_{k-1}\right]\right\|_{\mathscr{L}(\mathscr{H}_{K})}}{\lambda_{k}}\right.$
(D.54)
$\displaystyle\left.+\frac{\lambda_{k}-\lambda_{k+1}}{\lambda_{k}}\right)\|f_{\lambda}(k)-f^{\star}\|_{K}~{}\mathrm{a.s.},~{}\forall~{}k\in\mathbb{N}.$
(D.55)

###### Proof.

It follows from Assumption 1 and the definition of $f_{\lambda}(k)$ that

$\displaystyle\lambda_{k}f_{\lambda}(k)=\mathbb{E}\left[K_{x_{k}}\otimes
K_{x_{k}}|{\mathcal{F}}_{k-1}\right]f^{\star}-\mathbb{E}\left[K_{x_{k}}\otimes
K_{x_{k}}|{\mathcal{F}}_{k-1}\right]f_{\lambda}(k),~{}\forall~{}k\in\mathbb{N},$

from which we have

$\displaystyle~{}~{}~{}\left(\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right]+\lambda_{k+1}I\right)\left(f_{\lambda}(k+1)-f_{\lambda}(k)\right)$
$\displaystyle=\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right]f^{\star}-\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right]f_{\lambda}(k)-\lambda_{k+1}f_{\lambda}(k)$
$\displaystyle=\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right]f^{\star}-\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right]f_{\lambda}(k)-\frac{\lambda_{k+1}\lambda_{k}}{\lambda_{k}}f_{\lambda}(k)$
$\displaystyle=\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right]f^{\star}-\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right]f_{\lambda}(k)$
$\displaystyle~{}~{}~{}~{}-\frac{\lambda_{k+1}}{\lambda_{k}}\left(\mathbb{E}\left[K_{x_{k}}\otimes
K_{x_{k}}|{\mathcal{F}}_{k-1}\right]f^{\star}-\mathbb{E}\left[K_{x_{k}}\otimes
K_{x_{k}}|{\mathcal{F}}_{k-1}\right]f_{\lambda}(k)\right)$
$\displaystyle=\left(\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right]-\frac{\lambda_{k+1}}{\lambda_{k}}\mathbb{E}\left[K_{x_{k}}\otimes
K_{x_{k}}|{\mathcal{F}}_{k-1}\right]\right)(f^{\star}-f_{\lambda}(k))$
$\displaystyle=\frac{1}{\lambda_{k}}\left(\lambda_{k}\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right]-\lambda_{k+1}\mathbb{E}\left[K_{x_{k}}\otimes
K_{x_{k}}|{\mathcal{F}}_{k-1}\right]\right)(f^{\star}-f_{\lambda}(k))$
$\displaystyle=\frac{\lambda_{k}-\lambda_{k+1}}{\lambda_{k}}\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right](f^{\star}-f_{\lambda}(k))$
$\displaystyle~{}~{}~{}~{}+\frac{\lambda_{k+1}}{\lambda_{k}}\left(\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right]-\mathbb{E}\left[K_{x_{k}}\otimes
K_{x_{k}}|{\mathcal{F}}_{k-1}\right]\right)(f^{\star}-f_{\lambda}(k))~{}\mathrm{a.s.},~{}\forall~{}k\in\mathbb{N}.$

By multiplying $(\mathbb{E}[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}]+\lambda_{k+1}I)^{-1}$ on both sides of the
above equality, we get

$\displaystyle f_{\lambda}(k+1)-f_{\lambda}(k)$ (D.56)
$\displaystyle=\frac{\lambda_{k}-\lambda_{k+1}}{\lambda_{k}}\left(\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right]+\lambda_{k+1}I\right)^{-1}\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right](f^{\star}-f_{\lambda}(k))$ (D.57)
$\displaystyle+\frac{\lambda_{k+1}}{\lambda_{k}}\left(\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right]+\lambda_{k+1}I\right)^{-1}$ (D.58)
$\displaystyle\times\left(\mathbb{E}\left[K_{x_{k+1}}\otimes
K_{x_{k+1}}|{\mathcal{F}}_{k}\right]-\mathbb{E}\left[K_{x_{k}}\otimes
K_{x_{k}}|{\mathcal{F}}_{k-1}\right]\right)(f^{\star}-f_{\lambda}(k))~{}\mathrm{a.s.},~{}\forall~{}k\in\mathbb{N}.$
(D.59)

Noting that

$\displaystyle\begin{cases}\displaystyle\left\|\left(\mathbb{E}\left[K_{x_{k}}\otimes
K_{x_{k}}|{\mathcal{F}}_{k-1}\right]+\lambda_{k}I\right)^{-1}\mathbb{E}\left[K_{x_{k}}\otimes
K_{x_{k}}|{\mathcal{F}}_{k-1}\right]\right\|_{\mathscr{L}(\mathscr{H}_{K})}\leq
1~{}\mathrm{a.s.},~{}\forall~{}k\in\mathbb{N},\\\
\displaystyle\left\|\left(\mathbb{E}\left[K_{x_{k}}\otimes
K_{x_{k}}|{\mathcal{F}}_{k-1}\right]+\lambda_{k}I\right)^{-1}\right\|_{\mathscr{L}(\mathscr{H}_{K})}\leq\frac{1}{\lambda_{k}}~{}\mathrm{a.s.},~{}\forall~{}k\in\mathbb{N},\end{cases}$
(D.60)

then by (D.56)-(D.60), we obtain (D.54). ∎

## References

  * [1] V. N. Vapnik, _Adaptive and Learning Systems for Signal Processing, Communications, and Control_. John Wiley & Sons Inc., New York, 1998.
  * [2] E. Parzen, “An approach to time series analysis,” Ann. Math. Statist., vol. 32, no. 4, pp. 951-989, 1961.
  * [3] G. Wahba, _Spline Models for Observational Data_. PA, USA: SIAM, 1990.
  * [4] S. Smale and F. Cucker, “On the mathematical foundations of learning,” Bulletin of the American Mathematical Society, vol. 39, no. 1, pp. 1-49, 2001.
  * [5] S. Smale and F. Cucker, “Best choice for regularization parameters in learning theory: on the bias-variance problem”, Foundations of Computational Mathematics, vol. 2, pp. 413-418, 2002.
  * [6] S. Smale and D. -X. Zhou, “Learning theory estimates via integral operators and their approximations”, Constructive Approximation, vol. 26, no. 2, pp. 153-172, 2007.
  * [7] B. Bousselmi, J. -F. Dupuy and A. Karoui, “Reproducing kernels based schemes for nanparametric regression”, arXiv preprint: 2001. 11213, 2020.
  * [8] I. Steinwart, D. Hush, and C. Scovel, “Learning from dependent observations”, J. Multivariate Anal., vol. 100, no. 1, pp. 175-194, 2009.
  * [9] B. Yu, “Rates of convergence for empirical processes of stationary mixing sequences”, Ann. Probab., vol. 22, no. 1, pp. 94-116, 1994.
  * [10] R. Meir, “Nonparametric time series prediction through adaptive model selection”, Machine Learning, vol. 39, pp. 5-34, 2000.
  * [11] B. Zou, L. Q. Li, and Z. B. Xu, “The generalization performance of ERM algorithm with strongly mixing observations”, J. Mach. Learn. Res., vol. 75, no. 3, pp. 275-295, 2009.
  * [12] M. Mohri and A. Rostamizadeh, “Stability bounds for stationary $\phi$-mixing and $\beta$-mixing processes”, J. Mach. Learn. Res., vol. 11, pp. 789-814, 2010.
  * [13] B. Zou, Y. Y. Tang, Z. B. Xu, L. Q. Li, J. Xu, and Y. Lu, “The generalization performance of regularized regression algorithms based on markov sampling”, IEEE Trans. Cybernetics, vol. 44, no. 9, pp. 1497-1507, 2014.
  * [14] Z. W. Pan and Q. W. Xiao, “Least-square regularized regression with non-iid sampling”, J. Statist. Plann. Inference, vol. 139, no. 10, pp. 3579-3587, 2009.
  * [15] M. J. Zhang and H. W. Sun, “Regression learning with non-identically and non-independently sampling”, International Journal of Wavelets, Multiresolution and Information Processing, vol. 15, no. 1, 2017.
  * [16] A. Sancetta, “Estimation in reproducing kernel Hilbert spaces with dependent data”, IEEE Trans. Information Theory, vol. 67, no. 3, pp. 1782-1795, 2020.
  * [17] I. Ziemann and S. Tu, “Learning with little mixing”, Advances in Neural Information Processing Systems, vol. 35, pp. 4626-4637, 2022.
  * [18] A. Agarwal and J. C. Duchi, “The generalization ability of online algorithms for dependent data”, IEEE Trans. Information Theory, vol. 59, no. 1, pp. 573-587, 2013.
  * [19] J. Xu, Y. Y. Tang, B. Zou, Z. Xu, L. Li, and Y. Lu, “The generalization ability of online SVM classification based on Markov sampling”, IEEE Trans. Neural Networks and Learning Systems, vol. 26, no. 3, pp. 628-639, 2014.
  * [20] V. Kuznetsov and M. Mohri, “Time series prediction and online learning”, Proc. 29th Annual Conference on Learning Theory, PMLR, vol. 49, pp. 1190-1213, 2016.
  * [21] A. Godichon-Baggioni, N. Werge, and O. Wintenberger, “Learning from time-dependent streaming data with online stochastic algorithms”, arXiv preprint: 2205.12549, 2022.
  * [22] S. Smale and Y. Yao, “Online learning algorithms”, Found. Comput. Math., vol. 6, no. 2, pp. 145-170, 2006.
  * [23] Y. Yao, “A dynamic theory of learning”, Ph.D dissertation, Dept. Math., Univ. Calfornia, Berkeley, CA, USA, 2006.
  * [24] Y. Ying and M. Pontil, “Online gradient descent learning algorithms”, Found. Comput. Math., vol. 5, no. 5, pp. 561-596, 2008.
  * [25] P. Tarr$\grave{\text{e}}$s and Y. Yao, “Online learning as stochastic approximation of regularization paths”, IEEE Trans. Information Theory, vol. 60, no. 99, pp. 5716-5735, 2014.
  * [26] A. Dieuleveut and F. Bach, “Nonparametric stochastic approximation with large step-sezes”, The Annals of Statistics, vol. 44, no. 4, pp. 1363-1399, 2016.
  * [27] Y. Ying and D. -X. Zhou, “Unregularized online learning algorithms with general loss functions”, Applied and Computational Harmonic Analysis, vol. 42, no. 2, pp. 224-244, 2017.
  * [28] Z. C. Guo and L. Shi, “Fast and strong convergence of online learning algorithms”, Advances in Computational Mathematics, vol. 45, pp. 2745-2770, 2019.
  * [29] J. Lin and V. Cevher, “Optimal convergence for distributed learning with stochastic gradient methods and spectral algorithms”, Journal of Machine Learning Research, vol. 21, no. 147, pp. 1-63, 2020.
  * [30] X. Guo, Z. C. Guo, and L. Shi, “Capacity dependent analysis for functional online learning algorithms”, Applied and Computational Harmonic Analysis, vol. 67, 101567, 2023.
  * [31] Z. C. Guo, A. Christmann, and L. Shi, “Optimality of Robust Online Learning”, arXiv preprint: 2304.10060, 2023.
  * [32] S. Smale and D. -X. Zhou, “Online learning with Markov sampling”, Analysis and Applications, vol. 7, no. 1, pp. 87-113, 2009.
  * [33] T. Hu and D. -X. Zhou, “Online learning with samples drawn from non-identical distributions”, Journal of Machine Learning Research, vol. 10, no. 12, pp. 2873-2898, 2009.
  * [34] T. Li and X. Zhang, “Random inverse problems over graphs: decentralized online learning”, arXiv preprint: 2303.11789, 2023.
  * [35] M. Green and J. B. Moore, “Persistency of excitation in linear systems”, Syst. Control Lett., vol 7. no. 5, pp. 351-360, 1986.
  * [36] L. Guo, “Estimating time-varying parameters by Kalman filter based algorithm: Stability and convergence”, IEEE Trans. Automatic Control, vol. 35, no. 2, pp. 141-147, 1990.
  * [37] J. F. Zhang, L. Guo, and H. F. Chen, “$L_{p}$-stability of estimation errors of Kalman filter for tracking time-varying parameters”, Int. J. Adaptive Control and Signal Processing, vol. 5, pp. 155-174, 1991.
  * [38] L. Guo, “Stability of recursive stochastic tracking algorithms”, SIAM J. Control Optim., vol. 32, no. 5, pp. 1195-1225, 1994.
  * [39] L. Guo and L. Ljung, “Performance analysis of general tracking algorithms”, IEEE Trans. Automatic Control, vol. 40, no. 8, pp. 1388-1402, 1995.
  * [40] L. Guo and L. Ljung, and G. J. Wang, “Necessary and sufficient conditions for stability of LMS”, IEEE. Trans. Automatic Control, vol. 42, no. 6, pp. 761-770, 1997.
  * [41] S. Rosset and J. Zhu,“Piecewise linear regularized solution paths”, The Annals of Statistics, vol. 35, no. 3, pp. 1012-1030, 2007.
  * [42] H. W. Engl, M. Hanke, and A. Neubauer,“Regularization of inverse problems”, in Mathematics and its Applications. Boston, MA, USA: Kluwer, 1996.
  * [43] X. Zhang and T. Li, “Online learning in reproducing kernel Hilbert space with non-iid data”, Proc. 62nd IEEE Conference on Decision and Control (CDC), Singapore, 13-15, December, pp. 6610-6615, 2023.
  * [44] P. D. Lax, _Functional Analysis_ , John Wiley & Sons, Now York, 2002.
  * [45] D. X. Zhou, “Capacity of reproducing kernel spaces in learning theory”, IEEE Trans. Information Theory, vol. 49, no.7, pp. 1743-1752, 2003.
  * [46] C. Gu, _Smothing Spline ANOVA Models_ , Springer, 2002.
  * [47] S. Theodoridis, _Machine Learning: A Bayesian and Optimization Perspective_ , London, UK: Elsevier, 2015.


11institutetext: IBM Research, Almaden, San Jose, USA
11email<EMAIL_ADDRESS>22institutetext: Dept. of Electrical and
Computer Eng., Johns Hopkins University, Baltimore, USA

# mSPD-NN: A Geometrically Aware Neural Framework for Biomarker Discovery from
Functional Connectomics Manifolds

Niharika S. D’Souza * 11 Archana Venkataraman 22

###### Abstract

Connectomics has emerged as a powerful tool in neuroimaging and has spurred
recent advancements in statistical and machine learning methods for
connectivity data. Despite connectomes inhabiting a matrix manifold, most
analytical frameworks ignore the underlying data geometry. This is largely
because simple operations, such as mean estimation, do not have easily
computable closed-form solutions. We propose a geometrically aware neural
framework for connectomes, i.e., the mSPD-NN, designed to estimate the
geodesic mean of a collections of symmetric positive definite (SPD) matrices.
The mSPD-NN is comprised of bilinear fully connected layers with tied weights
and utilizes a novel loss function to optimize the matrix-normal equation
arising from Fréchet mean estimation. Via experiments on synthetic data, we
demonstrate the efficacy of our mSPD-NN against common alternatives for SPD
mean estimation, providing competitive performance in terms of scalability and
robustness to noise. We illustrate the real-world flexibility of the mSPD-NN
in multiple experiments on rs-fMRI data and demonstrate that it uncovers
stable biomarkers associated with subtle network differences among patients
with ADHD-ASD comorbidities and healthy controls.

###### Keywords:

Functional Connectomics SPD Manifolds Fréchet Mean Estimation Geometry-Aware
Neural Networks

##  1 Introduction

Resting state functional MRI (rs-fMRI) measures steady state patterns of co-
activation [11] (i.e., connectivity) as a proxy for communication between
brain regions. The ‘connectome’ is a whole-brain map of these connections,
often represented as a correlation or covariance matrix [16] or a network-
theoretic object such as adjacency matrix or graph kernel [10]. The rise of
connectomics has spurred many analytical frameworks for group-wise diagnostics
and biomarker discovery from this data. Early examples include statistical
comparisons of connectivity features [16], aggregate network theoretic
measures [10], and dimensionality reduction techniques [14, 8]. More recently,
the field has embraced deep neural networks to learn complex feature
representations from both the connectome and the original rs-fMRI time series
[2, 18, 7]. While these approaches have yielded valuable insights, they
largely ignore the underlying geometry of the connectivity data. Namely, under
a geometric lens, connectomes derived from rs-fMRI data lie on the manifold of
symmetric positive definite (SPD) matrices. A major computational bottleneck
for developing geometrically-aware generalizations [19, 1] is the estimation
of the geodesic mean on SPD manifolds. This is a far more challenging problem
than statistical estimation in Euclidean data spaces because extensions of
elementary operations such as addition, subtraction, and distances on the SPD
manifold entail significant computational overhead [17].

The most common approach for estimating the geodesic mean on the SPD manifold
is via gradient descent [20]. While this method is computationally efficient,
it is highly sensitive to the step size. To mitigate this issue, Riemannian
optimization methods [12], the majorization-maximization algorithm [25], and
fixed-point iterations [4] can be used. While these extensions have desirable
convergence properties, this comes at the cost of increased computational
complexity, meaning they do not scale well to higher input dimensionality and
larger numbers of samples [3]. In contrast, the work of [3] leverages the
approximate joint diagonalization [21] of matrices on the SPD manifold. While
this approach provides guaranteed convergence to a fixed point, the accuracy
and stability of the optimization is sensitive to the deviation of the data
from the assumed common principal component (CPC) generating process. Taken
together, existing methods for geodesic mean estimation on the SPD manifold
poorly balance accuracy, robustness and computational complexity, which makes
them difficult to fold into a larger analytical framework for connectomics
data.

We propose a novel end-to-end framework to estimate the geodesic mean of data
on the SPD manifold. Our method, the Geometric Neural Network (mSPD-NN),
leverages a matrix autoencoder formulation [9] that performs a series of bi-
linear transformations on the input SPD matrices. This strategy ensures that
the estimated mean remains on the manifold at each iteration. Our loss
function for training approximates the first order matrix-normal condition
arising from Fréchet mean estimation [17]. Using conventional backpropagation
via stochastic optimization, the mSPD-NN automatically learns to estimate the
geodesic mean of the input data. We demonstrate the robustness of our
framework using simulation studies and show that mSPD-NN can handle input
noise and high-dimensional data. Finally, we use the mSPD-NN for various
groupwise discrimination tasks (feature selection, classification, clustering)
on functional connectivity data and discover consistent biomarkers that
distinguish between patients diagnosed with ADHD-Autism comorbidities and
healthy controls.

##  2 Biomarker Discovery from Functional Connectomics Manifolds via the mSPD-
NN

Let matrices $\\{\mathbf{\Gamma}_{n}\\}^{N}_{n=1}\in\mathcal{M}$ be a
collection of $N$ functional connectomes belonging to the manifold
$\mathcal{M}$ of Symmetric Positive Definite (SPD) matrices of dimensionality
$P\times P$, i.e. $\mathcal{M}\in\mathcal{P}^{+}_{P}$ (and a real and smooth
Reimannian manifold). We can define an inner product that varies smoothly at
each vector $\mathcal{T}_{\mathbf{\Gamma}}(\mathcal{M})$ in the tangent space
defined at any point $\mathbf{\Gamma}\in\mathcal{M}$. Finally, a geodesic
denotes the shortest path joining any two points on the manifold along the
manifold surface.

Geodesic Mappings: The matrix exponential and the matrix logarithm maps allow
us to translate geodesics on the manifold back and forth to the local tangent
space at a reference point. The matrix exponential mapping translates a vector
$\mathbf{V}\in\mathcal{T}_{\mathbf{\Phi}}(\mathcal{M})$ in the tangent space
at $\mathbf{\Phi}\in\mathcal{M}$ to a point on the manifold
$\mathbf{\Gamma}\in\mathcal{M}$ via the geodesic emanating from
$\mathbf{\Phi}$. Conversely, the matrix logarithm map translates the geodesic
between $\mathbf{\Phi}\in\mathcal{M}$ to $\mathbf{\Gamma}\in\mathcal{M}$ back
to the tangent vector $\mathbf{V}\in\mathcal{T}_{\mathbf{\Phi}}(\mathcal{M})$.
Mathematically, these operations are parameterized as:

$\displaystyle\mathbf{\Gamma}=\textbf{Expm}_{\mathbf{\Phi}}(\mathbf{V})=\mathbf{\Phi}^{1/2}\textbf{expm}(\mathbf{\Phi}^{-1/2}\mathbf{V}\mathbf{\Phi}^{-1/2})\mathbf{\Phi}^{1/2}$
(1)
$\displaystyle\mathbf{V}=\textbf{Logm}_{\mathbf{\Phi}}(\mathbf{\Gamma})=\mathbf{\Phi}^{1/2}\textbf{logm}(\mathbf{\Phi}^{-1/2}\mathbf{\Gamma}\mathbf{\Phi}^{-1/2})\mathbf{\Phi}^{1/2}$
(2)

Here, $\textbf{expm}(\cdot)$ and $\textbf{logm}(\cdot)$ refer to the matrix
exponential and logarithm respectively, each requiring an eigenvalue
decomposition of the argument matrix, a point-wise transformation of the
eigenvalues, and a matrix reconstruction.

Distance Metric: Given two connectomes
$\mathbf{\Gamma}_{1},\mathbf{\Gamma}_{2}\in\mathcal{M}$, the Fisher
Information distance between them is the length of the geodesic connecting the
two points:

$\delta_{R}(\mathbf{\Gamma}_{1},\mathbf{\Gamma}_{2})={||{\textbf{logm}(\mathbf{\Gamma}^{-1}_{1}\mathbf{\Gamma}_{2})}||}_{F}={||{\textbf{logm}(\mathbf{\Gamma}^{-1}_{2}\mathbf{\Gamma}_{1})}||}_{F},$
(3)

where ${||{\cdot}||}_{F}$ denotes the Frobenius norm. The Reimannian norm of
$\mathbf{\Gamma}$ is the geodesic distance from the identity matrix
$\mathcal{I}$ i.e.
${||{\mathbf{\Gamma}}||}_{R}={||{\textbf{logm}(\mathbf{\Gamma})}||}_{F}$

Figure 1: The mSPD-NN architecture: The input is transformed by a cascade of
2D fully connected layers. The matrix logarithm function is used to obtain the
matrix normal form, which serves as the loss function for mSPD-NN during
training.

###  2.1 Geodesic Mean Estimation via the mSPD-NN:

The geodesic mean of $\\{\mathbf{\Gamma}_{n}\\}$ is defined as the matrix
$\mathbf{G}_{R}\in\mathcal{M}$ whose sum of squared geodesic distances (Eq.
(3)) to each element is minimal [17].

$\mathcal{G}_{R}(\\{\mathbf{\Gamma}_{n}\\})=\operatorname*{arg\,min}_{\mathbf{G}_{R}}\mathbf{L}(\mathbf{G}_{R})=\operatorname*{arg\,min}_{\mathbf{G}_{R}}\sum_{n}{||{\textbf{logm}(\mathbf{G}_{R}^{-1}\mathbf{\Gamma}_{n})}||}^{2}_{F}$
(4)

A pictorial illustration is provided in the green box in Fig 1. While Eq. (4)
does not have a closed-form solution for $N>2$, it is also is convex and
smooth with respect to the unknown quantity $\mathbf{G}_{R}(\cdot)$ [17]. To
estimate population means from the connectomes, mSPD-NN makes use of
Proposition 3.4 from [17].

Proposition 1: The geodesic mean $\mathbf{G}_{R}$ of a collection of $N$ SPD
matrices $\\{\mathbf{\Gamma}_{n}\\}$ is the unique symmetric positive-definite
solution to the nonlinear matrix equation
$\sum_{n}\textbf{logm}(\mathbf{G}^{-1/2}_{R}\mathbf{\Gamma}_{n}\mathbf{G}^{-1/2}_{R})=\mathbf{0}$.
$\mathbf{0}$ is a $P\times P$ matrix of all zeros.

Proof: The proof follows by computing the first order necessary (and here,
sufficient) condition for optimality for Eq. (4). First, we express the
derivative of a real-valued function of the form
$\mathbf{H}(\mathbf{S}(t))=\frac{1}{2}{||{\textbf{logm}(\mathbf{C}^{-1}\mathbf{S}(t))}||}_{F}^{2}$
with respect to $t$. In this expression, the argument
$\mathbf{S}(t)=\mathbf{\mathbf{G}_{R}}^{1/2}\textbf{expm}(t\mathbf{A})\mathbf{G_{R}}^{1/2}$
is the geodesic arising from $\mathbf{G}_{R}$ in the direction of
$\mathbf{\Delta}=\dot{\mathbf{S}}(\mathbf{0})=\mathbf{\mathbf{G}_{R}}^{1/2}\mathbf{A}\mathbf{G_{R}}^{1/2}$,
and the matrix $\mathbf{C}\in\mathcal{P}^{+}_{P}$ is a constant SPD matrix of
dimension $P$. By using the cyclic properties of the trace function and the
distributive equivalence of
$\textbf{logm}(\mathbf{A}^{-1}[\mathbf{B}]\mathbf{A})=\mathbf{A}^{-1}[\textbf{logm}(\mathbf{B})]\mathbf{A}$,
we obtain the following condition:

$\displaystyle\mathbf{H}(\mathbf{S}(t))=\frac{1}{2}{{||{\textbf{logm}(\mathbf{C}^{-1/2}\mathbf{S}(t)\mathbf{C}^{-1/2})}||}_{F}^{2}}$

By the symmetry of the term
$\textbf{logm}(\mathbf{C}^{-1/2}\mathbf{S}(t)\mathbf{C}^{-1/2})$ we have that:

$\displaystyle\therefore{\frac{d}{dt}{\mathbf{H}(\mathbf{S}(t))}}\Bigr{|}_{\begin{subarray}{c}t=0\end{subarray}}=\frac{1}{2}\frac{d}{dt}{\operatorname{Tr}\Big{(}[\textbf{logm}(\mathbf{C}^{-1/2}\mathbf{S}(t)\mathbf{C}^{-1/2})]^{2}\Big{)}}\Bigr{|}_{\begin{subarray}{c}t=0\end{subarray}}$
$\displaystyle\therefore{\frac{d}{dt}{\mathbf{H}(\mathbf{S}(t))}}\Bigr{|}_{\begin{subarray}{c}t=0\end{subarray}}=\operatorname{Tr}\Big{(}[\textbf{logm}(\mathbf{C}^{-1}\mathbf{G}_{R})\mathbf{G}_{R}^{-1}\mathbf{\Delta}]\Big{)}=\operatorname{Tr}[\mathbf{\Delta}\textbf{logm}(\mathbf{C}^{-1}\mathbf{G}_{R})\mathbf{G}_{R}^{-1}]$
$\displaystyle\therefore\nabla{\mathbf{H}}=\textbf{logm}(\mathbf{C}^{-1}\mathbf{G}_{R})\mathbf{G}^{-1}_{R}=\mathbf{G}^{-1}_{R}\textbf{logm}(\mathbf{G}_{R}\mathbf{C}^{-1})$

Notice that since $\nabla{\mathbf{H}}$ is symmetric, it belongs to the tangent
space $\mathcal{S}_{P}$ of $\mathcal{P}^{+}_{P}$. Therefore, we express the
gradient of $\mathbf{L}(\mathbf{G}_{R})$ defined in Eq. (4), as follows:

$\displaystyle\mathbf{L}(\mathbf{G}_{R})=\sum_{n}{||{\textbf{logm}(\mathbf{G}_{R}^{-1}\mathbf{\Gamma}_{n})}||}^{2}_{F}\
\ \
\implies\nabla\mathbf{L}(\mathbf{G}_{R})=\mathbf{G}^{-1}_{R}\sum_{n}\textbf{logm}({\mathbf{G}_{R}\mathbf{\Gamma}^{-1}_{n}})$
$\displaystyle\therefore\operatorname*{arg\,min}_{\mathbf{G}_{R}}\mathbf{L}(\mathbf{G}_{R})\implies\sum_{n}\textbf{logm}({\mathbf{G}_{R}\mathbf{\Gamma}^{-1}_{n}})=\sum_{n}\textbf{logm}({\mathbf{G}^{-1/2}_{R}\mathbf{\Gamma}_{n}\mathbf{G}^{-1/2}_{R}})=\mathbf{0}$

The final step uses the property that $\mathbf{L}(\mathbf{G}_{R})$ is a sum of
convex functions, with the first order stationary point is the necessary and
sufficient condition being the unique minima. Denoting
$\mathbf{G}^{-1/2}_{R}=\mathbf{V}\in\mathcal{P}_{P}^{+}$, the matrix
multiplications in the argument of the $\textbf{logm}(\cdot)$ term can be
efficiently expressed within the feed-forward operations of a neural network
with unknown parameters $\mathbf{V}$.

###  2.2 mSPD-NN Architecture

The mSPD-NN uses the form above to perform geodesic mean estimation. The
architecture is illustrated in Fig. 1. The encoder of the mSPD-NN is a 2D
fully-connected neural network (FC-NN) [5] layer
$\mathbf{\Psi}_{\text{enc}}(\cdot):\mathcal{P}_{P}^{+}\rightarrow\mathcal{P}_{P}^{+}$
that projects the input matrices $\mathbf{\Gamma}_{n}$ into a latent
representation. This mapping is computed as a cascade of two linear layers
with tied weights $\mathbf{W}\in\mathcal{R}^{P\times P}$, i.e.,
$\mathbf{\Psi}_{\text{enc}}(\mathbf{\Gamma}_{n})=\mathbf{W}\mathbf{\Gamma}_{n}\mathbf{W}^{T}$
The decoder $\mathbf{\Psi}_{dec}(\cdot)$ has the same architecture as the
encoder, but with transposed weights $\mathbf{W}^{T}$. The overall
transformation can be written as:

$\text{mSPD-
NN}(\mathbf{\Gamma}_{n})=\mathbf{\Psi}_{\text{dec}}(\mathbf{\Psi}_{\text{enc}}(\mathbf{\Gamma}_{n}))=\mathbf{W}\mathbf{W}^{T}(\mathbf{\Gamma}_{n})\mathbf{W}\mathbf{W}^{T}=\mathbf{V}(\mathbf{\Gamma}_{n})\mathbf{V}$
(5)

where $\mathbf{V}\in\mathcal{R}^{P\times P}$ and is symmetric and positive
definite by construction. We would like our loss function to minimize Eq. (4)
in order to estimate the first order stationary point as
$\mathbf{V}=\mathbf{G}^{-1/2}_{R}$, and therefore devise the following loss:

$\mathcal{L}(\cdot)=\frac{1}{P^{2}}{\Big{|}\Big{|}{\frac{1}{N}\sum_{n}\textbf{logm}\Big{[}{\mathbf{W}\mathbf{W}^{T}(\mathbf{\Gamma}_{n})\mathbf{W}\mathbf{W}^{T}}}\Big{]}\Big{|}\Big{|}}_{F}^{2}$
(6)

Formally, an error of $\mathcal{L}(\cdot)=0$ implies that the argument
satisfies the matrix normal equation exactly under the parameterization
$\mathbf{V}=\mathbf{W}\mathbf{W}^{T}=\mathbf{G}^{-1/2}_{R}$. Therefore, Eq.
(6) allows us to estimate the geodesic mean on the SPD manifold. We utilize
standard backpropagation to optimize Eq. (6). From an efficiency standpoint,
the mSPD-NN architecture maps onto a relatively shallow neural network.
Therefore, this module can be easily integrated into other deep learning
inference frameworks for example, for batch normalization on the SPD manifold.
This flexibility is the key advantage over classical methods, in which
integrating the geodesic mean estimation within a larger framework is not
straightforward. Finally, the extension of Eq. (6) to the estimation of a
weighted mean (with positive weights $\\{w_{n}\\}$) also follows naturally as
a multiplier in the summation.

Implementation Details: We train mSPD-NN for a maximum of $100$ epochs with an
initial learning rate of $0.001$ decayed by $0.8$ every $50$ epochs. The
tolerance criteria for the training loss is set at $1e^{-4}$. mSPD-NN
implemented in PyTorch (v1.5.1), Python 3.5 and experiments were run on an 4.9
GB Nvidia K80 GPU. We utilize the ADAM optimizer during training and a default
PyTorch initialization for the model weights. To ensure that $\mathbf{W}$ is
full rank, we add a small bias to the weights, i.e.,
$\tilde{\mathbf{W}}=\mathbf{W}+\lambda\mathcal{I}_{P}$ for regularization and
stability.

##  3 Evaluation and Results

###  3.1 Experiments on Synthetic Data

We evaluate the scalability, robustness, and fidelity of mSPD-NN using
simulated data. We compare the mSPD-NN against two popular mean estimation
algorithms, the first being the Riemannian gradient descent [20] on the
objective in Eq. (4) and the second being the Approximate Joint
Diagonalization Log Euclidean (ALE) mean estimation [3], which directly
leverages properties of the common principal components (CPC) data generating
process [21].

Our synthetic experiments are built off the CPC model [13]. In this case, each
input connectome $\mathbf{\Gamma}_{n}\in\mathcal{R}^{P\times P}$ is derived
from a set of components $\mathbf{B}\in\mathcal{R}^{P\times P}$ common to the
collection and a set of example specific (and strictly positive) weights
across the components $\mathbf{c}_{n}\in\mathcal{R}^{(+)P\times 1}$. Let the
diagonal matrix $\mathbf{C}_{n}$ be defined as
$\mathbf{C}_{n}=\mathbf{diag}(\mathbf{c}_{n})\in\mathcal{R}^{(+)P\times P}$.
From here, we have
$\mathbf{\Gamma}_{n}=\mathbf{B}\mathbf{C}_{n}\mathbf{B}^{T}$.

Evaluating Scalability: In the absence of corrupting noise, the theoretically
optimal geodesic mean of the examples $\\{\mathbf{\Gamma}_{n}\\}_{n=1}^{N}$
can be computed as: $\mathbf{G}^{*}_{R}=\mathbf{B}\
\mathbf{expm}\left[\frac{1}{N}\sum_{n=1}^{N}\mathbf{logm}(\mathbf{B}^{-1}\mathbf{\Gamma}_{n}\mathbf{B}^{-T})\right]\
\mathbf{B}^{T}$ [3]. We evaluate the scalability of each algorithm with
respect to the dataset dimensionality $P$ and the number of examples $N$ by
comparing its output to this theoretical optimum.

We randomly sample columns of the component matrix $\mathbf{B}$ from a
standard normal, i.e.,
$\mathbf{B}[:,j]\sim\mathcal{N}(\mathbf{0},\mathcal{I}_{P})\ \ \forall\ \
j\in\\{1,\dots,P\\}$, where $\mathcal{I}_{P}$ is an identity matrix of
dimension $P$. In parallel, we sample the component weights $\mathbf{c}_{nk}$
according to $\mathbf{c}_{nk}^{1/2}\sim\mathcal{N}(0,1)\ \ \forall\ \
k\in\\{1,\dots,P\\}$. To avoid degenerate behavior when the inputs are not
full-rank, we clip $\mathbf{c}_{nk}$ to a minimum value of $0.001$. We
consider two experimental scenarios. In Experiment 1, we fix the data
dimensionality at $P=30$ and sweep the dataset size as
$N\in\\{5,10,20,50,100,200\\}$. In Experiment 2, we fix the dataset size at
$N=20$ and sweep the dimensionality as $P\in\\{5,10,20,50,100,200\\}$. For
each parameter setting, we run all three estimation algorithms ten times using
different random initializations.

We score performance based on the correctness of the solution and the
execution time in seconds. Correctness is measured in two ways. First is the
final condition fit $\mathcal{L}(\mathbf{G}^{\text{est}}_{R})$ from Eq. (6),
which quantifies the deviation of the solution from the first order stationary
condition (i.e., $\mathcal{L}(\mathbf{G}^{\text{est}}_{R})=0$). Second is the
normalized squared Riemannian distance
$d_{\text{mean}}={d^{2}_{R}(\mathbf{G}^{\text{est}}_{R},\mathbf{G}^{*}_{R}})/{||{\mathbf{G}^{*}_{R}}||}^{2}_{R}$
between the solution and the theoretically optimal mean. Lower values of the
condition fit $\mathcal{L}(\mathbf{G}_{R})$ and deviation $d_{\text{mean}}$
imply a better quality solution.

Fig. 2 illustrates the performances of mSPD-NN, gradient descent and ALE mean
estimation algorithms. Figs. 2(a) and (d) plot the first-order condition fit
$\mathcal{L}(\mathbf{G}^{\text{est}}_{R})$ when varying the dataset size $N$
(Experiment 1) and the matrix dimensionality $P$ (Experiment 2), respectively.
Likewise, Figs. 2(b) and (e) plot the recovery performance for each
experiment. We observe that the first order condition fit for the mSPD-NN is
better than the ALE for all settings, and better than the gradient descent for
most settings. We note that the recovery performance of mSPD-NN is better than
the baselines in most cases while being a close approximation in the remaining
ones. Finally, Figs. 2(c) and (f) illustrate the time to convergence for each
algorithm. As seen, the performance of mSPD-NN scales with dataset size but is
nearly constant with respect to dimensionality. In all cases, it either beats
or is competitive with ALE.

Figure 2: Evaluating the estimates from mSPD-NN, gradient descent and ALE
according to (a) and (d) first-order condition fit (Eq. 6) (b) and (e)
deviation from the theoretical solution (c) and (f) execution time for varying
dataset size $N$ and data dimension $P$ respectively

Robustness to Noise: Going one step further, we evaluate the efficacy of the
mSPD-NN framework when there is deviation from the ideal CPC generating
process. In this case, we add rank-one structured noise to obtain the input
data:
$\mathbf{\Gamma}_{n}=\mathbf{B}\mathbf{C}_{n}\mathbf{B}^{T}+\frac{1}{P}\mathbf{x}_{n}\mathbf{x}_{n}^{T}$.
As before, the bases and coefficients are randomly sampled as
$\mathbf{B}[:,j]\sim\mathcal{N}(\mathbf{0},\mathcal{I}_{P})$ and
$\mathbf{c}_{nj}^{1/2}\sim\mathcal{N}(0,1)\ \ \forall\ \ j\in\\{1,\dots,P\\}$.
In a similar vein, the structured noise is generated as
$\mathbf{x}_{n}\sim\mathcal{N}(\mathbf{0},\sigma^{2}\mathcal{I}_{P})\in\mathcal{R}^{P\times
1}$, with $\sigma^{2}$ controlling the extent of the deviation. For this
experiment, we set $P=30,N=20$ and vary the noise over the range $[0.2-1]$ in
increments of $0.1$. One caveat in this setup is that the theoretically
optimal mean defined previously and cannot be used to evaluate performance.
Hence, we report only the first-order condition fit
$\mathcal{L}(\mathbf{G}_{R})$ We also calculate the pairwise concordance
$d_{\text{weights}}$ of the final mSPD-NN weights for different
initializations.

Figure 3: Performance of the mSPD-NN, gradient descent and ALE estimation
under increasing additive noise: (a) First order condition fit (Eq. 6) (b)
Pairwise distance between the recovered mSPD-NN solutions across random
initializations.

Fig. 3(a) illustrates the first-order condition fit
$\mathcal{L}(\mathbf{G}^{\text{est}}_{R})$ across all three methods for
increasing noise $\sigma$. As seen, $\mathcal{L}(\mathbf{G}^{est}_{R})$ for
the mSPD-NN is consistently lower than the corresponding value for the
gradient descent and ALE algorithm, suggesting improved performance despite
increasing corruption to the CPC process. The ALE algorithm is designed to
utilize the CPC structure within the generating process, but its poor
performance suggests that it is particularly susceptible to noise. Fig. 3(b)
plots the pairwise distances between the geodesic means estimated by mSPD-NN
across the 10 random initializations. As seen, mSPD-NN produces a consistent
solution, thus underscoring its robustness.

###  3.2 Experiments on Functional Connectomics Data

Dataset: To probe the efficacy of the mSPD-NN for representation learning on
real world matrix manifold data, we experiment on several groupwise
discrimination tasks (such as group-wise discrimination, classification and
clustering) on the publicly available CNI 2019 Challenge dataset [23]
consisting of preprocessed rs-fMRI time series, provided for $158$ subjects
diagnosed with Attention Deficit Hyperactivity Disorder (ADHD), $92$ subjects
with Autism Spectrum Disorder (ASD) with an ADHD comorbidity [15], and $257$
healthy controls. The scans were acquired on a Phillips 3T Achieva scanner
using a single shot, partially parallel, gradient-recalled EPI sequence with
TR/TE = $2500/30$ms, flip angle 70, voxel resolution = $3.05\times 3.15\times
3$mm, with a scan duration of either $128$ or $156$ time samples (TR). A
detailed description of the demographics and preprocessing can be found in
[23]. Connectomes are estimated via the Pearson’s correlation matrix,
regularized to be full-rank via two parcellations, the Automated Anatomical
Atlas (AAL) ($P=116$) and the Craddocks 200 atlas ($P=200$).

Groupwise Discrimination: We expect that FC differences between the ASD and
ADHD cohorts are harder to tease apart than differences between ADHD and
controls [23, 15]. We test this hypothesis by comparing the geodesic means
estimated via mSPD-NN for the three cohorts. For robustness, we perform
bootstrapped trials for mean estimation by sampling $25$ random subjects from
a given group (ADHD/ASD/Controls). We then compute the Riemannian distance
$d(\mathbf{G}_{R}(\\{\mathbf{\Gamma}_{g1}\\}),\mathbf{G}_{R}(\\{\mathbf{\Gamma}_{g2}\\}))$
between the mSPD-NN means associated with groups $g1$ and $g2$. A higher value
of $d(\cdot,\cdot)$ implies a better separation between the groups. We also
run a Wilcoxon signed rank test on the distribution of $d(\cdot,\cdot)$.

Fig. 4 illustrates the pairwise distances between the geodesic means of
cohorts $g1-g2$ across bootstrapped trials (t-SNE representations for the
group means are provided in Fig. 5(c)). As a sanity check, we note that the
mean estimates across samples within the same cohort (ADHD-ADHD) are closer
than those across cohorts (ADHD-controls, ASD-controls, ADHD-ASD). More
interestingly, we observe that ADHD-controls separation is consistently larger
than that of the ADHD-ASD groups for both parcellations. This result confirms
the hypothesis that the overlapping diagnosis for the two classes translates
to a reduced separability in the space of FC matrices and indicates that mSPD-
NN is able to robustly uncover population level differences in FC.

Figure 4: Groupwise discrimination between the FC matrices estimated via the
(a) AAL (b) Craddock’s 200 atlas, for the ADHD/ASD/Controls cohorts according
to pairwise distances between the mSPD-NN mean estimates. Results of pairwise
connectivity comparisons between group means for (c) ADHD-Controls (d) ADHD-
ASD groups for the AAL parcellation. The red connections are significant
differences ($p<0.001$).

Classification: Building on the observation that mSPD-NN provides reliable
group-separability, we adopt this framework for classification. Using the AAL
parcellation, we randomly sample $25$ subjects from each class for training,
and set aside the rest for evaluation with a $10\%/90\%$ validation/test
split. We estimate the geodesic mean for each group across the training
samples via $10$ bootstrapped trials, in which we sub-sample $80\%$ of the
training subjects from the respective group. Permutation testing is performed
on the mean estimates [24], and functional connections (i.e., entries of
$\mathbf{G}_{R}(\\{\mathbf{\Gamma}_{n}\\})$) that differ with an FDR-corrected
threshold of $p<0.001$ are retained for classification. Finally, a Random
Forest classifier is trained on the selected features to classify ADHD vs
Controls. The train-validation-test splits are repeated $10$ times to compute
confidence intervals.

We use classification accuracy and area under the receiver operating curve
(AU-ROC) as metrics for evaluation. The mSPD-NN feature selection plus Random
Forest approach provides an accuracy of $0.62\pm{0.031}$ and an AU-ROC of
$0.60\pm{0.04}$ for ADHD-Control classification on the test samples. We note
that this approach outperforms all but one method on the CNI challenge leader-
board [23]. Moreover, one focus of the challenge is to observe how models
trained on the ADHD vs Control discrimination task translate to ASD (with ADHD
comorbidity) vs Control discrimination in a transfer learning setup.
Accordingly, we apply the learned classifiers in each split to ASD vs Control
classification and obtain an accuracy of $0.54\pm{0.044}$ and an AU-ROC of
$0.53\pm{0.03}$. This result is on par with the best performing algorithm in
the CNI-TL challenge. The drop in accuracy and AU-ROC for the transfer
learning task is consistent with the performance profile of all the challenge
submissions. These results suggest that despite the comorbidity, connectivity
differences between the cohorts are subtle and hard to reliably capture.
Nonetheless, the mSPD-NN+RF framework is a first step to underscoring stable,
yet interpretable (see below) connectivity patterns that can discriminate
between diseased and healthy populations.

Qualitative Analysis: To better understand the group-level connectivity
differences, we plot the most consistently selected features (top 10 percent)
from the previous experiment (ADHD-control feature selection) in Fig. 4(c). We
utilize the BrainNetViewer Software for visualization. The blue circles are
the AAL nodes, while the solid lines denote edges between nodes. We observe
that the highlighted connections appear to cluster in the sensorimotor and
visual areas of the brain, along with a few temporal lobe contributions.
Altered sensorimotor and visual functioning has been previously reported among
children and young adults diagnosed with ADHD [6]. Adopting a similar
procedure, we additionally highlight differences among the ASD and ADHD
cohorts in Fig. 4(d). The selected connections concentrate around the pre-
frontal areas of the brain, which is believed to be associated with altered
social-emotional regulation in Autism [22]. We additionally provide an
extended version of the group connectivity difference results across trials in
Fig. 5 (a) ADHD vs Controls and (b) ADHD vs ASD. Across train-test-validation
splits, we observe that several connectivity differences appear fairly
consistently. Overall, the patterns highlighted via statistical comparisons on
the mSPD-NN estimates are both robust as well as in line with the
physiopathology of ADHD and ASD reported in the literature.

Figure 5: Pairwise differences between mSPD-NN group means for (a) ADHD-
Controls (b) ADHD-ASD groups across bootstrapped trials. Significant
differences marked in red ($p<0.001$). t-SNE plots for group means from
experiment on (c) Groupwise Discrimination using mSPD-NN (d) After data-driven
clustering via the mSPD-EM

Data-Driven Clustering: Finally, we evaluate the stability of the mapping
between the functional connectivity and diagnostic spaces via a geometric
clustering experiment. We use the geodesic mean estimates from the groupwise
discrimination experiment (generated using the ground truth Controls/ASD/ADHD
labels and mSPD-NN) as an initialization and track the shift in the diagnostic
assignments upon running an unsupervised Expectation-Maximization (EM)
algorithm. At each iteration of the mSPD-EM, the E-Step assigns cluster
memberships to a given subject according to the geodesic distance (Eq. (3))
from the cluster centroids, while the M-Step uses the mSPD-NN for recomputing
the centroids. Upon convergence, we evaluate the alignment between the
inferred clusters and diagnostic labels. To this end, we map each cluster to a
diagnostic label according to majority voting, and measure the cluster purity
(fraction of cluster members that are correctly assigned). mSPD-EM provides an
overall cluster purity of $0.59\pm 0.05$ (Controls), $0.52\pm 0.12$ (ADHD),
ASD $0.51\pm 0.09$ (ASD), indicating that there is considerable shift in the
assignment of diagnostic labels from ground truth. We also visualise the
cluster centroids using t-Stochastic Neighbor Embeddings (t-SNE) at
initialization and after convergence of the mSPD-EM in Fig. 5 (c) and (d)
respectively. We provide 3-D plots to better visualise the cluster separation.
Again, we observe that the diagnostic groups overlap considerably and are
challenging to separate in the functional connectivity space alone. One
possible explanation may be that the distinct neural phenotypes between the
disorders are being overwhelemed by other rs-fMRI signatures. Given the
migration of diagnostic assignments from the ground truth, the strict
diagnostic criteria used to separate the diseased and healthy cohorts group
may need to be more critically examined.

##  4 Conclusion

We have proposed a novel mSPD-NN framework to reliably estimate the geodesic
mean of a collection of functional connectivity matrices. Through extensive
simulation studies, we demonstrate that the mSPD-NN scales well to high-
dimensional data and can handle input noise when compared with current
iterative methods. By conducting a series of experiments on group-wise
discrimination, feature selection, classification, and clustering, we
demonstrate that the mSPD-NN is a reliable framework for discovering
consistent group differences between patients diagnosed with ADHD-Autism
comorbidities and controls. The mSPD-NN makes minimal assumptions about the
data and can potentially be a useful tool to advance data-scientific and
clinical research.

#### Acknowledgements

This work is supported by the National Science Foundation CAREER award 1845430
(PI Venkataraman), the National Institute of Health R01HD108790 (PI
Venkataraman) and R01EB029977 (PI Caffo).

## References

  * [1] Banerjee, M., et al.: Nonlinear regression on riemannian manifolds and its applications to neuro-image analysis. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 719–727. Springer (2015)
  * [2] Bessadok, A., Mahjoub, M.A., Rekik, I.: Graph neural networks in network neuroscience. IEEE Transactions on Pattern Analysis and Machine Intelligence (2022)
  * [3] Congedo, M., Afsari, B., Barachant, A., Moakher, M.: Approximate joint diagonalization and geometric mean of symmetric positive definite matrices. PloS one 10(4), e0121423 (2015)
  * [4] Congedo, M., Barachant, A., Koopaei, E.K.: Fixed point algorithms for estimating power means of positive definite matrices. IEEE Transactions on Signal Processing 65(9), 2211–2220 (2017)
  * [5] Dong, Z., et al.: Deep manifold learning of symmetric positive definite matrices with application to face recognition. In: Thirty-First AAAI Conference on Artificial Intelligence (2017)
  * [6] Duerden, E.G., Tannock, R., Dockstader, C.: Altered cortical morphology in sensorimotor processing regions in adolescents and adults with attention-deficit/hyperactivity disorder. Brain research 1445, 82–91 (2012)
  * [7] D’Souza, N.S., et al.: Integrating neural networks and dictionary learning for multidimensional clinical characterizations from functional connectomics data. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 709–717. Springer (2019)
  * [8] D’Souza, N.S., et al.: A joint network optimization framework to predict clinical severity from resting state functional mri data. NeuroImage 206, 116314 (2020)
  * [9] D’Souza, N.S., et al.: A matrix autoencoder framework to align the functional and structural connectivity manifolds as guided by behavioral phenotypes. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 625–636. Springer (2021)
  * [10] Fornito, A., Zalesky, A., Breakspear, M.: Graph analysis of the human connectome: promise, progress, and pitfalls. Neuroimage 80, 426–444 (2013)
  * [11] Fox, M.D., et al.: Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. Nat. Rev. Neuro. 8(9), 700 (2007)
  * [12] Jeuris, B.: Riemannian optimization for averaging positive definite matrices (2015)
  * [13] Jolliffe, I.T., Cadima, J.: Principal component analysis: a review and recent developments. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 374(2065), 20150202 (2016)
  * [14] Khosla, M., et al.: Machine learning in resting-state fmri analysis. Magnetic resonance imaging 64, 101–121 (2019)
  * [15] Leitner, Y.: The co-occurrence of autism and attention deficit hyperactivity disorder in children–what do we know? Frontiers in human neuroscience 8, 268 (2014)
  * [16] Lindquist, M.A.: The statistical analysis of fmri data. Statistical science 23(4), 439–464 (2008)
  * [17] Moakher, M.: A differential geometric approach to the geometric mean of symmetric positive-definite matrices. SIAM journal on matrix analysis and applications 26(3), 735–747 (2005)
  * [18] Nandakumar, N., et al.: A multi-task deep learning framework to localize the eloquent cortex in brain tumor patients using dynamic functional connectivity. In: MICCAI Workshop on Machine Learning in Clinical Neuroimaging, pp. 34–44. Springer (2020)
  * [19] Nguyen, X.S., et al.: A neural network based on spd manifold learning for skeleton-based hand gesture recognition. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 12036–12045 (2019)
  * [20] Pennec, X., Fillard, P., Ayache, N.: A riemannian framework for tensor computing. International Journal of computer vision 66(1), 41–66 (2006)
  * [21] Pham, D.T.: Joint approximate diagonalization of positive definite hermitian matrices. SIAM Journal on Matrix Analysis and Applications 22(4), 1136–1152 (2001)
  * [22] Pouw, L.B., et al.: The link between emotion regulation, social functioning, and depression in boys with asd. Research in Autism Spectrum Disorders 7(4), 549–556 (2013)
  * [23] Schirmer, M.D., et al.: Neuropsychiatric disease classification using functional connectomics-results of the connectomics in neuroimaging transfer learning challenge. Medical image analysis 70, 101972 (2021)
  * [24] Zalesky, A., Fornito, A., Bullmore, E.T.: Network-based statistic: identifying differences in brain networks. Neuroimage 53(4), 1197–1207 (2010)
  * [25] Zhang, T.: A majorization-minimization algorithm for the karcher mean of positive definite matrices. arXiv preprint arXiv:1312.4654 (2013)