We consider multivariate functional data \(\mathbf {X}(\mathbf {t})=({X}^{(1)}(t^{(1)}),\ldots , {X}^{(p)}(t^{(p)}))^\top \in \mathbb {R}^p~(p \ge 1)\) . Following the definitions from [1]} ([1]}), \(\mathbf {t}:=(t^{(1)},\ldots ,t^{(p)}) \in \mathcal {T}:=\mathcal {T}_1 \times \cdots \times \mathcal {T}_p\) , where \(\mathcal {T}_{j}~(j=1,\ldots ,p)\)  is a compact set in \(\mathbb {R}^{d_j}\) . \(X^{(j)}(t^{(j)}): \mathcal {T}_j \xrightarrow{}\mathbb {R}\)  is assumed to be in \({L}^2(\mathcal {T}_j)\) . Hence, \(\mathbf {t}\)  is a \(p\) -tuple of \(d_1\) ,..., \(d_p\) -dimensional vectors and not a scalar. In practice, instead of \(\mathbf {X}(\mathbf {t})\) , we observe \({\mathbf {Y}}(\mathbf {t})=({Y}^{(1)}(t^{(1)}),\ldots , {Y}^{(p)}(t^{(p)}))^\top \in \mathbb {R}^p\) , due to measurement errors \(\mathbf {\epsilon }_i=(\epsilon ^{(1)}_i,\ldots ,\epsilon ^{(p)}_i)^\top \stackrel{i.i.d.}{\sim }\mathcal {N}_p(\mathbf {0},\hbox{diag}\lbrace \sigma _1^2,\ldots ,\sigma _p^2\rbrace )\) . Taking the \(i\) th observation as an example, \({\mathbf {Y}}_i(\mathbf {t}_i)\)  is denoted as
\({\mathbf {Y}}_i(\mathbf {t}_i)=\mathbf {X}_i(\mathbf {t}_i)+\mathbf {\epsilon }_i=\mathbf {\mu }(\mathbf {t}_i)+\sum _{m=1}^{\infty }\rho _{i,m}\mathbf {\psi }_{i,m}(\mathbf {t}_i)+\mathbf {\epsilon }_i,~~i=1,\ldots ,n,\) 
