empty Latent Transport Models for
Multivariate Functional Data
[Research supported by NSF Grant DMS-2014626. Corresponding author<EMAIL_ADDRESS>
Cody Carroll and Hans-Georg Müller
Department of Statistics, University of California, Davis

July 2021

Multivariate functional data present theoretical and practical complications which are not found in  univariate functional data. One of these is a situation  where the component functions of multivariate functional data are positive and are subject to mutual time warping. That is, the component processes exhibit a similar shape but are subject to systematic phase variation across their time domains. We introduce  a novel model for multivariate functional data that incorporates such mutual time warping via nonlinear transport functions. This model allows for meaningful interpretation and is well suited to represent functional vector data.  The proposed approach combines a random amplitude factor for each component with population based registration across the components of a multivariate functional data vector and also includes a latent population function, which corresponds to  a common underlying trajectory as well as  subject-specific warping  component.
We also propose estimators for all components of the model. The proposed approach not only leads to a novel representation for multivariate functional data, but is also useful for  downstream analyses such as Fréchet regression.  Rates of convergence are established  when curves are fully observed or observed with measurement error. The usefulness of the model, interpretations and practical aspects are illustrated in simulations and with application to multivariate human growth curves as well as multivariate environmental pollution data.

KEY WORDS:  Functional data analysis, multivariate functional data,  component processes, time warping, cross-component registration, longitudinal studies.

§ INTRODUCTION

The increased prevalence of complicated and high-dimensional data in recent years has brought along with it the need for improved statistical techniques to guide practical data analyses. One such type of data with complex structure is functional data, in which observations are functions $X_i(t),~i =1, \dots, n, $ where the index $t$ (often representing time) ranges over some continuum  $\mathcal{T}\subset \mathbb{R}$.

Functional data analysis (FDA) has found important applications  in many fields of research (e.g. biology, ecology, economics) and  has spawned considerable methodological work as a subfield of statistics <cit.>. In particular, the analysis of  univariate functional data has driven the majority of developments in this area  such as functional principal component analysis <cit.>, regression <cit.>, and clustering <cit.>. Much less work exists for multivariate functional data, which  consist of samples of a finite dimensional vectors whose elements are random functions, while such data
arise in many applications.

Dimension reduction is a common approach to handle functional data with several components, with many studies focusing on extending univariate functional principal components analysis to the multivariate case <cit.>. Alternative approaches include  projections of the component functions
or decomposition into marginal component processes and their interactions <cit.>. Most methodological work was focused on traditional amplitude variation-based models for dimension reduction, though phase variation-based methods for multivariate functional data have also been investigated; <cit.> proposed a method for estimating multivariate structural means and <cit.> introduced a model for clustering multivariate functional data in the presence of phase variation. <cit.> combined the notions of dimension reduction and phase variability through a multivariate version of the shape-invariant model <cit.>, in which component processes share a common latent structure which is time-shifted across components. However the assumption of a rigid shift-warping framework imposes a major parametric constraint on the warping structure and often the class of models described with simple shifts is not rich enough to apply to real-world datasets. Our main contribution is a less-restrictive alternative, in which time characterization of individual-specific temporal effects and component-specific effects is achieved through a nonparametric deformation model.

A major motivation for this framework is
that the component functions of a multivariate data vector  often share a common structure which is subject to variation across modalities; the fundamental shape of growth curves is similar but not identical in timing patterns across body parts, for instance.
In addition to  a  latent shape  curve, a such a  shape-sharing model also requires characterization of individual- and component-level variation. While the idea of component-specific time effects has been previously explored through a time-shifting model, this was done within a restrictive parametric framework with limited practical applicability <cit.>. Allowing for more flexible and nuanced component effects substantially increases the viability of a component-warping  approach and allows for improved data fidelity when describing component-specific effects which inform the time-dynamics of a larger system at work. This motivates  a novel representation of multivariate functional data through a Latent Transport Model,  which combines  tools from time warping <cit.> and optimal transport theory <cit.>.

The organization of this paper is as follows. Section 2 discusses existing approaches to univariate curve registration and introduces the Latent Transport Model for component-warped multivariate functional data. We derive estimators of model components in Section 3 and illustrate the utility and performance of our methodology through data analysis and simulation studies in Sections 4 and 5, respectively. Lastly, asymptotic results are established in Section 6.

§ CURVE REGISTRATION AND THE LATENT TRANSPORT MODEL

§.§ The Univariate Curve Registration Problem

A classical univariate curve registration problem occurs when one has  a sample of curves $X_i(t),~ i = 1,\dots,n,$ observed on an interval $T$ which are realizations of a fixed template $\xi(t)$ that is randomly time-warped and is characterized by  monotone time-warping random functions $h_i(t)$, which correspond to  random homeomorphisms of $T$,

\begin{equation}
X_i(t) = (\xi \circ h_i)(t), \quad \text{for all} \quad t\in \mathcal{T},\quad i=1,\dots,n.
\label{eq.1}
\end{equation}

Curve registration aims at estimating  the time warping functions $h_i$, which are typically considered to be nuisance effects, in order to adjust for them before proceeding with further analysis, e.g., the estimation of $\xi(t)$. A version of this  problem was first considered in <cit.> with the  dynamic time warping algorithm.
A common approach is to align processes to some template function that  carries the main features  common across subjects. Versions of this are known as
landmark-based registration  <cit.>, pairwise curve alignment <cit.> or the Procrustes approach <cit.>, among many others; see

While the curve registration literature is varied and rich in methodology,
any reasonably performing method that satisfies some basic regularity and identifiability conditions can be used for our purposes; the regularity conditions are indispensable to obtain theoretical guarantees.
We emphasize that our emphasis  here is not on univariate warping,  but rather on extending curve registration methodology to  the multivariate case and to develop a  Latent Transport Model, which provides a new representation for multivariate functional data.

§.§ A Unified Model for Multivariate Time Dynamics

Let $\{X_j\}_{j=1}^p$ denote a generic set of random functions with positive component processes $X_j$ in $L_2(\mathcal{T})$ defined on an interval $\mathcal{T}$ that for simplicity we assume to be $[0,1]$.
For such processes we introduce a  Latent Transport Model (LTM). This model provides a representation of functional data that is particularly attractive when the functional forms of the component processes $X_j,~ j\in\{1,\dots,p\}$ have similar shapes.  The general  idea is to pool the information available for individual component processes to arrive at a useful representation of the multivariate process.
Denoting a random sample from a $p$-dimensional stochastic process by $\{\textbf{X}_i\}_{i=1}^n$, where $\bm{X}_i(t) = (X_{i1}(t),\dots,X_{ip}(t))^T$, at the core of this shared structure is a nonrandom latent curve $\lambda$, which characterizes the component curves through

\begin{equation}
X_{ij}(G^{-1}_{ij}(t)) = A_{ij} \lambda(t) , \quad  i = 1,\dots, n,~j = 1,\dots, p.
\label{eq.2}
\end{equation}

The random amplitude factors $A$ and random time distortion functions $G$ reflect differences in realized curves across components and individuals. Without loss of generality we assume $\underset{t \in \mathcal{T}}{\sup}|\lambda(t)|=||\lambda||_\infty = 1 $ since it is always possible to rescale the latent curve without changing the model by defining new amplitude factors $\tilde{A}_{ij} := A_{ij} ||\lambda||_\infty$ and $\tilde{\lambda}(t) = \lambda(t)/||\lambda||_\infty$.

The distortion functions $G$ are elements of $\mathcal{W}$, the convex space of all smooth, strictly increasing functions with common endpoints, i.e., $\mathcal{W}:=\{g:\mathcal{T}\to \mathcal{T}~ |~g\in C^2(\mathcal{T}), ~ g(T_1)=T_1, ~ g(T_2)=T_2,~ g~\text{is a strictly increasing homeomorphism}\}$.  The elements of this space represent random homeomorphisms of the time domain and capture the presence of non-linear phase variation. We further assume that  functions $G$ may be decomposed into the mixed-effects form,

\begin{equation}
G_{ij}(t) = (\Psi_j \circ H_i) (t)  , \quad  i = 1,\dots, n,~ j = 1,\dots, p,
\label{eq.3}
\end{equation}

where the deterministic functions $\Psi$ describe the component-based effects of time distortion and the random functions $H$ describe the subject-level phase variation. In this decomposition the mapping $\Psi_j$ conveys the internal time scale of the $j^{th}$ component, while $H_i$ carries the internal time scale of the $i^{th}$ subject. These mappings may be seen as time transports from standard clock time, $id(t)\equiv t$, to the  time of a given component or individual. As such we refer to the collection of functions $\varPsi = \{\Psi_j: j = 1,\dots, p\}$ as component-level transport functions and the collection of functions $\mathcal{H} = \{H_i: i = 1,\dots, n\}$ as subject-level transport functions.

For the random warping functions $H_i$
we assume that
$EH^{-1}_i(t) = t$ for $t \in \mathcal{T}$. This assumption is sometimes referred to as “standardizing" the registration procedure <cit.> and is commonly made so that there is no net variation.

The component transport functions are also assumed to be standardized, but because they are deterministic and not random, the assumption becomes  $\frac{1}{p}\sum_{j=1}^p\Psi^{-1}_j (t) = t$ for $t \in \mathcal{T}$. Together these conditions imply $E(\frac{1}{p}\sum_{j=1}^p G^{-1}_{ij} (t)) = t$.
Combining (<ref>) and (<ref>) yields the Latent Transport Model (LTM) for multivariate functional data,

\begin{equation}
X_{ij}(t) = A_{ij} \left(\lambda \circ \Psi_j \circ H_i\right) (t)  , \quad i = 1,\dots, n, ~ j = 1,\dots, p.
\label{eq.4}
\end{equation}

With component-warped versions of the latent curve  $\gamma_j = \lambda \circ \Psi_j$, the LTM can be equivalently expressed as

\begin{equation}
X_{ij}(t) = A_{ij} \left(\gamma_j \circ H_i\right) (t) , \quad i = 1,\dots, n, ~ j = 1,\dots, p.
\label{eq.5}
\end{equation}

In this form, the curves $\gamma_j(t)$ convey the “typical" time progression of the latent curve according to the $j^{th}$ component's system time, so we refer to them  as the $j^{th}$ component tempo functions. The component tempo functions can be viewed as the synchronized processes for each component after accounting for random subject-level time distortions. Figure <ref> provides a schematic of the data generating mechanism of the LTM for a simulated dataset.

Schematic of the Latent Transport Model, where $\lambda$  denotes the latent base curve (top-left), $\Psi$ denotes component transports (bottom-left), $\gamma$ denotes component tempos (top-center), $H$ denotes random subject-wise time distortion functions (bottom-right), and $X$  denotes the observed multivariate curve data (top-right) resulting from the complete data generating mechanism.

§.§ Cross-Component Transport Maps

Marginal Cross-Component Transports

To understand and quantify the relative timings between any pair of components, $j,k\in\{1,\dots,p\}$, it is useful to define their cross-component transport $T_{jk}$, which is the transport that, when applied to $j^{th}$ component, maps its tempo to that of the $k^{th}$ component, i.e.,

\begin{equation}
T_{jk} = \Psi^{-1}_{j}\circ\Psi_{k},
\end{equation}

so that $\gamma_j(T_{jk}) = \lambda \circ \Psi_{j} \circ \Psi^{-1}_{j} \circ \Psi_{k}  = \lambda \circ \Psi_{k} = \gamma_k$. Because on domains $[0,1]$ the component transports $\Psi_{k}$ can be represented as distribution functions and are closed under composition, the cross-component transports (XCT) may also be represented as  distribution functions.
For example, consider a pair of component processes, Component A and Component B, for which Component A tends to lag behind the latent curve, while the Component B precedes it. An example of this can be seen in the red and orange curves, respectively, in Figure <ref>. The corresponding red transport, $\Psi_{A}$, falls below the diagonal and conveys the lagged tempo, while the orange transport, $\Psi_{B}$ lies above the diagonal and expresses an accelerated time. The transport function $T_{AB}$ then sits above the diagonal and represents the time-acceleration needed to bring the red tempo in line with the orange component.

Subject-Level Cross-Component Transports

While the marginal XCTs describe the general time relations between components on a population level, we may also be interested in  how an individual's component processes relate to one another. This perspective may be especially useful when trying to understand  intercomponent dynamics which are mediated by covariate effects. Conceptually it is straightforward to extend the notion of cross-component transports to individuals by searching for the warping function $T^{(i)}_{jk}$ which brings the $i^{th}$ individual's $j^{th}$ component in line with the $k^{th}$. A natural definition under the LTM is then

\begin{equation}
T^{(i)}_{jk} = G^{-1}_{ij} \circ G_{ik},
\label{indiv_xct}
\end{equation}

since this choice gives $X_{ij}\circ T^{(i)}_{jk} \propto A_{ij} (\lambda \circ G_{ij} \circ G^{-1}_{ij} \circ G_{ik}) \propto (\lambda \circ G_{ik}) \propto X_{ik}$. In practice, this proportionality will become equality once random amplitude factors are dealt with during estimation.
The proposed cross-component transports  represent intercomponent time dynamics and do not rely on an absolute clock time. Statistics based on the XCT model can be used in downstream analyses like hypothesis testing and regression. Data illustrations where this comes into play are given in Section <ref>

§ MODEL ESTIMATION AND CURVE RECONSTRUCTION

§.§ Pairwise Warping

The model estimation procedure proposed here relies on solving several univariate warping problems of type $(1)$. Any of the warpingmethods described in Section 2 may be used for
practical implementation.
We use here the pairwise alignment method of <cit.>, which provides an explicit representation of warping functions and satisfies some properties required by our theory in order to derive convergence rates. This pairwise alignment  is  easily  implemented with
the R package  <cit.>.

Fix a component $j$ and to simplify notation write $U_i = X_{ij}/||X_{ij}||_\infty$ to denote a curve normalized by its maximum. For any two curves $U_i, U_{i'}, ~i,i'\in\{1,\dots,n\}$, the pairwise warping function $V_{i'i}: \mathcal{T}\to \mathcal{T}$ is the homeomorphism of the time domain which aligns $U_{i'}$ to $U_i$,

\begin{equation}
V_{i'i}(t) = H^{-1}_{i'}(H_{i}(t)), \quad\text{for all } t \in \mathcal{T}.
\end{equation}

Following <cit.>, we  assume that the functions $v$ can be represented  by a linear spline.  Let $t_\ell = \ell/(L+1),$ for $\ell = 1, \dots, L, $ denote $L$ equidistant knots over $\mathcal{T}$ with $t_0 = 0$ and $t_{L+1}= 1$.  For a warping function $v$, let  $\theta=\theta(v) = [v(t_1),\dots,v(t_\ell)]^T$ so that

\begin{equation}
v(t) = \theta^T\alpha(t), \quad\text{for all } t \in \mathcal{T},
\label{eq:spline}
\end{equation}

where $\alpha(t) = [\alpha_1(t) , \dots, \alpha_{L+1}(t)]^T$ and $\alpha_\ell(t)$ are linear basis functions defined as $\alpha_\ell(t) = \alpha^{(1)}_\ell(t) - \alpha^{(2)}_{\ell+1}(t)$ with $ \alpha^{(1)}_\ell(t) = \dfrac{t - t_{\ell-1}}{t_\ell- t_{\ell-1}} \textbf{1}_{[t_{\ell-1},t_\ell)}$ and $ \alpha^{(2)}_\ell(t) = \dfrac{t - t_{\ell}}{t_\ell- t_{\ell-1}} \textbf{1}_{[t_{\ell-1},t_\ell)}$, for $\ell = 1,..,L+1$.

The parameter space $\Theta$ of the coefficient vector is

\begin{equation*}
\Theta = \{\theta\in\mathbb{R}^{L+1}: 0<\theta_1<\dots<\theta_{L+1}=1\}
\end{equation*}

so that the monotonicity constraint is satisfied.
We then express the pairwise warping function between subjects $i$ and $i'$ as
$V_{i'i}(\cdot) = \theta^T\alpha(\cdot),$ with  $\theta\in\Theta$ and
the spline coefficient vector is obtained by solving the minimization problem

\begin{align}
\tilde{\theta}_{V_{i'i}} &= \underset{\theta\in\Theta}{\argmin} ~\mathcal{C}_{\eta_1}(\theta, U_i, U_{i'}), \nonumber \\
\quad\text{where~}\mathcal{C}_{\eta_1}(\theta, U_i, U_{i'}) &=  \int_\mathcal{T} d^2\left(U_{i'}(\theta^T\alpha(t)),U_i(t)\right) dt + {\eta_1}\int_\mathcal{T}(\theta^T\alpha(t)-t)^2dt.
\label{eq:objfun}
\end{align}

Here  $d$ is a distance  and $\eta_1$ is a regularization parameter, which penalizes against heavy warping and guards against the “pinching effect"

In the original pairwise warping framework of <cit.>, the default distance was taken as $d^2(f(t),g(t)) = \int (f(t)-g(t))^2 \, dt$, which we also adopt here for the pre-standardized curves as
defined right after   (<ref>).  We note that without the pre-standardization, this would correspond to using the semi-metric

$d^2_{\infty}(f(t),g(t)) =\left|\dfrac{f(t)}{~~||f||_\infty}-\dfrac{g(t)}{~~||g||_\infty}\right|^2$ in  (<ref>).
Once the $\tilde{\theta}_{V_{i'i}}$ have been obtained, we obtain representations of the pairwise subject warping functions

\begin{equation}
\tilde{V}_{i'i}(t) = \tilde{\theta}_{V_{i'i}}^T \alpha(t)\quad\text{for all } t \in \mathcal{T}.
\label{eq:pw}
\end{equation}

Since we assume $E(H^{-1}_i(t)) = t$ for all $t \in \mathcal{T}$, we have the identity  $E[V_{i'i}(t)|H_i] = E[H^{-1}_{i'}\left(H_{i}(t)\right)|H_i]\\ = H_i(t)$, the empirical version of which leads to the estimates

\begin{equation}
\tilde{H}_{i}(t) = \frac{1}{n}\sum_{i'=1}^n \tilde{V}_{i'i}(t)\quad\text{for all } t \in \mathcal{T}.
\label{pw:est}
\end{equation}

§.§ Component-wise Alignment

For a fixed component $j$, consider  the sample of univariate curves, $S_j:=\{X_{ij}\}_{i=1}^n$ that are contained in  the sample of observed curves of size $n$ with associated normailzed curves
$X^*_{ij} = X_{ij}/||X_{ij}||_\infty$. Estimates of $\gamma_j$ and $H_i$ for the $j^{th}$ component are suggested by the relation

\begin{align}
X^*_{ij}(t) &= (\lambda\circ\Psi_j \circ H_i)(t),
\label{eq_univ}
\end{align}

which falls within the warping framework of type (<ref>) with $\xi = \lambda\circ\Psi_j$, and $h_i = H_i$.
Replacing $X$ by $X^*$ in $(\ref{eq_univ})$ eliminates the random amplitude factors $A_{ij}$, because
$||X_{ij}||_\infty = A_{ij} ||\lambda\circ G_{ij}||_\infty = A_{ij}$; thus the normalized curves $X^*_{ij}(t) = (\lambda\circ\Psi_j\circ H_i)(t)$ do not depend on the $A_{ij}$.

Applying an estimation method like pairwise warping for each of the subcollections $S_1,\dots,S_p$, results in $p$ estimates of the subject-level warping function, $\tilde{H}^{(1)}_i(t),\dots,\tilde{H}^{(p)}_i(t)$. Taking the mean of the resulting $p$ warping functions then is a natural estimate for the subject-specific warp,

\begin{equation}
\hat{H}_i = p^{-1} \sum_{j=1}^p \tilde{H}^{(j)}_i, \quad i=1,\dots,n.
\end{equation}

Under this definition, the  overall penalty parameter associated with the warping function is $\eta_1= \underset{1\leq j \leq p}{\max}~\eta_{1j},$ where $\eta_{1j}=10^{-4}\times  \{n^{-1} \sum_{i=1}^n
\int_\mathcal{T} (X_{ij}(t)-\bar{X}_j(t))^2dt\}, ~j=1,\dots,p,$ is the default choice of penalty parameter for each of the $p$ registrations  <cit.>. A plug-in estimate of $\gamma_j$ is then obtained by averaging the component-aligned curves,

\begin{equation}
\hat{\gamma}_j = n^{-1} \sum_{i = 1}^n (X_{ij}\circ \hat{H}_i^{-1})/||X_{ij}||_\infty,\quad\text{for}~j=1,\dots,p.
\end{equation}

§.§ Global Alignment and Latent Curve Estimation

A central idea of the LTM is  that all univariate curves $X_{ij}$ contain information about the latent template function $\lambda$.
We therefore consider the full collection of curves,  $S=\cup_{j=1}^p S_j$ and for the estimation of  each individual $i$, select one of its component curves at random as a representative. Call this representative curve $Z_i$ and denote its normalized counterpart by $Z^*_i$. Selecting one of the components at random ensures that we have $P(Z_i = X_{ij}) = 1/p$ for all $i=1,\dots,n,~j=1,\dots,p$. The collection of curves $\{Z_i,~i=1,\dots,n\}$ can be thought of as realizations of $\lambda$ subject to some random distortion $D_i$, where $D_i=G_{ij}$ if the $j^{th}$ component curve is selected. Define $I_{ij}$ as the event that the curve $Z_i$ comes from the collection of $j^{th}$ component curves, $S_j$. Conditional on the event $I_{ij}$ (which happens with probability $1/p$ for all $i=1,\dots,n$), it follows that $D_{i} = G_{ij}= \Psi_j
\circ H_i$.  Then, on average

\begin{equation}
E[D^{-1}_{i}] = E\{E[D^{-1}_{i}|I_{ij}]\} = \sum_{j=1}^p E[H^{-1}_i\circ \Psi^{-1}_j] P(I_{ij}) =  p^{-1}\sum_{j=1}^p \Psi^{-1}_j = id,
\end{equation}

motivating  the warping problem

\begin{equation}
Z^*_i = \lambda \circ D_i, \quad\text{for $i=1,\dots,n$}.
\end{equation}

The critical implication of this relation is that if we expand our scope to the full collection $S$ and apply a traditional method like pairwise warping to obtain $\hat{D}_{i}$ for all $i=1,\dots,n$, the latent curve can be estimated by averaging the globally-aligned curves,

\begin{equation}
\hat{\lambda} = n^{-1} \sum_{i = 1}^{n} (Z_{i}\circ \hat{D}_{i}^{-1})/||Z_{i}||_\infty.
\end{equation}

The estimators of the component transports are motivated by recalling that

\begin{align*}
{\gamma}_j = \lambda \circ \Psi_j , \quad  j=1,\dots.,p.
\end{align*}

Using a spline representation similar to $(\ref{eq:spline})$, we write

\begin{equation}
\Psi_j(t) = \theta^T \alpha(t)
\label{eq:psispline}
\end{equation}

and estimate the component warps by solving the penalized minimization problem,

\begin{align}
\begin{split}
\tilde{\theta}_{\Psi_{j}} &= \underset{\theta\in\Theta}{\argmin} ~\mathcal{C}_{\eta_2}(\theta; \hat{\gamma}_j, \hat{\lambda}),\\
\mathcal{C}_{\eta_2}(\theta; \hat{\gamma}_j, \hat{\lambda})&=  \int_\mathcal{T} d^2\left(\hat{\gamma}_j,\hat{\lambda}(\theta^T\alpha(t))\right) dt+\eta_2\int_\mathcal{T}(\theta^T\alpha(t)-t)^2dt,
\end{split}
\label{eq:objfun2}
\end{align}

$\eta_{2}=10^{-4}\times  \{p^{-1} \sum_{j=1}^p
\int_\mathcal{T} (\hat{\gamma}_j(t)-\hat{\lambda}(t))^2dt\}$ as the default choice of  penalty parameter in line with <cit.>.

Finally, we estimate the component warps as

\begin{equation}
\hat{\Psi}_j(t) = \tilde{\theta}_{\Psi_j}^T \alpha(t).
\end{equation}

Note that under the assumption of fully observed curves without measurement error, the amplitude factors $A_{ij}=||X_{ij}||_\infty$ are known. Often in practice, this is not realistic, and the factors must be estimated by, e.g., $\hat{A}_{ij}=||\tilde{X}_{ij}||_\infty$ where $\tilde{X}$ denotes a smoothing estimate of a function $X$ that is observed with noise, as described in the following section.

§.§ Measurement Error

In practice, the functions $X_{ij}$ are often contaminated with measurement error and available only on a discrete grid. In this situation an initial step is to perform smoothing on the discrete and noisy observations before applying the model estimation method of Sections 3.2 and 3.3.
Let $\{\bm{X}_i, ~i=1,\dots,n\}$ be a random sample of a $p$-dimensional stochastic process $\bm{X}$ in $L_2(\mathcal{T})\times\cdots\times L_2(\mathcal{T})$, where  $\bm{X}_i(t) = (X_{i1}(t),\dots,X_{ip}(t))^T$. We assume that these processes are observed for $m$ equispaced points at the discrete times $t_{s}, ~s=1,\dots,m$. We write $\bm{Y}_{is} = (Y_{i1s},\dots,Y_{ips})^T$ to represent the $i^{th}$ subject's observation at time $t_{s}$, where
\begin{equation}
\bm{Y}_{is} = \bm{X}_i(t_{is}) +\bm{\epsilon}_{is},
\label{eq:measurementerror}
\end{equation}

with i.i.d. zero mean measurement errors $\bm{\epsilon}_{is} = (\epsilon_{i1s},\dots,\epsilon_{ips})^T$ that  are independent of $\bm{X}_i$

and have component-specific variances $\sigma^2_j<\infty$ for $j=1,\dots,p$.
Estimation of the smooth component curves, $\tilde{X}_{ij}$, may be performed by any smoothing technique, e.g. kernel methods, smoothing splines, or local polynomial fitting.

As an example we consider local linear smoothing. Let $K: \mathbb{R} \to \mathbb{R}$ represent a non-negative kernel function. Applying local linear smoothing with bandwidth $b$ leads to the estimates  $\tilde{X}_{ij}(t) = \hat{\beta}^{(0)}_{ij}$,  where for all $i=1,\dots,n,~j=1,\dots,p$
\begin{equation}
(\hat{\beta}^{(0)}_{ij},\hat{\beta}^{(1)}_{ij}) = \argmin_{\beta^{(0)}_{ij},\beta^{(1)}_{ij}} \sum_{s=1}^{m}K\left(\frac{t_{s}-t}{b}\right)(Y_{ijs}-\beta^{(0)}_{ij}-\beta^{(1)}_{ij}(t-t_{s}))^2,
\label{eq:lwls}
\end{equation}

Bandwidth selection may be data-adaptive <cit.>,  but a subjective choice often suffices in practice.
After the smoothing step, estimation may be carried out by substituting $\tilde{X}_{ij}$ for $X_{ij}$ and implementing the procedure described in Sections 3.2 and 3.3.

Once all model components are estimated,
plug-in estimates of the composite distortion functions and marginal and subject-level component transport functions follow immediately,

\begin{align}
\hat{G}_{ij} &= \hat{\Psi}_j\circ \hat{H}_i, \label{eq:mixedwarp_est}\\
\hat{T}_{jk} &= \hat{\Psi}_j^{-1}\circ\hat{\Psi}_k,\label{eq:xct_est}\\
\hat{T}^{(i)}_{jk} &= \hat{G}_{ij}^{-1}\circ\hat{G}_{ik}, \quad \quad i=1,\dots,n, \quad j,k=1,\dots,p. \label{eq:subjxct_est}
\end{align}

Additionally, fitted curves based on the LTM are obtained as

\begin{align}
\hat{X}_{ij}(t) &= \hat{A}_{ij} (\hat{\lambda}\circ\hat{G}_{ij})(t)
=\hat{A}_{ij} (\hat{\lambda}\circ\hat{\Psi}_{j}\circ \hat{H}_i)(t), \,\, i=1,\dots,n, \,\, j,k=1,\dots,p.
\label{eq:recon}
\end{align}

These fits can be viewed as dimension reduction provided through the LTM, as they require only $n+p+1$ estimated functions as opposed to $np$ curves in the original data. In Section 6, we derive asymptotic results for these  estimates.

§ DATA APPLICATIONS

§.§ Zürich Growth Study

From 1954 to 1978, a longitudinal study on human growth and development was conducted at the University Children’s Hospital in Zürich. The sitting heights, arm lengths, and leg lengths of a cohort of children were measured on a dense time grid such that the resulting data can be viewed as multivariate functional data. We focus on the timing of pubertal growth spurts, which usually occur between ages 9 and 18. It is standard in the growth curve literature to work with  the derivatives of the growth curves, i.e. the growth velocities, instead of the curves themselves <cit.>. The growth velocities have a peak during puberty, with the crest location representing the age when an individual is growing fastest. The timings and curvatures of these peaks are critical in informing growth patterns. Estimates of  growth velocities for the three modalities of sitting height, leg length and arm length are displayed in Fig. <ref>, obtained via local linear smoothing. There is a well known difference in the timing of the pubertal spurts between boys and girls that is evident from Fig. <ref>; girls start the pubertal spurt earlier than boys.  It is of further interest to study body shape changes as represented by the three growth modalities.
We demonstrate in the following that the  proposed representation provided by the LTM  facilitates such studies.

Growth velocities (in cm/year) during puberty for boys (blue) and girls (red). Scaled component tempo functions are marked for boys and girls with dashed and dotted lines, respectively.

Component tempos for boys and girls are a simple way to summarize these differences (Fig. <ref>, dashed and dotted lines, respectively).

Component tempos $\gamma$ (left) and transports $\Psi$ (right) for growth modalities. The dashed line represents the tempo and transport for the latent tempo, $\lambda$.

Turning our attention to the joint time dynamics of the $p=3$ modalities, aiming to study body shape changes, we restrict our analysis to the boys for the sake of brevity.
To compare growth patterns we study  the component tempos, which are displayed for each modality in the left panel of Fig. <ref>. The dynamics of joint development can be seen  by examining the order of peaks across modalities. In this case leg length is first, followed by arm length, while sitting height lags behind. The tempos obtain similar slopes during puberty, though leg length has the most gradual spurt and sitting height has the sharpest increase, This may reflect the fact that its lagged onset results in a smaller window between the onset of its growth spurt and the end of growth which is a consequence of the ossification of the epiphyses. We note that it is  not impossible for an individual to experience some minor growth past the age of 18, but in the Zürich study such cases were rare and so this complication was ignored. The component transports displayed in Fig. <ref> (right) further illustrate the nature of the  body part tempos relative to the baseline latent system time. Remarkably, the tempo of arm length is nearly identical to the latent curve. This suggests that arm length can be used a representative modality which mirrors a child's overall growth.

The cross-component transport matrix, which characterizes pairwise time relations between components.

The estimated cross-component transports, $\hat{T}_{jk},~j,k\in\{1,\dots,p\}$ as per (<ref>), are depicted in Figure <ref> and illustrate the intercomponent time dynamics for the three growth modalities. An XCT map can be interpreted as the transport required in order to accelerate or decelerate the tempo of component $j$ to that of component $k$. Taking the transport between sitting height and leg length, $\hat{T}_{23}$ as an example, we see that the XCT map falls above the identity, indicating that sitting height's tempo must be accelerate in order to synchronize with that of leg length. This matches the interpretation that individuals tend to experience growth spurts in leg lengths before sitting height, reflected in Fig. <ref>; see also

The cross-component transport map $T_{12}$ which expresses the sitting height's timing patterns relative to the leg length's as a baseline. The peak of pubertal growth rate for the leg occurs at approximately age $13.5$, while the maximum growth velocity for sitting height growth occurs at approximately $T_{12}(13.5) \approx 14.5$ years old.

The magnitude of the XCT map's deviation from the identity shows how dissimilar the two components are. For example, sitting height and leg length are the most distinct modalities of growth among those considered here, and their XCT map exhibits the most pronounced departure from the identity. An intuitive interpretation of the map is that $T_{jk}$ expresses the $k^{th}$ component's timing patterns relative to the $j^{th}$ component's as a baseline. For example, when the leg tempo is at time $t=13.5$, the comparable time point for the sitting height tempo is approximately at $T_{jk}(13.5) \approx 14.5$ as illustrated in Fig. $\ref{fig:sit_to_leg}$.

§.§ Air Pollutants in Sacramento, CA

The study of air pollutants has been a topic of interest for atmospheric scientists and environmentalists alike for several decades. In particular, increased ground-level ozone (O$_3$) concentrations have been shown to have harmful effects on human health.  Unlike many air pollutants, surface ozone is not directly emitted by sources of air pollution (e.g. road traffic); it is formed as a result of interactions between nitrogen oxides and volatile organic compounds in the presence of sunlight <cit.>. Because of this interaction, compounds such as nitrogen dioxide are known and important precursors of increased ozone concentrations <cit.>.

The California Environmental Protection Agency has monitored hourly air pollutant concentrations at several station locations since the 1980s. Here we consider the sample of weekday trajectories of ozone (O$_3$), and nitrogen oxides (NO$_x$) concentrations during the summer of 2005 in Sacramento (Fig. <ref>). Smooth trajectories were obtained from raw data using locally linear weighted least squares. <cit.> has previously investigated a similar dataset in the context of warped functional regression, where the primary aim was to model phase variation explicitly in order to relate the timing of peak concentrations of NO$_x$ to those of O$_3$. The chemistry of the compounds as well as a visual inspection of the curves suggests that NO$_x$ concentrations tend to peak around 8 a.m., reflecting standard morning commute hours and the impact of traffic emissions  on air quality, while  ozone levels peak around 2 to 3 p.m., indicating that the synthesis mechanism induces a lag of up to approximately 6 hours.

24-hour trajectories of NO$_x$ (left) and ozone (right),  concentrations in parts per billion (ppb) on a log scale.

It is of interest  if meteorological factors affect the rate of ozone synthesis. Individual component transports combined with Fréchet regression for distributions provide a natural framework to study this query <cit.>. Subject-specific transports from NO$_x$ concentrations to ozone concentrations, $T^{(i)}_{NO_x\to O_3}$, were calculated as per (<ref>) for each day. Global Fréchet regression was then applied through the model

\begin{align}
\begin{split}
\hat{m}_\oplus(x) &= \underset{T\in\mathcal{W}}{\argmin}~M_n(T,x),\\
M_n(T,x) &= n^{-1} \sum_{i=1}^n q_{in}d^2_W(T_i,T),
\end{split}
\end{align}

where $m_\oplus$ denotes the conditional Fréchet mean of the transport given the covariate $x$, the wind speed recorded a given day. Here, $d_W$ denotes the $2-$Wasserstein distance <cit.> and the weights are defined as $q_{in}=1+(x_i-\bar{x})(x-\bar{x})/\hat{s}^2_x$  <cit.>, where $\bar{x}$ and $\hat{s}^2_x$ represent the sample mean and variance of the observed wind speeds, respectively. The model was fit using the  package . This relies on the key observation  that the transports can be represented as distribution functions if the domain is normalized to $[0,1]$ <cit.>.

Fréchet regression of NO$_x$-to-O$_3$ cross-component transports onto daily max windspeeds in knots. Windier days correspond to more linear transport functions, which suggests  O$_3$ synthesis more closely follow NO$_x$ emission. Less windy days are associated with more pronounced lags between the pollutants.

Figure $\ref{fig:freg}$ displays the observed transports and the fits obtained from Fréchet regression using windspeed as a predictor. The rainbow gradient corresponds to windspeeds ranging from 3 to 10 knots and their associated fitted transports are overlaid the original data. The regression fits suggest that days with lower windspeeds correspond with transports which are further from the diagonal, indicating an augmented lag between peak concentrations of NO$_x$ and ozone. On the other hand, days with high wind speeds have fitted transports very near the diagonal which suggests that windier settings accelerate the synthesis process. Intuitively this is reasonable, as  more wind will result in a higher rate of collisions of the particles, and thus quicker production of ozone after peak NO$_x$ emission. The Fréchet $R^2_\oplus$ value was $0.44$, which suggests that wind speed explains a considerable amount of variation in the observed transports.

§ SIMULATION STUDY

In this section we illustrate the latent transport model through a simulated $p~=~4-$ dimensional dataset which exhibits intercomponent warping of a base latent curve. Component curves are simulated on a grid $\mathcal{T} = [0,1]$ traversed by  increments of $0.05$ as

\begin{align}
X_{ij}(t_s) &=  A_{ij} (\lambda \circ \Psi_j \circ H_i \circ R_{ij})(t_s) + \epsilon_{ijs},
\end{align}

where the latent curve is defined as $\lambda(t) = \lambda_0(t)/||\lambda_0||_\infty$, for $\lambda_0(t) = 20+15t^2-5\cos(4\pi t)+3\sin(\pi t^2)$, the random amplitudes are distributed as $A_{ij}\overset{iid}{\sim} \mathcal{N}(100,\sigma^2_A)$, the component distortion functions $\Psi_j,~j=1,\dots,4$ are mixtures of Beta distributions, and $H_i$ are random distortions of the unit interval. Specifically, we set $\Psi_j(t) = \vartheta B_t(a_j, b_j) + (1-\vartheta) t$ for $j=1,2$, where $B_t$ denotes the regularized incomplete Beta function, and $a = (2,1)^T$, $b = (2,1/2)^T$,  and $\vartheta = 0.5$. Then we define $\Psi^{-1}_{j+2}(t) = 2t - \Psi^{-1}_{j}(t), ~j=1,2,$  so the net identity warp constraint $p^{-1}\sum_{j=1}^p \Psi^{-1}_j(t) = t$ is satisfied. The sources of random phase variation are introduced by the subject-level warping functions defined by  $H^{-1}_i(t) = \dfrac{\exp(tw_i)-1}{\exp(w_i)-1}$ where $w_i\overset{iid}{\sim}\mathcal{N}(0, \sigma^2_W)$ and a random nuisance distortion function defined by $R^{-1}_{ij}(t) = \dfrac{\exp(t d_{ij})-1}{\exp(d_{ij})-1}$ where $d_{ij}\overset{iid}{\sim}\mathcal{N}(0, \sigma^2_D)$. Finally, measurement error is also added to contaminate the raw observations through the random variables $\epsilon_{ijs}\overset{iid}{\sim}\mathcal{N}(0, \sigma_E^2)$. The levels of these distortions was varied by choosing  $\sigma_W=0,.5,1,$ $\sigma_D=0,.5,1,$ and $\sigma_E=0,1,5,10$. Figure <ref> displays an example of simulated data and component warping functions for the case where $\sigma_W=0.5$, $\sigma_D=0$, and $\sigma_E=1$.

Estimated (dashed) vs. true (solid) component transport functions (left) and cross-component transport functions (right) under the latent transport model at the noise setting with
$\sigma_W=0.5$, $\sigma_D=0.5$, and $\sigma_E=1$.

The component warping functions were chosen to induce both distortions which accelerate (or delay) the latent trajectory uniformly and more complicated distortions which alternate between stretching and shrinking the time domain. The former induces a simpler warping structure as the warping functions remain on either side of the diagonal. In some cases, these kinds of warps may be approximated adequately with simple shift parameters. On the other hand, a component transport which oscillating across the diagonal corresponds to a nuanced warping effect that requires added flexibility beyond that of shift-warping. The LTM satisfies this need for additional flexibility. We also provide details on the finite sample performance of our methods on the fits of both the latent curves and the realizations of the full process.

Pre-smoothing of curves was performed using local linear regression with the Epanechnikov kernel and the number of knots was set at $K=4$. We  note that the choice of this  hyperparameter did not drastically change the quality of estimates in practice unless it was set at a blatantly low value $(K\leq2)$, which then resulted in lack of fit, or raised excessively high $(K\geq7)$, which distorted the estimates of warping functions.

Estimates were obtained from fitting the model on a sample of $n=50$ $4-$dimensional processes at each noise level. This process was repeated for $B=250$ Monte Carlo runs and the results are summarized in terms of performance measures, including integrated squared errors for the estimated latent curve, component warps, and fitted processes,

\begin{align}
\begin{split}
\textit{LISE} = \int_0^1 \{\hat{\lambda}(t) -\lambda(t)\}^2dt, &~~
\textit{HMISE} = \frac{1}{n}\sum_{i=1}^n \int_0^1 \{\hat{H}_{i}(t)-H_{i}(t)\}^2dt,~  \text{and} \\
% PSISE_j &= \int_0^1 \{\hat{\Psi}_j(t) -\Psi_j(t)\}^2dt, \quad \text{and}\\
% \quad TMISE_{jk} &=  \int_0^1 \{\hat{T}_{jk}(t)-T_{jk}(t)\}^2dt, \quad \text{and}  \\
\textit{XMISE} = \frac{1}{np}\sum_{j=1}^p \sum_{i=1}^n &\int_0^1 \{\hat{X}_{ij}(t)-X_{ij}(t)\}^2dt, \quad \text{for}\quad j,k=1,\dots,p.
\label{eq:perform}
\end{split}
\end{align}

Tables $\ref{tab:sim_lambda} - \ref{tab:sim_fit}$ depict these quantities based on samples generated at each noise level.
For ease of reading, the quantities LISE and HMISE are magnified by a factor of $10^2$.
Table $\ref{tab:sim_lambda}$ shows that estimation of the latent curve $\lambda$ generally degrades with increased warping and nuisance distortion levels. Interestingly the effect of these distortions is comparable between levels when warping and nuisance  distortions are 0 and 0.5, but jumps when warping or nuisance levels  increase to 1, which reflects a severe amount of phase variation in terms of the warping distortion level.  Effects on the integrated squared error for the latent curve are similar for both sources of distortion, $\sigma_W$ and $\sigma_D$. Increases in measurement error on the other hand do not seem to affect the quality of the estimate of $\lambda$ much. A potential explanation for this is the fact that the estimator of $\lambda$ is a global average, so additive measurement errors tend to cancel each other out.

Integrated squared error of the latent curve estimates under simulated noise settings, magnified by factor of 100.

2cLISE $\times~ 10^2$   4cMeasurement Error

Warping Distortion      Nuisance Distortion    $\sigma_{E}=0$    $\sigma_{E}=1$     $\sigma_{E}=5$     $\sigma_{E}=10$

3*$\sigma_{W}=0$    $\sigma_{D}=0.0$     0.01     0.01     0.01     0.01

$\sigma_{D}=0.5$     0.00     0.00     0.01     0.01

$\sigma_{D}=1.0$    0.09     0.08     0.06     0.08

3*$\sigma_{W}=0.5$    $\sigma_{D}=0.0$    0.01     0.01     0.01     0.01

$\sigma_{D}=0.5$     0.01     0.01     0.01     0.02

$\sigma_{D}=1.0$     0.10     0.09     0.07     0.09

3*$\sigma_{W}=1$    $\sigma_{D}=0.0$    0.11     0.10     0.08     0.10

$\sigma_{D}=0.5$    0.13     0.11     0.09     0.11

$\sigma_{D}=1.0$     0.32     0.31     0.28     0.31

As for the estimates of the subject-specific warping functions $H$, we see a similar trend across both kinds of distortion (Table $\ref{tab:sim_warp}$). As $\sigma_W$ increases, the warping function targets have a greater likelihood to be more extreme warps, which are harder to estimate.  In addition the bias imposed by regularization is greater for these extreme cases. As expected, greater nuisance distortion also degrades performance.  Measurement error is associated with generally  worse performance, in contrast to its effect on the estimation of $\lambda$. Since the warping functions $H$ are subject-specific, and not averaged over the sample, they do not enjoy the robustness to measurement error that the latent curve displays.

Mean integrated squared error (magnified by $10^2$) of the subject-level warping estimates $\hat{H}_i(t)$ under simulated noise settings.

2cHMISE $\times~ 10^2$   4cMeasurement Error

Warping Distortion      Nuisance Distortion    $\sigma_{E}=0$    $\sigma_{E}=1$     $\sigma_{E}=5$     $\sigma_{E}=10$

3*$\sigma_{W}=0$    $\sigma_{D}=0.0$    0.00     0.00     0.01     0.02

$\sigma_{D}=0.5$    0.02     0.02     0.02     0.03

$\sigma_{D}=1.0$    0.15     0.15     0.12     0.10

3*$\sigma_{W}=0.5$    $\sigma_{D}=0.0$    0.00     0.00     0.01     0.02

$\sigma_{D}=0.5$     0.02     0.02     0.02     0.03

$\sigma_{D}=1.0$     0.15     0.15     0.12     0.11

3*$\sigma_{W}=1$    $\sigma_{D}=0.0$    0.12     0.12     0.13     0.21

$\sigma_{D}=0.5$    0.14     0.13     0.15     0.22

$\sigma_{D}=1.0$     0.29     0.28     0.31     0.38

Finally, performance measures for the sample of fitted curves are displayed in Table $\ref{tab:sim_fit}$. Mean integrated squared errors for fitted curves $X_{ij}$ can be thought of as a composite measure of the errors in estimates of the latent curve $\lambda$, component transports $\Psi_j$, individual warps $H_i$, and amplitude factors $A_{ij}$. Note that the amplitude factors vary substantially between components
and therefore the errors for fitted curves are much larger than that of the latent curve and warping functions.
The various distortions are seen to have a slightly different effect on XMISE as compared to  LISE and HMISE. Again, larger warping and nuisance distortions result in worse performance, however in the case of fitted curves, there is a difference between warping and nuisance distortion. Increasing the level of warping distortion $\sigma_W$ has less of an impact on XMISE than does increasing that of nuisance distortion $\sigma_D$. For example, in the case of no measurement error, increasing $\sigma_D$ from 0 to 0.5 results in a roughly 175-fold increase in XMISE, while a similar change in $\sigma_W$
increases the error by a factor of 6.5. A similar pattern occurs over all levels of measurement error. This observation suggests that increased nuisance distortion degrades the estimation of the component transports $\Psi_j$, as this effect is not present in the estimation of $\lambda$ or $H$, and the amplitude factors are invariant to distortion effects of all kinds. Measurement errors have a similar effect as for HMISE, which again may be explained by the fact that fitted curves are subject-specific.

Mean integrated squared error of the fitted curve estimates $\hat{X}_{ij}(t)$ under simulated noise settings. Note the increase in scale of the XMISE when compared to LISE and HMISE is due to the fact that curves are scaled up by the amplitude factor which varies randomly around the value 100.

2cXMISE   4cMeasurement Error

Warping Distortion      Nuisance Distortion    $\sigma_{E}=0$    $\sigma_{E}=1$     $\sigma_{E}=5$     $\sigma_{E}=10$

3*$\sigma_{W}=0$    $\sigma_{D}=0.0$    0.04     0.50     13.09     57.19

$\sigma_{D}=0.5$     6.97     7.35     20.05     63.65

$\sigma_{D}=1.0$    80.98     81.50     92.78     140.97

3*$\sigma_{W}=0.5$    $\sigma_{D}=0.0$     0.26     0.74     13.32     57.37

$\sigma_{D}=0.5$     7.11     7.54     20.14     63.81

$\sigma_{D}=1.0$     80.75     81.84     93.75     143.37

3*$\sigma_{W}=1$    $\sigma_{D}=0.0$    17.35     19.22     34.28     91.97

$\sigma_{D}=0.5$    25.48     25.20     42.68     98.38

$\sigma_{D}=1.0$    104.60     102.80     120.65     181.99

§ THEORETICAL RESULTS

Our results focus on convergence of the components of the XCT model described in (4) as the number of curves $n$ and the number of observations per curve $m$ tend to infinity. We require the following assumptions on (L) the components of the Latent Transport Model  and (S) the smoothing methodology in the presence of discrete observations with measurement error.

* The latent curve $\lambda(t) \in C^2(D)$ is a bounded function. For any non-degenerate interval $\mathcal{T}_0\subset \mathcal{T}, ~0<\int_{\mathcal{T}_0}\lambda'(t)^2dt < \infty$.

* For $j=1,\dots,p$,  $\underset{1\leq i \leq n}{\sup} ~A_{ij} = \mathcal{O}_P(1)$ and $\underset{1\leq i \leq n}{\sup} ~A^{-1}_{ij} = \mathcal{O}_P(1)$.

Assumption (L1) bounds the latent curves and its derivatives and ensures there are no flat stretches, ensuring identifiability of the model components.
(L2) ensures that the ranges of the random processes are bounded away from zero and infinity with high probability; this condition is needed for the uniform convergence of the smoothing estimate.

* The time points $t_1,\dots,t_m$, where the sequence $m$ is assumed to depend on the sample size $n$, $m=m(n)$,  correspond to a dense regular design with smooth design density $f$ with $\underset{t \in \mathcal{T}}{\inf }f(t)>0$ that generates the time points according to $t_s = F^{-1}(\frac{s-1}{m-1}), ~s=1,\dots,m,$ where $F^{-1}$ denotes the quantile function associated with $f$. The second derivative $f''$ is bounded, $\underset{t \in \mathcal{T}^\circ}{\sup}|f''(t)|<\infty$.

* The kernel function $K$ is a probability density function with support $[-1,1]$, symmetric around zero, and uniformly continuous on its support,  with $\int_{-1}^1 K^2(u)du<\infty$.

* For each $j=1,\dots,p$, the sequences $m = m(n)$ and  $b = b(n)$ satisfy (1) $0<b< \infty$, and (2) $m\to\infty$, $b\to0$, and $mb^2(\log b)^{-1}\to\infty$ as $n\to\infty$.

These assumptions guarantee the consistent estimation of $n$ curves simultaneously, as shown in the following Proposition. We observe that (S2) is for example satisfied if the bandwidth sequence is chosen such that $b=b(n) \sim m(n)^{-1/6}$.
We next establish the rate of convergence for the smoothed curves and also the normalized versions to the true underlying processes.

Under assumptions (S0$-$S2), if $E||X^{(\nu)}(t)||_\infty^2<\infty, ~\nu=0,1,2$,  we have the uniform convergence

\begin{equation}
\underset{t \in \mathcal{T}}{\sup}|\tilde{X}_{ij}(t)-X_{ij}(t)| = \mathcal{O}_P(m^{-1/3}).
\end{equation}

The rate also extends to the standardized versions $X_{ij}^* = X_{ij}/||X_{ij}||_\infty$,

\begin{equation}
\underset{t \in \mathcal{T}}{\sup}\left|\frac{\tilde{X}_{ij}(t)}{||\tilde{X}_{ij}||_\infty}-\frac{X_{ij}(t)}{||X_{ij}||_\infty}\right| = \mathcal{O}_P(m^{-1/3}).
\end{equation}

The estimators of the latent curve and component transports involve averages of the smoothing estimates over the sample of curves as $n\to\infty$. The corresponding rates of convergence will thus rely on the uniform summability of the difference between the smoothed and true curves over $n$ and we then have a uniform rate of $\tau_m=m^{-(1-\delta)/3}$ for an arbitrarily small $\delta>0$ in lieu of the above rate $m^{-1/3}$; see Lemma 1 in the Appendix.

The proposed estimators also rely on the mechanics of the pairwise warping methods, whose convergence properties have been established in a general form in <cit.> and <cit.>. Lemma 2 in the Appendix states these rates in the specific framework of the Latent Transport Model.
We are now in a position  to state our main result, which establishes rates of convergence for the estimators of the components of the Latent Transport Model as follows.

Under assumptions  (L1), (L2), and (S0$-$S2), with $\tau_m=m^{-(1-\delta)/3}$ for an arbitrarily small $\delta>0$ and penalty parameters $\eta_1, \eta_2$ as in (<ref>), (<ref>),  we have for all $i=1,\dots,n,~j=1,\dots,p$,

a. $\underset{t \in \mathcal{T}}{\sup}|\hat{H}_{i}(t) - H(t)| = \mathcal{O}_P(n^{-1/2}) + \mathcal{O}_P(\tau_m^{1/2})+ \mathcal{O}(\eta_1^{1/2})$,

b. $\underset{t \in \mathcal{T}}{\sup}|\hat{\gamma}_j(t) - \gamma_j(t)| = \mathcal{O}_P(n^{-1/2}) +  \mathcal{O}_P(\tau_m^{1/2})+ \mathcal{O}(\eta_1^{1/2}) $

c. $\underset{t \in \mathcal{T}}{\sup}|\hat{\lambda}(t) - \lambda(t)| = \mathcal{O}_P(n^{-1/2})+  \mathcal{O}_P(\tau_m^{1/2})+ \mathcal{O}(\eta_1^{1/2})$,

d. $\underset{t \in \mathcal{T}}{\sup}|\hat{\Psi}_j(t) - \Psi_j(t)| = \mathcal{O}_P(n^{-1/2}) +  \mathcal{O}_P(\tau_m^{1/2})+ \mathcal{O}(\max(\eta_1,\eta_2)^{1/2})$,

e. $\underset{t \in \mathcal{T}}{\sup}|\hat{G}_{ij}(t) - G_{ij}(t)| = \mathcal{O}_P(n^{-1/2}) +  \mathcal{O}_P(\tau_m^{1/2})+ \mathcal{O}(\max(\eta_1,\eta_2)^{1/2})$, and

f. $|\hat{A}_{ij} - A_{ij}| = \mathcal{O}_P(m^{-1/6}).$

The three terms in the rates correspond, in order, to (1) the parametric rate achieved through the standard central limit theorem, (2) the smoothing rate which is dependent on the number of observations per curve $m$, and (3) a rate due to the well-known bias introduced by the penalty parameters used in the regularization steps. Additionally, if we suppose that $m$ is bounded below by a multiple of $n^{3(1-\delta)^{-1}}$, then the rates corresponding to the smoothing steps are bounded above by $n^{-1/2}$. If we take the penalty parameters to be $\eta_1\sim\eta_2=\mathcal{O}(n^{-1})$, a $n^{-1/2}$ rate of convergence can be achieved for each of the estimators in Theorem 1 $a.$-$e.$ Otherwise if $m\sim n^{\Delta(1-\delta)^{-1}}$, for any $\Delta<3$,  the convergence is limited by the smoothing step and achieves the rate of $n^{-\Delta/6}$. This line of reasoning gives the following result for curves which are observed fully or  on sufficiently dense designs.

If the penalty parameters satisfy $\eta_1\sim\eta_2=\mathcal{O}(n^{-1})$ and the random trajectories are fully observed without error or the trajectories are recorded with at least a multiple of $m\sim n^{\Delta(1-\delta)^{-1}}$  observations per curve, with $\Delta>3$, then under the assumptions of Theorem 1, we have for all $i=1,\dots,n,~j=1,\dots,p$,

a. $\underset{t \in \mathcal{T}}{\sup}|\hat{H}_{i}(t) - H(t)| = \mathcal{O}_P(n^{-1/2})$,

b. $\underset{t \in \mathcal{T}}{\sup}|\hat{\gamma}_j(t) - \gamma_j(t)| = \mathcal{O}_P(n^{-1/2}) $

c. $\underset{t \in \mathcal{T}}{\sup}|\hat{\lambda}(t) - \lambda(t)| = \mathcal{O}_P(n^{-1/2})$,

d. $\underset{t \in \mathcal{T}}{\sup}|\hat{\Psi}_j(t) - \Psi_j(t)| = \mathcal{O}_P(n^{-1/2})$,

e. $\underset{t \in \mathcal{T}}{\sup}|\hat{G}_{ij}(t) - G_{ij}(t)| = \mathcal{O}_P(n^{-1/2})$, and

f. $|\hat{A}_{ij} - A_{ij}| = \mathcal{O}_P(n^{-1/2}).$

The asymptotic results for the cross-component transports then follow immediately from the rates established in Theorem 1.

Under assumptions of Theorem 1 for  $i=1,\dots,n,~1\leq j,k\leq p$,

$\underset{t \in \mathcal{T}}{\sup}|\hat{T}_{jk}(t)- T_{jk}(t)| = \mathcal{O}_P(n^{-1/2}) +  \mathcal{O}_P(\tau_m^{1/2})+ \mathcal{O}( \max(\eta_1,\eta_2)^{1/2}),$  and

b. $\underset{t \in \mathcal{T}}{\sup}|\hat{T}^{(i)}_{jk}(t)- T^{(i)}_{jk}(t)| = \mathcal{O}_P(n^{-1/2}) +  \mathcal{O}_P(\tau_m^{1/2})+ \mathcal{O}( \max(\eta_1,\eta_2)^{1/2}).$

A similar corollary for cross-component transports follows in the case of fully observed curves or dense enough designs.

If the penalty parameters satisfy $\eta_1\sim\eta_2=\mathcal{O}(n^{-1})$ and  the random trajectories are fully observed without error or are recorded with at least a multiple of $m\sim n^{\Delta(1-\delta)^{-1}}$  observations per curve, with $\Delta>3$, then under the assumptions of Theorem 1, we have for  $i=1,\dots,n,~1\leq j,k\leq p$,

a.	$\underset{t \in \mathcal{T}}{\sup}|\hat{T}_{jk}(t)- T_{jk}(t)| = \mathcal{O}_P(n^{-1/2}),$  and

b. $\underset{t \in \mathcal{T}}{\sup}|\hat{T}^{(i)}_{jk}(t)- T^{(i)}_{jk}(t)| = \mathcal{O}_P(n^{-1/2}).$

Corollaries 1 and 2 suggest that, on dense enough measurement schedules, parametric rates of convergence are achievable for the components of the LTM. The collection of cross-component transports also exhibit the following algebraic structure.

For any cycle of components indexed by the sequence,
$$\pi_1\to \pi_2\to \pi_3\to\dots \to \pi_L\to \pi_1,$$
with arbitrary length $L$ and $\pi_1,\dots,\pi_L\in \{1,\dots, p\}$, their respective cross-component transports satisfy

$$T_{\pi_1\pi_2} \circ T_{\pi_2\pi_3} \circ \dots \circ T_{\pi_L\pi_1} = id.$$

This result ensures that the system of cross-component transport maps is free of internal inconsistencies. For example, if for three components $A$, $B$, and $C$, the pairwise transports $T_{AB}$ and $T_{BC}$ suggest that Component $A$ tends to precede Component $B$ which itself tends to precede Component $C$, this implies that the transport $T_{AC}$ must indicate that Component $A$ tends to precede Component $C$. Furthermore, mapping a component tempo through other components and then back to itself will result in the original component tempo, unchanged.
Next we consider the convergence rates of reconstructed curves as per (<ref>), putting all model components together.

Under assumptions of Theorem 1 for  $i=1,\dots,n,~j=1,\dots ,p$,

$$\underset{t \in \mathcal{T}}{\sup}|\hat{X}_{ij}(t)-X_{ij}(t)|= \mathcal{O}_P(n^{-1/2}) +  \mathcal{O}_P(\tau_m^{1/2})+ \mathcal{O}( \max(\eta_1,\eta_2)^{1/2}).$$

Again a parametric rate is achievable on dense enough designs.

If the penalty parameters satisfy $\eta_1\sim\eta_2=\mathcal{O}(n^{-1})$ and the random trajectories are fully observed without error or the trajectories are recorded with at least a multiple of $m\sim n^{\Delta(1-\delta)^{-1}}$  observations per curve, with $\Delta>3$, then under the assumptions of Theorem 1, we have for $i=1,\dots,n,~j=1,\dots,p$,

$$\underset{t \in \mathcal{T}}{\sup}|\hat{X}_{ij}(t)-X_{ij}(t)|= \mathcal{O}_P(n^{-1/2}).$$

§ CONCLUDING REMARKS

In this paper we have introduce the Latent Transport Model (LTM) to represent and decompose multivariate functional data and to quantify their inter-component time dynamics. This model is an attractive alternative to current functional principal component based models  and provides a simple representation for multivariate functional data that exhibit  subject-specific and inter-component warping: Information concerning their time dynamics is synthesized and compressed into two fixed effect terms (the latent curve and a collection of component-level warping functions) and two random effect terms (a random amplitude vector and a collection of subject-level warping functions). This representation requires the estimation of only one random warping function and amplitude vector per subject, in addition to $p+1$ deterministic functions overall.

In some cases these components may be reduced even further. For example, when subject-level warping is negligible or part of a pre-processing step, a special case of the model arises in which time-dynamics are fully characterized by the $p+1$ fixed effect curves and one random scalar per component. Alternatively, if subject-level time warping is present but further dimension reduction is desired, transformation of warps by the LQD (or an alternative) transformation <cit.>
permits a Karhunen-Loève expansion in $\mathcal{L}^2-$space. Applying the LTM and truncating this expansion at an appropriate number of eigenfunctions, say $K_0$, creates a representation of multivariate functional data using only $p+K_0$ random scalars, as opposed to, for example, a standard FPCA representation which requires $p\times K_0$ variables.

The LTM serves both as an extension of existing univariate functional warping methods, as well as a model class  for multivariate functional data analysis and registration. Future directions of note include implementing the LTM for alternative alignment algorithms besides pairwise warping, harnessing the flexibility of cross-component transport maps for imputation of components in partially observed multivariate functional data, or relaxing structural assumptions to allow for more flexible functional relationships between different latent curves for distinct subsets of components. Spatiotemporal applications are also promising for the LTM, in which the vector function components are indexed by location, as long as repeated observations are available.  In such a situation, component warping functions may reveal time trends across geographic regions. In areas like spatiotemporal data analysis where the number of vector components is high, downstream application of dimension reduction techniques like multidimensional scaling may prove useful in comparing and understanding the mutual time dynamics across a large number of components.

§ APPENDIX

§.§ Intermediate Lemmata

The following two Lemmata provide rates of convergence for the smoothing procedure and the pairwise warping estimates which are needed in the proofs of Theorems 1 and 2.

Under the conditions of Proposition 3, for any $\delta\in(0,1)$ if $b\sim m^{-(1-\delta)/6}$ for all $i=1,\dots,n,$ and $\underset{n\to\infty}{\limsup}~nm^{-\delta/2}<\infty$, one has   for all $j=1,\dots,p$,

\begin{align}
\frac{1}{n}\sum_{i=1}^n\underset{t \in \mathcal{T}}{\sup}|\tilde{X}_{ij}(t)-X_{ij}(t)| &= \mathcal{O}_P(\tau_m)\label{lemma1a},\\
\frac{1}{n}\sum_{i=1}^n\underset{t \in \mathcal{T}}{\sup}\left|\frac{\tilde{X}_{ij}(t)}{||\tilde{X}_{ij}||_\infty}-\frac{X_{ij}(t)}{||X_{ij}||_\infty}\right| &= \mathcal{O}_P(\tau_m)\label{lemma1b},
\end{align}

where $\tau_m=m^{-(1-\delta)/3}$.

For all $j = 1,\dots,p$, under assumptions  (L1), (L2), and (S0-S2), with bandwidth $b(m) \sim m^{-1/6}$ and $\eta\to0$, the spline coefficients of the warping functions $\tilde{\theta}^{(j)}_{V_{i'i}}$ satisfy

\begin{equation*}
||\tilde{\theta}^{(j)}_{V_{i'i}} - \theta^{(j)}_{V_{i'i}}|| =  \mathcal{O}_P(m^{-1/6} ) + \mathcal{O}(\eta_1^{1/2}),
\end{equation*}
and the corresponding estimates of the pairwise warping function $\tilde{V}^{(j)}_{i'i}$ in (<ref>) satisfy

\begin{equation*}
\underset{t \in \mathcal{T}}{\sup}|\tilde{V}_{i'i}^{(j)}(t)- V_{i'i}^{(j)}(t)| =  \mathcal{O}_P(m^{-1/6}) + \mathcal{O}(\eta_1^{1/2}) .
\end{equation*}

§.§ Proofs of Theoretical Results

Throughout we use $C, C',C''$ etc. to denote generic constants.

Proof of Proposition  1.

Let $(\Omega, \mathscr{F}, P)$
be the probability space on which the observed data $Y_{ijs} = X_{ij}(t_s) + \epsilon_{ijs}$, $i=1,\dots,n,~j=1,\dots,p,~s=1,\dots,m,$ are defined, where $\Omega$ is the sample space, $\mathscr{F}$ is the $\sigma$-algebra of events, and $P: \mathscr{F}\to[0,1]$ is a probability measure. First we observe that because of the independence between $(A_{ij},H_i)$, the sources of amplitude and phase variation,  and $\epsilon_{ijs}$, the measurement errors, the probability space $(\Omega, \mathscr{F}, P)$ can be considered as a product space of two probability spaces, $(\Omega_1, \mathscr{F}_1, P_1)$ where the amplitude factors, warping functions, and hence the multivariate processes, $\mathbf{X}$,
are defined, and $(\Omega_2, \mathscr{F}_2, P_2)$ where the measurement errors are defined. Fixing an element $\omega_1\in \Omega_1$ corresponds to a realization of the amplitude factor $A$, the
warping function $H$, and the multivariate process $\mathbf{X}$. We note that the errors $\{\epsilon_{ijs}\}_{s=1}^m$
are independent in $(\Omega_2, \mathscr{F}_2, P_2)$. Hence, we may substitute the index $i$ in the subscripts of $A_{ij}, H_i, X_{ij}, Y_{ijs}$ and $\epsilon_{ijs}$ by $\omega_1$; specifically, we use $A_{\omega_1,j}, H_{\omega_1}, X_{\omega_1,j}, Y_{\omega_1,js}$ and $\epsilon_{\omega_1,js}$,
respectively, in the following.
For any fixed $\omega_1\in\Omega_1$, $\{\epsilon_{\omega_1,js}\}_{s=1}^m$ are i.i.d. realizations of $\epsilon_{\omega_1,j}$, with $Y_{\omega_1,j}(t_s) =X_{\omega_1,j}(t_s) +  \epsilon_{\omega_1,js}$.

A second observation is that

assumption (L2) guarantees that for arbitrarily small $\rho>0$, there exist constants $0<c_\rho < C_\rho<\infty$ and an event $Q_\rho:=\{c_\rho<\underset{1\leq j \leq p}{\sup}\underset{1\leq i \leq n}{\sup}~ A_{ij}(t)<C_\rho\}$ such that $P(Q_\rho)>1-\rho$. When the event  $Q_\rho$ holds, the range of $X_{ij}(t)$ is the bounded set $R^* = [c_\rho,C_\rho]$ for all $j=1,\dots,p$. We may then ignore the event $Q_\rho^C$ which happens with arbitrarily small probability.

A third observation is that for  any $t \in \mathcal{T}$, the unobserved true process $X_{\omega_1,j}$ satisfies

\begin{align}
\begin{split}
X_{\omega_1,j}(t) &= E(Y_{\omega_1,j}(t)|X_{\omega_1,j}) = \underset{x\in R^*}{\argmin} ~M^{(j)}_{\omega_1}(x, t)\\
M^{(j)}_{\omega_1}(x,t) &= E_{\Omega_2}\left[(Y_{\omega_1,j}(t)-x)^2|X_{\omega_1,j}\right]
\label{eq:xmin}
\end{split}
\end{align}

and  the target of local linear estimation is

\begin{align}
\begin{split}
\beta_{\omega_1,j}(t) &= (\beta^{(0)}_{\omega_1,j}(t), \beta^{(1)}_{\omega_1,j}(t))^T,
\end{split}
\end{align}

which minimizes

\begin{align}
\begin{split}
\check{M}^{(j)}_{\omega_1,b}(\beta_{j},t) &= E_{\Omega_2}\left[(Y_{\omega_1,j}-Z(t)\beta_{j})^TW(Y_{\omega_1,j}-Z(t)\beta_{j})|X_{\omega_1,j}\right],
\label{eq:xlocal}
\end{split}
\end{align}

where $Z(t) = \begin{bmatrix}
1& t_0-t \\
1& t_1-t\\
\vdots\\
1& t_m-t
\end{bmatrix}$,
$W = \text{diag}(K_b(t_s-t))$ for$~s=1,\dots,m$, $K_b(\cdot)=b^{-1}K(\cdot/b)$, and $\check{X}_{\omega_1,j}(t) = \beta^{(0)}_{\omega_1,j}(t)$.
Finally,  the local linear estimates can be expressed as

\begin{align}
\begin{split}
\hat{\beta}_{\omega_1,j} &= (\hat{\beta}^{(0)}_{\omega_1,j}, \hat{\beta}^{(1)}_{\omega_1,j})^T,
\end{split}
\end{align}

\begin{align}
\begin{split}
\tilde{M}^{(j)}_{\omega_1,m}(\beta_{\omega_1,j},t) &= (Y_{\omega_1,j}-Z(t)\beta_{\omega_1,j})^TW(Y_{\omega_1,j}-Z(t)\beta_{\omega_1,j}),
\label{eq:xsmooth}
\end{split}
\end{align}

with the final smoothing estimate taken as  $\tilde{X}_{\omega_1,j}(t)=\hat{\beta}^{(0)}_{\omega_1,j}$.

We also require the following conditions:

* The marginal density $f$ is bounded away from zero on $\mathcal{T}$ and the second derivative is bounded, $\underset{t \in \mathcal{T}^\circ}{\sup}|f''(t)|<\infty$.

* For all $\omega_1\in \Omega_1$ and $t \in \mathcal{T}$, the minimizers $X_{\omega_1,j}(t)$, $\beta_{\omega_1,j}$ and $\hat{\beta}_{\omega_1,j}$ exist and are unique, the last $P_{\Omega_2}$-almost surely. In addition for any $\epsilon>0$,

\begin{align*}
\underset{\omega_1\in \Omega_1, t \in \mathcal{T}}{\inf}& ~\underset{|X_{\omega_1,j}(t)-x|>\epsilon}{\inf} \left( M^{(j)}_{\omega_1}(x,t)- M^{(j)}_{\omega_1}(X_{\omega_1,j}(t),t)\right)>0,\\
\underset{b\to0}{\liminf}\underset{\omega_1\in \Omega_1, t \in \mathcal{T}}{\inf}&~\underset{||\beta_{\omega_1,j} - \beta_j||>\epsilon}{\inf} \left( \check{M}^{(j)}_{\omega_1,b}(\beta_j,t)- \check{M}^{(j)}_{\omega_1,b}(\beta_{\omega_1,j},t)\right)>0.
\end{align*}

* There exist constants $r_{1j}, r_{2j}>0, ~C_{1j},C_{2j}>0$, such that for $j=1,\dots,p$,

\begin{equation*}
\underset{\omega_1\in \Omega_1,t \in \mathcal{T}}{\inf}~ \underset{|X_{\omega_1,j}(t)-x|<r_{1j}}{\inf} ~\left[ M^{(j)}_{\omega_1}(x, t)- M^{(j)}_{\omega_1}(X_{\omega_1,j}(t), t)  - C_{1j}| x-X_{\omega_1,j}(t)|^2\right]>0,
\end{equation*}

\begin{equation*}
\hspace{-0.5cm}\underset{b\to 0}{\liminf }~\underset{\omega_1\in \Omega_1,t \in \mathcal{T}}{\inf}~  \underset{||\beta_{\omega_1,j} - \beta||<r_{2j}}{\inf} ~\left[ \check{M}^{(j)}_{\omega_1,b}(\beta_j, t)- \check{M}^{(j)}_{\omega_1,b}(X_{\omega_1,j}(t), t)  - C_{2j}||\beta_j- \beta_{\omega_1,j}||^2 \right]>0.
\end{equation*}

* Let $B_r(X_{\omega_1,j}(t))\subset R_C$ be a ball of radius $r$ centered at $X_{\omega_1,j}(t)$ and denote its covering number using balls of radius $\epsilon$ by  $N(r\epsilon, B_r(X_{\omega_1,j}(t)),||\cdot||_\infty)$. Then

\begin{align*}
\underset{r>0}{\sup}        \underset{\omega_1\in \Omega_1}{\sup} \int_{0}^1 \underset{t \in \mathcal{T}}{\sup} \sqrt{1+\log N(r\epsilon, B_r(X_{\omega_1,j}(t)),||\cdot||_\infty)} d\epsilon < \infty .
\end{align*}

Note that  (i) is satisfied due to  assumption (S0). For (ii) and (iii),
observe that $M^{(j)}_{\omega_1}(x,t)$ and  $\check{M}^{(j)}_{\omega_1,b}(x,t)$ are concave for all $t \in \mathcal{T}$ and   $\tilde{M}^{(j)}_{\omega_1,m}(x,t)$  is $P_{\Omega_2}$-almost surely concave for all $t \in \mathcal{T}$, so their respective minimizers exist and are unique. Next we verify that the minimum difference in objective functions outside of an $\epsilon-$neighborhood of the minimizer is positive. By independence of $X_{\omega_1,j}$ and $\varepsilon_{\omega_1,j}$, we have

\begin{align*}
M^{(j)}_{\omega_1}(x,t) &= E_{\Omega_2}[(Y_{\omega_1,j,b}(t)-x)^2|X_{\omega_1,j}(t)]\\
&= E_{\Omega_2}[(X_{\omega_1,j}(t)-x+\varepsilon_{\omega_1,j})^2|X_{\omega_1,j,b}(t)]\\
&= E_{\Omega_2}[(X_{\omega_1,j}(t)-x)^2|X_{\omega_1,j}(t)]+E_{\Omega_2}[\varepsilon_{\omega_1,j}^2|X_{\omega_1,j}(t)]\\
&= (X_{\omega_1,j}(t)-x)^2 + \sigma^2_j
\end{align*}

\begin{align*}
M^{(j)}_{\omega_1}(X_{\omega_1,j}(t),t) &= E_{\Omega_2}[(Y_{\omega_1,j}(t)-X_{\omega_1,j}(t))^2|X_{\omega_1,j}(t)]\\
&= E_{\Omega_2}[\varepsilon_{\omega_1,j}^2|X_{\omega_1,j}(t)]= \sigma^2_j.
\end{align*}

Together these imply

\begin{align}
M^{(j)}_{\omega_1}(x,t)- M^{(j)}_{\omega_1}(X_{\omega_1,j}(t),t) &= (X_{\omega_1,j}(t)-x)^2
\label{curv1}
\end{align}

and so condition (ii) is satisfied as

\begin{align*}
% \underset{\omega_1\in \Omega_1, t \in \mathcal{T}}{\inf} ~\underset{|X_{\omega_1,j}(t)-x|>\varepsilon}{\inf} \left( M^{(j)}_{\omega_1}(x,t)- M^{(j)}_{\omega_1}(X_{\omega_1,j}(t),t)\right)
% &=
\underset{\omega_1\in \Omega_1, t \in \mathcal{T}}{\inf} ~\underset{|X_{\omega_1,j}(t)-x|>\epsilon}{\inf} (X_{\omega_1,j}(t)-x)^2> \epsilon^2 >0.
\end{align*}

Similarly, for the localized target

\begin{align*}
\check{M}^{(j)}_{\omega_1,b}(\beta_j,t) &= E_{\Omega_2}\left[(Y_{\omega_1,j}-Z(t)\beta_{j})^TW(Y_{\omega_1,j}-Z(t)\beta_{j})|X_{\omega_1,j}\right]\\
&=E_{\Omega_2}\left[\sum_{s=1}^m K(\frac{t_s-t}{b})(Y_{\omega_1,js}-\beta_0-\beta_1(t-t_s))^2|X_{\omega_1,j}\right]\\
&=E_{\Omega_2}\left[\sum_{s=1}^m K(\frac{t_s-t}{b})(X_{\omega_1,js}-\beta_0-\beta_1(t-t_s))^2|X_{\omega_1,j}\right] \\
& \quad\quad\quad\quad+ E_{\Omega_2}\left[\sum_{s=1}^m K(\frac{t_s-t}{b})\varepsilon_{\omega_1,js}^2|X_{\omega_1,j}\right] \\
&=(X_{\omega_1,j}-Z(t)\beta_j)^TW(X_{\omega_1,j}-Z(t)\beta_j) + \sigma_j^2\sum_{s=1}^m K(\frac{t_s-t}{b}),
\end{align*}

\begin{align*}
\check{M}^{(j)}_{\omega_1,b}(\beta_{\omega_1,j},t) &= E_{\Omega_2}\left[(Y_{\omega_1,j}-Z(t)\beta_{\omega_1,j})^TW(Y_{\omega_1,j}-Z(t)\beta_{\omega_1,j})|X_{\omega_1,j}\right]\\
% &=E_{\Omega_2}\left[\sum_{s=1}^m K(\frac{t_s-t}{b})(Y_{\omega_1,js}-\beta_0-\beta_1(t-t_s))^2|X_{\omega_1,j}\right]\\
&=E_{\Omega_2}\left[\sum_{s=1}^m K(\frac{t_s-t}{b})(X_{\omega_1,js}-\beta^{(0)}_{\omega_1,j}-\beta^{(1)}_{\omega_1,j}(t-t_s))^2|X_{\omega_1,j}\right] \\
& \quad\quad\quad\quad+ E_{\Omega_2}\left[\sum_{s=1}^m K(\frac{t_s-t}{b})\varepsilon_{\omega_1,js}^2|X_{\omega_1,j}\right] \\
&=(X_{\omega_1,j}-Z(t)\beta_{\omega_1,j})^TW(X_{\omega_1,j}-Z(t)\beta_{\omega_1,j}) + \sigma_j^2\sum_{s=1}^m  K(\frac{t_s-t}{b}).
\end{align*}

Together these lead to

\begin{align}
\check{M}^{(j)}_{\omega_1,b}(\beta,t)- \check{M}^{(j)}_{\omega_1,b}(\beta_{\omega_1,j},t)&=[Z(t)(\beta_{\omega_1,j}-\beta_j)]^TW[Z(t)(\beta_{\omega_1,j}-\beta_j)]\\
&= ||W^{1/2}Z(t)(\beta_{\omega_1,j}-\beta)||^2,
\label{curv2}
\end{align}

and thus

\begin{align*}
\underset{b\to0}{\liminf}&\underset{\omega_1\in \Omega_1, t \in \mathcal{T}}{\inf}\underset{||\beta_{\omega_1,j} - \beta||>\epsilon}{\inf}\check{M}^{(j)}_{\omega_1,b}(\beta,t)- \check{M}^{(j)}_{\omega_1,b}(\beta_{\omega_1,j},t)=\\
& \underset{b\to0}{\liminf}\underset{\omega_1\in \Omega_1, t \in \mathcal{T}}{\inf}\underset{||\beta_{\omega_1,j} - \beta||>\epsilon}{\inf} ||W^{1/2}Z(t)(\beta_{\omega_1,j}-\beta_j)||^2 > \underset{b\to0}{\liminf}||W^{1/2}Z(t)||^2 \epsilon^2 >0.
\end{align*}

so condition (ii) is met. Equations $(\ref{curv1})$ and $(\ref{curv2})$ also imply that the inequality in (iii) holds for $C_{1j}=1$, and $C_{2j}=||W^{1/2}Z(t)||^2$.

Finally, condition (iv) is met since by Lemma 2.7.8 of <cit.>, there exists some constant $C>0$,

\begin{align*}
\log N(r\epsilon, B_r(X_{\omega_1,j}(t)),||\cdot||_\infty)\leq C\epsilon^{-1}
\end{align*}

where $C$ depends only on the bounded set of $\mathbb{R}$ on which $X_{\omega_1,j}$ takes values and is thus uniform over $\omega_1\in\Omega_1$. Then, for any $r>0$, we may bound the integral by

\begin{align*}
\underset{\omega_1\in \Omega_1}{\sup} &\int_{0}^1 \underset{t \in \mathcal{T}}{\sup} \sqrt{1+\log N(r\epsilon, B_r(X_{\omega_1,j}(t)),||\cdot||_\infty)} d\epsilon\\
&= \underset{\omega_1\in \Omega_1}{\sup} \int_{0}^1 \sqrt{1+C\epsilon^{-1}} d\epsilon \leq 1+2\sqrt{C}<\infty,
\end{align*}

whence (iv) holds. This concludes the proof.

Proof of Proposition 2.

Note that

\begin{align*}
\underset{t \in \mathcal{T}}{\sup}\left|\frac{\tilde{X}_{ij}(t)}{||\tilde{X}_{ij}||_\infty} -  \frac{X_{ij}(t)}{||X_{ij}||_\infty}\right|
&\leq\underset{t \in \mathcal{T}}{\sup}\left|\frac{\tilde{X}_{ij}(t)}{||\tilde{X}_{ij}||_\infty} - \frac{\tilde{X}_{ij}(t)}{||X_{ij}||_\infty}\right| + \underset{t \in \mathcal{T}}{\sup}\left|\frac{\tilde{X}_{ij}(t)}{||X_{ij}||_\infty} -  \frac{X_{ij}(t)}{||X_{ij}||_\infty}\right|\\
&:= I + II,
\end{align*}

where the difference in sup-norms can be bounded by

\begin{align*}
\left| ||\tilde{X}_{ij}||_\infty-|| X_{ij}||_\infty \right|  \leq ||\tilde{X}_{ij}-X_{ij}||_\infty
\end{align*}

$$|\hat{A}_{ij}- A_{ij}| = \left| ||\tilde{X}_{ij}||_\infty- ||X_{ij}||_\infty \right|= \mathcal{O}_P\left(m^{-1/3}\right).$$

Using the fact that $A_{ij}^{-1}=1/||X_{ij}(t)||_\infty = \mathcal{O}_P(1)$,

\begin{align*}
I  = \underset{t \in \mathcal{T}}{\sup}\left|
\frac{(||X_{ij}||_\infty - ||\tilde{X}_{ij}||_\infty) \tilde{X}_{ij}(t)}{||X_{ij}||_\infty ||\tilde{X}_{ij}||_\infty}\right|
\leq ~ \frac{\left|||X_{ij}||_\infty - ||\tilde{X}_{ij}||_\infty\right|}{||X_{ij}||_\infty}
= \mathcal{O}_P\left(m^{-1/3}\right),
\end{align*}

\begin{align*}
II=\underset{t \in \mathcal{T}}{\sup}\left|\frac{\tilde{X}_{ij}(t)}{||X_{ij}||_\infty} -  \frac{X_{ij}(t)}{||X_{ij}||_\infty}\right|
\leq ||\tilde{X}_{ij}(t)-  X_{ij}(t)||_\infty
= \mathcal{O}_P\left(m^{-1/3}\right).
\end{align*}

Proof of Lemma 1.

Under  notations as  in the proof of Proposition 1, where $\tilde{X}_{ij}$ is as per $(\ref{eq:xsmooth})$ and denotes the smoothing estimator  and $\check{X}_{ij}$ is as per $(\ref{eq:xlocal})$ and denotes the local target of $\tilde{X}_{ij}$, we bound the sup-norm between the smoothed function and the true curve by

\begin{equation*}
n^{-1}\sum_{i=1}^n\underset{t \in \mathcal{T}}{\sup} |X_{ij}(t)-\tilde{X}_{ij}(t)|\leq n^{-1}\sum_{i=1}^n\underset{t \in \mathcal{T}}{\sup} |X_{ij}(t)-\check{X}_{ij}(t)|+n^{-1}\sum_{i=1}^n\underset{t \in \mathcal{T}}{\sup} |
\tilde{X}_{ij}(t)-\check{X}_{ij}(t)|.
\end{equation*}

For the first term, corresponding to the bias part, observing that

\begin{equation*}
\underset{t \in \mathcal{T}}{\sup} |X_{ij}(t)-\check{X}_{ij}(t)|\leq \underset{1\leq i \leq n}{\sup} ~\underset{t \in \mathcal{T}}{\sup} |X_{ij}(t)-\check{X}_{ij}(t)|,
\end{equation*}

we apply Theorem 2 of <cit.> to conclude

\begin{align}
\underset{1\leq i\leq n}{\sup}~ \underset{t \in \mathcal{T}}{\sup} |X_{ij}(t)-\check{X}_{ij}(t)| = \mathcal{O}(b^2),
\label{bias}
\end{align}

and thus

\begin{align*}
n^{-1}\sum_{i=1}^n \underset{t \in \mathcal{T}}{\sup} |X_{ij}(t)-\check{X}_{ij}(t)| =
\underset{1\leq i \leq n}{\sup}~\underset{t \in \mathcal{T}}{\sup} |X_{ij}(t)-\check{X}_{ij}(t)| = \mathcal{O}(b^2).
\end{align*}

For the second term, corresponding to the stochastic part, letting $a_{m}= b\sqrt{m}$,

\begin{align*}
\underset{n\to\infty}{\limsup}  ~&P\left(n^{-1}\sum_{i=1}^n  \underset{t \in \mathcal{T}}{\sup} |\check{X}_{ij}(t) -\tilde{X}_{ij}(t)|>C
a_{m} m^{-\delta/2})^{-1}
\right)\\
\leq
&~\underset{n\to\infty}{\limsup}  \sum_{i=1}^n P\left( \underset{t \in \mathcal{T}}{\sup} |\check{X}_{ij}(t) -\tilde{X}_{ij}(t)|>C
a_{m} m^{-\delta/2})^{-1}
\right)\\
\leq &~\underset{n\to\infty}{\limsup}\sum_{i=1}^n \underset{\omega_1\in \Omega_1}{\sup} P_{\Omega_2}\left(a_{m}m^{-\delta/2} \underset{t \in \mathcal{T}}{\sup} |\check{X}_{ij}(t) -\tilde{X}_{ij}(t)|>C
\right),
\end{align*}

so that the stochastic term will converge at the rate $a_m^{-1} m^{\delta/2}= b^{-1}m^{(\delta-1)/2}$ if we show that

\begin{align}
\underset{n\to\infty}{\limsup}\sum_{i=1}^n \underset{\omega_1\in \Omega_1}{\sup}& P_{\Omega_2}\left(a_{m}m^{-\delta/2} \underset{t \in \mathcal{T}}{\sup} |\check{X}_{ij}(t) -\tilde{X}_{ij}(t)|>C
\right) \to 0,\quad  \text{as} \quad C\to\infty.
\label{var2}
\end{align}

Adopting arguments in <cit.> with $\beta_1=\beta_2=2$ results in the following bound for any $q\in \mathbb{N}^+$,

\begin{align}
\underset{n\to\infty}{\limsup}\sum_{i=1}^n \underset{\omega_1\in \Omega_1}{\sup}& P_{\Omega_2}\left(a_{m}m^{-\delta/2} \underset{t \in \mathcal{T}}{\sup} |\check{X}_{ij}(t) -\tilde{X}_{ij}(t)|>2^{q}
\right) \nonumber \\
&\leq 4C \sum_{k>q} 2^{-k} ~\underset{n\to\infty}{\limsup}~ m^{-\delta/2} +
4C \sum_{k>q} 2^{-k}~ \underset{n\to\infty}{\limsup} ~nm^{-\delta/2}.
\end{align}

Thus the probability tends to 0 as $q\to\infty$ uniformly over $\omega_1\in\Omega_1$ if $\underset{n\to\infty}{\limsup}~nm^{-\delta/2}<\infty$, whence $(\ref{var2})$ follows.
Then with $b\sim m^{-(1-\delta)/6}$, we have  $b^2\sim b^{-1}m^{(\delta-1)/2}\sim m^{-(1-\delta)/3}$ and (<ref>) follows.  Following the same arguments,  we infer (<ref>).

Proof of Lemma 2.

Without loss of generality, consider a fixed component $j$.  Define the unpenalized cost function
\begin{equation*}
\mathcal{C}_{\gamma_j}(\theta,H_{i'},H_i) =  \int_\mathcal{T} d^2\left[(\gamma_j\circ H_{i'})[\theta^T\alpha(t)],(\gamma_j\circ H_{i})(t)\right]dt ,\quad j=1,\dots,p,
\end{equation*}

where in our case $d=d(f,g)=(f-g)(t)$. Then the arguments of Theorem 3 of <cit.> imply that, for some positive constants $C,C',C''<\infty,$

\begin{align*}
||\tilde{\theta}^{(j)}_{V_{i'i}} - \theta^{(j)}_{V_{i'i}}|| &\leq C
\underset{\theta\in\Theta}{\sup}|\mathcal{C}_{\eta}(\theta,H_i',H_i) - \mathcal{C}_{\gamma_j}(\theta,H_i',H_i)|^{1/2}\\
&\leq C' \left[\underset{t \in \mathcal{T}}{\sup}|\tilde{X}_{ij}(t)-X_{ij}(t)|^{1/2} + \underset{t \in \mathcal{T}}{\sup}|\tilde{X}_{i'j}(t)-X_{i'j}(t)|^{1/2}\right] + C''\eta_1^{1/2}\\
&= \mathcal{O}_P(m^{-1/6})+\mathcal{O}(\eta_1^{1/2}).
\end{align*}

Since $\sup_{t \in \mathcal{T}}|\alpha_\ell(t)|\leq1$ for all $\ell = 1,\dots,L+1$, it follows that

\begin{align}
\underset{t \in \mathcal{T}}{\sup}|\tilde{V}_{i'i}^{(j)}(t)- V_{i'i}^{(j)}(t)| &=  \underset{t \in \mathcal{T}}{\sup}|(\tilde{\theta}^{(j)}_{V_{i'i}} - \theta^{(j)}_{V_{i'i}})^T\alpha(t)| \leq ||\tilde{\theta}^{(j)}_{V_{i'i}} - \theta^{(j)}_{V_{i'i}}||  \nonumber \\
&= \mathcal{O}_P(m^{-1/6})+\mathcal{O}(\eta_1^{1/2}).
\end{align}

Proof of Theorem 1.

A. From the definition of the $\hat{H}_i$,

\begin{align*}
\underset{t \in \mathcal{T}}{\sup}|\hat{H}_i(t) - H_i(t)| &=
\underset{t \in \mathcal{T}}{\sup}|p^{-1}\sum_{j=1}^p (\tilde{H}^{(j)}_i(t) - H_i(t))|\\
&\leq p^{-1}\sum_{j=1}^p \underset{t \in \mathcal{T}}{\sup}|\tilde{H}^{(j)}_i(t) - H_i(t)|\\
&\leq   \sup_{1\leq j\leq p} \underset{t \in \mathcal{T}}{\sup}|\tilde{H}^{(j)}_i(t) - H_i(t)|.
\end{align*}

Then for each $j=1,\dots,p$,

\begin{align*}
\underset{t \in \mathcal{T}}{\sup}|\tilde{H}^{(j)}_i(t) - H_i(t)| &\leq \frac{1}{n}\sum_{i'=1}^n \underset{t \in \mathcal{T}}{\sup}|\tilde{V}^{(j)}_{i'i}(t) - V^{(j)}_{i'i}(t)| +  \underset{t \in \mathcal{T}}{\sup}\left|\frac{1}{n}\sum_{i'=1}^n V^{(j)}_{i'i}(t) - H_i(t)\right|.
\end{align*}

For the first term, an application of Lemmas 1 and 2 gives,

\begin{align*}
\frac{1}{n}\sum_{i'=1}^n \underset{t \in \mathcal{T}}{\sup}&|\tilde{V}^{(j)}_{i'i}(t) - V^{(j)}_{i'i}(t)| \leq \frac{1}{n}\sum_{i'=1}^n ||\tilde{\theta}^{(j)}_{i'i} - \theta^{(j)}_{i'i}||\\
&\leq C \bigg[\underset{t \in \mathcal{T}}{\sup}|\tilde{X}_{ij}(t)-X_{ij}(t)|^{1/2} +\frac{1}{n}\sum_{i'=1}^n\underset{t \in \mathcal{T}}{\sup}|\tilde{X}_{i'j}(t)-X_{i'j}(t)|^{1/2}+\eta_1^{1/2} \bigg]\\
&= \mathcal{O}_P(m^{-(1-\delta)/6})+\mathcal{O}(\eta_1^{1/2}).
\end{align*}

For the second term, by Theorem 2.7.5 of <cit.>,

\begin{align*}
\quad\quad\quad\quad \underset{t \in \mathcal{T}}{\sup}\left|\frac{1}{n}\sum_{i'=1}^n V^{(j)}_{i'i}(t) - H_i(t)\right|
&= \underset{t \in \mathcal{T}}{\sup}\left|\frac{1}{n}\sum_{i'=1}^n H^{-1}_{i'}(t) - t\right| \\
&= \underset{t \in \mathcal{T}}{\sup}\left|\frac{1}{n}\sum_{i'=1}^n H^{-1}_{i'}(t) - EH^{-1}_{i'}(t)\right| = \mathcal{O}_P(n^{-1/2}),
\end{align*}

and therefore

\begin{align*}
\underset{t \in \mathcal{T}}{\sup}|\hat{H}_i(t) - H_i(t)|
&\leq  p \sup_{1\leq j\leq p} \underset{t \in \mathcal{T}}{\sup}|\tilde{H}^{(j)}_i(t) - H_i(t)|\\
&=  \mathcal{O}_P(n^{-1/2})+ \mathcal{O}_P(m^{-(1-\delta)/6})+ \mathcal{O}(\eta_1^{1/2}) , ~i = 1,\dots,n.
\end{align*}

B. For all $j = 1,\dots,p$, we have

\begin{align*}
\underset{t \in \mathcal{T}}{\sup} |\hat{\gamma}_j(t) - {\gamma}_j(t)| &=     \underset{t \in \mathcal{T}}{\sup} \left| \frac{1}{n}\sum_{i=1}^n \frac{\tilde{X}_{ij}\circ\hat{H}^{-1}_i(t)}{||\tilde{X}_{ij}||_\infty} - \gamma(t)\right|\\
\underset{t \in \mathcal{T}}{\sup} \left| \frac{1}{n}\sum_{i=1}^n  \left[\frac{\tilde{X}_{ij}\circ\hat{H}^{-1}_i(t)}{||\tilde{X}_{ij}||_\infty} - \frac{X_{ij}\circ\hat{H}^{-1}_i(t)}{||X_{ij}||_\infty}\right] \right| \\
&\hspace{0.5cm} +
\underset{t \in \mathcal{T}}{\sup} \left| \frac{1}{n}\sum_{i=1}^n \frac{X_{ij}\circ\hat{H}^{-1}_i(t)}{||X_{ij}||_\infty} - \gamma(t)\right|,
\end{align*}

\begin{align*}
\underset{t \in \mathcal{T}}{\sup} \left| \frac{1}{n}\sum_{i=1}^n  \left[\frac{\tilde{X}_{ij}\circ\hat{H}^{-1}_i(t)}{||\tilde{X}_{ij}||_\infty} - \frac{X_{ij}\circ\hat{H}^{-1}_i(t)}{||X_{ij}||_\infty}\right] \right|&\leq
\frac{1}{n}\sum_{i=1}^n  \underset{t \in \mathcal{T}}{\sup} \left|\frac{\tilde{X}_{ij}\circ\hat{H}^{-1}_i(t)}{||\tilde{X}_{ij}||_\infty} - \frac{X_{ij}\circ\hat{H}^{-1}_i(t)}{||X_{ij}||_\infty}
\right|\\
&= \frac{1}{n}\sum_{i=1}^n  \underset{t \in \mathcal{T}}{\sup} \left|\frac{\tilde{X}_{ij}(t)}{||\tilde{X}_{ij}||_\infty} - \frac{X_{ij}(t)}{||X_{ij}||_\infty}
\right|\\
\end{align*}

by Lemma 1. For the second term,

\begin{align*}
\underset{t \in \mathcal{T}}{\sup} \left| \frac{1}{n}\sum_{i=1}^n \frac{X_{ij}\circ\hat{H}^{-1}_i(t)}{||X_{ij}||_\infty} - \gamma_j(t)\right|&=\underset{t \in \mathcal{T}}{\sup} \left| \frac{1}{n}\sum_{i=1}^n
\left[
\frac{A_{ij}\lambda\circ\Psi_j\circ H_i\circ\hat{H}^{-1}_i(t)}{A_{ij}} - (\lambda\circ\Psi_j)(t)
\right]\right|\\
&=\underset{t \in \mathcal{T}}{\sup} \left| \frac{1}{n}\sum_{i=1}^n  \left[(\lambda\circ\Psi_j)(H_i(t)) - (\lambda\circ\Psi_j)(\hat{H}_i(t))
\right]\right|\\
& \leq \frac{C}{n}\sum_{i=1}^n~
\underset{t \in \mathcal{T}}{\sup}|\hat{H}_i(t)-H_i(t)| \\
&\leq  C' \underset{1\leq j \leq p}{\sup}~ \frac{1}{n}\sum_{i=1}^n~
\underset{t \in \mathcal{T}}{\sup}|\tilde{H}^{(j)}_i(t)-H_i(t)|.
\end{align*}

\begin{align*}
\frac{1}{n}\sum_{i=1}^n~
\underset{t \in \mathcal{T}}{\sup}|\tilde{H}^{(j)}_i(t)-H_i(t)|&\leq  \frac{1}{n}\sum_{i=1}^n\bigg[ \frac{1}{n}\sum_{i'=1}^n \underset{t \in \mathcal{T}}{\sup}|\tilde{V}^{(j)}_{i'i}(t) - V^{(j)}_{i'i}(t)|
\\
&\hspace{2.5cm}+\underset{t \in \mathcal{T}}{\sup}\bigg|\frac{1}{n}\sum_{i'=1}^n V^{(j)}_{i'i}(t) - H_i(t)\bigg|\bigg]\\
&= \frac{1}{n^2}  \sum_{i=1}^n\sum_{i'=1}^n ||\tilde{\theta}^{(j)}_{i'i} - \theta^{(j)}_{i'i} || + \underset{t \in \mathcal{T}}{\sup}\left|\frac{1}{n}\sum_{i'=1}^n H^{-1}_{i'}(t) - t\right| \\
&\leq C\frac{1}{n}\sum_{i=1}^n \underset{t \in \mathcal{T}}{\sup}|X_{ij}(t)-\tilde{X}_{ij}(t)|^{1/2} \\
&\quad + C\frac{1}{n}\sum_{i'=1}^n \underset{t \in \mathcal{T}}{\sup}|X_{i'j}(t)-\tilde{X}_{i'j}(t)|^{1/2} +C\eta_1^{1/2}\\
&\quad  + \underset{t \in \mathcal{T}}{\sup}\left|\frac{1}{n}\sum_{i'=1}^n H^{-1}_{i'}(t) - t\right|\\
&=\mathcal{O}_P(n^{-1/2}) + \mathcal{O}_P(m^{-(1-\delta)/6}) + \mathcal{O}(\eta_1^{1/2}),
\end{align*}

where the final rates come from Lemma 1 and Theorem 2.7.5 of <cit.>, we may conclude

\begin{align*}
\underset{t \in \mathcal{T}}{\sup} |\hat{\gamma}_j(t) - {\gamma}_j(t)| &=\mathcal{O}_P(m^{-(1-\delta)/3})+ \mathcal{O}_P(n^{-1/2})+ \mathcal{O}_P(m^{-(1-\delta)/6})+\mathcal{O}(\eta_1^{1/2})\\
&= \mathcal{O}_P(n^{-1/2})+ \mathcal{O}_P(m^{-(1-\delta)/6})+\mathcal{O}(\eta_1^{1/2}).
\end{align*}

C.   The argument for the latent curve mirrors the proof of part B.,

\begin{align*}
&\underset{t \in \mathcal{T}}{\sup} |\hat{\lambda}(t) - {\lambda}(t)| = \frac{1}{n}\sum_{i = 1}^{n} \left| \frac{\tilde{Z}_{i}\circ\hat{D}^{-1}_{i}(t)}{||\tilde{Z}_{i}||_\infty} - \lambda(t) \right|\\
\underset{t \in \mathcal{T}}{\sup} \left| \frac{1}{n}\sum_{i = 1}^{n} \left[\frac{\tilde{Z}_{i}\circ\hat{D}^{-1}_{i}(t)}{||\tilde{Z}_{i}||_\infty} - \frac{Z_{i}\circ\hat{D}^{-1}_{i}(t)}{||Z_{i}||_\infty}\right] \right| +
\underset{t \in \mathcal{T}}{\sup} \left| \frac{1}{n}\sum_{i = 1}^{n} \frac{Z_{i}\circ\hat{D}^{-1}_{i}(t)}{||Z_{i}||_\infty} - \lambda(t)\right|.
\end{align*}

Defining $\bm{1}_{ij}$ as the indicator that event $I_{ij}$ occurs (that is, the $j^{th}$ curve is selected as the representative curve for the $i^{th}$ observation), the first term relies on the uniform rate of the smoother,

\begin{align*}
\underset{t \in \mathcal{T}}{\sup} \left| \frac{1}{n}\sum_{i = 1}^{n} \left[\frac{\tilde{Z}_{i}\circ\hat{D}^{-1}_{i}(t)}{||\tilde{Z}_{i}||_\infty} - \frac{Z_{i}\circ\hat{D}^{-1}_{i}(t)}{||Z_{i}||_\infty}\right] \right|&\leq
\frac{1}{n}\sum_{i = 1}^{n} \underset{t \in \mathcal{T}}{\sup} \left|\frac{\tilde{Z}_{i}\circ\hat{D}^{-1}_{i}(t)}{||\tilde{Z}_{i}||_\infty} - \frac{Z_{i}\circ\hat{D}^{-1}_{i}(t)}{||Z_{i}||_\infty}
\right|\\
\frac{1}{n} \sum_{j = 1}^{p} \sum_{i = 1}^{n}  \bm{1}_{ij}\left(\underset{t \in \mathcal{T}}{\sup} \left|\frac{\tilde{X}_{ij}(t)}{||\tilde{X}_{ij}||_\infty} - \frac{X_{ij}(t)}{||X_{ij}||_\infty}
\right|\right)\\
\frac{p}{n}\underset{1\leq j \leq p}{\sup} \sum_{i = 1}^{n} \underset{t \in \mathcal{T}}{\sup} \left|\frac{\tilde{X}_{ij}(t)}{||\tilde{X}_{ij}||_\infty} - \frac{X_{ij}(t)}{||X_{ij}||_\infty}
\right|\\
\end{align*}

where the rate is implied by Lemma 1.
For the  second term, sing the Lipschitz continuity of $\lambda$ and applying Lemma 2 and arguments used in the proof of Theorem 1,

\begin{align*}
\underset{t \in \mathcal{T}}{\sup} \left| \frac{1}{n}\sum_{i = 1}^{n} \frac{Z_{i}\circ\hat{D}^{-1}_{i}(t)}{||X_{ij}||_\infty} - \lambda(t)\right|&=\underset{t \in \mathcal{T}}{\sup} \left| \frac{1}{n}\sum_{i = 1}^{n}
\left[
\frac{A_{ij}\lambda\circ D_{i}\circ\hat{D}^{-1}_{i}(t)}{A_{ij}} - \lambda(t)
\right]\right|\\
&=\underset{t \in \mathcal{T}}{\sup} \left| \frac{1}{n}\sum_{i = 1}^{n} \left[\lambda(D_{i}(t)) - \lambda(\hat{D}_{i}(t))
\right]\right|\\
& \leq \underset{1\leq j \leq p}{\sup}~ \frac{C}{n}\sum_{i = 1}^{n}
\underset{t \in \mathcal{T}}{\sup}|\hat{D}_{i}(t)-D_{i}(t)| \quad\text{}\\
&= \mathcal{O}_P(n^{-1/2})+ \mathcal{O}_P(m^{-(1-\delta)/6})+\mathcal{O}(\eta_1^{1/2}),
\end{align*}

\begin{align*}
\underset{t \in \mathcal{T}}{\sup} |\hat{\lambda}(t) - \lambda(t)| %=\mathcal{O}_P(m^{-(1-\delta)/3})+ \mathcal{O}_P(n^{-1/2})+ \mathcal{O}_P(m^{-(1-\delta)/6})+\mathcal{O}(\eta_1^{1/2})
= \mathcal{O}_P(n^{-1/2})+ \mathcal{O}_P(m^{-(1-\delta)/6})+\mathcal{O}(\eta_1^{1/2}).
\end{align*}

D.  Using  the unpenalized objective function for the component warps

\begin{equation}
\mathcal{C}^*(\theta;\gamma_j,\lambda) = \int_\mathcal{T} (\gamma_j(t)-\lambda(\theta^T\alpha(t))^2dt,
\end{equation}

arguments similar to those in  Theorem 3 of <cit.> imply
that there exists a constant $C_0>0$ such that

\begin{align*}
||\tilde{\theta}_{\Psi_{j}} - \theta_{\Psi_{j}}||^2 &\leq C_0~
\left[ \mathcal{C}^*(\tilde{\theta}_{\Psi_{j}}; \gamma_j, \lambda)- \mathcal{C}^*(\theta_{\Psi_{j}}; \gamma_j, \lambda)- \mathcal{C}_{\eta_2}(\tilde{\theta}_{\Psi_{j}};\hat{\gamma}_j,\hat{\lambda}) + \mathcal{C}_{\eta_2}(\theta_{\Psi_{j}};\hat{\gamma}_j,\hat{\lambda}) \right].
\end{align*}

This leads to

\begin{align*}
\mathcal{C}^*(\tilde{\theta}_{\Psi_{j}}; \gamma_j, \lambda)- &\mathcal{C}^*(\theta_{\Psi_{j}}; \gamma_j, \lambda)- \mathcal{C}_{\eta_2}(\tilde{\theta}_{\Psi_{j}};\hat{\gamma}_j,\hat{\lambda}) +\mathcal{C}_{\eta_2}(\theta_{\Psi_{j}};\hat{\gamma}_j,\hat{\lambda})\\ 	&=  \int_\mathcal{T} \left( ( \gamma_j(t)- \lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t)))^2 -
(\hat{\gamma}_j(t)- \hat{\lambda}(\tilde{\theta}_{\Psi_{j}}^T\alpha(t)))^2
\right)dt\\
&\hspace{1cm}+  \int_\mathcal{T} \left((\gamma_j(t)- \lambda(\theta_{\Psi_{j}}^T\alpha(t)))^2 -
(\hat{\gamma}_j(t)- \hat{\lambda}(\theta_{\Psi_{j}}^T\alpha(t)))^2
\right)dt\\
&\hspace{1cm}+\eta_2 \left( \int_\mathcal{T} (\tilde{\theta}_{\Psi_{j}}^T\alpha(t))-t)^2 -
\right)dt \\
&:= I + II + III,
\end{align*}

for $I$, we have

\begin{align*}
I&=\int_\mathcal{T} \left( (\gamma_j(t)- \lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t)))^2 -
(\hat{\gamma}_j(t)- \hat{\lambda}(\tilde{\theta}_{\Psi_{j}}^T\alpha(t)))^2
\right)dt\\
&=\int_\mathcal{T} \bigg[ \left( \gamma_j(t)-\hat{\gamma}_j(t)+
\hat{\lambda}(\tilde{\theta}_{\Psi_{j}}^T\alpha(t)) - \lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t))
\right)\\
& \hspace{3.5cm}
\times \left( \gamma_j(t)+\hat{\gamma}_j(t)-  \lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t)) -
\hat{\lambda}(\tilde{\theta}_{\Psi_{j}}^T\alpha(t))
\right) \bigg] dt\\
&\leq \left(\int_\mathcal{T}  \left| \gamma_j(t)-\hat{\gamma}_j(t)+
\hat{\lambda}(\tilde{\theta}_{\Psi_{j}}^T\alpha(t)) - \lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t))
\right|^2dt\right)^{1/2} \\
& \hspace{3.5cm}
\times\left(\int_\mathcal{T}  \left| \gamma_j(t)+\hat{\gamma}_j(t)-  \lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t)) -
\hat{\lambda}(\tilde{\theta}_{\Psi_{j}}^T\alpha(t))
\right|^2dt\right)^{1/2}.
\end{align*}

For the first term in the product,

\begin{align*}
&\left(\int_\mathcal{T}  \left| \gamma_j(t)-\hat{\gamma}_j(t)+
\hat{\lambda}(\tilde{\theta}_{\Psi_{j}}^T\alpha(t)) - \lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t))
\right|^2dt\right)^{1/2}\\
&\leq\bigg(2\int_\mathcal{T}  \left| \gamma_j(t)-\hat{\gamma}_j(t)\right|^2 dt + 2\int_\mathcal{T}  \left| \hat{\lambda}(\tilde{\theta}_{\Psi_{j}}^T\alpha(t)) - \lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t)) \right|^2 dt \bigg)^{1/2}\\
&\leq \bigg(C_1 ||\hat{\gamma}_j-\gamma_j||^2_\infty + C_2  ||\hat{\lambda}-\lambda||^2_\infty \bigg)^{1/2}\\
&\leq C_3\left(||\hat{\gamma}_j-\gamma_j||_\infty + ||\hat{\lambda}-\lambda||_\infty\right),
\end{align*}

and for the second term, using

$\gamma_j(t) = \lambda(\Psi_j(t)) = \lambda(\theta_{\Psi_j}^T\alpha(t))$,

\begin{align*}
&\left(\int_\mathcal{T}  \left| \gamma_j(t)+\hat{\gamma}_j(t)-  \lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t)) -
\hat{\lambda}(\tilde{\theta}_{\Psi_{j}}^T\alpha(t))
\right|^2dt\right)^{1/2}\\
&= \left(\int_\mathcal{T}  \left| \hat{\gamma}_j(t)- \gamma_j(t) + \lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t)) -
\hat{\lambda}(\tilde{\theta}_{\Psi_{j}}^T\alpha(t)) +2(\gamma_j(t)-\lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t)) )
\right|^2dt\right)^{1/2}\\
&\leq \bigg(2\int_\mathcal{T}  \left|\hat{\gamma}_j(t)- \gamma_j(t) + \lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t)) -
\hat{\lambda}(\tilde{\theta}_{\Psi_{j}}^T\alpha(t))\right|^2 dt \\
+8\int_\mathcal{T}|\lambda(\theta_{\Psi_j}^T\alpha(t))-\lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t))|^2 dt\bigg)^{1/2}.
% 		&\leq \int_\mathcal{T}  \left| \hat{\gamma}_j(t)-\gamma_j(t)\right|dt + \int_\mathcal{T}  \left|
% 		 \hat{\lambda}(\tilde{\theta}_{\Psi_{j}}^T\alpha(t))-\lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t))\right| dt \\
% 		&\hspace{0.5cm}+ 2 \int_\mathcal{T}  \left|\lambda(\Psi_j(t))-\lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t))
% 		\right|dt\\
% 		&\leq C_3 ||\hat{\gamma}_j-\gamma_j||_\infty + C_4 ||\hat{\lambda}-\lambda||_\infty+2\int_\mathcal{T}  \left|\lambda(\theta_{\Psi_{j}}^T\alpha(t))-\lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t))\right|dt\\
% 		&\leq C_3 ||\hat{\gamma}_j-\gamma_j||_\infty + C_4 ||\hat{\lambda}-\lambda||_\infty +C_5 ||\theta_{\Psi_{j}}^T\alpha(t)-\tilde{\theta}_{\Psi_{j}}^T\alpha(t)||_\infty\\
% 		&\leq C_3 ||\hat{\gamma}_j-\gamma_j||_\infty + C_4 ||\hat{\lambda}-\lambda||_\infty +C_5 ||\tilde{\theta}_{\Psi_{j}}-\theta_{\Psi_{j}}||
\end{align*}

Using the  Lipschitz continuity of $\lambda$  and $\sup_t |\alpha_\ell (t)| \leq 1$ for all $\ell=1,...,L+1$,

\begin{align*}
\int_\mathcal{T}|\lambda(\theta_{\Psi_j}^T\alpha(t))-\lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t))|^2
\leq C_4||\theta_{\Psi_{j}}^T\alpha(t)-\tilde{\theta}_{\Psi_{j}}^T\alpha(t)||^2_\infty
\leq C_5  ||\theta_{\Psi_{j}}-\tilde{\theta}_{\Psi_{j}}||^2,
\end{align*}

\begin{align*}
\hspace{-5cm}    \left(\int_\mathcal{T}  \left| \gamma_j(t)+\hat{\gamma}_j(t)-  \lambda(\tilde{\theta}_{\Psi_{j}}^T\alpha(t)) -
\hat{\lambda}(\tilde{\theta}_{\Psi_{j}}^T\alpha(t))
\right|^2dt\right)^{1/2}\\
\leq C_6\big(||\hat{\gamma}_j-\gamma_j||_\infty + ||\hat{\lambda}-\lambda||_\infty + ||\tilde{\theta}_{\Psi_{j}}-\theta_{\Psi_{j}}||\big),
\end{align*}

\begin{align*}
I&\leq C_3\left[ ||\hat{\gamma}_j-\gamma_j||_\infty +  ||\hat{\lambda}-\lambda||_\infty\right] \\
&\hspace{3.5cm} \times C_6\left[ ||\hat{\gamma}_j-\gamma_j||_\infty + ||\hat{\lambda}-\lambda||_\infty + ||\theta_{\Psi_{j}}-\tilde{\theta}_{\Psi_{j}}|| \right]\\
&\leq  C_7\left( ||\hat{\gamma}_j-\gamma_j||_\infty + ||\hat{\lambda}-\lambda||_\infty\right)^2 \\
&\hspace{3.5cm}+ C_{8} \left( ||\hat{\gamma}_j-\gamma_j||_\infty + ||\hat{\lambda}-\lambda||_\infty\right) \times ||\tilde{\theta}_{\Psi_{j}}-\theta_{\Psi_{j}}||.
\end{align*}

Next, $II$ can be simplified. Again noting that $\gamma_j(t)= \lambda(\Psi_j(t)) = \lambda(\theta_{\Psi_{j}}^T\alpha(t))$,

\begin{align*}
&\int_\mathcal{T} \left((\gamma_j(t)- \lambda(\theta_{\Psi_{j}}^T\alpha(t)))^2 -
(\hat{\gamma}_j(t)- \hat{\lambda}(\theta_{\Psi_{j}}^T\alpha(t)))^2
\right)dt = \int_\mathcal{T}
\left(\hat{\gamma}_j(t)- \hat{\lambda}(\theta_{\Psi_{j}}^T\alpha(t))\right)^2
\end{align*}

\begin{align*}
II&= \int_\mathcal{T}
\left((\hat{\gamma}_j(t)-\lambda(\theta_{\Psi_{j}}^T\alpha(t))+\lambda(\theta_{\Psi_{j}}^T\alpha(t))- \hat{\lambda}(\theta_{\Psi_{j}}^T\alpha(t))\right)^2
% 		&=  \int_\mathcal{T}
% 		\left(\hat{\gamma}_j(t)-\lambda(\theta_{\Psi_{j}}^T\alpha(t))\right)^2dt + \int_\mathcal{T}
% 		\left(\lambda(\theta_{\Psi_{j}}^T\alpha(t))- \hat{\lambda}(\theta_{\Psi_{j}}^T\alpha(t))\right)^2 dt\\
% 		&\hspace{3.5cm}+
% 		2 \int_\mathcal{T}
% 		\left(\hat{\gamma}_j(t)-\lambda(\theta_{\Psi_{j}}^T\alpha(t))\right)\times\left(\lambda(\theta_{\Psi_{j}}^T\alpha(t))- \hat{\lambda}(\theta_{\Psi_{j}}^T\alpha(t))\right) dt\\
&\leq 2\int_\mathcal{T}
\left(\hat{\gamma}_j(t)-\gamma_j(t)\right)^2dt + 2\int_\mathcal{T}
\left(\lambda(\theta_{\Psi_{j}}^T\alpha(t))- \hat{\lambda}(\theta_{\Psi_{j}}^T\alpha(t))\right)^2 dt\\
% 		&\hspace{3.5cm}+
% 		2 \left(\int_\mathcal{T} \left|
% 		\hat{\gamma}_j(t)-\gamma_j(t)\right|^2dt\times\int_\mathcal{T}\left|\lambda(\theta_{\Psi_{j}}^T\alpha(t))- \hat{\lambda}(\theta_{\Psi_{j}}^T\alpha(t)\right|^2 dt\right)^{1/2}\\
&\leq C_{9}||\hat{\gamma}_j-\gamma_j||^2_\infty + C_{10}||\hat{\lambda}-\lambda||^2_\infty\\
&\leq C_{11} \left(||\hat{\gamma}_j-\gamma_j||_\infty+||\hat{\lambda}-\lambda||_\infty\right)^2.
\end{align*}

Lastly $III$ is bounded as follows:

\begin{align*}
\eta_2 \bigg(& \int_\mathcal{T} (\tilde{\theta}_{\Psi_{j}}^T\alpha(t))-t)^2 -
\bigg)dt \\
&\hspace{0.5cm}\leq \eta_2 \left( \int_\mathcal{T} (\tilde{\theta}_{\Psi_{j}}^T\alpha(t))-t)^2 +
\right)dt\\
&\hspace{0.5cm}\leq 2\eta_2\int_\mathcal{T}t^2dt\\
&\hspace{0.5cm}\leq 2\eta_2 |\mathcal{T}|^3/3.
\end{align*}

Then combining these bounds, we have

\begin{align*}
\mathcal{C}^*(\tilde{\theta}_{\Psi_{j}}; \gamma_j, \lambda)- &\mathcal{C}^*(\theta_{\Psi_{j}}; \gamma_j, \lambda)- \mathcal{C}_{\eta_2}(\tilde{\theta}_{\Psi_{j}};\hat{\gamma}_j,\hat{\lambda}) + \mathcal{C}_{\eta_2}(\theta_{\Psi_{j}};\hat{\gamma}_j,\hat{\lambda})\\
&\leq I + II + III\\
&\leq (C_7 +C_{11})\left( ||\hat{\gamma}_j-\gamma_j||_\infty + ||\hat{\lambda}-\lambda||_\infty\right)^2\\
&\hspace{.5cm}+ C_{8} \left( ||\hat{\gamma}_j-\gamma_j||_\infty + ||\hat{\lambda}-\lambda||_\infty\right) \times ||\tilde{\theta}_{\Psi_{j}}-\theta_{\Psi_{j}}|| + 2\eta_2 |\mathcal{T}|^3/3\end{align*}
which implies

\begin{align*}
||\tilde{\theta}_{\Psi_{j}}-\theta_{\Psi_{j}}||^2 - &C_0^2 C_{8} \left( ||\hat{\gamma}_j-\gamma_j||_\infty + ||\hat{\lambda}-\lambda||_\infty\right) \times ||\tilde{\theta}_{\Psi_{j}}-\theta_{\Psi_{j}}||\\
\hspace{0.5cm}&\leq C_0 \bigg[(C_7 +C_{11})\left( ||\hat{\gamma}_j-\gamma_j||_\infty + ||\hat{\lambda}-\lambda||_\infty\right)^2 + 2\eta_2 |\mathcal{T}|^3/3\bigg],
\end{align*}

By completing the square

\begin{align*}
\bigg| ||\tilde{\theta}_{\Psi_{j}}-\theta_{\Psi_{j}}|| - &\frac{C_0 C_{8}}{2}\left( ||\hat{\gamma}_j-\gamma_j||_\infty + ||\hat{\lambda}-\lambda||_\infty \right) \bigg|^2\\
\hspace{0.5cm}&\leq C_0 \bigg[(C_7 +C_{11})\left( ||\hat{\gamma}_j-\gamma_j||_\infty + ||\hat{\lambda}-\lambda||_\infty\right)^2 + 2\eta_2 |\mathcal{T}|/3\bigg]\\
&\hspace{0.5cm}+ \frac{C_0^2 C_{7}^2 }{4} \left( ||\hat{\gamma}_j-\gamma_j||_\infty + ||\hat{\lambda}-\lambda||_\infty\right)^2\\
&\leq C_{12}\left( ||\hat{\gamma}_j-\gamma_j||_\infty + ||\hat{\lambda}-\lambda||_\infty\right)^2 + C_0\eta_2 |\mathcal{T}|^3/3,
\end{align*}

and taking the square root, we have

\begin{align*}
\bigg| ||\tilde{\theta}_{\Psi_{j}}-\theta_{\Psi_{j}}|| - &C^*\left( ||\hat{\gamma}_j-\gamma_j||_\infty + ||\hat{\lambda}-\lambda||_\infty \right) \bigg|\\
&\leq C_{13}\left( ||\hat{\gamma}_j-\gamma_j||_\infty + ||\hat{\lambda}-\lambda||_\infty\right) + C_{14}\eta_2^{1/2}.
\end{align*}

Then by the proofs B and C of Theorem 1,
$$C_{13}\left( ||\hat{\gamma}_j-\gamma_j||_\infty + ||\hat{\lambda}-\lambda||_\infty\right) + C_{14}\eta_2^{1/2} = \mathcal{O}_P(n^{-1/2}) + \mathcal{O}_P(\tau_m^{1/2}) +\mathcal{O}(\max(\eta_1,\eta_2)^{1/2}),$$

$$C^*\left( ||\hat{\gamma}_j-\gamma_j||_\infty + ||\hat{\lambda}-\lambda||_\infty\right)=\mathcal{O}_P(n^{-1/2}) + \mathcal{O}_P(\tau_m^{1/2}) +\mathcal{O}(\max(\eta_1,\eta_2)^{1/2}),$$

leading to

\begin{align*}
||\tilde{\theta}_{\Psi_{j}}-\theta_{\Psi_{j}}|| = \mathcal{O}_P(n^{-1/2}) + \mathcal{O}_P(\tau_m^{1/2}) +\mathcal{O}(\max(\eta_1,\eta_2)^{1/2}),
\end{align*}

and with $\sup_t |\alpha_\ell (t)| \leq 1$ for all $\ell=1,...,L+1$,

\begin{align*}
\underset{t\in \mathcal{T}}{\sup} |\hat{\Psi}_j(t)-\Psi_j(t)|&= \underset{t\in \mathcal{T}}{\sup} |(\tilde{\theta}_{\Psi_j}(t)-\theta_{\Psi_j})^T\alpha(t)| \\
&\leq ||\tilde{\theta}_{\Psi_j}-\theta_{\Psi_j}||\\
&=\mathcal{O}_P(n^{-1/2}) + \mathcal{O}_P(\tau_m^{1/2}) +\mathcal{O}(\max(\eta_1,\eta_2)^{1/2}).
\end{align*}

E.  For the overall warping functions $G_{ij}$, we again break the difference into two parts, the rates of which have already been derived, and   use the Lipschitz continuity of $\Psi_j$,

\begin{align*}
\underset{t \in \mathcal{T}}{\sup} |\hat{G}_{ij}(t)-G_{ij}(t)|&= \underset{t \in \mathcal{T}}{\sup} |(\hat{\Psi}_{j}\circ\hat{H}_i)(t)-(\Psi_j\circ \hat{H}_i)(t)+(\Psi_j\circ \hat{H}_i)(t)-(\Psi_j\circ H_i)(t)|\\
&\leq \underset{t \in \mathcal{T}}{\sup} |(\hat{\Psi}_{j}\circ\hat{H}_i)(t)-(\Psi_j\circ \hat{H}_i)(t)|+\underset{t \in \mathcal{T}}{\sup}|(\Psi_j\circ \hat{H}_i)(t)-(\Psi_j\circ H_i)(t)|\\
&\leq \underset{t \in \mathcal{T}}{\sup} |\hat{\Psi}_{j}(t)-\Psi_j(t)|+C\underset{t \in \mathcal{T}}{\sup}| \hat{H}_i(t)-H_i(t)|\\
&= \mathcal{O}_P(n^{-1/2})+ \mathcal{O}_P(m^{-(1-\delta)/6})+\mathcal{O}(\max(\eta_1,\eta_2)^{1/2}).
\end{align*}

F.  The convergence rates for the amplitude factors derived in the proof of Proposition 1.

Proof of Theorem 2.

A.  Bounding the maximum difference between the estimated XCT and its target by

\begin{align*}
\underset{t \in \mathcal{T}}{\sup}|\hat{T}_{jk}(t)- T_{jk}(t)|&= \underset{t \in \mathcal{T}}{\sup}\left|(\hat{\Psi}_j^{-1}\circ \hat{\Psi}_k)(t)- (\Psi_j^{-1}\circ \Psi_k)(t)\right|\\
&\leq \underset{t \in \mathcal{T}}{\sup}\left|(\hat{\Psi}_j^{-1}\circ \hat{\Psi}_k)(t)- (\Psi_j^{-1}\circ \hat{\Psi}_k)(t)\right| \\
&\hspace{0.5cm}+ \underset{t \in \mathcal{T}}{\sup}\left|(\Psi_j^{-1}\circ \hat{\Psi}_k)(t) - (\Psi_j^{-1}\circ \Psi_k)(t)\right|,
\end{align*}

the rate of the first term can be obtained in a straightforward manner by observing,

\begin{align*}
\underset{t \in \mathcal{T}}{\sup}\left|(\hat{\Psi}_j^{-1}\circ \hat{\Psi}_k)(t)- (\Psi_j^{-1}\circ \hat{\Psi}_k)(t)\right| & = \underset{t \in \mathcal{T}}{\sup}\left|(\hat{\Psi}_j^{-1}\circ \hat{\Psi}_k)(\hat{\Psi}^{-1}_k(t))- (\Psi_j^{-1}\circ \hat{\Psi}_k)(\hat{\Psi}^{-1}_k(t))\right|\\
&= \underset{t \in \mathcal{T}}{\sup}\left|\hat{\Psi}_j^{-1}(t)- \Psi_j^{-1}(t)\right|\\
&= \mathcal{O}_P(n^{-1/2}) +  \mathcal{O}_P(m^{-(1-\delta)/6})+ \mathcal{O}( \max(\eta_1,\eta_2)^{1/2}).
\end{align*}

where the rates are as in part D. of the proof of Theorem 1.
For the second term, we use the Lipschitz continuity of $\Psi_j^{-1}$,

\begin{align*}
\underset{t \in \mathcal{T}}{\sup}\left|(\Psi_j^{-1}\circ \hat{\Psi}_k)(t) - (\Psi_j^{-1}\circ \Psi_k)(t)\right|&\leq C\underset{t \in \mathcal{T}}{\sup}\left|\hat{\Psi}_k(t) - \Psi_k(t)\right|\\
&= \mathcal{O}_P(n^{-1/2}) +  \mathcal{O}_P(m^{-(1-\delta)/6})+ \mathcal{O}(\max(\eta_1,\eta_2)^{1/2}),
\end{align*}

again by part D. of the proof of Theorem 1. Then for all $j,k=1,\dots,p,$

\begin{align*}
\underset{t \in \mathcal{T}}{\sup}|\hat{T}_{jk}(t)- T_{jk}(t)|
&\leq \underset{t \in \mathcal{T}}{\sup}\left|(\hat{\Psi}_j^{-1}\circ \hat{\Psi}_k)(t)- (\Psi_j^{-1}\circ \hat{\Psi}_k)(t)\right|\\
& \hspace{0.5cm}+ \underset{t \in \mathcal{T}}{\sup}\left|(\Psi_j^{-1}\circ \hat{\Psi}_k)(t) - (\Psi_j^{-1}\circ \Psi_k)(t)\right|\\
&= \mathcal{O}_P(n^{-1/2}) + \mathcal{O}_P(m^{-(1-\delta)/6})+ \mathcal{O}(\max(\eta_1,\eta_2)^{1/2}),
\end{align*}

B. The argument for the subject-specific cross-component transports follows the same argument as above, noting that
$G_{ij}=\Psi_j\circ H_i$ is Lipschitz continuous,

\begin{align*}
\underset{t \in \mathcal{T}}{\sup}|\hat{T}^{(i)}_{jk}(t)- T^{(i)}_{jk}(t)|&= \underset{t \in \mathcal{T}}{\sup}\bigg|(\hat{G}^{-1}_{ij}\circ \hat{G}_{ik})(t)-(G^{-1}_{ij}\circ \hat{G}_{ik})(t)~\\
&\hspace{3.8cm}+(G^{-1}_{ij}\circ \hat{G}_{ik})(t)-(G^{-1}_{ij}\circ G_{ik})(t)\bigg|\\
%&\leq \underset{t \in \mathcal{T}}{\sup}\left|(\hat{G}^{-1}_{ij}\circ \hat{G}_{ik})(t)-(G^{-1}_{ij}\circ \hat{G}_{ik})(t)\right|\\
%&\hspace{1cm}+\underset{t \in \mathcal{T}}{\sup}\left|(G^{-1}_{ij}\circ \hat{G}_{ik})(t)-(G^{-1}_{ij}\circ G_{ik})(t)\right|\\
&\leq \underset{t \in \mathcal{T}}{\sup}\left|\hat{G}^{-1}_{ij}(t)-G_{ij}(t)\right| + C\underset{t \in \mathcal{T}}{\sup}\left|\hat{G}_{ik}(t)-G_{ik}(t)\right|\\
&=\mathcal{O}_P(n^{-1/2}) + \mathcal{O}_P(m^{-(1-\delta)/6})+ \mathcal{O}(\max(\eta_1,\eta_2)^{1/2}),
\end{align*}

for all $i=1,\dots,n,~j,k=1,\dots,p,$ by Theorem 1E.
This completes the proof.

Proof of Theorem 3.

The reconstructed curves depends on the rates of convergence for the $A_{ij}$, $\lambda$, and $G_{ij}$ terms, which are established in Theorem 1. Now

\begin{align*}
\underset{t \in \mathcal{T}}{\sup}|\hat{X}_{ij}(t)-&X_{ij}(t)|= \underset{t \in \mathcal{T}}{\sup}|\hat{A}_{ij}(\hat{\lambda}\circ \hat{G}_{ij})(t)-A_{ij}(\lambda\circ G_{ij})(t)|\\
&\leq\underset{t \in \mathcal{T}}{\sup}|\hat{A}_{ij}(\hat{\lambda}\circ \hat{G}_{ij})(t) -\hat{A}_{ij}(\lambda\circ G_{ij})(t)| \\
&\hspace{0.5cm}+\underset{t \in \mathcal{T}}{\sup}|\hat{A}_{ij}(\lambda\circ G_{ij})(t)-A_{ij}(\lambda\circ G_{ij})(t)|\\
&=\hat{A}_{ij}\underset{t \in \mathcal{T}}{\sup}|(\hat{\lambda}\circ \hat{G}_{ij})(t) -(\lambda\circ G_{ij})(t)|+\underset{t \in \mathcal{T}}{\sup}|(\hat{A}_{ij}-A_{ij})(\lambda\circ G_{ij})(t)|\\
&:=I + II,
\end{align*}

$$I \leq\hat{A}_{ij}\left(\underset{t \in \mathcal{T}}{\sup}|(\hat{\lambda}-\lambda)(t)| + C\underset{t \in \mathcal{T}}{\sup}| \hat{G}_{ij}(t)- G_{ij}(t)|\right)$$

\begin{align*}
&\hat{A}_{ij}\left(\underset{t \in \mathcal{T}}{\sup}|(\hat{\lambda}-\lambda)(t)| + C\underset{t \in \mathcal{T}}{\sup}| \hat{G}_{ij}(t)- G_{ij}(t)|\right)\hspace{7cm}\\
&\hspace{2cm}\leq A_{ij}\left(\underset{t \in \mathcal{T}}{\sup}|(\hat{\lambda}-\lambda)(t)| + C\underset{t \in \mathcal{T}}{\sup}| \hat{G}_{ij}(t)- G_{ij}(t)|\right)+\\
&\hspace{4cm}|\hat{A}_{ij}-A_{ij}|\left(\underset{t \in \mathcal{T}}{\sup}|(\hat{\lambda}-\lambda)(t)| + C\underset{t \in \mathcal{T}}{\sup}| \hat{G}_{ij}(t)- G_{ij}(t)|\right)\\
&\hspace{2cm}=\mathcal{O}_P(1)\left[\mathcal{O}_P(n^{-1/2}) + \mathcal{O}_P(m^{-(1-\delta)/6})+ \mathcal{O}(\max(\eta_1,\eta_2)^{1/2})\right]+\\
&\hspace{4cm}\mathcal{O}_P(m^{-1/6})\left[\mathcal{O}_P(n^{-1/2}) + \mathcal{O}_P(m^{-(1-\delta)/6})+ \mathcal{O}(\max(\eta_1,\eta_2)^{1/2})\right]\\
&\hspace{2cm}=\mathcal{O}_P(n^{-1/2}) + \mathcal{O}_P(m^{-(1-\delta)/6})+ \mathcal{O}(\max(\eta_1,\eta_2)^{1/2})
\end{align*}

leads to

\begin{align*}
II %\underset{t \in \mathcal{T}}{\sup}|(\hat{A}_{ij}-A_{ij})(\lambda\circ G_{ij})(t)|&
\leq|\hat{A}_{ij}-A_{ij}|~\underset{t \in \mathcal{T}}{\sup}|(\lambda\circ G_{ij})(t)|
\end{align*}

Combining these rates gives

\begin{align*}
\underset{t \in \mathcal{T}}{\sup}|\hat{X}_{ij}(t)-X_{ij}(t)|%&=\mathcal{O}_P(n^{-1/2}) + \mathcal{O}_P(m^{-(1-\delta)/6})+ \mathcal{O}(\max(\eta_1,\eta_2)^{1/2})+\mathcal{O}_P(m^{-1/6})\
=\mathcal{O}_P(n^{-1/2}) + \mathcal{O}_P(m^{-(1-\delta)/6})+ \mathcal{O}(\max(\eta_1,\eta_2)^{1/2}),
\end{align*}

which completes the proof.

Proof of Corollaries 1-3.

The arguments for all Corollaries 1-3 are almost identical so we show Corollary 1 only.
If $\eta_1\sim\eta_2=\mathcal{O}(n^{-1})$ and $m(n)=\mathcal{O}(n^{\Delta(1-\delta)^{-1}})$, with $\Delta>3$, then

\begin{align*}
\underset{t \in \mathcal{T}}{\sup}|\hat{H}_{i}(t) - H(t)| &= \mathcal{O}_P(n^{-1/2}) + \mathcal{O}_P(\tau_m^{1/2})+ \mathcal{O}(\eta_1^{1/2})\\
&=\mathcal{O}_P(n^{-1/2}) + \mathcal{O}_P(n^{-\Delta/6}) + \mathcal{O}_P(n^{-1/2}).
\end{align*}

The result follows by observing that
$\Delta/6>1/2$ if $\Delta>3$ and thus,

\begin{align*}
\underset{t \in \mathcal{T}}{\sup}|\hat{H}_{i}(t) - H(t)|
\end{align*}

Following previous arguments we also have

and the result follows.